Timescale analyses of fluctuations in coexisting populations of a native and invasive tree squirrel

Abstract Competition from invasive species is an increasing threat to biodiversity. In Southern California, the western gray squirrel (Sciurus griseus, WGS) is facing competition from the fox squirrel (Sciurus niger, FS), an invasive congener. We used spectral methods to analyze 140 consecutive monthly censuses of WGS and FS within a 11.3 ha section of the California Botanic Garden. Variation in the numbers for both species and their synchrony was distributed across long timescales (>15 months). After filtering out annual changes, concurrent mean monthly temperatures from nearby Ontario Airport yielded a spectrum with a large semi‐annual peak and significant spectral power at long timescales (>28 months). The cospectrum between WGS numbers and temperature revealed a significant negative correlation at long timescales (>35 months). Cospectra also revealed significant negative correlations with temperature at a six‐month timescale for both WGS and FS. Simulations from a model of two competing species indicate that the risk of extinction for the weaker competitor increases quickly as environmental noise shifts from short to long timescales. We analyzed the timescales of fluctuations in detrended mean annual temperatures for the time period 1915–2014 from 1218 locations across the continental USA. In the last two decades, significant shifts from short to long timescales have occurred, from <3 years to 4–6 years. Our results indicate that (i) population fluctuations in co‐occurring native and invasive tree squirrels are synchronous, occur over long timescales, and may be driven by fluctuations in environmental conditions; (ii) long timescale population fluctuations increase the risk of extinction in competing species, especially for the inferior competitor; and (iii) the timescales of interannual environmental fluctuations may be increasing from recent historical values. These results have broad implications for the impact of climate change on the maintenance of biodiversity.


