Ecology

ModelTo investigate the evolution of cyclically competing species with intraspecific interaction which sensitively plays to the territory awareness, we employ the spatial RPS model11,19,20,23. At the microscopic level, the model can be demonstrated on a lattice system, and for convenience, we consider a square lattice of size N with periodic boundary conditions where all sites have von Neumann neighbors. Each site can be occupied by an individual from one of the three species (referred to as A, B, and C, respectively) or left empty(E), and thus the system describes a limited carrying capacity. In addition, to explore the effect of territory awareness on intraspecific interaction, we assume that the given lattice is divided into two areas of the same size which may possibly realize different territorial ranges. Here we simply divide the two regions into the top and bottom halves of the given square lattice. To reflect the territorial awareness on intraspecific interaction, we distribute population in each group into two sub-networks randomly, and denote species (X_1) for the top and (X_2) for the bottom ((X in {A, B, C})) to distinguish the emergence of intraspecific interaction between individuals who lie on different domains. The distribution of all species with respect to the separation of the domain is illustrated in Fig. 1a.Figure 1Schematic diagrams of network structure and the invasion rules among species. (a) Each circle represents a node, and individuals of species A, B, and C are evenly and randomly distributed on each node. To realize territorial awareness, the lattice is divided into two regions of equal size: the top and the bottom where the dashed line indicates the regional boundary. Two genera of the same species are distributed in different regions, and different color markers represent different species types. Nodes without color markers are empty nodes. (b) Interspecific interaction among three species A, B, and C (indicated by three boxes) occurs cyclically with a rate (p_1). A box of each group describes the intraspecific interaction between individuals who belong to different territories where the interaction is regulated by territorial consciousness. Here intraspecific interaction in each group occurs with a rate (k_i cdot p_2) ((i in {A, B, C})).Full size imageUnder the given assumption for the lattice, all the interactions between individuals occur within nearest neighboring sites by the following set of rules (see Fig. 1b):$$begin{aligned}&A_{i}+B_{j}overset{p_{1}}{longrightarrow }A_{i}+E,quad B_{i}+C_{j}overset{p_{1}}{longrightarrow }B_{i}+E,quad C_{i}+A_{j}overset{p_{1}}{longrightarrow }C_{i}+E, end{aligned}$$
(1)
$$begin{aligned}&A_{i}+A_{j}overset{k_{A} cdot p_{2}}{longrightarrow }A_{i}+E,mathrm{or},E+A_{j}, quad B_{i}+B_{j}overset{k_{B} cdot p_{2}}{longrightarrow }B_{i}+E,mathrm{or},E+B_{j}, quad C_{i}+C_{j}overset{k_{C} cdot p_{2}}{longrightarrow }C_{i}+E,mathrm{or},E+C_{j},quad ine j, end{aligned}$$
(2)
$$begin{aligned}&A_{i}+Eoverset{r}{longrightarrow }A_{i}+A_{i},quad B_{i}+Eoverset{r}{longrightarrow }B_{i}+B_{i},quad C_{i}+Eoverset{r}{longrightarrow }C_{i}+C_{i}, end{aligned}$$
(3)
$$begin{aligned}&A_{i}+otimes overset{m}{longrightarrow }otimes +A_{i},quad B_{i}+otimes overset{m}{longrightarrow }otimes +B_{i},quad C_{i}+otimes overset{m}{longrightarrow }otimes +C_{i}, end{aligned}$$
(4)
where (i, j=1, 2). The mark (otimes) stands for any species or empty sites. Relation (1) describes interspecific interaction among three species which occurs cyclically with a rate (p_1): (A_{i}) dominates (B_{i}), (B_{i}) dominates (C_{i}), and (C_{i}) dominates (A_{i}) ((i=1, 2)). The defeated individual dies and the site becomes an empty site. Relation (2) demonstrates the intraspecific interaction which will sensitively depend on territorial awareness. Since we assume the intraspecific interaction is related to the territorial consciousness, the rate in each species may be defined by (k_{A} cdot p_{2}), (k_{B} cdot p_{2}), (k_{C} cdot p_{2}) for species A, B, C, respectively, where (p_{2}) is the given rate of interaction, k is the the sensitive parameter to territorial awareness. Similar to previous works, the result of intraspecific interaction eventually results in a death of one individual at random with a 1/2 chance. Relation (3) stands for the reproduction with a rate r which is allowed when an empty site in neighbors is selected, and migration defined by an exchange between two neighboring sites is denoted in Relation (4). Based on the theory of random walks38, it occurs with a rate (m=2MN) where M and N indicate individuals’ mobility and a system size, respectively, as usual to previous works. Thus, an actual time step is defined when each individuals has interacted with others once on average, i.e., N pairwise interactions will occur in one actual time step unit. In order to make an unbiased comparison with previous works15,19,20,21 and for the convenience of interpretations, we assume parameters as (p_{1}=p_{2}=r=1) and (k_{A}=k_{B}=k_{C}=k) (see the Methods for the meanings of specific parameters) in our simulations. Three species are divided into two types to distinguish distributions on different regions: (A_{1,2}), (B_{1,2}), and (C_{1,2}), and randomly distributed initially on a square lattice of size (N= 300 times 300). In addition, in all our simulations, species coexistence refers to the coexistence of (A_i), (B_j), and (C_k) for any combination of (i,,j,,k in left{ 1,,2 right}).Biodiversity under territorial awarenessWe first consider the effect of territorial awareness on species biodiversity. In general, it is well-known that, the spatial RPS game exhibits a transition of survival states from coexistence to extinction (which is presented by the uniform state) as individuals’ mobility increases. The phase transition occurs when M exceed a certain value, referred to as a critical mobility (M_c = (4.5 pm 0.5) times 10^{-4}), which is identified in Ref.11. To address the effect of territorial awareness, we consider two different mobility values (M=1 times 10^{-5}) and (M=1 times 10^{-3}) which eventually yield different survival states: coexistence and extinction, respectively, for different sensitivity parameter k.In general, the total simulation time T in classic spatial RPS games is considered as (T = N) which can yield the extinction for the critical mobility (M_c)11. In this regard, using the time (T=N) may yield different results for species evolution and corresponding survival states due to stochastic events, and such behaviors may be induced by the choice of mobility. In our simulations, since we consider two different mobility values where the one (Fig. 2a–c) is quite lower and the other (Fig. 2d–f) is higher than (M_c), we thus consider different simulation times at (M=1 times 10^{-5}) and (M=1 times 10^{-3}): more than 490, 000 and 180, 000 steps, respectively, to obtain robust features on species survival states. The time dependent evolution of densities are illustrated in Fig. 2 where the top and bottom panels are obtained from simulations with the first 250, 000 and 140, 000 steps, respectively.Figure 2Time dependent evolution of densities in the system for different M and k. Top and bottom panels are obtained with (M=1times 10^{-5}) and (M=1times 10^{-3}), respectively, and the sensitivity parameter k is given by (k=5), 10, and 20 from the left to right in each row. (a–c) Regardless of the choice k, the low mobility still leads species coexistence as usual. (d–f) At high mobility regimes, the system also always exhibit the extinction state.Full size imageEven if different k are considered, the panels in Fig. 2 show features similar to previous works11,14,15,16,18,20: coexistence and extinction for tops and bottoms, respectively. At the low mobility (M=1times 10^{-5}) as shown in Fig. 2a–c, even if the individuals located in different domains in each group disappear, the spatial RPS game eventually exhibits coexistence as k increases since some of individuals in A, B, C are survived. For instance, in our simulations, coexistence can be presented by survival of species (A_1), (B_1), (C_1) (Fig. 2a–b) or (A_2), (B_2), (C_1) (Fig. 2c). Since the typical waiting time for extinction is exponentially increasing to the size N at low mobility11, there will be extinction and eventually only one species will dominate the system after extremely long times. Thus, within the finite time steps, one type of each species will disappear slowly with the increase of k and the system exhibits coexistence.On the other hand, the high mobility (M=1times 10^{-3}) leads the extinction and only one species dominate the whole domain. As shown in Fig. 2d–f, the extinction that is defined by the two types of the species disappear occurs and the only one species finally dominate the system [e.g., (C_2), (B_2), and (C_1) for (k=5), 10, and 20, respectively]. In this case, the increase in k has little effect on the disappearance of one of the species, but has a tendency to accelerate the complete extinction of the second species. Take Fig. 2d for example, when species (A_2), (B_1) and (B_2) became extinct, (A_1), (C_1) and (C_2) is left in the system, the density of (A_1) in the system had an absolute advantage, while (C_1) and (C_2) had intraspecific interaction. Since the intensity of intraspecific interaction sensitive to territorial awareness was greater than that of interspecific interaction, the interaction was mainly intraspecific interaction between (C_1) and (C_2), then (C_1) was defeated into extinction, (C_2) preyed on the only specie (A_2) remaining in the system and eventually occupied the whole system. As k value affects the intraspecific interaction intensity, it determines the waiting time for the extinction of two species in the system. For example, we found that the larger the k value is, the shorter the waiting time for the extinction of two species is, as shown in Fig. 2e–f. But this is only the observation result of a single simulation. Due to the randomness of the simulation, this phenomenon needs further verification, so we give specific results about the effect of k value on the average extinction time in the next section. Due to stochastic events during Monte-Carlo simulations, the combination of survival species for coexistence and extinction at the final step can be different, but the such states at two mobility regimes will be still maintained. Fig. 2 may impose the follows: territorial awareness on intraspecific interaction can eventually yield similar feature to previous works in a broad aspect, but the composition of the surviving species type for each state may vary.Average extinction time versus territorial awarenessWhile survival states in both cases are consistent with previous works on the effects of species migration in Fig. 2, we found an interesting feature that the evolutionary time when some type of species disappear is changed depending on k. To be concrete, at (M=1 times 10^{-5}), we found that one type of each species (A_1), (B_1), (C_1) (Fig. 2a) will eventually coexist while their companion species (A_2), (B_2), (C_2) are extinct as t exceeds (t approx 50{,}000). As k is increasing, the time point when one genus of each species disappears shows an increasing pattern as presented in Fig. 2b,c. The opposite trend can be captured at high mobility (M=1 times 10^{-3}), that is, the increase of k seems to shorten the evolution time of two species extinction in the system. Based on these observations, we may assume that the critical time for such disappearance phenomena has a certain relationship with k and the relation may differ to the choice of M.To answer the issue, we measure the average extinction time T. In classic RPS games, traditionally, the extinction state on spatially extended systems has been identified by the uniform state that only one species dominates whole domain11,14,15,16,18,20. As shown in Fig. 2a–c which ultimately describe a coexistence state in a finite time, however, any one of type in each species disappeared and the time associated with the phenomena is changed by the strength of k. In a slightly different aspect to the classic meaning of extinction, we here define the average extinction time T with respect to the regime of mobility: (a) the evolutionary time when one genus of each species disappears for low mobility and (b) the time when two of the three species disappear completely for high mobility. In this consideration, for both given cases of M in Fig. 2, the average extinction time T in each k is measured from 30 independent realizations and presented in Fig. 3.Figure 3The average extinction time T as a function of the territorial sensitive parameter k. (a) Two cases of fixed mobility in territorial sensitive intraspecific interaction. For low mobility (M=1times 10^{-5}), the time T which is measured by detecting the time when one genus of each species disappears tends to increase with the increase of k, i.e., the high sensitivity of territorial awareness has the effect of delaying the waiting time for extinction. Similarly, it can be seen that at high mobility (M=1times 10^{-3}), an increase in k value will also delay the waiting time for extinction, but the effect is much more gentle. (b) At low mobility value (M=1times 10^{-5}), traditional intraspecific interaction (i.e., intraspecific interaction among all individuals of the same species, regardless of territorial residence) was compared with territorial sensitive intraspecific interaction. Here, for the traditional case, k represents intraspecific interaction intensity, and the running time of the simulation is 810, 000 steps. In the case that the final steady state has not occurred before the end of the simulation, we take the maximum time step ((t=810{,}000)) as the extinction time T value, which causes the blue line to become gentle when (k >14). Compared with the traditional situation, the intraspecific interaction affected by territorial awareness significantly reduced the average extinction time, that is, accelerated the damage of species diversity in the system. The results were averaged from 30 independent simulations, and error bars (using standard errors, which defined as the sample standard deviation divided by the square root of the number of samples) are shown in the figure.Full size imageAs shown in Fig. 3a, we find clearly that the average extinction time is obviously affected by the strength of sensitivity coefficient k, especially, when the mobility is low. When species has no consciousness on territories ((k=0)), the system becomes exactly the classic RPS model11 since intraspecific interaction is undefined, and the waiting time T generally tends to increase exponentially to the choice of M. However, our simulation shows the T is approximately measured at (T=110{,}000) at (k=0). Traditionally, it is well known that the average waiting time for extinction in the classic RPS game is taken (T=N) near the critical mobility regime ((M approx M_c)), and the coexistence duration is exponentially increasing as M decreases from (M_c). Within this knowledge, our simulation results may seem inconsistent with the general concept of extinction time. In our model, however, the definition of extinction is different at the low mobility regime, and the change into a single RPS system as one genus of each individual disappears may have a similar meaning to the previous definition of extinction in some sense, the above result can be said to be reasonable.The important point is actually addressed for (k >0). In this case, species can allow intraspecific interaction and the strength of intraspecific interaction is also increasing since the territorial awareness is intensified. As a result, it is found that the average extinction time T shows a tendency to gradually increase with the increase of k at (M=1 times 10^{-5}). In addition, this trend can also be observed at (M = 1 times 10^{-3}), but it is more gradual. To investigate whether the tendency to prolong the waiting time for extinction time at low migration rates is caused by territorial awareness or the existence of intraspecific interaction, we compared traditional intraspecific interaction (i.e., intraspecific interaction among all individuals of the same species, regardless of territorial residence, which equivalent to removing the condition (ine j) from Relation (2)) with territorial-sensitive intraspecific interaction in our model, the results are shown in Fig. 3b. We found that in the presence of intraspecific interaction, the average extinction time increased with the intensity of intraspecific interaction. Specifically, the stronger the intraspecific interaction, the slower the loss of species diversity. However, compared with the traditional situation, intraspecific interaction influenced by territorial consciousness controlled the delay of extinction to a certain extent. Even if our simulations have been carried on for two specific M, it is obvious that the territorial awareness can affect the average extinction time, and we suggest that a strong sense of territoriality can also delay species extinction and lead to long-term coexistence of systems at low mobility regimes, although the introduction of territoriality leads to faster damage to species diversity than is traditionally the case, while it does not affect significantly on the extinction time and the biodiversity (which eventually appears as extinction) at high mobility regimes.Evolution of the interface between territoriesFrom the investigation on the average extinction time in Fig. 3, we know that the territorial awareness can affect not only species survival but also the maintenance period of survival states. Here, we may wonder why the territorial awareness can affect the waiting time to extinction. In order to investigate such an issue, we observe evolution of the spatial system, in particular invasion between species near the border on two territories, i.e., the evolution of the interface. To capture the phenomena in detail, we consider pattern formations associated with the given two mobility values at the initial state of the evolution (e.g. (t=1000)) which are represented in Fig. 4.Figure 4The typical snapshots of evolution on patterns at (t=1000) for different k: 0.5 for (a) and (e), 2 for (b) and (f), 10 for (c) and (g), and 20 for (d) and (h), where the mobility is considered as (M = 1times 10^{-5}) for tops and (M=1 times 10^{-3}) for bottoms. Different colors correspond to different species types, as shown in Figs. 1 and 2, with white indicating vacancy. As k increases, the invasion among species between two territories occurs more gradually, and such phenomena are clearly observed for the high mobility as shown in the panel (h).Full size imageThe top and bottom panels in Fig. 4 exhibit spatial patterns for the low and high mobility values, respectively. For (M=1 times 10^{-5}), when the value k is small such as (k=0.5) (see Fig. 4a), interspecific interaction can occur more frequently than intraspecific interaction among all pairwise reactions (1)–(4). The system can exhibit similar pattern formations to the classic RPS game11. Three species, even if they are distinguished into six subgroups, are spirally entangled with clearly exhibiting spiral waves which are appeared in both two territories. Since the given lattice has periodic boundaries, species in both territories can migrate to the other region each other, but such migration is weak because the normalized probability for migration (Relation (4)) is small at the low mobility. Thus, when the system exhibits coexistence, it may be possible to predict that the top and bottom territories present dominance of species (X_1) and (X_2) ((X in {A, B, C})), respectively, while our simulations only present spatial patterns at the first 1,000 steps which may be too short to lead phase transitions.We also find that the spiral-wave patterns are getting to fuzzy as k increases. In particular, such fuzzy patterns are conspicuous near the boundary between the two territories at the large k (see Fig. 4c,d). The increase of k directly means the intensification of intraspecific interaction, and according to the setting on the initial distribution of population, intraspecific interaction will have many chances to occur in the vicinity of the boundary than near the top and bottom periodic boundaries. Frequent intraspecific interaction can provide as many chances to allow reproduction as possible, and high intraspecific interaction rate can dominate on pairwise invasions than interspecific interaction.In the vicinity of the border between two territories, the occurrence of intraspecific interaction is observed more prominently at (M=1 times 10^{-3}), and such features are clear as k increases. To be concrete, compared with figures among Fig. 4e–h, we found that empty sites are produced near the border and their presence is clear for high strength k such as 10 and 20 (Fig. 4g–h). In this case, the two domains appear to be more clearly divided and each domain is dominated by a single RPS system. Each single RPS system shows extinction state (only one genus survives) at high mobility, and eventually shows extinction state through interspecific or intraspecific interaction depending on the type of surviving genus. This is in good agreement with the results we obtained in Fig. 2. However, it can be seen from Fig. 3 that the time for each domain system to reach extinction at high mobility is very short compared to that at low mobility, and this has no relation with the degree of territorial awareness in interspecific interaction.From our simulations, we find that: the relationship between territorial awareness and the average extinction time is particularly prominent at the low mobility, and the likelihood of intraspecific interaction is relatively high near territorial boundaries. Under these considerations, we may expect a new relationship between the delay of the extinction time and boundary of two territories. To uncover this veil, we try to quantify a width for occurrence of intraspecific interaction near the border between two area with respect to the territorial awareness. Specifically, we give each node on a two-dimensional grid a coordinate, defined by its row and column position. For each column (j=1,ldots ,L), calculate the interface width, defined as I:$$begin{aligned} I_{j}= {left{ begin{array}{ll} P(1,j)-P(2,j),&{}quad text {if } P(1,j) >P(2,j)\ 0,&{}quad text {if } P(1,j) More

1.Gibson, L. et al. Primary forests are irreplaceable for sustaining tropical biodiversity. Nature 478(7369), 378–381. https://doi.org/10.1038/nature10425 (2011).CAS 
Article 
ADS 

Google Scholar 
2.Newbold, T. et al. A global model of the response of tropical and sub-tropical forest biodiversity to anthropogenic pressures. Proc. R. Soc. B. 281, 20141371. https://doi.org/10.1098/rspb.2014.1371 (2014).Article 
PubMed 
PubMed Central 

Google Scholar 
3.Hansen, M. C. et al. The fate of tropical forest fragments. Sci. Adv. 6(11), eaax8574. https://doi.org/10.1126/sciadv.aax8574 (2020).Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
4.Peres, C. A. et al. Biodiversity conservation in human-modified Amazonian Forest landscapes. Biol. Conserv. 143, 2314–2327. https://doi.org/10.1016/j.biocon.2010.01.021 (2010).Article 

Google Scholar 
5.PRODES INPE. Monitoring Deforestation of the Brazilian Amazon Forest by Satellite. TerraBrasilis (inpe.br) (accessed in october 2020, 2020).6.Barlow, J. et al. Quantifying the biodiversity value of tropical primary, secondary, and plantation forests. Proc. Natl. Acad. Sci. 104, 18555–18560. https://doi.org/10.1073/pnas.0703333104 (2007).Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
7.Peres, C. A., Barlow, J. & Laurance, W. F. Detecting anthropogenic disturbance in tropical forests. Trends Ecol. Evol. 21, 227–229. https://doi.org/10.1016/j.tree.2006.03.007 (2006).Article 
PubMed 

Google Scholar 
8.Arroyo-Rodríguez, V. et al. Designing optimal human-modified landscapes for forest biodiversity conservation. Ecol. Lett. 23, 1404–1420. https://doi.org/10.1111/ele.13535 (2020).Article 
PubMed 

Google Scholar 
9.Gardner, T. A. et al. Prospects for tropical forest biodiversity in a human-modified world. Ecol. Lett. 12, 1–21. https://doi.org/10.1111/j.1461-0248.2009.01294.x (2009).Article 

Google Scholar 
10.Hardwick, S. R. et al. The relationship between leaf area index and microclimate in tropical forest and oil palm plantation: Forest disturbance drives changes in microclimate. Agric. For. Meteorol. 201, 187–195. https://doi.org/10.1016/j.agrformet.2014.11.010 (2015).Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
11.Sambuichi, R. H. et al. Cabruca agroforests in southern Bahia, Brazil: Tree component, management practices and tree species conservation. Biodivers. Conserv. 21, 1055–1077. https://doi.org/10.1007/s10531-012-0240-3 (2012).Article 

Google Scholar 
12.Devictor, V., Julliard, R. & Jiguet, F. Distribution of specialist and generalist species along spatial gradients of habitat disturbance and fragmentation. Oikos 117, 507–514. https://doi.org/10.1111/j.0030-1299.2008.16215.x (2008).Article 

Google Scholar 
13.Banks-Leite, C. Using ecological thresholds to evaluate the costs and benefits of set-asides in a biodiversity hotspot. Science 345, 1041–1045. https://doi.org/10.1126/science.1255768 (2014).CAS 
Article 
PubMed 
ADS 

Google Scholar 
14.Newbold, T. et al. Global patterns of terrestrial assemblage turnover within and among land uses. Ecography 39, 1151–1163. https://doi.org/10.1111/ecog.01932 (2016).Article 

Google Scholar 
15.Paglia, A. P. et al. Annotated checklist of Brazilian mammals. Occas. Pap. Conserv. Int. 6, 1–82 (2012).
Google Scholar 
16.Dirzo, R. et al. Defaunation in the anthropocene. Science 345, 401–406. https://doi.org/10.1126/science.1251817 (2014).CAS 
Article 
PubMed 
ADS 

Google Scholar 
17.Estrada, A. et al. Impending extinction crisis of the world’s primates: Why primates matter. Sci. Adv. 3, e1600946. https://doi.org/10.1126/sciadv.1600946 (2017).Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
18.Newbold, T. et al. Global effects of land use on local terrestrial biodiversity. Nature 520, 45–50. https://doi.org/10.1038/nature14324 (2015).CAS 
Article 
PubMed 
ADS 

Google Scholar 
19.Phillips, H. R., Newbold, T. & Purvis, A. Land-use effects on local biodiversity in tropical forests vary between continents. Biodivers. Conserv. 26, 2251–2270. https://doi.org/10.1007/s10531-017-1356-2 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
20.Teixeira, D. F., Guillera-Arroita, G., Hilário, R. R., Fonseca, C. & Rosalino, L. M. Influence of life-history traits on the occurrence of carnivores within exotic Eucalyptus plantations. Divers. Distrib. 26, 1071–1082. https://doi.org/10.1111/ddi.13114 (2020).Article 

Google Scholar 
21.Asner, G. P. et al. Selective logging in the Brazilian Amazon. Science 310, 480–482. https://doi.org/10.1126/science.1118051 (2005).CAS 
Article 
PubMed 
ADS 

Google Scholar 
22.Robinson, J. G. & Redford, K. H. Body size, diet, and population density of neotropical forest mammals. Am. Nat. 128, 665–680. https://doi.org/10.1086/284596 (1986).Article 

Google Scholar 
23.Cardillo, M. et al. Multiple causes of high extinction risk in large mammal species. Science 309, 1239–1241. https://doi.org/10.1890/05-0112 (2005).CAS 
Article 
PubMed 
ADS 

Google Scholar 
24.Almeida-Maués, P.C.R. Efeitos antropogênicos sobre a diversidade de mamíferos de médio e grande porte na Amazônia Oriental. PhD. Thesis, Graduate Program in Ecology, Federal University of Pará, Belém, Pará, Brazil (2019).25.Parry, L., Barlow, J. & Peres, C. A. Large-vertebrate assemblages of primary and secondary forests in the Brazilian Amazon. J. Trop. Ecol. 23, 653–662. https://doi.org/10.1017/S0266467407004506 (2007).Article 

Google Scholar 
26.Mendes-Oliveira, A. C. et al. Oil palm monoculture induces drastic erosion of an Amazonian forest mammal fauna. PLoS ONE 12, e0187650. https://doi.org/10.1371/journal.pone.0187650 (2017).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
27.Coelho, M., Juen, L. & Mendes-Oliveira, A. C. The role of remnants of Amazon savanna for the conservation of Neotropical mammal communities in eucalyptus plantations. Biodivers. Conserv. 23, 3171–3184. https://doi.org/10.1007/s10531-014-0772-9 (2014).Article 

Google Scholar 
28.Bicknell, J. E., Struebig, M. J. & Davies, Z. G. Reconciling timber extraction with biodiversity conservation in tropical forests using reduced-impact logging. J. Appl. Ecol. 52, 379–388. https://doi.org/10.1111/1365-2664.12391 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
29.Chazdon, R. L. et al. The potential for species conservation in tropical secondary forests. Conserv. Biol. 23, 1406–1417. https://doi.org/10.1111/j.1523-1739.2009.01338.x (2009).Article 
PubMed 

Google Scholar 
30.Koh, L. P. & Wilcove, D. S. Is oil palm agriculture really destroying tropical biodiversity?. Conserv. Lett. 1, 60–64. https://doi.org/10.1111/j.1755-263X.2008.00011.x (2008).Article 

Google Scholar 
31.Putz, F. E. & Pinard, M. A. Reduced-impact logging as a carbon-offset method. Conserv. Biol. 7, 755–757. https://doi.org/10.1046/j.1523-1739.1993.7407551.x (1993).Article 

Google Scholar 
32.Pinard, M. A. & Putz, F. E. Retaining forest biomass by reducing logging damage. Biotropica 28, 278–295. https://doi.org/10.2307/2389193 (1996).Article 

Google Scholar 
33.Prudente, B. S., Pompeu, P. S., Juen, L. & Montag, L. F. A. Effects of reduced-impact logging on physical habitat and fish assemblages in streams of Eastern Amazonia. Freshw. Biol. 62, 303–316. https://doi.org/10.1111/fwb.12868 (2017).Article 

Google Scholar 
34.Kanowski, J., Catterall, C. P. & Wardell-Johnson, G. W. Consequences of broadscale timber plantations for biodiversity in cleared rainforest landscapes of tropical and subtropical Australia. For. Ecol. Manage. 208, 359–372. https://doi.org/10.1016/j.foreco.2005.01.018 (2005).Article 

Google Scholar 
35.Correa, F. S., Juen, L., Rodrigues, L. C., Silva-Filho, H. F. & Santos-Costa, M. C. Effects of oil palm plantations on anuran diversity in the eastern Amazon. Anim. Biol. 65, 321–335. https://doi.org/10.1163/15707563-00002481 (2015).Article 

Google Scholar 
36.Peres, C. A. & Cunha, A. A. Line-Transect Censuses of Large-Bodied Tropical Forest Vertebrates: A Handbook (Wildlife Conservation Society, 2011).
Google Scholar 
37.Chao, A. & Jost, L. Coverage-based rarefaction and extrapolation: Standardizing samples by completeness rather than size. Ecology 93, 2533–2547. https://doi.org/10.1890/11-1952.1 (2012).Article 
PubMed 

Google Scholar 
38.Oksanen, J. F. et al. vegan: Community Ecology Package. R package version 2.5–6. https://CRAN.R-project.org/package=vegan (2019).39.Ceballos, G., Ehrlich, P. R. & Dirzo, R. Biological annihilation via the ongoing sixth mass extinction signaled by vertebrate population losses and declines. Proc. Natl. Acad. Sci. 114, 6089–6096. https://doi.org/10.1073/pnas.1704949114 (2017).CAS 
Article 

Google Scholar 
40.Kricher, J. Tropical Ecology 632 (Princeton University Press, 2011).
Google Scholar 
41.Edwards, D. P. et al. Reduced-impact logging and biodiversity conservation: A case study from Borneo. Ecol. Appl. 22, 561–571. https://doi.org/10.1890/11-1362.1 (2012).Article 
PubMed 

Google Scholar 
42.Melo, F. P. L., Arroyo-Rodríguez, V., Fahrig, L., Martínez-Ramos, M. & Tabarelli, M. On the hope for biodiversity friendly tropical landscapes. Trends Ecol. Evol. 28, 462–468. https://doi.org/10.1016/j.tree.2013.01.001 (2013).Article 
PubMed 

Google Scholar 
43.Benton, T. G., Vickery, J. A. & Wilson, J. D. Farmland biodiversity: Is habitat heterogeneity the key?. Trends Ecol. Evol. 18, 182–188. https://doi.org/10.1016/S0169-5347(03)00011-9 (2003).Article 

Google Scholar 
44.Almeida-Rocha, J. M., Peres, C. A. & Oliveira, L. C. Primate responses to anthropogenic habitat disturbance: A pantropical meta-analysis. Biol. Conserv. 215, 30–38. https://doi.org/10.1016/j.biocon.2017.08.018 (2017).Article 

Google Scholar 
45.Palmeirim, A. F., Vieira, M. V. & Peres, C. A. Herpetofaunal responses to anthropogenic forest habitat modification across the neotropics: Insights from partitioning β-diversity. Biodivers. Conserv. 26, 2877–2891. https://doi.org/10.1007/s10531-017-1394-9 (2017).Article 

Google Scholar 
46.Christie, A. P. et al. Quantifying and addressing the prevalence and bias of study designs in the environmental and social sciences. Nat. Commun. 11, 6377. https://doi.org/10.1038/s41467-020-20142-y (2020).CAS 
Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
47.Whitworth, A. et al. Human disturbance impacts on rainforest mammals are most notable in the canopy, especially for larger-bodied species. Divers. Distrib. 25, 1166–1178. https://doi.org/10.1111/ddi.12930 (2019).Article 

Google Scholar 
48.Johns, A. D. & Skorupa, J. P. Responses of rain-forest primates to habitat disturbance: A review. Int. J. Primatol. 8, 157–191. https://doi.org/10.1007/BF02735162 (1987).Article 

Google Scholar 
49.Wearn, O. R. et al. Mammalian species abundance across a gradient of tropical land-use intensity: A hierarchical multi-species modelling approach. Biol. Conserv. 212, 162–171. https://doi.org/10.1016/j.biocon.2017.05.007 (2017).Article 

Google Scholar 
50.Benchimol, M. & Peres, C. A. Determinants of population persistence and abundance of terrestrial and arboreal vertebrates stranded in tropical forest land-bridge islands. Conserv. Biol. 35(3), 870–883. https://doi.org/10.1111/cobi.13619 (2020).Article 
PubMed 

Google Scholar 
51.Gittleman, J. L. & Harvey, P. H. Carnivore home-range size, metabolic needs and Ecology. Behav. Ecol. Sociobiol. 10(1), 57–63. https://doi.org/10.1007/BF00296396 (1982).Article 

Google Scholar 
52.Edwards, D. P., Tobias, J. A., Sheil, D., Meijaard, E. & Laurance, W. F. Maintaining ecosystem function and services in logged tropical forests. Trends Ecol. Evol. 29, 511–520. https://doi.org/10.1016/j.tree.2014.07.003 (2014).Article 
PubMed 

Google Scholar 
53.Mollinari, M. M., Peres, C. A. & Edwards, D. P. Rapid recovery of thermal environment after selective logging in the Amazon. Agric. For. Meteorol. 278, 107637. https://doi.org/10.1016/j.agrformet.2019.107637 (2019).Article 
ADS 

Google Scholar 
54.Azevedo-Ramos, C., de Carvalho, O. & de Amaral, B. D. Short-term effects of reduced-impact logging on eastern Amazon fauna. For. Ecol. Manag. 232, 26–35. https://doi.org/10.1016/j.foreco.2006.05.025 (2006).Article 

Google Scholar 
55.Bicknell, J. E. & Peres, C. A. Vertebrate population responses to reduced-impact logging in a neotropical forest. For. Ecol. Manage. 259, 2267–2275. https://doi.org/10.1016/j.foreco.2010.02.027 (2010).Article 

Google Scholar 
56.Laufer, J., Michalski, F. & Peres, C. A. Effects of reduced-impact logging on medium and large-bodied forest vertebrates in eastern Amazonia. Biota Neotrop. 15, e20140131. https://doi.org/10.1590/1676-06032015013114 (2015).Article 

Google Scholar 
57.Carvalho Jr, E. A. R., Mendonça, E. N., Martins, A. & Haugaasen, T. Effects of illegal logging on Amazonian medium and large-sized terrestrial vertebrates. For. Ecol. Manage. 466, 118105. https://doi.org/10.1016/j.foreco.2020.118105 (2020).Article 

Google Scholar 
58.Kuussaari, M. et al. Extinction debt: A challenge for biodiversity conservation. Trends Ecol. Evol. 24, 564–571. https://doi.org/10.1016/j.tree.2009.04.011 (2009).Article 
PubMed 

Google Scholar 
59.Richardson, V. A. & Peres, C. A. Temporal decay in timber species composition and value in Amazonian logging concessions. PLoS ONE 11, e0159035. https://doi.org/10.1371/journal.pone.0159035 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
60.Chazdon, R. L. Second Growth: The Promise of Tropical Forest Regeneration in an Age of Deforestation (University of Chicago Press, 2014).Book 

Google Scholar 
61.Acevedo-Charry, O. & Aide, T. M. Recovery of amphibian, reptile, bird and mammal diversity during secondary forest succession in the tropics. Oikos 128, 1065–1078. https://doi.org/10.1111/oik.06252 (2019).Article 

Google Scholar 
62.Sodhi, N. S. et al. Conserving Southeast Asian forest biodiversity in human-modified landscapes. Biol. Conserv. 143, 2375–2384. https://doi.org/10.1016/j.biocon.2009.12.029 (2010).Article 

Google Scholar 
63.Dunn, R. R. Recovery of faunal communities during tropical forest regeneration. Conserv. Biol. 18, 302–309. https://doi.org/10.1111/J.1523-1739.2004.00151.X (2004).Article 

Google Scholar 
64.Luskin, M. S. & Potts, M. D. Microclimate and habitat heterogeneity through the oil palm lifecycle. Basic Appl. Ecol. 12, 540–551. https://doi.org/10.1016/j.baae.2011.06.004 (2011).Article 

Google Scholar 
65.Fitzherbert, E. B. et al. How will oil palm expansion affect biodiversity?. Trends Ecol. Evol. 23(10), 538–545. https://doi.org/10.1016/j.tree.2008.06.012 (2008).Article 
PubMed 

Google Scholar 
66.Martello, F. et al. Homogenization and impoverishment of taxonomic and functional diversity of ants in Eucalyptus plantations. Sci. Rep. 8, 3266. https://doi.org/10.1038/s41598-018-20823-1 (2018).CAS 
Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
67.da Rocha, P. L. B. What is the value of eucalyptus monocultures for the biodiversity of the Atlantic Forest? A multitaxa study in southern Bahia, Brazil. J. For. Res. 24, 263–272. https://doi.org/10.1007/s11676-012-0311-z (2013).Article 

Google Scholar 
68.Martin, P. S., Gheler-Costa, C., Lopes, P. C., Rosalino, L. M. & Verdade, L. M. Terrestrial non-volant small mammals in agro-silvicultural landscapes of Southeastern Brazil. For. Ecol. Manag. 282, 185–195. https://doi.org/10.1016/j.foreco.2012.07.002 (2012).Article 

Google Scholar 
69.Fayle, T. M. et al. Oil palm expansion into rain forest greatly reduces ant biodiversity in canopy, epiphytes and leaf-litter. Basic Appl. Ecol. 11, 337–345. https://doi.org/10.1016/j.baae.2009.12.009 (2010).Article 

Google Scholar 
70.Koh, L. P. Can oil palm plantations be made more hospitable for forest butterflies and birds?. J. Appl. Ecol. 45, 1002–1009. https://doi.org/10.1007/s10531-009-9760-x (2008).Article 

Google Scholar 
71.Martins, C. A. & Júnior, A. P. P. Production of biodiesel: Source strategies and efficiency in the Brazilian energy matrix. Energy Sour. Part A Recov. Util. Environ. Eff. 38, 277–285. https://doi.org/10.1080/15567036.2012.716139 (2016).CAS 
Article 

Google Scholar 
72.Peres, C. A. Why we need megareserves in Amazonia. Cons. Biol. 19, 728–733. https://doi.org/10.1111/j.1523-1739.2005.00691.x (2005).Article 

Google Scholar  More

Sampling site and sediment collectionPaleolimnological analysis by using sediment core samples were applied to reconstruct historical variations of Daphnia eDNA concentrations in Lake Biwa (Fig. 1b). Lake Biwa is the largest monomictic and mesotrophic lake in Japan. In this lake, during the last several decades the industrial revolution, multiple stressors of human origins impacted this ecosystem and the resident biological communities34,35,58. In our study, four 30-cm-long gravity core samples (namely LB1, LB2, LB4, and LB7; Fig. 1a–e) were collected on 17 August 2017 at the anchoring site of Hasu, a Center of Ecological Research boat from Kyoto University (Supplementary Fig. S2). A gravity corer with an inner diameter of 10.9 cm and a length of 30 cm was used to obtain the core samples. LB7 core was analyzed for chronological and reconstruction of temporal variation in Daphnia remain abundance (Fig. 1a,b). LB2 and LB1, LB4 cores were analyzed for reconstruction of temporal variation in sedimentary Daphnia DNA concentrations and resting egg production, respectively (Fig. 1b). Additionally, two 30-cm-long gravity core samples (namely IM1 and IM8; Fig. 1c,f) were collected at a pelagic site in the northern basin of Lake Biwa in August 2019 (Supplementary Fig. S2). The collected cores were sectioned at intervals of 1-cm thickness using a vertical extruder with a cutting apparatus, except for core number IM8, which was sectioned at 5-cm intervals (Fig. 1f). During the sectioning process, several millimeters of outer edge in each layer disturbed during the splitting process were carefully removed from the entire samples using a knife. After sectioning, each sliced sample were homogenized by shaking and then, all subsamples were taken from each homogenized sample. The pipes, knives, and cutting apparatus were cleaned with 0.6% sodium hypochlorite, tap water, and Milli-Q water to avoid DNA cross-contamination. Each sliced sample was transferred to lightproof bags and frozen at − 80 °C until further analysis.To examine the contamination due to core splitting, DNA extraction, and qPCR analysis, control water samples were inserted at a depth of 14.5–29.5 cm in the sediment cores, and the water samples for core IM1 were used as the negative control (Fig. 1c).Chronology of sediment coresSediment chronology was performed for the LB7 core based on the constant rate of supply (CRS) method of 210Pb dating59 and verified using the 137Cs peak traced in the period 1962 to 196360. Details of the chronological method have been reported elsewhere61. Briefly, dried samples were sealed in holders for a month to allow 222Rn and its short-lived decay product (214Pb) to equilibrate, which were determined by gamma counting using a germanium detector (GXM25P; EG & G ORTEC, Tokyo, Japan) equipped with a multi-channel analyzer (MCA7700; SEIKO EG & G, Tokyo, Japan) at the Center for Marine Environmental Studies, Ehime University. The activity of supported 210Pb was estimated by measuring the activity of 214Pb, whereas that of 210Pbexcess was determined according to the difference between the total and the supported 210Pb (210Pbexcess = 210Pbtotal − 214Pb). The age and age error of the remaining cores (LB1, LB2, and LB4) were indirectly estimated using stratigraphic correlations between the cores based on chronological controls in chlorophyll pigments and magnetic susceptibilities of the chronological LB7 core61. To compare these proxies, the marked peak or trough layers were used as reference layers (Supplementary Fig. S3).DNA extraction and purificationDNA extraction in the sediment samples was performed according to methods described in previous studies45,62. In brief, 9 g of each sediment sample was incubated at 94 °C for 50 min in a 9 mL alkaline solution comprising 6 mL of 0.33 M sodium hydroxide and a 3 mL Tris–EDTA buffer (pH 6.7). After centrifugation at 10,000×g for 60 min, 7.5 mL of the supernatant of the alkalized mixture was neutralized with 7.5 mL of 1 M Tris–HCl (pH 6.7). After adding 1.5 mL of 3 M sodium acetate (pH 5.2) and 30 mL absolute ethanol, the solution was preserved at − 20 °C for more than 1 h and then centrifuged at 10,000×g for 60 min. The pellet was transferred into a power bead tube that was installed in a fecal-soil DNA extraction kit (Power Soil DNA Isolation Kit, Qiagen, Germany). The ‘Experienced User Protocol 3 to 22’ of the Power Soil DNA Isolation Kit was followed. Finally, 200 μL of the DNA solution was obtained and stored at − 20 °C until qPCR analysis.12S rRNA gene primer-prove development for Daphnia geleata and Daphnia pulicaria
As the primer–probe for Daphnia galeata and D. pulicaria in qPCR analysis were not purchased by a company, thus we developed them for the two Daphnia species (see Supplementary Table S1). We preliminary obtained the mitochondrial 12S, 16S and COI gene of Daphnia genus from the National Center for Biotechnology Information (NCBI, http://www.ncbi.nlm.nih.gov/) and compared among them. From the preliminary results, we decided to use 12S because of the variability of sequences among Daphnia genus. Then we obtained the 12S sequences of Daphnia genus and other inhabiting plankton species in Lake Biwa, including Copepoda. We designed the primer–probe using Primer3plus (https://www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi). The reference sequences for the targeted gene regions are queried for potential amplicons between 50–150 bp using NCBI primer blast. The specificities of the primers and probes were then assessed in silico with homologous sequences from other Daphnia species in Japan using NCBI targeting 154 bp of the mitochondrial 12S rRNA gene. Once suitable amplicons are found the respective primers and probes are tested against template DNA originating from the species of D. galeata and D. pulicaria to verify amplification. During the in silico screening for specificity, we performed Primer-BLAST (http://www.ncbi.nlm.nih.gov/tools/primer-blast/). We checked all species from Japan of the order Daphnia. Using the D. galeata primer-set, we did not detect any Daphnia species. However, the D. pulicaria primer can amplify D. pulex DNA, as these species are known to have very similar sequences47. In Lake Biwa, another subgenus Daphnia (D. pulex group) different from D. pulicaria was temporally found around during the 1920s, although thereafter it was never reported47. Thus, the D. pulicaria primer may temporally detect another subgenus Daphnia (D. pulex group). However, their appearance time do not overlap, therefore we used the primer for our measurement to detect D. pulicaria during the last several decades.Quantitative PCRThe DNA samples were quantified by real-time TaqMan quantitative PCR using the PikoReal Real-Time PCR system (Thermo Fisher Scientific, Waltham, MA, USA). The primer–probe sets for the two Daphnia species were used for qPCR (Supplementary Table S1). The TaqMan reaction contained 900 nM of each forward and reverse primer, 125 nM TaqMan-Probe, 5 μL qPCR master mix (TaqPath; Thermo Fisher Scientific), and 2.0 μL sedimentary DNA solution. The final volume of the PCR was 10 μL after adding distilled water (DW). The qPCR conditions were as follows: 50 °C for 2 min and 95 °C for 10 min, followed by 50 cycles of 95 °C for 15 s and 60 °C for 60 s. We used a dilution series of 10,000, 1000, 100, and 10 copies per PCR reaction (n = 4) for the standard curve using the target DNA cloned into a plasmid. The R2 values of the standard curves ranged from 0.988 to 0.996 (PCR efficiencies = 93.1–102.0%). The quantitative data of the DNA copies (copies g−1 dry sed.) were reported by mean values ± standard deviation, which were calculated from DNA copies µL−1 PCR reaction with four replicates including zero (i.e., no detection). We also performed four replicates for each sample and an NTC (n = 4). No positives were detected from the NTC and the negative control of DNA extraction, confirming that there was no cross-contamination in any of the DNA measurements.To confirm primer specificity, an in vivo test for the primer/probe set was performed using the extracted DNA (10 pg per PCR reaction, n = 4) of D. galeata and D. pulicaria. In addition, qPCR amplicons were sequenced directly from a positive PCR from each site (n = 21) after treatment with ExoSAP-IT (USB Corporation, Cleveland, OH, USA). Sequences were determined using a commercial sequencing service (Eurofins Genomics, Tokyo, Japan).Inhibitor testSpike tests were performed for the LB2 core sample to evaluate the PCR inhibition effect of several substances and minerals in the sediment samples (Fig. 1c). For the spike test, 1 µL plasmid, including the internal positive control (IPC, 207-bp, Nippon Gene Co. Ltd., Tokyo, Japan;100 copies per PCR reaction), was added to the PCR template with 1.6 µL DNA-free DW. We used the primer and probe sets for IPC as follows:

IPC1-5′: CCGAGCTTACAAGGCAGGTT

IPC1-3′: TGGCTCGTACACCAGCATACTAG

IPC1-Taq: [FAM] TAGCTTCAAGCATCTGGCTGTCGGC [TAMRA].

To measure the relative degree of PCR inhibition in the samples, the Ct shift was compared between the samples and controls with the same number of known target DNA copies. The presence of PCR inhibitors was evaluated as ΔCt = Ct sample − Ct positive control. ΔCt ≥ 3 cycles was considered evidence of inhibition63 because the presence of PCR inhibitors will delay the Ct with a given quantity of template DNA.
Daphnia abundance and resting egg production as potential sources of Daphnia DNA archived in sedimentsTo unveil the potential source of sedimentary DNA of Daphnia, we reconstructed the historical variation in Daphnia abundance by counting remains of the post abdominal claw for LB7 core. There are two dominant Daphnia species: D. galeata Sars (Hyalodaphnia) and D. pulicaria Forbes (Daphnia)47,61, which have different post-abdominal claw characteristics64 and are known to be preserved in centuries-old sediments65. The post-abdominal claw remains were counted for core LB7 from the surface to a depth layer of 21.5 cm and additionally 23.5 cm, 25.5 cm, and 29.5 cm, totaling 25 samples, though each layer was expressed as mid-depth; e.g., 0.5 cm for the 0–1 cm depth layer. The enumeration method was based on a simplified standard method65 as previously reported29.Daphnia resting eggs enveloped by thickened carapaces, referred to as ephippial cases, and these ephippia can be preserved in sediments for decades to centuries29,30,33. In Lake Biwa, Daphnia species in Lake Biwa are distinguished on the basis of the size of the ephippium, with a boundary length of approximately 860 μm between them61. We collected ephippia from the surface to a depth layer of 29.5 cm for cores LB1 and LB4 (except for several layers of the LB1 core), totaling 56 samples (Supplementary Table S4). A detailed method for collecting ephippia is described in a previous study61. The total number of collected ephippia with an almost perfect shape, namely complete formation, or with a partial body constituting more than half of the original shape, namely incomplete formation, are shown in Supplementary Table S4. In our study, at least 16 ephippia in each sample were measured by photographs taken by a digital camera, excluding those from the samples in which fewer than 16 complete ephippia were detected (Supplementary Fig. S4). Species identification was then performed based on length.To determine whether the Daphnia sedimentary DNA concentrations were regulated by DNA derived directly from Daphnia remains or ephippia included in the analytical sediment, we divided the sediment sample into two fractions to exclude the remains and ephippia (Fig. 1d). The minimum size of Daphnia remains in this lake was approximately 55 μm (Tsugeki et al., in preparation). The analytical sediments for DNA extraction were divided into particles  38 μm using 38-μm mesh sieves on three-layer samples (specifically, LB2-5; 4.5 cm, LB2-7; 6.5 cm, and LB2-17; 16.5 cm expressed in middle depth of each sample) for core sample LB2, whose layers were known to include abundant ephippia and Daphnia remains. Furthermore, to test the possibility of the vertical movement of Daphnia sedimentary DNA through pore waters, we examined the sedimentary DNA concentration in pore water and its residual sediment by qPCR analysis (Fig. 1e). All DNA extractions were evaluated for sediment with and without sieves, and pore waters and the associated residual sediment samples were evaluated according to previous studies45.Measurement of DNA concentration in sediment ephippiaTo determine the potential source of sedimentary Daphnia DNA, we quantified the DNA concentration extracted from several ephippia obtained from the 0–5 cm and 5–10 cm layers of core IM8 using qPCR analysis (Fig. 1f). We selected 34 and 23 ephippia for D. galeata and D. pulicaria, respectively. We then measured the ephippial lengths and determined whether they contained resting eggs using a microscope. Among the selected ephippia, the well-preserved 17 ephippia with almost complete formation were set aside and grouped into 6 samples together in two or three ephippia for DNA analysis (Supplementary Table S5). Grouping was performed because of the low DNA concentrations typically associated with individual ephippium61.Possible factors regulating sedimentary Daphnia DNATo explore potential factors regulating temporal variation in sedimentary DNA concentrations, we analyzed chlorophyll pigments and algal remains. Sedimentary pigments of chlorophyll a were investigated for the LB 2 core, and algal remains were investigated for the LB7 core (Fig. 1a). Details of the method used for chlorophyll-a and algal remains are described in previous study61. In short, the concentrations of chlorophyll-a and phaeopigments were calculated according to the method66 and the diatom remains were analyzed according to the simplified method67. Green algae, Micrasterias hardyi, Staurastrum dorsidentiferum, S. arctiscon, S. limneticum, S. pingue, and Pediastrum biwae, were enumerated in a Sedgewick–Rafter chamber, following the method of zooplankton enumeration.Data analysisRegression models along with the standardized major axis method were used to determine the relationship between the sedimentary DNA concentration obtained from qPCR analysis and abundance or resting egg production in the sediment layers. Since qPCR (LB2), remains (LB7), and ephippia (LB1, LB4) analyses were performed on different cores, the chronological age of each analytical sample differed slightly. Therefore, prior to performing the statistical analysis, the sedimentary DNA (LB2) and ephippia data (LB1, LB4) in each chronological age were converted to annual data by linear interpolation and averaged for the year corresponding to the period in each sample of the chronology core (LB7). This conversion was possible because the time resolution at 1-cm intervals represented several years, depending on the sediment depth29,61. We employed the Gaussian type II model because our preliminary evaluation showed higher R2 values for type II regression models with a Gaussian distribution than for those with a logarithmic distribution, in all cases. All statistical analyses were performed using R ver. 4.0.3 (R Core Team 2020) with the package “smatr” ver. 3.4-8 for type II regressions. The significant criteria of all analyses were set as α = 0.05. In addition, to explore the potential environmental factors driving temporal variation in sedimentary DNA concentrations, we performed Pearson’s correlation analysis among the sedimentary DNA concentrations, chlorophyll a concentration, and algal remains using the SPSS version 20.0 statistical package. More

Changepoint detection methodAlthough most of the reflectance time series used in the BinSeg–Normal–Mean and BinSeg–Normal–MeanVar algorithms had a normal distribution, several lagoons had distributions that were skewed or did not follow a normal distribution (Fig. S1). However, results suggested that the accuracy of the detected changepoints were not sensitive to the normality assumption or distributional characteristics.The BinSeg-Normal-Mean algorithm had the highest performance (81% of the 340 validation sites) in detecting the correct year of swine waste lagoon construction, followed by BinSeg-Normal-MeanVar (77%). The two algorithms did not detect the same year of construction for 19 waste lagoons; of these 19, the BinSeg-Normal-Mean detected the correct year for 84% of them, while the BinSeg-Normal-MeanVar detected the correct year for only 16%. Therefore, the BinSeg-Normal-MeanVar algorithm was abandoned given it did not provide additional useful information relative to the BinSeg-Normal-Mean algorithm.Despite good performance, the BinSeg-Normal-Mean algorithm consistently detected a changepoint during the period of record for all sites included in the 10% validation set (n = 340 swine waste lagoons). However, 58 of the 340 swine waste lagoons were constructed prior to 1986, before the period of record suitable for detecting an accurate changepoint. Changepoints before 1986 either (1) detected the correct construction year, or (2) incorrectly detected a changepoint due to artifact signals identified on the images taken in 1984, probably associated with the initial satellite commissioning. In the latter circumstance, if the algorithm detected a changepoint due to this signal, it meant that no land-use change was detected after 1986. Therefore, these waste lagoons were estimated as having been constructed before 1986. In some conditions, when a large number of images was available for the year 1985 and 1986, the algorithm was able to detect the changepoint occurring for the years 1985 or 1986. Further, the BinSeg-Normal-Mean algorithm detected a false year of construction for swine waste lagoons for which the mean of the segment after the changepoint (S2) had a greater average than the segment before the changepoint (S1).To increase algorithm performance, we developed a workflow to address some of the aforementioned caveats (Fig. 4). In this workflow, the BinSeg-Normal-Mean algorithm is applied to a B4 reflectance time series at location j. If the BinSeg-Normal-Mean changepoint is identified for a time in or prior to 1986 (Fig. 4a,i,b,i) we assume that the lagoon was constructed in or prior to 1986. Similarly, a lagoon is assumed to be constructed in or prior to 1986 if a BinSeg-Normal-Mean changepoint is identified after 1986 and the mean of S2 is greater than the mean of S1 (Fig. 4a,ii,b,ii). If a changepoint occurred after 1986 and the mean of S1 was greater than S2, then the changepoint was estimated as having occurred between 1987 and 2010 (Fig. 4a,iii,b,iii).Figure 4Changepoint detection algorithm for determining the year of construction of swine waste lagoons. Panel (a) summarizes the algorithm workflow, while panel (b) illustrates specific examples corresponding to each step (i–iii) in the workflow.Full size imageThe performance of the workflow was evaluated using the validation set composed of 10% of the total number of swine waste lagoons (n = 340). With the new approach, 94% of the swine waste lagoon construction years (+ /- one year) were accurately retrieved. A tolerance of + /− 1 year was chosen to account for a lack of images in some years due to issues with image quality (e.g. high cloud cover) (e.g., Fig. 5a), or because construction spanned at least a year (e.g., Fig. 5b). The changepoint detection workflow incorrectly estimated the construction years for 19 of the 340 swine waste lagoons in the validation set; the differences between the observed and predicted years of construction of these lagoons ranged from 2 to 26 years with a median of 8 years.Figure 5Examples of limitations to the changepoint detection algorithm. In some cases, an insufficient number of high-quality Landsat 5 images were available to capture the year of construction of an individual swine waste lagoon (a), resulting in errors of + /− 1 year. In other cases, the changepoint algorithms detected the start of the construction of the swine waste lagoon but the swine waste lagoon was not fully operational until later years due to prolonged construction timelines (b).Full size imageBy visually inspecting historical Google Earth images for each of the lagoon sites for which the model incorrectly estimated construction year, we identified that model errors were associated with swine waste lagoon expansion, pixel transitions to land-use classes other than swine waste lagoons, or issues with pixels being partly covered by clouds or incompletely covered by the lagoon (i.e., narrow and small waste lagoons that do not entirely cover a pixel).Estimating swine waste lagoon construction yearsUsing the newly developed algorithm (Fig. 4), construction years were estimated for each swine waste lagoon in the NC Coastal Plain (Fig. 6); the years of construction for each swine waste lagoon are included in the supplementary material. Most swine waste lagoons were built in the early 90s and prior to the moratorium of 1997. More specifically, 80% of the swine waste lagoons (n = 2,736) were built between 1987 and 1997. Sixteen percent of the swine waste lagoons were constructed in or prior to 1986. A large decrease in the construction of swine waste lagoons occurred after the moratorium of 1997, with only 3.7% of swine waste lagoons being constructed after the moratorium. These results suggest that the 1997 moratorium did not completely halt the construction of lagoons, but dramatically slowed the rate of expansion.Figure 6Spatiotemporal distribution of swine waste lagoon construction (+/- 1 year) across the HUC6 watersheds. This figure was produced using QGIS version QGIS 3.18.3 (https://www.qgis.org/).Full size imageWith regards to hydrological boundaries (Fig. 7a–h), the Cape Fear River watershed had the highest number of swine waste lagoons (i.e., 56%; Fig. 7b), followed by the Neuse River (i.e., 23%; Fig. 7d), the Lower Pee Dee River (i.e., 9%; Fig. 7c) watersheds. The Albemarle-Chowan (Fig. 7a), Onslow Bay (Fig. 7e), Pamlico (Fig. 7f), Roanoke (Fig. 7g), and Upper Pee Dee (Fig. 7h) watersheds all had less than 9% of the total lagoons within the study area.Figure 7Year of construction of the swine waste lagoons (+ /− 1 year) for the HUC6 watersheds. The y-axis scales are unequal between the plots to improve readability. The dashed red lines correspond to the establishment of the moratorium in 1997.Full size imageResults suggested that the Cape Fear River watershed was the center of the historical growth of the swine industry, where over 300 swine waste lagoons were built prior to 1987. The Cape Fear River watershed experienced a steady increase in the number of swine waste lagoons from 1987 to 1990, with an average of 46 swine waste lagoons being built annually. However, after 1991, the pace of swine waste lagoon construction increased dramatically with an average of 192 swine waste lagoons built annually between 1991 and 1997. The highest construction rate occurred in 1994, with 242 swine waste lagoons built. However, after the 1997 moratorium, the construction rate decreased dramatically; in 1997, 153 swine waste lagoons were constructed, and this number dropped to 23 in 1998. After 1998, the annual average number of swine waste lagoons constructed plunged to 5. Although the swine waste lagoon construction rate fell considerably after the 1997 moratorium, the decrease had already started in 1995. The same pattern was observed for the Neuse, Pamlico, Albemarle-Pamlico, and Onslow Bay watersheds.The spatiotemporal distribution of swine waste lagoons at the HUC12 watershed scale emphasized the historical clustering of the swine industry in the NC Coastal Plain. After the moratorium, swine waste lagoons were present within 436 HUC12 watersheds. However, before 1986, they were spread across only 197 HUC12 watersheds (Fig. 8). Before 1986, the density of waste lagoons was relatively low with an average of 3.38 swine waste lagoons per 100 km2 and a maximum of 15.13 swine waste lagoons per 100 km2 (i.e., Clayroot Swamp-Swift Creek watershed) (Fig. 8). In the 90s, swine waste lagoon construction expanded and continued to intensify in the region. After the moratorium of 1997, the average density of waste lagoons per HUC12 watersheds was 10 per 100 km2 with a maximum of 78 waste lagoons per 100 km2 identified in the Maxwell Creek-Stocking Head Creek basin. After 1997, 16 of 436 HUC12 watersheds had a swine waste lagoon density greater than 40 per 100 km2 (Fig. 8).Figure 8Cumulative swine waste lagoon density per 100 km2 reported at the HUC12 watershed scale; HUC6 watersheds shown in gray for reference. This figure was produced using QGIS version QGIS 3.18.3 (https://www.qgis.org/).Full size imageSpatiotemporal distribution of swine waste lagoons in relation to water resourcesDistance of swine waste lagoon sites to the nearest water feature (i.e., reservoir, canal/ditch, lake/pond, stream/river, estuary) were assessed using the NHD. The analysis revealed that over 150 swine waste lagoons were misclassified by the NHD and were documented in the NHD as lake/pond (n = 102) or swamp/marsh (n = 46). Further, we observed that some NHD water features were misclassified as other non-water features (e.g., forest, pasture), and most of these misclassifications were for polygons with an area less than 0.05 km2. Therefore, NHD water features with areas less than 0.05 km2 were removed from subsequent analyses. Distances between swine waste lagoons and waterways were computed from the NHD without features with areas less than 0.05 km2. The new analysis revealed that 3 swine waste lagoons remained misclassified as lake/pond (n = 1) and swamp/marsh (n = 2). Canal/Ditch, lake/pond, stream/river, and swamp/marsh were identified as the NHD features that were most commonly near swine waste lagoons (Fig. 9). Two swine waste lagoons were near a reservoir in which one was identified as a treatment-sewage pond by the NHD.Figure 9Nearest water features distance to swine waste lagoons.Full size imageThe average and median distance of all swine waste lagoons (including those built early and late in the period of record) to the nearest water features were 234 and 177 m, respectively. Further, 92% of the swine waste lagoons were less than 500 m from the nearest waterways. The Mann–Kendall results revealed a significant upward trend over time of swine waste lagoon distances to the nearest water features (alpha = 0.05, p-value = 0.01). A slight increase over time of swine waste lagoon distances to the nearest water feature is also documented in Table 1.Table 1 Temporal average and median of nearest distance (m) of swine waste lagoons to water features. NA indicated that the water feature was not the closest waterway to any of the studied swine waste lagoons for the time period.Full size table More

1.MacLachlan, N. J. & Guthrie, A. J. Re-emergence of bluetongue, African horse sickness, and other Orbivirus diseases. Vet. Res. 41, 35 (2010).Article 

Google Scholar 
2.Zientara, S., Weyer, C. T. & Lecollinet, S. African horse sickness. OIE Revue Sci. Tech. 34, 315–327 (2015).CAS 
Article 

Google Scholar 
3.Ayelet, G. et al. Outbreak investigation and molecular characterization of African horse sickness virus circulating in selected areas of Ethiopia. Acta Trop. 127, 91–96 (2013).Article 

Google Scholar 
4.Diarra, M. et al. Spatial distribution modelling of Culicoides (Diptera: Ceratopogonidae) biting midges, potential vectors of African horse sickness and bluetongue viruses in Senegal. Parasit. Vectors 11, 1–15 (2018).Article 

Google Scholar 
5.Karamalla, S. T. et al. Sero-epidemioloical survey on African horse sickness virus among horses in Khartoum State, Central Sudan. BMC Vet. Res. 14, 1–6 (2018).Article 

Google Scholar 
6.Escobar, L. E. Ecological Niche modeling: An introduction for veterinarians and epidemiologists. Front. Vet. Sci. 7, 519059. https://doi.org/10.3389/fvets.2020.519059 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
7.Okely, M., Anan, R., Gad-Allah, S. & Samy, A. M. Mapping the environmental suitability of etiological agent and tick vectors of Crimean-Congo hemorrhagic fever. Acta Trop. 203, 105319 (2020).CAS 
Article 

Google Scholar 
8.Chavy, A. et al. Ecological niche modelling for predicting the risk of cutaneous leishmaniasis in the Neotropical moist forest biome. PLoS Negl. Trop. Diseases 13, e0007629 (2019).Article 

Google Scholar 
9.Sloyer, K. E. et al. Ecological niche modeling the potential geographic distribution of four Culicoides species of veterinary significance in Florida, USA. PLoS ONE 14, e0206648 (2019).CAS 
Article 

Google Scholar 
10.Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).Article 

Google Scholar 
11.Cao, Z., Jin, Y., Shen, T., Xu, F. & Li, Y. Risk factors and distribution for peste des petits ruminants (PPR) in Mainland China. Small Rumin. Res. 162, 12–16 (2018).Article 

Google Scholar 
12.Naimi, B. & Araújo, M. B. sdm: a reproducible and extensible R platform for species distribution modelling. Ecography 39, 368–375 (2016).Article 

Google Scholar 
13.Naimi, B., Hamm, N. A. S., Groen, T. A., Skidmore, A. K. & Toxopeus, A. G. Where is positional uncertainty a problem for species distribution modelling. undefined 37, 191–203 (2014).
Google Scholar 
14.R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, (2020).15.Thuiller, W., Lafourcade, B., Engler, R. & Araújo, M. B. BIOMOD—a platform for ensemble forecasting of species distributions. Ecography 32, 369–373 (2009).Article 

Google Scholar 
16.Araújo, M. B. & New, M. Ensemble forecasting of species distributions. Trends Ecol. Evol. 22, 42–47 (2007).Article 

Google Scholar 
17.Uusitalo, R. et al. Predicting spatial patterns of sindbis virus (Sinv) infection risk in finland using vector, host and environmental data. Int. J. Environ. Res. Public Health 18, 7064 (2021).Article 

Google Scholar 
18.Raffini, F. et al. From nucleotides to satellite imagery: Approaches to identify and manage the invasive pathogen Xylella fastidiosa and its insect vectors in Europe. Sustainability (Switzerland) 12, 4508 (2020).CAS 
Article 

Google Scholar 
19.Phillips, S. B., Aneja, V. P., Kang, D. & Arya, S. P. Maximum entropy modeling of species geographic distributions. Ecol. Model. 190, 231–259 (2006).Article 

Google Scholar 
20.Barbet-Massin, M., Jiguet, F., Albert, C. H. & Thuiller, W. Selecting pseudo-absences for species distribution models: how, where and how many?. Methods Ecol. Evol. 3, 327–338 (2012).Article 

Google Scholar 
21.Hernández-Urcera, J., Murillo, F. J., Regueira, M., Cabanellas-Reboredo, M. & Planas, M. Preferential habitats prediction in syngnathids using species distribution models. Marine Environ. Res. 172, 105488 (2021).Article 

Google Scholar 
22.Smeraldo, S. et al. Generalists yet different: distributional responses to climate change may vary in opportunistic bat species sharing similar ecological traits. Mammal Rev. 51, 571–584 (2021).Article 

Google Scholar 
23.Thomson, A. M. et al. RCP4.5: A pathway for stabilization of radiative forcing by 2100. Clim. Change 109, 77–94 (2011).ADS 
CAS 
Article 

Google Scholar 
24.QGIS Development Team. QGIS Geographic Information System. Open-Source Geospatial Foundation Project. (2020).25.Ramirez-Reyes, C. et al. Embracing ensemble species distribution models to inform at-risk species status assessments. J. Fish Wildl. Manag. 12, 98–111 (2021).Article 

Google Scholar 
26.Stephenson, F. et al. Presence-only habitat suitability models for vulnerable marine ecosystem indicator taxa in the South Pacific have reached their predictive limit. ICES J. Mar. Sci. 78, 2830–2843 (2021).Article 

Google Scholar 
27.Zhu, G., Fan, J. & Peterson, A. T. Cautions in weighting individual ecological niche models in ensemble forecasting. Ecol. Modelling 448, 109502 (2021).Article 

Google Scholar 
28.Leta, S. et al. Modeling the global distribution of Culicoides imicola: an Ensemble approach. Sci. Rep. 9, 1–9 (2019).ADS 
CAS 
Article 

Google Scholar 
29.Onyango, M. G. et al. Delineation of the population genetic structure of Culicoides imicola in East and South Africa. Parasit. Vectors 8, 660 (2015).Article 

Google Scholar 
30.Carpenter, S., Mellor, P. S., Fall, A. G., Garros, C. & Venter, G. J. African horse sickness virus: history. Transm. Curr. Status. 62, 343–358. https://doi.org/10.1146/annurev-ento-031616-035010 (2017).CAS 
Article 

Google Scholar 
31.Carpenter, S., Mellor, P. S., Fall, A. G., Garros, C. & Venter, G. J. African Horse Sickness Virus: History, Transmission, and Current Status. Annu. Rev. Entomol. 62, 343–358 (2017).CAS 
Article 

Google Scholar 
32.Fall, M. et al. Culicoides (Diptera: Ceratopogonidae) midges, the vectors of African horse sickness virus—a host/vector contact study in the Niayes area of Senegal. Parasit. Vectors 8, 1–13 (2015).Article 

Google Scholar 
33.Mellor, P. S. Epizootiology and vectors of African horse sickness virus. Comp. Immunol. Microbiol. Infect. Dis. 17, 287–296 (1994).CAS 
Article 

Google Scholar 
34.Wu, X., Lu, Y., Zhou, S., Chen, L. & Xu, B. Impact of climate change on human infectious diseases: Empirical evidence and human adaptation. Environ. Int. 86, 14–23 (2016).Article 

Google Scholar 
35.Nosrat, C. et al. Impact of recent climate extremes on mosquito-borne disease transmission in Kenya. PLOS Negl. Trop. Diseases 15, e0009182 (2021).CAS 
Article 

Google Scholar 
36.Abiodun, G. J., Maharaj, R., Witbooi, P. & Okosun, K. O. Modelling the influence of temperature and rainfall on the population dynamics of Anopheles arabiensis. Malar. J. 15, 1–15 (2016).Article 

Google Scholar  More

1.Holway, D. A. & Suarez, A. V. Animal behavior: An essential component of invasion biology. TREE 14, 328–330 (1999).CAS 
PubMed 

Google Scholar 
2.Chapple, D. G., Simmonds, S. M. & Wong, B. B. M. Can behavioral and personality traits influence the success of unintentional species introductions? Trends Ecol. Evol. 27, 57–64 (2012).PubMed 

Google Scholar 
3.Weis, J. & Sol, D. Behaviour and the Invasion Process. in Biological Invasions and Animal Behaviour 5–116 (Cambridge University Press, 2016).4.Cote, J., Fogarty, S., Weinersmith, K., Brodin, T. & Sih, A. Personality traits and dispersal tendency in the invasive mosquitofish (Gambusia affinis). Proc. R. Soc. B Biol. Sci. 277, 1571–1579 (2010).
Google Scholar 
5.Myles-Gonzalez, E., Burness, G., Yavno, S., Rooke, A. & Fox, M. G. To boldly go where no goby has gone before: Boldness, dispersal tendency, and metabolism at the invasion front. Behav. Ecol. 26, 1083–1090 (2015).
Google Scholar 
6.Mutascio, H. E., Pittman, S. E. & Zollner, P. A. Investigating movement behavior of invasive Burmese pythons on a shy–bold continuum using individual-based modeling. Perspect. Ecol. Conserv. 15, 25–31 (2017).
Google Scholar 
7.Chuang, A. Living Life on the Edge: The Role of Invasion Processes in Shaping Personalities in a Non-Native Spider Species (The University of Tennessee, Knoxville, 2019). https://doi.org/10.1017/CBO9781107415324.004.Book 

Google Scholar 
8.Blackburn, T. M. et al. A proposed unified framework for biological invasions. Trends Ecol. Evol. 26, 333–339 (2011).PubMed 

Google Scholar 
9.Pintor, L. M., Sih, A. & Kerby, J. L. Behavioral correlations provide a mechanism for explaining high invader densities and increased impacts on native prey. Ecology 90, 581–587 (2009).PubMed 

Google Scholar 
10.Petren, K. & Case, T. J. An experimental demonstration of exploitation competition in an ongoing invasion. Ecology 77, 118–132 (1996).
Google Scholar 
11.Wright, T. F., Eberhard, J. R., Hobson, E. A., Avery, M. L. & Russello, M. A. Behavioral flexibility and species invasions: The adaptive flexibility hypothesis. Ethol. Ecol. Evol. 22, 393–404 (2010).
Google Scholar 
12.Dick, J. T. A. Role of behaviour in biological invasions and species distributions; lessons from interactions between the invasive Gammarus pulex and the native G. duebeni (Crustacea: Amphipoda). Contrib. Zool. 77, 91–98 (2008).
Google Scholar 
13.Dick, J. T. A. et al. Invader Relative Impact Potential: A new metric to understand and predict the ecological impacts of existing, emerging and future invasive alien species. J. Appl. Ecol. 54, 1259–1267 (2017).
Google Scholar 
14.Dick, J. T. A., Elwood, R. W. & Montgomery, W. I. The behavioural basis of a species replacement: differential aggresssion and predation between the introduced Gammarus pulex and the native G. duebeni celticus (Amphipoda). Behav. Ecol. Sociobiol. 37, 393–398 (1995).
Google Scholar 
15.Dick, J. T. A. et al. Ecological impacts of an invasive predator explained and predicted by comparative functional responses. Biol. Invasions 15, 837–846 (2013).
Google Scholar 
16.Dick, J. T. A. et al. Advancing impact prediction and hypothesis testing in invasion ecology using a comparative functional response approach. Biol. Invasions 16, 735–753 (2014).
Google Scholar 
17.Iacarella, J. C., Dick, J. T. A. & Ricciardi, A. A spatio-temporal contrast of the predatory impact of an invasive freshwater crustacean. Divers. Distrib. 21, 803–812 (2015).
Google Scholar 
18.Toscano, B. J. & Griffen, B. D. Trait-mediated functional responses: Predator behavioural type mediates prey consumption. J. Anim. Ecol. 83, 1469–1477 (2014).PubMed 

Google Scholar 
19.MacCrimmon, H. R. World distribution of rainbow trout (Salmo gairdneri): further observations. J. Fish. Res. Board Canada 28, 663–704 (1971).
Google Scholar 
20.MacCrimmon, H. R., Marshall, T. L. & Gots, B. L. World distribution of brown trout, Salmo trutta: further observations. J. Fish. Res. Board Canada 27, 811–818 (1970).
Google Scholar 
21.Crawford, S. S. & Muir, A. M. Global introductions of salmon and trout in the genus Oncorhynchus: 1870–2007. Rev. Fish Biol. Fish. 18, 313–344 (2008).
Google Scholar 
22.Crowl, T. A., Townsend, C. R. & Mcintosh, A. R. The impact of introduced brown and rainbow trout on native fish: The case of Australasia. Rev. Fish Biol. Fish. 241, 217–241 (1992).
Google Scholar 
23.Hasegawa, K. Invasions of rainbow trout and brown trout in Japan: A comparison of invasiveness and impact on native species. Ecol. Freshw. Fish 29, 419–428 (2020).
Google Scholar 
24.Cambray, J. A. The global impact of alien trout species—A review; with reference to their impact in South Africa. African J. Aquat. Sci. 28, 61–67 (2003).
Google Scholar 
25.Dunham, J. B., Wheeler, A. & Rosenberger, A. Assessing the consequences of nonnative trout in headwater ecosystems in western North America. Fisheries 29, 37–41 (2004).
Google Scholar 
26.Fausch, K. D., Taniguchi, Y., Nakano, S., Grossman, G. D. & Townsend, C. R. Flood disturbance regimes influence rainbow trout invasion success among five holarctic regions. Ecol. Appl. 11, 1438–1455 (2001).
Google Scholar 
27.Anderson, R. M. & Nehring, R. B. Effects of a catch-and-release regulation on a wild trout population in Colorado and its acceptance by Anglers. North Am. J. Fish. Manag. 4, 257–265 (1984).
Google Scholar 
28.Young, K. A. et al. A trial of two trouts: Comparing the impacts of rainbow and brown trout on a native galaxiid. Anim. Conserv. 13, 399–410 (2010).
Google Scholar 
29.Conrad, J. L., Weinersmith, K. L., Brodin, T., Saltz, J. B. & Sih, A. Behavioural syndromes in fishes: A review with implications for ecology and fisheries management. J. Fish Biol. 78, 395–435 (2011).CAS 
PubMed 

Google Scholar 
30.Mowles, S. L., Cotton, P. A. & Briffa, M. Consistent crustaceans: The identification of stable behavioural syndromes in hermit crabs. Behav. Ecol. Sociobiol. 66, 1087–1094 (2012).
Google Scholar 
31.Sih, A., Bell, A. & Johnson, J. C. Behavioral syndromes: An ecological and evolutionary overview. Trends Ecol. Evol. 19, 372–378 (2004).PubMed 

Google Scholar 
32.Bell, A. M. Behavioural differences between individuals and two populations of stickleback (Gasterosteus aculeatus). J. Evol. Biol. 18, 464–473 (2005).CAS 
PubMed 

Google Scholar 
33.Bourne, G. R. & Sammons, A. J. Boldness, aggression and exploration: evidence for a behavioural syndrome in male pentamorphic livebearing fish, Poecilia parae. AACL Bioflux 1, 39–50 (2008).
Google Scholar 
34.Lukas, J. et al. Consistent behavioral syndrome across seasons in an invasive freshwater fish. Front. Ecol. Evol. 8, 466 (2021).ADS 

Google Scholar 
35.Gjedrem, T., Gjøen, H. M. & Gjerde, B. Genetic origin of Norwegian farmed Atlantic salmon. Aquaculture 98, 41–50 (1991).
Google Scholar 
36.Huntingford, F. & Adams, C. Behavioural syndromes in farmed fish: Implications for production and welfare. Behaviour 142, 1207–1221 (2005).
Google Scholar 
37.Alvarez, D. & Nicieza, A. G. Predator avoidance behaviour in wild and hatchery-reared brown trout : The role of experience and domestication. J. Fish Biol. 63, 1565–1577. https://doi.org/10.1046/j.1095-8649.2003.00267.x (2003).Article 

Google Scholar 
38.Geffroy, B. et al. Evolutionary dynamics in the anthropocene: Life history and intensity of human contact shape antipredator responses. PLoS Biol. 18, 1–17 (2020).
Google Scholar 
39.Lincoln, R. F. & Scott, A. P. Production of all-female triploid rainbow trout. Aquaculture 30, 375–380 (1983).
Google Scholar 
40.Maxime, V. The physiology of triploid fish: Current knowledge and comparisons with diploid fish. Fish Fish. 9, 67–78 (2008).
Google Scholar 
41.Chatterji, R., Longley, D., Sandford, D., Roberts, D. & Stubbing, D. Performance of stocked triploid and diploid brown trout and their effects on wild brown trout in UK rivers. (2008).42.Benfey, T. J. The physiology and behavior of triploid fishes. Rev. Fish. Sci. 7, 39–67 (1999).
Google Scholar 
43.Carter, C. G. et al. Food consumption, feeding behaviour, and growth of triploid and diploid Atlantic salmon, Salmo salar L., parr.. Can. J. Zool. 72, 609–617 (1994).
Google Scholar 
44.Weber, G. M., Hostuttler, M. A., Cleveland, B. M. & Leeds, T. D. Growth performance comparison of intercross-triploid, induced triploid, and diploid rainbow trout. Aquaculture 433, 85–93 (2014).
Google Scholar 
45.Øverli, Ø., Pottinger, T. G., Carrick, T. R., Øverli, E. & Winberg, S. Differences in behaviour between rainbow trout selected for high- and low-stress responsiveness. J. Exp. Biol. 205, 391–395 (2002).PubMed 

Google Scholar 
46.Sadoul, B., Leguen, I., Colson, V., Friggens, N. C. & Prunet, P. A multivariate analysis using physiology and behavior to characterize robustness in two isogenic lines of rainbow trout exposed to a confinement stress. Physiol. Behav. 140, 139–147 (2015).CAS 
PubMed 

Google Scholar 
47.Adriaenssens, B. & Johnsson, J. I. Learning and context-specific exploration behaviour in hatchery and wild brown trout. Appl. Anim. Behav. Sci. 132, 90–99 (2011).
Google Scholar 
48.Näslund, J. & Johnsson, J. I. State-dependent behavior and alternative behavioral strategies in brown trout (Salmo trutta L.) fry. Behav. Ecol. Sociobiol. 70, 2111–2125 (2016).PubMed 
PubMed Central 

Google Scholar 
49.Mortensen, E. Density-dependent mortality of trout fry (Salmo trutta L.) and its relationship to the management of small streams. J. Fish Biol. 11, 613–617 (1977).
Google Scholar 
50.Armstrong, J. D. & Nislow, K. H. Critical habitat during the transition from maternal provisioning in freshwater fish, with emphasis on Atlantic salmon (Salmo salar) and brown trout (Salmo trutta). J. Zool. 269, 403–413 (2006).
Google Scholar 
51.Walsh, R. N. & Cummins, R. A. The open-field test: A critical review. Psychol. Bull. 83, 482–504 (1976).CAS 
PubMed 

Google Scholar 
52.Adriaenssens, B. & Johnsson, J. I. Shy trout grow faster: Exploring links between personality and fitness-related traits in the wild. Behav. Ecol. 22, 135–143 (2010).
Google Scholar 
53.Sneddon, L. U. The bold and the shy: Individual differences in rainbow trout. J. Fish Biol. 62, 971–975 (2003).
Google Scholar 
54.Adriaenssens, B. Individual variation in behaviour: personality and performance of brown trout in the wild (University of Gothenburg, 2010).55.Elias, A., Thrower, F. & Nichols, K. M. Rainbow trout personality: Individual behavioural variation in juvenile Oncorhynchus mykiss. Behaviour 155, 205–230 (2018).
Google Scholar 
56.Dick, J. T. A. et al. Functional responses can unify invasion ecology. Biol. Invasions 19, 1667–1672 (2017).
Google Scholar 
57.Sloman, K. A., Metcalfe, N. B., Taylor, A. C. & Gilmour, K. M. Plasma cortisol concentrations before and after social stress in rainbow trout and brown trout. Physiol. Biochem. Zool. 74, 383–389 (2001).CAS 
PubMed 

Google Scholar 
58.Sadoul, B., Blumstein, D. T., Alfonso, S. & Geffroy, B. Human protection drives the emergence of a new coping style in animals. PLoS Biol. 19, 1–11 (2021).
Google Scholar 
59.Campbell, J. M., Carter, P. A., Wheeler, P. A. & Thorgaard, G. H. Aggressive behavior, brain size and domestication in clonal rainbow trout lines. Behav. Genet. 45, 245–254 (2015).PubMed 

Google Scholar 
60.Berejikian, B. A., Mathews, S. B. & Quinn, T. P. Effects of hatchery and wild ancestry and rearing environments on the development of agonistic behavior in steelhead trout (Oncorhynchus mykiss) fry. Can. J. Fish. Aquat. Sci. 53, 2004–2014 (1996).
Google Scholar 
61.Laverty, C. et al. Assessing the ecological impacts of invasive species based on their functional responses and abundances. Biol. Invasions 19, 1653–1665 (2017).
Google Scholar 
62.Alexander, M. E., Dick, J. T. A., Weyl, O. L. F., Robinson, T. B. & Richardson, D. M. Existing and emerging high impact invasive species are characterized by higher functional responses than natives. Biol. Lett. 10, 20130946 (2014).PubMed 
PubMed Central 

Google Scholar 
63.Dickey, J. W. E., Cuthbert, R. N., Steffen, G. T., Dick, J. T. A. & Briski, E. Sea freshening may drive the ecological impacts of emerging and existing invasive non-native species. Divers. Distrib. 27, 144–156 (2021).
Google Scholar 
64.Sadler, J., Pankhurst, P. M. & King, H. R. High prevalence of skeletal deformity and reduced gill surface area in triploid Atlantic salmon (Salmo salar L.). Aquaculture 198, 369–386 (2001).
Google Scholar 
65.Benfey, T. J. & Biron, M. Acute stress response in triploid rainbow trout (Oncorhynchus mykiss) and brook trout (Salvelinus fontinalis). Aquaculture 184, 167–176 (2000).CAS 

Google Scholar 
66.Sadler, J., Pankhurst, N. W., Pankhurst, P. M. & King, H. Physiological stress responses to confinement in diploid and triploid Atlantic salmon. J. Fish Biol. 56, 506–518 (2000).
Google Scholar 
67.Berrebi, P., Splendiani, A., Palm, S. & Berna, R. Genetic diversity of domestic brown trout stocks in Europe. Aquaculture 544, 737043 (2021).CAS 

Google Scholar 
68.Gross, R., Lulla, P. & Paaver, T. Genetic variability and differentiation of rainbow trout (Oncorhynchus mykiss) strains in northern and Eastern Europe. Aquaculture 272, 139–146 (2007).
Google Scholar 
69.Whelan, K. Assessing and mitigating the impact of a major rainbow trout escape on the wild salmon and trout populations of the Mourne river system, Northern Ireland. (2017).70.Shelton, J. et al. Temperature mediates the impact of non-native rainbow trout on native freshwater fishes in South Africa’s Cape Fold Ecoregion. Biol. Invasions 20, 2927–2944 (2018).
Google Scholar 
71.Michelangeli, M. et al. Sex-dependent personality in two invasive species of mosquitofish. Biol. Invasions 22, 1353–1364 (2020).
Google Scholar 
72.Friard, O. & Gamba, M. BORIS: A free, versatile open-source event-logging software for video/audio coding and live observations. Methods Ecol. Evol. 7, 1325–1330 (2016).
Google Scholar 
73.R Core Team. R: A language and environment for statistical computing. (2018).74.RStudio Team. RStudio Team (2020). RStudio: Integrated Development for R. RStudio, PBC, Boston, MA. http://www.rstudio.com/. 2019 (2020).75.Zuur, A. F., Ieno, E. N., Walker, N. J., Saveliev, A. A. & Smith, G. M. Mixed effects models and extensions in ecology with R. Springer https://doi.org/10.1086/648138 (2008).Article 
MATH 

Google Scholar 
76.Bates, D., Mächler, M., Bolker, B. M. & Walker, S. C. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 18637 (2015).
Google Scholar 
77.Wickham, H., François, R., Henry, L. & Müller, K. dplyr: A Grammar of Data Manipulation. R package version. Media https://doi.org/10.1007/978-0-387-98141-3 (2019).Article 

Google Scholar 
78.Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer-Verlag, 2016).MATH 

Google Scholar 
79.Barton, K. MuMIn: Multi-Model Inference. 2020 (2020).80.Lenth, R., Singmann, H., Love, J., Buerkner, P. & Herve, M. emmeans: estimated marginal means, aka least-squares means. R package version 1.5.2-1 (2020).81.Pritchard, D. frair: tools for functional response analysis. R package version 0.0.100 (2017).82.Juliano, S. A. Predation and functional response curves. in Design and Analysis of Ecological Experiments (eds. Scheiner, S. & Gurevitch, J.) Chapter 10 (2001).83.Rogers, D. Random search and insect population models. J. Anim. Ecol. 41, 369–383 (1972).
Google Scholar 
84.Bolker, B. M. Rogers random predator equation: extensions and estimation by numerical integration. 1–20 (2012). More

1.Angel, S. et al. The dimensions of global urban expansion: Estimates and projections for all countries, 2000–2050. Prog. Plan. 75, 53–107 (2011).
Google Scholar 
2.Grimm, N. B. et al. Global change and the ecology of cities. Science 319, 756–760 (2008).ADS 
CAS 
PubMed 

Google Scholar 
3.McKinney, M. L. Urbanization as a major cause of biotic homogenization. Biol. Conserv. 127, 247–260 (2006).
Google Scholar 
4.Groffman, P. M. et al. Ecological homogenization of urban USA. Front. Ecol. Environ. 12, 74–81 (2014).
Google Scholar 
5.Bolnick, D. I. et al. (Non)Parallel evolution. Annu. Rev. Ecol. Evol. Syst. 49, 303–330 (2018).
Google Scholar 
6.Donihue, C. M. & Lambert, M. R. Adaptive evolution in urban ecosystems. Ambio 44, 194–203 (2015).PubMed 

Google Scholar 
7.Johnson, M. T. J. & Munshi-South, J. Evolution of life in urban environments. Science 358, eaam8327 (2017).
Google Scholar 
8.Rivkin, L. R. et al. A roadmap for urban evolutionary ecology. Evol. Appl. 12, 384–398 (2019).PubMed 

Google Scholar 
9.Santangelo, J. S. et al. Urban environments as a framework to study parallel evolution. In Urban Evolutionary Biology (eds Szulkin, M. et al.) (Oxford University Press, 2020).
Google Scholar 
10.Cosentino, B. J., Moore, J.-D., Karraker, N. E., Ouellet, M. & Gibbs, J. P. Evolutionary response to global change: Climate and land use interact to shape color polymorphism in a woodland salamander. Ecol. Evol. 7, 5426–5434 (2017).PubMed 
PubMed Central 

Google Scholar 
11.Koprowski, J. L., Munroe, K. E. & Edelman, A. J. Gray not grey: Ecology of Sciurus carolinensis in their native range in North America. In Grey Squirrels: Ecology and Management of an Invasive Species in Europe (eds Shuttleworth, C. M. et al.) (European Squirrel Initiative, 2016).
Google Scholar 
12.McRobie, H., Thomas, A. & Kelly, J. The genetic basis of melanism in the gray squirrel (Sciurus carolinensis). J. Hered. 100, 709–714 (2009).CAS 
PubMed 

Google Scholar 
13.Gibbs, J. P., Buff, M. F. & Cosentino, B. J. The biological system: Urban wildlife, adaptation and evolution: Urbanization as a driver of contemporary evolution in gray squirrels (Sciurus carolinensis). In Understanding Urban Ecology (eds Hall, M. A. & Balogh, S.) (Springer, 2019).
Google Scholar 
14.Lehtinen, R. M. et al. Dispatches form the neighborhood watch: Using citizen science and field survey data to document color morph frequency in space and time. Ecol. Evol. 10, 1526–1538 (2020).PubMed 
PubMed Central 

Google Scholar 
15.Perlut, N. G. Long-distance dispersal by eastern gray squirrels in suburban habitats. Northeast. Nat. 27, 195–200 (2020).
Google Scholar 
16.Goheen, J. R., Swihart, R. K., Gehring, T. M. & Miller, M. S. Forces structuring tree squirrel communities in landscapes fragmented by agriculture: Species differences in perceptions of forest connectivity and carrying capacity. Oikos 102, 95–103 (2003).
Google Scholar 
17.Ducharme, M. B., Larochelle, J. & Richard, D. Thermogenic capacity in gray and black morphs of the gray squirrel, Sciurus carolinensis. Physiol. Zool. 62, 1273–1292 (1989).
Google Scholar 
18.Linnen, C. R. & Hoekstra, H. E. Measuring natural selection on genotypes and phenotypes in the wild. Cold Spring Harb. Symp. Quant. Biol. 74, 155–168 (2010).PubMed Central 

Google Scholar 
19.Campbell-Staton, S. C. et al. Parallel selection on thermal physiology facilitates repeated adaptation of city lizards to urban heat islands. Nat. Ecol. Evol. 4, 652–658 (2020).PubMed 

Google Scholar 
20.Reid, N. M. et al. The genomic landscape of rapid repeated evolutionary adaptation to toxic pollution in wild fish. Science 354, 1305–1308 (2016).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
21.Bowers, M. A. & Breland, B. Foraging of gray squirrels on an urban-rural gradient: Use of the GUD to assess anthropogenic impact. Ecol. Appl. 6, 1135–1142 (1996).
Google Scholar 
22.McCleery, R. A., Lopez, R. R., Silvy, N. J. & Gallant, D. L. Fox squirrel survival in urban and rural environments. J. Wildl. Manage. 72, 133–137 (2008).
Google Scholar 
23.Benson, E. The urbanization of the eastern gray squirrel in the United States. J. Am. Hist. 100, 691–710 (2013).
Google Scholar 
24.Leveau, L. United colours of the city: A review about urbanization impact on animal colours. Austral Ecol. 46, 670–679 (2021).
Google Scholar 
25.Ducrest, A.-L., Keller, L. & Roulin, A. Pleiotropy in the melanocortin system, coloration, and behavioural syndromes. Trends Ecol. Evol. 23, 502–510 (2008).PubMed 

Google Scholar 
26.Stothart, M. R. & Newman, A. E. M. Shades of grey: Host phenotype dependent effect of urbanization on the bacterial microbiome of a wild mammal. Anim. Microbiome. 3, 46 (2021).PubMed 
PubMed Central 

Google Scholar 
27.Vasemägi, A. The adaptive hypothesis of clinal variation revisited: Single-locus clines as a result of spatially restricted gene flow. Genetics 173, 2411–2414 (2006).PubMed 
PubMed Central 

Google Scholar 
28.Merrick, M. J., Evans, K. L. & Bertolino, S. Urban grey squirrel ecology, associated impacts, and management challenges. In Grey Squirrels: Ecology and Management of an Invasive Species in Europe (eds Shuttleworth, C. M. et al.) (European Squirrel Initiative, 2016).
Google Scholar 
29.Chipman, R., Slate, D., Rupprecht, C. & Mendoza, M. Downside risk of wildlife translocation. In Towards the Elimination of Rabies in Eurasia (eds Dodet, B. et al.) (Dev. Biol Basel, Karger, 2008).
Google Scholar 
30.Allen, D. L. Michigan Fox Squirrel Management (Michigan Department of Conservation, 1943).
Google Scholar 
31.Schorger, A. W. Squirrels in early Wisconsin. Trans. Wis. Acad. Sci. Arts Lett. 39, 195–247 (1949).
Google Scholar 
32.Robertson, G. I. Distribution of Color Morphs of Sciurus carolinensis in Eastern North America (University of Western Ontario, 1973).
Google Scholar 
33.MacCleery, D. W. American Forests: A History of Resiliency and Recovery (Forest History Society, 2011).
Google Scholar 
34.Foster, D. R. et al. Wildlands and Woodlands: A Vision for the New England Landscape (Harvard University Press, 2010).
Google Scholar 
35.Thompson, R. T., Carpenter, D. N., Cogbill, C. V. & Foster, D. R. Four centuries of change in northeastern United States forests. PLoS ONE 8(9), e72540 (2013).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
36.Lambert, M. R. et al. Adaptive evolution in cities: Progress and misconceptions. Trends Ecol. Evol. 36, 239–257 (2021).PubMed 

Google Scholar 
37.Farquhar, D. N. Some Aspects of Thermoregulation as Related to the Geographic Distribution of the Northern Melanic Phase of the Grey Squirrel (York University, 1974).
Google Scholar 
38.Innes, S. & Lavigne, D. M. Comparative energetics of coat colour polymorphs in the eastern gray squirrel Sciurus carolinensis. Can. J. Zool. 57, 585–592 (1979).
Google Scholar 
39.Santangelo, J. S. et al. Predicting the strength of urban-rural clines in a Mendelian polymorphism along a latitudinal gradient. Evol. Lett. 4, 212–225 (2020).PubMed 
PubMed Central 

Google Scholar 
40.Fidino, M. et al. Landscape-scale differences among cities alter common species’ responses to urbanization. Ecol. Appl. 31, e02253 (2021).PubMed 

Google Scholar 
41.Dickinson, J. L., Zuckerberg, B. & Bonter, D. N. Citizen science as an ecological research tool: Challenges and benefits. Annu. Rev. Ecol. Evol. Syst. 41, 149–172 (2010).
Google Scholar 
42.Alberti, M. Global urban signatures of phenotypic change in animal and plant populations. Proc. Natl. Acad. Sci. U.S.A. 114, 8951–8956 (2017).CAS 
PubMed 
PubMed Central 

Google Scholar 
43.United States Census Bureau. 2019 TIGER/Line Shapefiles (machine-readable data files) https://www2.census.gov/geo/tiger/TIGER2019/UAC/ (2019).44.XX. Statistics Canada. Population Centre Boundary File, Census year 2016 https://www150.statcan.gc.ca/n1/en/catalogue/92-166-X (2017).45.Aiello-Lammens, M. E. et al. spThin: An R package for spatial thinning of species occurrence records for use in ecological niche models. Ecography 38, 541–545 (2015).
Google Scholar 
46.R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. (2020).47.Brown de Colstoun, E. C. et al. Documentation for the Global Man-made Impervious Surface (GMIS) Dataset from Landsat (NASA Socioeconomic Data and Applications Center, 2017).
Google Scholar 
48.Steele, M. A. & Koprowski, J. L. North American Tree Squirrels (Smithsonian Books, 2001).
Google Scholar 
49.Hansen, M. C. et al. High-resolution global maps of 21st-century forest cover change. Science 342, 850–853 (2013).ADS 
CAS 
PubMed 

Google Scholar 
50.Fick, S. E. & Hijmans, R. J. WorldClim 2: New 1km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).
Google Scholar 
51.Hijmans, R. L. raster: Geographic data analysis and modeling. R package version 3.3–13. https://CRAN.R-project.org/package=raster (2020).52.Baston, D. exactextractr: Fast extraction from raster datasets using polygons. R package version 0.5.1. https://CRAN.R-project.org/package=exactextractr (2020).53.Harrison, X. A. et al. A brief introduction to mixed effects modelling and multi-model inference in ecology. PeerJ 6, e4794 (2018).PubMed 
PubMed Central 

Google Scholar 
54.Bates, D., Maechler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).
Google Scholar 
55.Gelman, A. & Su, Y. arm: Data analysis using regression and multilevel/hierarchical models. R package version 1.11–2. https://CRAN.R-project.org/package=arm (2020).56.Gelman, A. & Hill, J. Data Analysis Using Regression and Multilevel/Hierarchical Models (Cambridge University Press, 2007).
Google Scholar 
57.Crase, B., Liedloff, A. C. & Wintle, B. A. A new method for dealing with residual spatial autocorrelation in species distribution models. Ecography 35, 879–888 (2012).
Google Scholar 
58.Bivand, R. S. & Wong, D. W. S. Comparing implementations of global and local indicators of spatial association. TEST 27, 716–748 (2018).MathSciNet 
MATH 