| INTRODUC TI ON
Competition from non-native, invasive species is an increasing threat to the biodiversity of native species in a globalized world. Invasive species are often considered one of the most important threats to ecological function and a top driver of species extinctions (Dueñas et al., 2021;Flory & Lockwood, 2020). The presence of invasive species can alter animal communities, trigger trophic cascades, displace native species, and even lead to hybridizations with similar or related species (Doody et al., 2017;Huxel, 1999). The ability to be more competitive over limited resources is one of the characteristics that enables invasive species to be successful. In addition, they are often characterized by having life history traits with colonizer characteristics as follows: short generation times, high reproduction rates, and fast growth rates (Sakai et al., 2001). With this competitive edge, they can invade and displace native species.
An example where a native species is threatened in some habitats by competition from an invasive species occurs in Southern California, where the western gray squirrel (Sciurus griseus, WGS,  (Carraway & Verts, 1994;Escobar-Flores et al., 2011). Populations of WGSs have been declining in areas of Washington, Oregon, and California (Cooper, 2013;Cooper & Muchlinski, 2015;Muchlinski et al., 2009;Stuart, 2012). In Washington, they are listed as a statethreatened species (Linders & Stinson, 2007), while in Oregon they are an Oregon Conservation Strategy Species (Oregon Department of Fish & Wildlife, 2016). While there have been only a few studies regarding populations of WGSs in California, there is a noticeable trend in the decline of these squirrels in areas below an elevation of 457 m (Cooper, 2013;Cooper & Muchlinski, 2015). As of now, the WGS does not have special conservation status in California.
Fox squirrels have dispersed from original points of introduction through natural dispersal and through intentional movement of animals by humans (Frey & Campbell, 1997;Geluso, 2004;King et al., 2010). Since the original introduction to Los Angeles County (Becker & Kimball, 1947), the FS has expanded its range at a rate of 1.60 to 3.00 km/year in heavily suburbanized areas of Southern California (Garcia & Muchlinski, 2017). Although the FS has generally remained restricted to areas of human habitation, with continued range expansion the FS has become sympatric in some isolated suburban habitat fragments and in certain foothill areas with the native WGS (Hoefler & Harris, 1990).
FSs may compete with native WGSs for resources such as nesting sites and food, and the FS has replaced the WGS within certain habitats in Southern California (Cooper & Muchlinski, 2015;Muchlinski et al., 2009). Los Angeles County can be considered an ideal location for invasion by the FS given the mild Mediterranean climate and year-round food supply offered by exotic plant species, accompanied by the absence of the native WGS throughout much of the Los Angeles Basin. The FS is both morphologically, ecologically, and behaviorally similar to this native species; thus, these overlaps in form, function, activity, and presence provide a situation where interactions between the two species can be studied (Ortiz, 2021).
Many factors can influence population persistence, but one that has received comparatively less attention is the timescale of environmental and population fluctuations. In the current study, we investigated the effects of timescales using three approaches. First, we applied spectral analyses to examine the timescale distribution of the variance and covariance of WGS, FS, and weather time series. Second, we conducted model simulations to examine the implications of changes in the timescale of environmental fluctuations on the coexistence of a native species facing competition from an invasive species. Third, we used spectral and wavelet methods to examine the long-term changes in the timescale distribution of climate data.
By analogy with the spectrum of visible light, time series fluctuations that occur over long timescales are referred to as having a red spectrum and those occurring over short timescales as having a blue spectrum (Lawton, 1988). These are distinguished from white noise random fluctuations, which have no serial autocorrelations. In general, theoretical analyses from single-species discrete-time unstructured population models suggest that the response to colored environmental noise depends on the type of population dynamics.
Deterministic models with stable equilibria can exhibit undercompensatory dynamics, where the population approaches equilibrium monotonically, or overcompensatory dynamics, where the approach to equilibrium exhibits damped oscillations. Many studies have shown that reddened environmental spectra increase extinction risk for undercompensatory populations and blue spectra increase extinction risk for overcompensatory populations (Danielian, 2016; García-Carreras & Reuman, 2011; Mustin et al., 2013;Petchey et al., 1997;Ripa & Heino, 1999;Ripa & Lundberg, 1996;Ruokolainen et al., 2009;Schwager et al., 2006). While some studies reached different conclusions (Heino, 1998;Heino et al., 2000), Schwager et al. (2006 showed that these contrasting results depend on the modeling details and a consideration of the likelihood of catastrophic events. In their simulations of three competing species, Ruokolainen and Fowler (2008) found that extinction risk increased with reddened environmental noise when species responded independently to the environment but decreased when there was a strong correlation between species-specific responses. On the empirical side, Pimm and Redfearn (1988) looked at 100 time series from insects, birds, and mammals and found that the variance of the population fluctuations increased with the window of time used in the calculation, suggesting that these populations have red spectra.
García-Carreras and Reuman (2011) analyzed the dynamics of 147 animal populations and climate data for the population locations and found a positive correlation between the biotic and climatic spectral exponents (a measure of spectral color), with most spectra being red-shifted. Inchausti and Halley (2003) directly examined the relationship between population variability and quasi-extinction time (measured as the time required to observe a 90% decline of population abundance) for a large set of data comprised of 554 populations for 123 animal species that were censused for more than 30 years.
The results showed that the quasi-extinction time was shorter for populations having higher temporal variability and redder dynamics.
In a laboratory microcosm experiment, Fey and Wieczynski (2017) looked at how the autocorrelation in thermal warming affected the ability of a non-native cladoceran, Daphnia lumholtzi, to establish itself in the presence of a native congener, D. pulex. The non-native species was able to attain significantly higher population densities in the treatment with an autocorrelated warming regime relative to the treatment with uncorrelated warming but the same mean temperature and the unwarmed control. Although, in all three treatments, D. lumholtzi went extinct by the end of the experiment, their results demonstrated that the timescale of environmental fluctuations can impact the ability of an invasive species to establish itself in the presence of a native competitor.
Spectral methods are a powerful tool for characterizing the timescales of fluctuations in a time series (Brillinger, 2001 The cross-spectrum is a complex-valued function of frequency.
The real part is the cospectrum, which describes the distribution of the in-phase covariance between the time series at different frequencies, and the imaginary part is the quadrature spectrum, which is a phase-shifted covariance. The sum of the cospectral powers across frequencies is proportional to the total covariance of the two time series. The cospectrum can also be viewed as the distribution of the correlation coefficient across frequencies.
Since frequency, f, is the inverse of the period, the spectral and cospectral power provide information on the variance and correlation, respectively, at the timescale 1/f. The color of a power spectrum can be characterized using a spectral exponent (García-Carreras & Reuman, 2011;Vasseur & Yodzis, 2004). If S f is the power of the spectrum at frequency f, the spectral exponent can be computed as the slope of a least-squares linear regression of log(S f ) versus log(f). Negative spectral exponents are characteristic of spectra dominated by long timescale variation (red spectra), and positive values are indicative of short timescale variation (blue spectra). White noise spectra will have a spectral exponent of zero. When applied to environmental and population time series, spectral color allows one to better assess the risk of ecological extinction.
Wavelet analyses have been used in ecology to identify changes in the spectral distributions of population and environmental fluctuations over time (Cazelles et al., 2008

| Spectral analyses of census data
We used spectral methods to analyze the monthly census data.
We used fast Fourier transforms to compute the raw spectra and cross-spectrum of the bivariate time series. Computations were conducted using the spec.pgram algorithm from R modified to run in MATLAB. No trends were removed from the data prior to analysis. Since raw spectra and cross-spectra are usually jagged, we applied two iterations of a window-averaging smoothing Daniell kernel with spans of 5 and 7 data points, modified with clipped windows at the endpoints to preserve the number of data values. We divided the spectral powers by their sum across frequencies. This yielded a normalized spectral power plot for each species, which shows the distribution of variation across timescales. We used the real part of the cross-spectrum to obtain a smoothed cospectral power plot for the covariance between the two species. We normalized the cospectrum so that its sum equals the correlation coefficient.
We conducted computations to detect significant (p < .05) peaks or valleys in the observed spectra for the null hypothesis that there is no frequency dependence in the variance and covariance of the time series fluctuations (i.e., independent "white noise" time series).
We shuffled the temporal order of the bivariate time series by generating a random permutation of the integers 1 through n, where n = 140 is the number of monthly observations. We then used the permutation to reorder the bivariate monthly censuses of the two species. Next, we computed two smoothed normalized spectra and a smoothed normalized cospectrum in the same way as the correctly ordered data. We repeated this random reshuffling process 2000 times. For the spectra, which must be nonzero, we used the 95th per-

| Analyses of weather data
We obtained weather data for Ontario Airport (ONT) from the Climate Data Online website of NOAA's National Centers for Environmental Information (https://www.ncdc.noaa.gov/cdo-web/).
ONT is located about 12 km from the CBG and should be an accurate representation of the temperature profile of the study site. We focused on the reported "average monthly temperature," which is computed by averaging the daily maximum and minimum temperatures for each month. We avoided rainfall totals because many months have zeros, which is a problem for spectral analyses, and much of the vegetation in the CBG is irrigated. We obtained a temperature time series for the same months as the census data and applied the same spectral methods to obtain a smoothed normalized power spectrum.
Since annual seasonal changes dominated the temperature time series, we used the MATLAB "band-stop" function to attenuate cyclic components with periodicities in the range of 9-15 months.
This produced a filtered time series with annual effects removed.
We then produced a smoothed normalized spectrum for the filtered temperature time series. We also generated smoothed normalized cospectra between the filtered temperature time series and both the WGS and FS census time series. Using the methods described above, we obtained 95% confidence intervals for these spectra and cospectra.

| Model simulations
We conducted model simulations to obtain a better understanding of the implications of timescale-specific environmental variation on the dynamics of two competing species. We used the following discrete-time version of the Lotka-Volterra competition equation: where r 1 and r 2 are the intrinsic rates of population increase, K 1 and K 2 are the carrying capacities, and α and β are the competition coefficients for the two species. The variables ε 1 (t) and ε 2 (t) represent random environmental noise with a mean of zero and variance of .5. We used the coefficients σ 1 and σ 2 to scale the magnitude of the noise. For the purposes of discussion, species 1 will represent a native species and species 2 will represent an invasive species.
We introduced frequency-specific biases into the noise variables using an algorithm devised by Chambers (1995). This method generates a multivariate random time series based on any specified theoretical spectral matrix that is a function of frequency. The diagonal elements of that matrix are the theoretical spectra (frequency decompositions of the variances), and the off-diagonal elements are theoretical crossspectra (complex numbers). The real parts of the cross-spectra are the theoretical cospectra (frequency decompositions of the covariances), and the complex parts are the quadrature spectra (frequency-specific phase shifts). For the model (1), we used identical spectra that were linear functions of frequency for the two species. High-frequency- between the random variables ε 1 (t) and ε 2 (t), we used a cospectrum function that was equal to a constant fraction, 0.9, of the spectrum.
This resulted in a frequency-specific correlation of 0.9 across all frequencies. A high correlation was used since it was assumed that the native and invasive species are ecologically similar and occupy the same habitat. The quadrature spectrum was set to zero (no frequencyspecific phase shifts). To summarize, the timescales of the random environmental noise were varied from short (blue) to uniform (white) to long (red) with a frequency-independent correlation in the effects of the noise on the growth of the two species.
In addition to the spectral frequency of the environmental noise, the simulation protocol also involved varying the competition coefficient α, which represents the intensity of the competitive effects of the invasive species on the native species. We set the value of the α to .25 (weak competition), .50 (moderate competition), and .75 (strong competition). We kept the competitive effects of the native species on the invasive species at a value of β = .25 while increasing α because we are interested in situations of concern to conservationists where an invasive species outcompetes the native species.
We chose values for the intrinsic rate of increase r 1 = r 2 = .3 that are appropriate for tree squirrels of the genus Sciurus (Appendix A). The remaining model parameters had constant values of K 1 = K 2 = 50, and σ 1 = σ 2 = 0.75. For the assessment of extinction risk, when the population density of a species fell below 5% of its carrying capacity, we set it to zero. For simulations not involving extinction risk, the threshold was set to zero. We ran each simulation for 100 time steps.
For every set of parameter values and environmental noise color, we conducted 2000 replicate simulations. For blue, white, and red environmental noise, we computed smoothed normalized power spectra and cospectrum of the species and averaged these over replicates to see how the timescale for population fluctuations is affected by different colors of noise. To investigate gradual shifts in the effects of frequency-biased environmental noise on the population spectra and probability of extinction, we chose a slope for the environmental spectra, varying the slopes in 101 gradual increments, beginning at blue noise (slope = 2) and ending at red noise (slope = −2). For each choice of the environmental spectra, we simulated the population trajectories of the two species, estimated the unsmoothed normalized population spectra, computed the two spectral exponents, and averaged them.
In the same way, we computed average spectral exponents for the environmental noise realizations. We repeated these computations for each of the 2000 replicate simulations and computed an overall average for the spectral exponents of the populations and noise. To investigate the effects of frequency-biased environmental noise on ecological persistence, the number of instances where each population went extinct was divided by 2000 to yield estimates of the extinction risks for both the native and the invasive species.

| Analyses of climate data
We obtained climate data from the U.S. Historical Climatology Network (USHCN), which is freely available online (https://www. ncei.noaa.gov/produ cts/land-based -stati on/us-histo rical -clima tolog y-network). We used version 2.5 of the monthly temperature records, which contains long-term data from 1218 stations across the continental United States. Menne et al. (2009)  Although it would be tempting to analyze the changes in the spectral exponents using a repeated measures ANOVA, with stations as the subjects, spatial autocorrelations exist among stations that are in the same geographical proximity, inflating the Type I error rates. A solution to this problem was suggested by Clifford et al. (1989) and modified by Dutilleul (1993), which yields an "effective sample size" based on the spatial structure of the data. It is appropriate for paired observations distributed in space. We used the software package SAM (Spatial Analysis in Macroecology; Rangel et al., 2006) to compute effective sample sizes for the following three sets of paired data:  vs. ,  vs. , and  vs. . We conducted paired sample t-tests for the spectral exponents from these three paired data sets and adjusted the standard errors for the test statistics and degrees of freedom for the statistical significance values using the effective sample sizes. We then applied a Bonferroni correction to account for the multiple comparisons.
We also conducted a mean field wavelet analysis on the 100-year time series of mean annual temperatures. For each station, we detrended the time series using a quadratic polynomial and computed the standardized residuals. Next, we used the MATLAB continuous wavelet transform function "cwt" to compute wavelet powers for the residual time series using the analytic Morse filter (Olhede & Walden, 2002) with the default values of 3 for the symmetry parameter and 60 for the time-bandwidth product. Lastly, we averaged the wavelet powers across all stations for each time-frequency combination. We chose the Morse wavelet because it is useful for analyzing signals with time-varying amplitude and frequency. We investigated varying the symmetry and time-bandwidth product parameters, but the results were not much different from what was obtained using the default values. We also used a Morlet wavelet which has equal variance in time and frequency, but, again, the results were like the Morse wavelet with default parameters. We experimented with cubic and quartic polynomials for detrending, but these gave mean field wavelets that were much like the one obtained with a quadratic function. Our mean field wavelet differs from the "wavelet mean field" defined by Sheppard et al. (2016), which was designed as an average measure of the time-and timescale synchrony for time series from different locations.
To identify wavelet powers that were statistically significant, we used the surrogate time series approach (Schreiber & Schmitz, 2000). We took a random permutation of the mean annual temperature time series for all stations in tandem and computed a mean field wavelet as described above. We repeated this process 2000 times and computed the upper 95th percentile of the wavelet powers for each combination of time and frequency. This provided a set of critical values for identifying "hot spots" on the mean field wavelet under the null hypothesis of no timescale dependence in the fluctuations of the mean annual temperature residuals around the trends. The estimated Pearson correlation coefficient in animal numbers is R = .581 which is statistically significant from zero (p = 5.20 × 10 −14 ).

| Census data
The total variation in the numbers for each species and their synchrony seems to be distributed across different timescales. Long intervals can be seen where the numbers of both species are elevated and depressed ( Figure 2). Superimposed on this long timescale variation are random short timescale fluctuations. We quantify this timescale component of variation with the spectral analyses in the next section. The spectrum for the FS shows a small peak at 6 months, but that peak is not statistically significant. Since the total variation remains constant across frequencies for the confidence bands from the randomly ordered data, the larger variation in WGS and FS at long timescales is compensated for by smaller variation at timescales of about 4 months or less.

| Population spectra and cospectrum
The smoothed normalized cospectrum (Figure 3c) shows how the total correlation in population numbers between the two species is distributed across timescales. Covariance between WGS and FS is significantly biased toward long timescales, with a smaller nonsignificant peak at a timescale of around 6 months. The cospectrum crosses the upper significance threshold at timescale of about 18 months. The total correlation between the numbers of the WGS and FS is R = .581. Using the unsmoothed cospectrum, we can partition this total correlation by timescale intervals: R 1 = .409 for >12 months, R 2 = .162 for 4-12 months, and R 3 = .010 for ≤4 months, where R = R 1 + R 2 + R 3 . Thus, 70% of the total correlation occurs at timescales exceeding one year. We can infer that population synchrony for these two species occurs mostly at long timescales. Ontario Airport (ONT), which is 12 km from the study site. As one would expect, there is a strong seasonal component to these temperatures. Figure 4b shows the smoothed normalized spectrum for the mean monthly temperatures, which is dominated by a strong peak for the annual cycle. Since the squirrel spectra show no indication of an annual cycle (Figure 3), we applied a band-stop filter to remove the annual cycle and plotted the resulting time series (Figure 4a, dashed line). The smoothed normalized spectrum for the filtered mean monthly temperatures appears in Figure 4c. There is a peak at low frequencies which crosses the upper threshold for statistical significance at a timescale of approximately 28 months and reaches a minimum at a timescale of 12 months. There is also a large spectral peak at 6 months.
There is a negative correlation between the filtered temperature time series and the squirrel census data. For the WGS, the correlation is statistically significant (R = −.194, p = .022) and, as indicated by the smoothed normalized cospectrum (Figure 5a), is distributed at long timescales (>35 months) and at a timescale of 6 months. The correlation between the filtered temperature time series and the FS census data is also negative, but not statistically significant overall (R = −.146, p = .085). The cospectrum between the filtered temperature time series and FS census data shows a large significant peak at a 6-month timescale (Figure 5b). These results suggest that the distribution of variation in the squirrels' population fluctuations may be driven, in part, by fluctuations in weather and climate outside of the annual seasonal cycle.

| Simulation results
Our analyses of the simulations of the Lotka-Volterra competition model (1)     If the non-native species has a competitive advantage, the influence of reddened environmental spectra on the persistence of the native species becomes more pronounced. Figure 6d shows how increasing the competition coefficient for the invading species to α = .50 and α = .75 increases the likelihood that the native species will be lost, while slightly lowering the extinction risk for the invasive species. For, α = .75, a reddening of the environmental spectrum quickly elevates the probability of extinction for the native species from a value of about 3% for the bluest environmental noise to a value which asymptotes at about 94% for the reddest environmental noise (Figure 6d). This suggests the possibility of a synergy between the effects of reddening environmental noise and competition from invasive species for the risk of extinction for native populations.

| Climate data
We know that human impact on   1915-1939, 1940-1964, 1965-1989, and 1990-2015, respectively. It appears that there was a shift from 1990 to 2014 from blue-shifted spectra to red-shifted spectra. The significance values for the changes be- Our spectral analyses of the WGS and FS census data suggest that most of the variation in animal numbers occurs on timescales that exceed 15 months. In the case of the FS, there is also evidence for variation on a six-month timescale. This timescale-specific variation may be due to changes in resource abundance, the timing and frequency of reproduction, and reproductive output.
Changes in population numbers on a long timescale could be due to variation in the supply of food resources on multi-year, highly variable timescales. For example, acorns provide a valuable source of food for tree squirrels (Steele & Yi, 2020), but a very large (>600 g/ m 2 ) mast crop was only produced in one of the nine years in which we measured relative acorn production (Appendix B). We observed the production of a very large mast crop within our study area in the fall of 2013 (Table A2) -1939 1940-1964 1965-1989 1990-2014 (a) (c) period without a modest to large-sized acorn crop corresponded to relatively low census counts for both species (≤ 20 animals). A modest acorn crop produced in the fall of 2020 again corresponded to an increase in census counts for both species during the summer and fall of 2020. Acorns are present in the trees for a prolonged period before they appear in significant quantities on the ground, so this food resource is also available to the animals prior to the fall of the year which may account for the high census counts in the summers prior to our acorn crop sampling periods.
Acorn production by coastal live oaks (Quercus agrifolia), a common tree species within our study area, is influenced by the amount of rain in the one or two years prior to the year in which acorns are produced (Koenig et al., 1996). Mann and Gleick (2015), as well as Diffenbaugh et al. (2015), documented that an increase in ambient temperatures has been accompanied by a decrease in rainfall within California. A plot of temperature and precipitation anomalies over the period of 1895 through November 2014 showed the 3-year period ending in 2014 was by far the hottest and driest on record in California (Mann & Gleick, 2015). Diffenbaugh et al. (2015) documented that although there has not been a large change in the probability of either negative or moderately negative precipitation anomalies in recent decades, the occurrence of drought years has been greater in the two decades prior to their study than in the preceding century. In addition, the probability that precipitation deficits co-occur with warm conditions and the probability that precipita-

Wavelet power
Two distinct periods of potential reproduction for the FS in Southern California were documented by King's (2004) study of 135 L submitted to three wildlife rehabilitation centers during 2002.
Approximately, 60% of litter production documented by King (2004) was associated with the months of February, March, and April, with the largest number of litters born in March. A second pulse of litter production occurred during the months of August, September, and October with the largest number of litters born in September, six months after the largest pulse of litters born during the spring.
Although production of litters by the FS on a semi-annual basis is possible, thus leading to an increase in observed population size on a semi-annual basis, the number of juvenile/subadult animals observed during census counts in this study varied widely among years ( Figure A1).
The WGS appears to exhibit a yearly pattern of reproduction different than the FS. Most research documents breeding activity in late fall and early winter months with birth of most litters in spring and early summer months (Carraway & Verts, 1994;King, 2004).
A few pregnant females were observed in June, July, August, and September (Fletcher, 1963), and lactating females have been observed as late as October in Californian (Swift, 1977). However, no definite records of multiple pregnancies not attributable to intrauterine loss of the first litter are available (Bailey, 1936;Fletcher, 1963;Jameson & Peeters, 1988;Swift, 1977). The difference in reproductive patterns between the FS and the WGS could bring about the presence of a 6-month cycle in abundance of the FS and the absence of a similar 6-month cycle in the WGS. The difference in reproductive patterns could also give a competitive advantage to the FS in certain habitats through higher natality in years of good resource production.

Muchlinski et al. (2012) produced a Habitat Suitability Model
(HSM) for the WGS and the FS which allowed short-term and longerterm coexistence habitats to be identified using a linear combination of three habitat variables: percent canopy cover, percent of deciduous trees, and average height of ground cover. Habitats with a low percentage of canopy cover, a high percentage of deciduous trees, and a low height of ground cover were classified as short-term coexistence habitats. Locations with a high percentage of canopy coverage, a low percentage of deciduous trees, and a low height of ground cover were classified as longer-term coexistence sites. (Sites with a high height of ground cover, a high percentage canopy cover, and a low percentage deciduous trees were identified as "exclusion habitats" where only the WGS is found, but the FS exists in adjacent habitats.) For example, Muchlinski et al. (2009) King, 2004;King et al., 2010). The study area at CBG has been classified as a longer-term coexistence habitat by Muchlinski et al. (2012). How long coexistence can continue in longer-term coexistence habitats is unknown. Many longer-term coexistence sites are fragments of habitat where the FS, but not the WGS, exists in surrounding habitats. The WGS is also subject to loss of genetic diversity in these habitat fragments as described by DeMarco et al. (2020).
The predictions of the competition model presented in section 3.4 can be interpreted in terms of the HSM developed by Muchlinski et al. (2012). The HSM implies that the competitive effects of the FS on the WGS are high in a short-term coexistence site such as California State Polytechnic University, Pomona, and other former lowland coexistence sites (Cooper & Muchlinski, 2015). In terms of the competition model presented in Figure 6, Appendix D, we used least squares to fit a model with annual and semi-annual harmonics to the unfiltered mean monthly temperature data in Figure 4a and show that a model that includes both annual and semi-annual cycles provides a significantly better fit to the data than a model based on the annual cycle alone. Figure 4c shows that the filtered temperature time series has the same period, phase, and approximate amplitude as the theoretical semi-annual cycle. One can also see the effects of long timescale variability in the way the filtered data meanders above and below the semi-annual cycle. The cospectra of Figure 5 indicate a significant negative correlation between the squirrel census data and the ONT filtered weather data at a timescale of six months. While we cannot demonstrate a direct causal mechanism for this correlation, this observation could motivate further research.
It was not possible to specify an estimated peak value for the timescale of low-frequency variation in the squirrel numbers or mean monthly temperatures. The 140 months of the time series represent less than 12 years. In spectral analyses, estimates of longperiod, low-frequency cycles are less precise since they cannot be as readily observed as short-period, high-frequency oscillations. In the estimated spectra of Figures 3 and 4, the spectral power continues to increase as the frequency decreases. However, the record of mean annual temperatures for Ontario Airport extends back to 1999, providing a 22-year time series. In Appendix D, we show that a significant peak in the spectrum of annual temperatures occurs on a timescale of about 7 years, which is consistent with the wavelet analysis in Figure 9. If annual changes in environmental conditions are driving the long-term variation in squirrel numbers, which seems to be the case for the WGS (Figure 5a), this estimate could also represent the timescale of those fluctuations.
We presented simulation results in section 3.4 that were designed to explore the effects of changes in the timescale of environmental noise on the outcome of competition between ecologically similar native and non-native species. We showed that an increase in the timescale of environmental noise reddens the spectrum of population fluctuations and decreases the likelihood of coexistence, especially when the non-native is a better competitor. This result differs from one of the findings of Ruokolainen and Fowler (2008), who concluded that extinction risk decreased with a reddening of environmental noise when, like in our model, there was a strong correlation in the species response to the environmental fluctuations. However, their simulation protocols differed from ours in several ways. First, they looked at a community of three competing species. Second, their environmental noise was generated using an autoregressive process and was added to the carrying capacity for each species. Third, and most importantly, in their models the intrinsic rate of increase for each species was set to r i = 1.8, whereas in our model we chose r 1 = r 2 = .3, which is more consistent with the reproductive capabilities of tree squirrels (Appendix A). In deterministic models of the type used in our simulations and those of Ruokolainen and Fowler (2008), values of 1 > r > 2 lead to overcompensating dynamics, where the approach to equilibrium exhibits damped oscillations. Although Ruokolainen and Fowler (2008) claimed that their results are qualitatively similar for values of r < 1, previous work with single-species models (e.g., Danielian, 2016) has shown that extinction risk increases with a reddening of environmental noise when the deterministic model, like ours, has undercompensating dynamics (r i < 1), but decreases with a reddening of environmental noise when the deterministic model has overcompensating dynamics (1 < r i < 2). Since the apparent contradiction between our results and those of Ruokolainen and Fowler (2008) occurred when there was a strong synchrony among the three species due to the high correlation of the effects of environmental noise, their result is consistent with what one would expect from a single-species model. When we repeated our simulations using values of r 1 = r 2 = 1.8, we observed a result consistent with Ruokolainen and Fowler (2008).
Our simulation results were limited in their scope. They were

A PPE N D I X A I NTR I N S I C R ATE S O F I N CR E A S E FO R S Q U I R R E L S
We used the Euler-Lotka equation to compute the intrinsic rate of increase, r, for the Sciurus species in table 1 of Wood et al. (2007). For this equation, r is the solution to where l x is the survival probability to age x and b x is the number of offspring born to an individual of age x. In their population viability analysis of tree squirrels, Wood et al. (2007) give values for the percentage of females breeding per season, P, the number of breedings per year, β, the litter size per female, L, and the mortality rates for age zero, m 0 , and individuals one year and older, m (Table A1) Taking logarithms, we get.
A root-finding algorithm was used to numerically solve equation A2 for r. Wood et al. (2007) provided three sets of the life history values described for the following species: Sciurus aberti, S. carolinensis, S. granatensis, S. niger (FS), and S. vulgaris. The three sets of parameters are labeled as "optimistic," "average," and "pessimistic." Using the values labeled as average, we computed the values of r, respectively, for the five species (Table A1). The average of these five values is 0.32. Therefore, we used a value of r 1 = r 2 = .3 in our simulations. This is close to the value of 0.29 obtained for FS (Table A1).

A PPE N D I X B N U M B E R S O F ACO R N S
Estimates of acorn production were made in years 2012 through 2020 by measuring during the middle of October the mass of acorns in a one square meter plot under each of six oak trees within the study area. The same trees and the same plots were used each year for estimation of acorn production to assess variability among years and relationship to changes in abundance of the tree squirrels. The data appear in Table A2.  Wood et al. (2007,

A PPE N D I X C N U M B E R S O F J U V E N I LE A N D S U BA D U LT S Q U I R R E L S
Starting in January 2015, observers at our field site started distinguishing among adults, subadults, and juveniles. The latter two categories can serve as a proxy for reproduction. Figure A1 shows the numbers of juveniles and subadults for the WGS and the FS. For both species, juveniles and subadults are seen more often in spring and fall ( Figure A1, shaded areas), especially in the case of the FS. However, there is a great deal of year-to-year variation (long timescales) in the numbers.

A PPE N D I X D A N N UA L A N D S E M I -A N N UA L TE M PE R ATU R E C YCLE S
Meteorologists and climatologists have used harmonic analysis to identify seasonal cycles in atmospheric temperature (White & Wallace, 1978). In addition to a strong annual cycle, a semi-annual cycle, which can vary by year and location, has been identified (North et al., 2021;White & Wallace, 1978). In their equation (1) and ϕ i is the phase shift (in months) of the annual (i = 1) and semi-annual (i = 2) component cycles.
We fit equation (A3) to the mean monthly temperature data from Ontario Airport (ONT) (Figure 4a) using the method of nonlinear least squares. The parameter estimates and their standard errors appear in Table A3. The phase shifts are relative to the month of September. The temperature data and the fitted function appear in panel A of Figure A2. The fitted function provides a good description of the observed temperature time series.
The two component cycles are shown separately in panel B of Figure A2. The annual cycle is the larger of the two and is obtained by setting A 2 = 0 in equation (A3). It peaks in July-August and has its trough in January-February. The smaller semi-annual cycle, obtained by setting A 1 = 0 in equation (A3), peaks twice per year in February-March and August-September and has its troughs in November-December and May-June. The amplitude of the semi-annual cycle is 21% the size of the amplitude of the annual cycle.
In panel (c) of Figure A2, we plotted the temperature time series obtained after filtering out the annual cycle (Figure 4a, dashed line) with the semi-annual cycle from panel (b) of Figure A2. The theoretical and filtered cycles have the same periodicity, phase, and approximate amplitude.
Moreover, one can see the effects of long timescale variability in the way the filtered data meanders above and below the semi-annual cycle.
This gives us confidence that the band-stop filter used to attenuate the annual cycle worked well.
We fit equation (A3) to the observed temperature time series with A 2 = 0, so that only the annual cycle was present. The parameter estimates and their standard errors are also presented in Table A3. We computed the following Akaike Information Criterion (AIC) for both fitted models: (A3) x t = a 0 + A 1 cos 2 t + 1 12 + A 2 cos 4 t + 2 12 , AIC = nln SS E n + 2k, TA B L E A 2 Mass of acorns (g) in a 1 m 2 plot under each tree in the study area We computed 95% significance threshold by generating 2000 random permutations of the residuals, obtaining a smoothed normalized spectrum for each, and using the 95th percentiles of these surrogate spectra for the 95% limits. Figure A3 shows the ONT temperature data and spectrum. The fitted quadratic polynomial (dashed line in Figure A3, panel (a)) is nearly linear and suggests a trend of increasing temperatures.
Although the significance band is wide due to the short length of the time series, the spectrum for the residuals shows a significant peak at a timescale of about 7 years ( Figure A3, panel (b)). This corresponds to the timescale range of the local peak in mean wavelet power in the lower right corner of Figure 9.