Google Scholar 
59.Bardos, D. C., Guillera-Arroita, G. & Wintle, B. A. Valid auto-models for spatially autocorrelated occupancy and abundance data. Methods Ecol. Evol. 6, 1137–1149 (2015).
Google Scholar  More

Time-to-event models for wild animals generally model exposure of individuals to natural conditions that may affect the risk of mortality and disappearance. Most models neglect to consider seasons of high human activity that may affect such risks, or interactions between endpoint hazards (reflected in incidences) that may illuminate ecology. For many large carnivores, which suffer from low natural mortality yet are also subject to high risk of anthropogenic mortality and poaching, seasons of anthropogenic activity may be as important as natural ones in mediating cause-specific mortality and disappearance.Importantly, such anthropogenic seasons of higher mortality need not be specific to the animals being studied, especially if the species is controversial and much mortality illegal: our anthropogenic seasons consist of state hunting and hounding seasons for species other than wolves (i.e., deer or bear hunting, and hounding; not wolf hunting), but that mediate human activity on the landscape during those seasons. Our results support the hypothesis that increases in poaching risk during hunting seasons may be attributable to the surge of individuals with inclination to poach on the landscape14,18,29. Alternatively, it could also suggest enhanced criminal activity of a few poachers during the same periods. We temper this increase in poaching risk by establishing snow cover as a major environmental factor strongly associated with poaching. Moreover, our time-to-event analyses illuminate how to evaluate the effects that such anthropogenic seasons may have on risk of mortality and disappearance of monitored animals throughout their lifetime, and how considering such seasons may elucidate the mechanisms behind anthropogenic mortality and disappearance.Additionally, our analysis period precedes and completely excludes any established public wolf hunting seasons. Hence, our modeled anthropogenic seasons represent the periods of most relevant anthropogenic activity for wolves, as hypothesized by other studies14,29,33 and suggested by social science studies on inclinations to poach self-reported by both deer hunters and bear hunters, as well as acceptance of poaching by hunters and farmers30,31,32.Our analyses show increases in the hazard of disappearances of collared wolves (LTF) relative to the baseline period (which excludes environmental and anthropogenic risks) for all seasons. The highest hazard of LTF occurs during the snow season, whereas increases in hazard are lower (and similar) for the two seasons that included hounding and hunting. LTF may experience changes in hazard due to changes in the hazard of any/all of its components: migration, collar failure, or cryptic poaching.Constant and steep increases in LTF hazard throughout a wolf’s lifetime suggests mechanisms other than migration regulating LTF hazard, given migration for adults is most frequent by yearlings and younger adults, around 1.5 to 2.2 years34,35,36. Moreover, only migration out of state would end monitoring, not routine extraterritorial movements of radio-collared wolves. That our seasonal LTF curves depict the cumulative hazards more than doubling beyond those t generally associated with dispersal (~ t  More