Ecology

In this multi-scale assessment of temporal variation in bird-window collisions, our predictions related to the diel timing of collisions were only partly supported. We predicted more casualties would occur during morning than other times of day, which should have resulted in our detection of the greatest number of fatal collisions on midday surveys and more non-fatal collisions during morning surveys than midday and evening surveys. However, greatest numbers of both fatal and non-fatal collisions were observed on morning surveys, indicating that more collisions occurred during overnight and early morning periods than mid-to-late morning and afternoon. As predicted, this diel pattern was consistent across seasons. Our predictions about monthly and seasonal patterns were also only partly supported. Unexpectedly, total collision mortality was highest in the spring migration month of May, and avian residency status interacted with season such that roughly equivalent high numbers of resident and migrant collisions occurred in spring, more resident than migrant collisions occurred in summer (and overall from Apr to Oct), and more migrant than resident collisions occurred in fall. Further, unlike many past studies, collision mortality of migrants was roughly equivalent in spring and fall migration.Diel collision patternsWe observed more fatal and non-fatal collisions in the morning than in midday and evening surveys combined, even though we included fewer morning surveys in diel analyses. These differences are likely conservative in that an even greater proportion of fatal collisions than we observed likely occurred overnight and in the early morning, but went undetected because of bias associated with observer detection and scavenger removal of carcasses. Concurrent work in our study area20 found that relatively inexperienced volunteers detected a slightly lower proportion of carcasses (0.69) than experienced surveyors (0.76). Because roughly 10% of morning surveys were conducted by less-experienced volunteers and all midday or evening surveys were conducted by full-time technicians or the authors, we likely missed more carcasses on morning surveys. Additionally, most scavenging events (68%) were at night20, so bird carcasses from overnight collisions were the least likely to persist until the subsequent survey. Thus, underestimation of fatal collisions in the preceding interval was almost certainly greater for morning surveys than for midday and evening surveys.A prevailing hypothesis for why daytime bird-window collisions occur is that birds making local (e.g., foraging) movements fail to perceive a barrier when flying toward objects either on the other side of glass or reflected on a glassy surface8,15. Under this hypothesis, daytime collisions for both residents and migrating birds at stopover locations would be expected to occur most frequently when birds are most active, which is typically near dawn regardless of season. Our finding of the greatest number of collision fatalities on morning surveys circumstantially supports the above hypothesis and expectation, as do past descriptive studies of diel variation in bird-window collisions26,27,38,39. However, our study design did not allow differentiation between nocturnal and early morning collisions, and a nontrivial proportion of carcasses detected on morning surveys likely represented collisions from the preceding overnight period. Nighttime collisions may occur at any structural component not easily detectable at night (i.e., they are not limited to glass surfaces), and can be exacerbated by artificial light emission that attracts and disorients migrating birds36,37,44,45. Nonetheless, the observation of more non-fatal collisions (including directly witnessed collisions) during morning than midday or evening surveys does strongly suggest that morning carcass counts included many collisions that occurred near or shortly after dawn.A potential limitation of our study regarding time-of-day analyses is the longer interval between evening and morning surveys than between other survey periods. If collisions occurred uniformly or randomly in time, we would expect to find more bird carcasses during morning surveys due to the longer preceding time interval. However, as described above, the early morning peak observed for non-fatal collisions (Fig. 1), which are less persistent than carcasses and therefore do not accumulate over time, suggests that collisions do not have a uniform or random temporal distribution and that a real peak in collision frequency occurs sometime shortly prior to when we conducted morning surveys. Another possible bias is that we conducted surveys during fixed time periods rather than adjusting them to seasonally varying sunrise and sunset times. This limitation would have most strongly affected our summer results, potentially inflating summer morning counts as a result of collisions for both early morning and late evening (the two periods when birds are most active) being grouped into the same survey period. However, this limitation was unlikely to greatly influence our conclusions about diel collision patterns because: (1) these patterns were fairly consistent across seasons, suggesting a relatively small influence of varying sunrise and sunset times, and (2) sunrise and sunset times vary by only a few hours over the course of the entire year whereas the broad time periods for which we compared collisions (overnight, morning, afternoon) consist of longer lengths of time. Further research is needed to identify the exact timing of collisions, including during overnight periods, and this could be accomplished with collision surveys conducted more frequently during the day and night or remote detection methods, such as video cameras, motion-triggered still cameras, microphones that record sounds of impact, and glass-mounted pressure sensors that detect vibrations from collision impact.We predicted diel collision patterns would not vary seasonally because the pattern of birds being most active near sunrise and sunset is relatively consistent across seasons. Hager and Cosentino46 provide excellent guidelines for conducting bird-window collision surveys, but their recommendation to conduct surveys in mid-to-late afternoon is based on summer monitoring that found mortality to peak between late morning and early afternoon in Illinois, USA29. We suspect differences in diel patterns between that study and ours relate to geographic variation and/or seasonal sampling coverage, as our larger sample of surveys included spring and fall migration in addition to summer. Although many collision-prone species migrate nocturnally, the diel collision peak for migrants could still occur in early morning because nocturnally-migrating birds often set-down into stopover habitats during early morning47,48 and may be most susceptible to collisions at this time. There could be subtle seasonal variation in diel collision patterns that we failed to detect; however, the majority of collisions across seasons appear to occur near or before dawn (see also26,27,38,39). In combination with the previous study showing that scavenging peaks overnight20, we suggest that conducting daily collision surveys in the morning could result in the least biased mortality estimates, especially in urbanized areas where humans (e.g., cleaning crews) may remove carcasses in the early morning. As noted by Hager and Cosentino46, further research is needed to identify how the optimal survey time is influenced by factors such as geography, the bird community, animal scavengers, and removal of bird carcasses by humans.Monthly and seasonal collision patternsWe expected more collisions in fall than other seasons because bird populations in North America are larger after summer breeding and include higher proportions of young birds that have less experience with flight, migration, and human structures. Also, numerous studies have found the greatest window collision mortality in fall, a pattern driven largely by migrant birds 12,18,23,27,31,40,49. Contrary to expectation, both raw counts and bias-adjusted estimates of collision fatalities were highest in the spring migration month of May and higher overall in spring than fall. This pattern resulted from a high number of both migrant and resident bird collisions. In fact, we found that roughly equal numbers of resident and migrant collisions occurred in spring. When considering migrating individuals only, we found roughly similar numbers of collisions in spring and fall migration, a finding that also is unexpected in light of past studies. Notably, two other studies that found that a large proportion of total collisions consisted of resident birds30,50 also documented a seasonal pattern of total collisions that was less skewed toward fall. An explanation for the relatively large amount of total spring mortality, and for our finding that migrant mortality was roughly similar in spring and fall, is that some long-distance migrants follow elliptical migration paths where migration routes in fall are farther east than in spring51,52, such that in central North America, numbers of some species are greatest during spring migration. This explanation is supported by our observation of some elliptical migrants colliding during spring but not fall (e.g., Swainson’s Thrush [Catharus ustulatus]). Our study adds further nuance to the understanding of seasonal variation in bird collisions and exemplifies the need to study bird-window collisions in a wider array of geographic contexts to allow region-relevant management recommendations.Our predictions regarding avian residency status were only partly supported; more migrants than non-migrants were indeed killed during fall migration. In spring, however, roughly equal numbers of each group collided, and far more residents than migrants collided in summer (and overall from Apr–Oct). This result does not account for varying abundance (overall and seasonally) of each species group, so it does not necessarily imply that resident species are more vulnerable to, or at greater risk of, window collisions relative to their abundance or period of residency in our study area. However, the finding of a high proportion of non-migrant casualties was still unexpected given that previous similar studies have almost universally reported higher mortality among migrants15,25,26,39,49,53,54, although most sampled during migratory periods only. Even with our individual-based residency designations, we may have slightly underestimated migrant mortality because all individuals of some migratory species were classified as unknown due to overlapping resident and migratory periods. However, even if all unknown individuals were migrating, there were too few birds in this category (22 of 341 [7%] total carcasses) to change our conclusions regarding the migrant-resident comparisons. Anecdotally, many spring and summer collision fatalities were recently fledged juveniles (n = 26 [25%] in May; n = 17 [30%] in June), clearly indicating that some collision victims were indeed not migrating, and therefore, that the high number of resident collisions is not an artifact of our classification system. Moreover, we did not observe collisions of any migrant individuals during June, even though a few species migrate through our study area in small numbers during this period (e.g., shorebirds [order Charadriiformes] and some tyrant flycatchers [family Tyrannidae])55. It is possible, however, that some resident individuals were undergoing post-breeding dispersal at the time of collision, as evidenced by a small late-June peak of Tufted Titmouse (Baeolophus bicolor) and Carolina Chickadee (Poecile carolinensis) collision victims with brood patches (TJO unpublished data).Although our sampling captured the peak months of spring and fall migration in our study region, we undoubtedly missed some early-spring migrants before April and late-fall migrants after October. However, greater than 10 years of near-daily collision surveys at one of the most collision-prone buildings in our study (TJO unpublished data) suggests this number of missed collisions was relatively small. Specifically, total collisions at this building (including residents and migrants) dropped from an average of  > 8 birds in October to  More

Global Health Program, Smithsonian Conservation Biology Institute, Washington DC, USAJames M. Hassell & Dawn ZimmermanDepartment of Epidemiology of Microbial Disease, Yale School of Public Health, New Haven, CT, USAJames M. Hassell & Dawn ZimmermanCentre for Biodiversity & Environment Research (CBER), Department of Genetics, Evolution and Environment, University College London, London, UKTim Newbold & Lydia H. V. FranklinosDepartment of Ecology & Evolutionary Biology, Princeton University, Princeton, NJ, USAAndrew P. DobsonSanta Fe Institute, Santa Fe, NM, USAAndrew P. DobsonWalter Reed Biosystematics Unit (WRBU), Smithsonian Institution Museum Support Center, Suitland, MD, USAYvonne-Marie LintonDepartment of Entomology, Smithsonian National Museum of Natural History, Washington DC, USAYvonne-Marie LintonWalter Reed Army Institute of Research (WRAIR), Silver Spring, MD, USAYvonne-Marie LintonMarine Disease Ecology Laboratory, Smithsonian Environmental Research Center, Edgewater, MD, USAKatrina M. Pagenkopp Lohan More

Plant growthPrimary root length and rosette size were estimated. Control root length, measured until 21 DAG, was lightly minor in ddc vs WT (Fig. 1A). Cd differentially inhibited root growth in the two samples: at 21 DAG, 25 and 50 µM Cd-treated roots were 1.2 and 2.2 fold shorter than control roots in ddc, while in the WT Cd-treated roots were 1.8 and 2.8 fold shorter than control ones (Fig. 1A). Consequently, at 21 DAG root of Cd-treated samples was longer in ddc vs WT, particularly at the lowest Cd concentration.Figure 1(A) Primary root length (B) Picture of rosette leaf series and (C) rosette leaf area (cm2) of WT and ddc plants of A. thaliana, germinated and grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl). Root length was monitored up to 21 days after germination (DAG) every two days from germination. The results represent the mean value (± SD) of three independent biological replicates (n = 45). Asterisks indicate significant pairwise differences using Student’s t-test (*P ≤ 0.05; ** P ≤ 0.01; *** P ≤ 0.001), performed between ddc vs WT subjected to the same treatment. Bars, 0.5 cm.Full size imageRosette size was estimated at 21 DAG, corresponding to the period necessary for its full development11, by evaluating leaf number and area. Control plants of both ddc and WT exhibited a complete leaf series, although most leaves resulted smaller in ddc (Fig. 1B,C). Cd affected rosette development reducing leaf number and area, less in ddc than WT, resulting into a higher leaf area and/or number in ddc under both Cd concentrations (Fig. 1B,C).Gene expression profileRNA-Seq analysis provided an overview of gene expression profile of Cd-treated and control plants of both ddc and WT. The following comparisons were performed: ddc vs WT under control (Ctrl) conditions (ddc vs WT-Ctrl) and 25 and 50 µM Cd treatment (ddc vs WT-25 µM Cd; ddc vs WT-50 µM Cd); 25/50 µM Cd-treated vs Ctrl in ddc (25 µM Cd vs Ctrl-ddc; 50 µM Cd vs Ctrl-ddc); 25/50 µM Cd-treated vs Ctrl in the WT (25 µM Cd vs Ctrl-WT; 50 µM Cd vs Ctrl-WT).After DEGs identification (see Supplementary Fig. S1 online) 14 of them were analysed through qRT-PCR to validate transcriptomic analysis (see Supplementary Fig. S2 online). Results were fully consistent with RNA-seq data. Gene Enrichment analysis was also performed, evidencing that Cd strongly impacted on transcriptome in both ddc and WT, but in a largely different way (see Supplementary Figs. S3–S9 online). Notwithstanding, a common aspect was that in both ddc and WT the genetic pathways (GPs) more impacted by Cd dealt with photosynthesis, stress responses and hormone biosynthesis and signalling.Expression pattern of genetic pathways related to hormonesIn view of hormones pivotal role in plant development and stress response and considering the assessed epigenetic control on their action and signalling12, in this work we analysed in depth how the expression pattern (EP) of hormone-related GPs was modulated in ddc vs WT under Cd stress. The most relevant differences are discussed.AuxinsUnder control conditions, GPs related to auxin biosynthesis showed comparable EP in ddc and WT and no DEGs were detected (Fig. 2A). 25 µM Cd induced significant changes only in ddc resulting into: i) TAA1 and YUC5 downregulation along indole-3-pyruvic acid (IPA) pathway; ii) CYP71A13 and NIT2 overexpression along indole-3-acetaldoxime (IAOX) auxiliary pathway, while CYP79B3 was downexpressed (Fig. 2A). Differently, 50 µM Cd induced similar changes in ddc and WT consisting in: (i) a downexpression of YUC2 along IPA pathway in both samples and YUC5 and YUC9 in ddc and WT, respectively; (ii) overexpression of CYP71A12, CYP71A13, NIT2, NIT4 and downexpression of CYP71A16 along IAOX pathway in ddc and WT (Fig. 2A).Figure 2Genes differentially expressed (DEGs) along the pathway of (A) auxin biosynthesis, auxin conjugation, (B) indole-3-acetyl-amino acid biosynthesis, (C) methyl-indole-3-acetate interconversion and (D) auxin signalling in ddc and WT plants identified through a transcriptomic approach. For each comparison, the log2(fold change) of the analysed DEGs was shown in orange and in blue for the upregulated and downregulated genes, respectively. Plants were grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl).Full size imageAuxin level and homeostasis also depend on its oxidative degradation, conjugation and methylation13. Under control conditions, GPs related to auxin conjugation and methylation showed comparable EPs in ddc and WT and no DEGs were detected (Fig. 2B,C), but were differentially impacted by Cd, mainly at 25 µM concentration. Namely, at 25 µM Cd several genes related to auxin conjugation (GH3.3, GH3.17, YDK1) and methylation (MES7, MES17) were downregulated in ddc, whereas in WT only MES18, involved in methyl-indole-3-acetate production, was downregulated (Fig. 2B,C). At 25 µM Cd most of the above genes were downregulated in ddc vs WT (Fig. 2B,C), while at 50 µM Cd only two genes, working in auxin methylation were differentially modulated in ddc and WT (Fig. 2B,C).Under control conditions, GP related to auxin signalling exhibited similar EP in ddc and WT (Fig. 2D). Differences were induced by Cd. Namely, at 25 µM Cd, AUX/IAA family genes, which acts in signalling repression, were globally downregulated in ddc, while in WT the repressor IAA34 was overexpressed (Fig. 2D). Differently, 50 µM Cd effects on ddc and WT were quite similar, dealing with AUX/IAA family genes downregulation (Fig. 2D). In ddc vs WT comparisons, at 25 µM Cd above genes were downregulated, while no differences were found at 50 µM Cd (Fig. 2D). Moreover, and somehow unexpectedly, following 50 µM Cd it was observed in both ddc and WT a downregulation of several SAURs members, belonging to a large family of auxin responsive genes, which in turn can also have an impact on auxin pathway14 (Fig. 2D). Interestingly, such effect was more pronounced in ddc than WT. However, it must be mentioned that, although most of them are induced by auxin, several other hormones and co-factors acts upstream SAUR genes, regulating their activity in response to both endogenous stimuli and environmental cues14.In summary, in both ddc and WT, Cd induced: i) a downregulation of IPA pathway, which is the main auxin biosynthetic pathway15 and a simultaneous upregulation of IAOX auxiliary biosynthetic pathway; ii) an enhancement of hormone signalling. However, in the WT such effects occurred only at 50 μM Cd. Moreover, in ddc Cd also induced a downregulation of GPs related to auxin conjugation.CytokininsIn all comparisons, GPs related to CKs biosynthesis showed similar EPs, unless for the downregulation in 50 µM Cd-treated WT vs Ctrl of IPT5, encoding rate-limiting enzyme along the pathway16 (Fig. 3A).Figure 3Genes differentially expressed (DEGs) along the pathway of (A) trans-zeatin biosynthesis, (B) CKs degradation, (C) CKs N7- and N9-glucoside biosynthesis, (D) CKs O-glycosylation and (E) CKs signalling in ddc and WT plants identified through a transcriptomic approach. For each comparison, the log2(fold change) of the analysed DEGs was shown in orange and in blue for the upregulated and downregulated genes, respectively. Plants were grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl).Full size imageMajor differences were observed for GPs related to CKs catabolism and conjugation, occurring through cleavage by oxidation and glycosylation, respectively16. Under control conditions, these GPs also exhibited similar EPs in ddc and WT (Fig. 3B-D). At both 25 and 50 µM Cd a downregulation of CKX5 and CKX6, encoding cytokinin-oxidases, occurred only in ddc (Fig. 3B). Differently, Cd impact on GPs related to CKs N-glycosylation was almost comparable in ddc and WT, resulting into the overexpression of two different genes working in N7- and N9-glycosylation pathways at 25 µM Cd, and one gene at the higher concentration (Fig. 3C).Note that at 25 µM Cd the above genes were both overexpressed in ddc vs WT, while no differences were found at 50 µM Cd (Fig. 3C). Cd impact on GP related to cytokinin O-glycosylation was major, especially in ddc, involving at 25 µM Cd the overexpression of seven genes along this pathway compared to Ctrl (Fig. 3D), and one gene in WT 25 µM Cd vs Ctrl (Fig. 3D). By contrast, at 50 μM Cd the expression pattern along this pathway was similar in ddc and WT, being characterised by the upregulation of the same seven genes above mentioned and the downregulation of AT5G38010 (Fig. 3D). Finally, at 25 μM Cd, five genes along these pathways resulted upregulated in ddc vs WT while a similar EP occurred in ddc and WT under 50 μM Cd treatment (Fig. 3D).Concerning GP involved in CKs signalling, under control conditions A-ARRs, encoding negative regulators of CKs signalling17, were downregulated and signalling was likely enhanced in ddc vs WT (Fig. 3E). 25 µM Cd induced a downregulation of ARR11 A-type ARRs and B-ARR family ARR10 transcription factors, which control primary plant response to CKs, only in ddc (Fig. 3E). Whereas, at 50 µM Cd both ddc and WT showed A-ARR downregulation, supposedly leading to pathway upregulation (Fig. 3E). At 25 µM Cd AHP1, encoding positive regulators of CKs signalling18, was downregulated in ddc vs WT, while ARR17 was overexpressed, suggesting that signalling was downregulated also in ddc vs WT (Fig. 3E). No differences occurred between ddc and WT at 50 µM Cd (Fig. 3E).In summary, transcriptomic analysis evidenced that GP related to the biosynthesis of trans-zeatin, the most relevant CK, was negatively affected by Cd only in the WT at 50 μM Cd. In response to Cd, GPs related to CKs inactivation were enhanced in both ddc and WT, but in ddc a downregulation of GP related to CKs cleavage also occurred. Finally, hormone signalling was differentially modulated by Cd in relation to both the sample (ddc vs WT) and heavy metal concentration, resulting into a downregulation at 25 µM Cd only in ddc and an enhancement in ddc and WT at 50 µM Cd.GibberellinsUnder control conditions, GPs related GAs biosynthesis showed similar EPs in ddc vs WT (Fig. 4A). 25 μM Cd induced in ddc: i) a downregulation of GA2 encoding the ent-kaurene synthase, a pivotal enzyme along the early GAs biosynthetic pathways to synthetize GA12; ii) a downregulation of GA4, a key gene of GAs biosynthesis, along which bioactive GAs are synthetized19 (Fig. 4A). No Cd-induced modulation was observed in the WT (Fig. 4A). On the contrary, at 50 μM Cd both ddc and WT showed a downregulation of GA5 (Fig. 4A). Finally, in ddc vs WT the only difference dealt with GA2 downregulation at 25 µM Cd (Fig. 4A).Figure 4Genes differentially expressed (DEGs) along the pathway of (A) GAs biosynthesis, (B) GAs inactivation and (C) GAs signalling in ddc and WT plants identified through a transcriptomic approach. For each comparison, the log2(fold change) of the analysed DEGs was shown in orange and in blue for the upregulated and downregulated genes, respectively. Plants were grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl).Full size imageGPs controlling GAs inactivation also showed a comparable transcriptional pattern in ddc and WT under control conditions, and no DEGs were detected (Fig. 4B). 25 μM Cd induced a downregulation of DAO2 and AOP1, encoding GA2ox enzymes, only in ddc (Fig. 4B). Accordingly, in ddc vs WT these genes were downregulated only at the lowest Cd concentration (Fig. 4B).Under control conditions, also GP related to GAs signalling was not differentially modulated in ddc vs WT (Fig. 4C). A downregulation of genes encoding DELLA proteins, which act as repressors20, was induced only in ddc by 25 µM Cd (Fig. 4C) and in both ddc and WT at 50 µM Cd (Fig. 4C), suggesting an enhancement of hormone signalling. Consequently, DELLA-codifying genes resulted downregulated also in ddc vs WT only at the lowest Cd concentration (Fig. 4C).In summary, in ddc the lowest Cd treatment negatively affected GPs related to GAs biosynthesis but, at the same time, hormone signalling resulted enhanced. In the WT similar effects were observed only at 50 µM Cd.Jasmonic acidUnder control conditions, six genes along the GP related to JA biosynthesis were downregulated in ddc vs WT (Fig. 5A). 25 µM Cd induced in ddc a downregulation of this GP, except for LOX4 upregulation (Fig. 5A) and a downregulation involving eight genes in WT (Fig. 5A). At 50 μM, Cd effects were limited to LOX5 downregulation and LOX4 and OPR1 upregulation in ddc and LOX5 and AOS downregulation in WT (Fig. 5A). No Cd–induced differences were found in ddc vs WT (Fig. 5A).Figure 5Genes differentially expressed (DEGs) along the pathway of (A) JA biosynthesis, (B) JA signalling, (C) ABA biosynthesis, (D) ABA degradation, (E) ABA glucose ester biosynthesis and (F) ABA signalling in ddc and WT plants identified through a transcriptomic approach. For each comparison, the log2(fold change) of the analysed DEGs was shown in orange and in blue for the upregulated and downregulated genes, respectively. Plants were grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl).Full size imageConcerning the JA signalling-related GP, in control conditions JAZ5 gene, encoding a protein acting as repressor21, was downregulated in ddc vs WT, highlighting a signalling enhancement (Fig. 5B). Interestingly, JAZ10 and JAZ9 were differentially impacted by 25 µM Cd in ddc and WT resulting upregulated and downregulated, respectively. (Fig. 5B). At 50 µM Cd, JAZs were overexpressed in both ddc and WT (Fig. 5B). No differences were detected in ddc vs WT exposed to Cd (Fig. 5B).Globally, Cd negatively impacted on the GP related to JA biosynthesis especially in WT. Under Cd treatment hormone signalling was downregulated more in ddc than in WT, whatever concentration was applied.Abscisic acidUnder control conditions, ABA biosynthesis-related GP showed comparable EP in ddc and WT, and no DEGs were detected (Fig. 5C). The only significant Cd effect dealt with NCED3 downregulation both in ddc and WT, regardless of applied concentration (Fig. 5C). Under control conditions, also the GPs related to ABA catabolism showed a comparable EP in ddc and WT and no DEGs were detected (Fig. 5D), but Cd differentially impacted on CYP genes, involved in phaseic acid degradative production22. Namely, at 25 μM Cd, CYP707A3 was downregulated only in ddc (Fig. 5D). Moreover, also CYP707A2 appeared downregulated in ddc vs WT (Fig. 5D). At 50 µM Cd it was observed an upregulation of both CYP707A2 and CYP707A4 in ddc, and of only CYP707A4 in WT (Fig. 5D).Under control conditions, GP related to ABA inactivation through glucose conjugation showed similar EP in ddc vs WT (Fig. 5E). 25 μM Cd determined AT4G15260 upregulation and UGT71C3 downregulation only in ddc (Fig. 5E). At 50 μM, Cd equally impacted on ddc and WT, resulting into UGT71C1 and UGT2 downregulation, AT5G49690, UGT71B5, UGT71B6 upregulation and, limited to WT, AT4G15260 upregulation (Fig. 5E). No differences were highlighted in ddc vs WT exposed to Cd (Fig. 5E).Under control conditions, the GP related to hormone signalling also presented a comparable EP in ddc and WT (Fig. 5F). In ddc, 25 µM Cd impact on this GP appeared rather complex, resulting in an upregulation of PYL3, encoding ABA receptor, and a downregulation of PP2Cs (PP2CA and HAI1) encoding negative regulators of ABA signalling23. Moreover, ABI5, codifying a key transcription factor in ABA signalling24 belonging to AREBs/ABFs family, was upregulated. However, SnRK2.7 gene, codifying a protein which activate the AREBs/ABFs transcription factors24, was downregulated. Based on the prominent role of SnRK2s in plant response to ABA, it is likely that at 25 µM Cd ABA signalling was downregulated in ddc (Fig. 5F). Instead, 25 µM Cd determined in WT the upregulation of PYL6 and the downregulation of PP2Cs, suggesting an enhancement of ABA signalling (Fig. 5F). At 50 µM Cd, PYL6 was upregulated and SnRK2.7 downregulated in both ddc and WT, while HAI1 was downregulated only in ddc (Fig. 5F). When comparing ddc vs WT, at 25 µM Cd only SnRK2.7 was downregulated, while no differences occurred at 50 μM Cd (Fig. 5F).In summary, Cd determined a slight downregulation of GP related to ABA biosynthesis in both ddc and WT regardless of its concentration. ABA catabolic pathway was lightly downregulated in ddc at 25 µM Cd but upregulated in both samples at 50 μM Cd. At the transcriptomic level, ABA signalling featured as enhanced in WT and downregulated in ddc regardless of Cd concentration.EthyleneAlong GP related to ethylene biosynthesis, under control conditions ACS8 and ACS11 were upregulated in ddc vs WT (Fig. 6A). 25 µM Cd determined ACS8 and ACO5 downregulation in ddc and ACS4 upregulation in the WT (Fig. 6A). 50 μM Cd induced ACS7 upregulation and ACO5 downregulation in both ddc and WT and ACS2 and ACS11 overexpression only in WT (Fig. 6A). Finally, the only Cd-induced difference in ddc vs WT dealt with ACO1 downregulation at 25 µM Cd (Fig. 6A).Figure 6Genes differentially expressed (DEGs) along the pathway of (A) ethylene biosynthesis, (B) ethylene signalling and (C) SA signalling in ddc and WT plants identified through a transcriptomic approach. For each comparison, the log2(fold change) of the analysed DEGs was shown in orange and in blue for the upregulated and downregulated genes, respectively. Plants were grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl).Full size imageGP related to ethylene signalling showed similar EP in ddc and WT both under control conditions (Fig. 6B) and at 25 µM Cd, except for the upregulation of ETR2, encoding ethylene receptor25 in WT (Fig. 6B). At 50 µM Cd, both ddc and WT exhibited ETR2 and ERF1 overexpression suggesting an upregulation of ethylene signalling (Fig. 6B); no differences occurred in ddc vs WT (Fig. 6B).In summary, in control conditions GP related to ethylene biosynthesis was upregulated in ddc vs WT. Cd determined a downregulation and upregulation of this GP in ddc and WT, respectively. Concerning hormone signalling, at the highest Cd concentration in both ddc and WT an upregulation of this GP occurred.Salicylic acidRegarding SA, only the GP related to signalling resulted differentially expressed. Under control conditions, the GP related to SA signalling showed similar EP in ddc and WT. However, PR1, a useful molecular marker for the systemic acquired resistance (SAR) in response to pathogens26, was downregulated in ddc vs WT (Fig. 6C). 25 μM Cd induced a downregulation of genes codifying TGA10 transcription factor only in the ddc (Fig. 6C). Whereas, 50 μM Cd induced a downregulation of PRB1 in both ddc and WT and of TGA8 only in ddc (Fig. 6C). In ddc vs WT, differences were found only at 25 µM Cd, with the downregulation of TGA10 and CAPE3 (Fig. 6C).Altogether, these results evidenced a Cd-induced downregulation of this GP, likely resulting in an impairment of hormone signalling in both WT and ddc, but in the latter this effect already occurred at the lowest Cd concentration.Phytohormone levelBased on the major effects induced by 25 µM Cd treatment, hormone quantification was carried out on plants exposed to this concentration, compared to untreated control plants.Under control conditions, IAA amount was higher in ddc than WT, although not significantly. After Cd treatment, a decreasing trend was observed only in WT, resulting into a significant lower level as compared to ddc (Fig. 7A).Figure 7(A) IAA, (B–G) CKs, (H–R) GAs, (S) JA, (T) ABA, (U) SA and (V) SAG amount in A. thaliana ddc mutant and WT plants grown in Ctrl conditions and treated with 25 µM Cd estimated by GC–MS. The results represent the mean value (± SD) of three independent biological replicates. Statistical analysis was performed by using two-way ANOVA with Tukey post hoc test (P ≤ 0.05) after Shapiro–Wilk normality test. Means with the same letter are not significantly different at P ≤ 0.05.Full size imageConcerning CKs, both biological active (tZ) and inactive conjugate (tZR, cZR, tZOG, cZOG, iPR) forms were analysed (Fig. 7B–G). Under control conditions, all analysed CKs were present in ddc, but CKs conjugate forms and above all O-glycosylated exhibited the highest levels (Fig. 7B–G). By contrast, in WT tZ was not detectable and all the other CKs forms exhibited a lower level compared to the mutant, which appeared significant for tZR and cZOG (Fig. 7B–G). Following Cd treatment, CKs levels increased in ddc, except for tZ decrease. In WT, the unique Cd effect dealt with tZR increase and tZOG decrease. Consequently, under Cd treatment the level of all CKs forms remained higher (from 0.25 to 3 times) in ddc than in WT (Fig. 7B–G).Concerning GAs, precursors (GA9, GA19, GA20), biologically active forms (GA1, GA3, GA4, GA7) and catabolites (GA8, GA34, GA29, GA51) were analysed (Fig. 7H–R). Under control conditions, GA19 amount was significantly higher in ddc vs WT, while the amount of GA20, the other in serie precursor of hydroxylated forms, was comparable between the samples (Fig. 7H,I). Following Cd treatment, GA19 amount significantly increased in WT while a slight downtrend occurred in ddc, leading to comparable values in the two samples. The same trend was observed, but at less extent, for GA20 (Fig. 7H,I). Consistently, under control conditions also the amount of the active hydroxylated forms GA1 and GA3 was higher in WT than in ddc (Fig. 7J,K). Following Cd treatment, a decrease of their amount was detected only in WT, globally leading to a higher level of these GAs in ddc mutant compared to the WT (Fig. 7J,K). In addition, in ddc mutant also the related catabolites GA8 and GA29 were globally lower than in WT, under both control conditions and Cd treatment (Fig. 7L,M).Differences were observed also for GA9, precursor of non-hydroxylated GAs: under Cd treatment its amount decreased in the WT and was instead induced in ddc mutant, resulting in a quite comparable value between the two samples (Fig. 7N). Consistently, the amount of active non-hydroxylated forms, GA4 and GA7, increased under Cd treatment only in ddc mutant; also in this case, at the end of heavy metal treatment, comparable values were detected in ddc and WT (Fig. 7O,P) In agreement with these results, following Cd treatment, the amount of catabolites GA51 and GA34 did not change in the WT, whereas in ddc it increased and decreased, respectively (Fig. 7Q,R).As evident in Fig. 7S–V, differences were reported also for JA, ABA, SA and its predominant inactive conjugate, SA 2-O-β-D-glucoside (SAG). Under control conditions both JA and ABA amount was significantly lower in ddc vs WT and significantly decreased following Cd treatment only in the WT (Fig. 7S,T). Notwithstanding, under such condition the ABA amount remained lower in ddc than in WT while JA values were comparable in the two samples due its light, but not significant, increase in ddc (Fig. 7S,T). By contrast, under control conditions both SA and SAG amounts were significantly higher in ddc than in WT (Fig. 7U,V). Following Cd treatment, their amounts significantly decreased more in ddc than in WT, leading to an opposite condition (Fig. 7U,V).Testing of the involvement of SUPPRESSOR OF DRM1 DRM2 CMT3 (SDC) gene in ddc response to CdFinally, we planned to inquire on the possible involvement of SDC gene in the response of ddc triple mutant to Cd exposure. Indeed, it has been reported that in ddc mutant the misexpression of such gene, which encodes a F-Box protein, is ultimately responsible of the developmental phenotypes of ddc, such as curled leaves and reduced growth, as evidenced by its reversion in the drm1 drm2 cmt3 sdc quadruple mutant27. Note that in the WT SDC is silenced, being methylated in all its sequence contexts because of the redundant action of DRM2 and CMT3 enzymes. By contrast, in ddc, where DRM2 and CMT3 expression is silenced, the loss of non-CG methylation in the promoter region of SDC F-box gene determines its overexpression27.According to the above mentioned data27, we firstly verified that under control conditions SDC resulted silent in the WT and overexpressed in ddc also in our transcriptomic analysis (data not shown, complete raw transcriptomic data are available at NCBI SRA under the BioProject accessionPRJNA641242;https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA641242). Moreover, our data also showed that at the transcriptomic level SDC is not modulated by Cd since its expression level did not significantly change in ddc nor in WT whatever heavy metal concentration was applied.Thereafter, we tested the involvement of SDC in the growth response of ddc mutant under Cd exposure, by monitoring primary root length of Arabidopsis thaliana WT, ddc and sdc plants grown under the following conditions: (i) on growth medium without Cd as control (Ctrl) (ii) on a medium supplemented with 25/50 µM Cd; (iii) limited to the WT and sdc mutant plants, on a medium supplemented with 25/50 µM Cd plus 15 µM 5-Azacytidine (5-Aza), an inhibitor of DNA methylation applied in order to mimic the hypomethylated state of ddc mutant.Under control conditions, all three samples showed a similar root length. However, at 21 DAG, root was lightly shorter in ddc vs WT, while sdc displayed an intermediate length (Fig. 8A). Under both Cd treatments, roots were averagely longer in ddc than in WT. Again, sdc roots exhibited an intermediate length, more similar to WT than ddc (Fig. 8 B,C). Interestingly, WT and sdc plants treated with Cd plus 5-Aza had longer roots than the plants treated only with Cd, and quite comparable to ddc roots exposed to Cd (Fig. 8 D,E).Figure 8Primary root length of WT, sdc and ddc plants of A. thaliana, germinated and grown in long day condition (A) on growth medium without Cd as control (Ctrl), (B,C) on a medium supplemented with 25/50 µM Cd, (D,E) limited to WT and sdc plants, on a medium supplemented with 25/50 µM Cd plus 15 µM 5-Azacytidine (5-Aza). Root length was monitored up to 21 days after germination (DAG) every two days from germination. The results represent the mean value (± SD) of three independent biological replicates (n = 45). Statistical analysis was performed between samples at the same growth stage, by using two-way ANOVA with Tukey post hoc test (P ≤ 0.05) after Shapiro–Wilk normality test. Means with the same letter are not significantly different at P ≤ 0.05.Full size image More

Genotyping control and datasets creationsAfter filtering the initial raw dataset of 1,071 goats for poor quality genotype and related animals we obtained the dataset used for the haplotype sharing analysis, which included 980 animals and 48,396 SNPs. This dataset was further pruned for linkage disequilibrium (LD; r2  More

1.Maslo, B. & Fefferman, N. H. A case study of bats and white-nose syndrome demonstrating how to model population viability with evolutionary effects. Conserv. Biol. 29, 1176–1185 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
2.Morris, W. F. & Doak, D. F. Quantitative Conservation Biology (Sinauer, Sunderland, 2002).
Google Scholar 
3.Liebhold, A. & Bascompte, J. The Allee effect, stochastic dynamics and the eradication of alien species. Ecol. Lett. 6, 133–140 (2003).Article 

Google Scholar 
4.Stephens, P. A. & Sutherland, W. J. Consequences of the Allee effect for behaviour, ecology and conservation. Trends Ecol. Evol. 14, 401–405 (1999).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
5.Nunney, L. & Elam, D. R. Estimating the effective population size of conserved populations. Conserv. Biol. 8, 175–184 (1994).Article 

Google Scholar 
6.Lande, R. & Barrowclough, G. Effective population size, genetic variation, and their use in population. Viable populations for conservation, 87 (1987).7.Frankham, R. Effective population size/adult population size ratios in wildlife: A review. Genet. Res. 66, 95–107 (1995).Article 

Google Scholar 
8.Tallmon, D. A., Luikart, G. & Waples, R. S. The alluring simplicity and complex reality of genetic rescue. Trends Ecol. Evol. 19, 489–496 (2004).PubMed 
Article 
PubMed Central 

Google Scholar 
9.Whiteley, A. R., Fitzpatrick, S. W., Funk, W. C. & Tallmon, D. A. Genetic rescue to the rescue. Trends Ecol. Evol. 30, 42–49 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
10.Jiao, J., Gilchrist, M. A. & Fefferman, N. H. The impact of host metapopulation structure on short-term evolutionary rescue in the face of a novel pathogenic threat. Glob. Ecol. Conserv. 23, 01174 (2020).
Google Scholar 
11.Hanski, I. Metapopulation Ecology (Oxford University Press, Oxford, 1999).
Google Scholar 
12.Mortier, F., Jacob, S., Vandegehuchte, M. L. & Bonte, D. Habitat choice stabilizes metapopulation dynamics by enabling ecological specialization. Oikos 128, 529–539 (2019).Article 

Google Scholar 
13.Jiao, J., Riotte-Lambert, L., Pilyugin, S. S., Gil, M. A. & Osenberg, C. W. Mobility and its sensitivity to fitness differences determine consumer–resource distributions. R. Soc. Open Sci. 7, 200247 (2020).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
14.Anderson, S. C., Moore, J. W., McClure, M. M., Dulvy, N. K. & Cooper, A. B. Portfolio conservation of metapopulations under climate change. Ecol. Appl. 25, 559–572 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
15.Case, T. J. Invasion resistance, species build-up and community collapse in metapopulation models with interspecies competition. Biol. J. Lin. Soc. 42, 239–266 (1991).Article 

Google Scholar 
16.Gyllenberg, M. & Hanski, I. Habitat deterioration, habitat destruction, and metapopulation persistence in a heterogenous landscape. Theor. Popul. Biol. 52, 198–215 (1997).CAS 
PubMed 
MATH 
Article 
PubMed Central 

Google Scholar 
17.Jiao, J., Pilyugin, S. S. & Osenberg, C. W. Random movement of predators can eliminate trophic cascades in marine protected areas. Ecosphere 7, e01421 (2016).Article 

Google Scholar 
18.Nee, S. & May, R. M. Dynamics of metapopulations: Habitat destruction and competitive coexistence. J. Anim. Ecol. 61, 37–40 (1992).Article 

Google Scholar 
19.Ying, Y., Chen, Y., Lin, L. & Gao, T. Risks of ignoring fish population spatial structure in fisheries management. Can. J. Fish. Aquat. Sci. 68, 2101–2120 (2011).Article 

Google Scholar 
20.Hess, G. Disease in metapopulation models: Implications for conservation. Ecology 77, 1617–1632 (1996).Article 

Google Scholar 
21.Daszak, P., Cunningham, A. A. & Hyatt, A. D. Emerging infectious diseases of wildlife: Threats to biodiversity and human health. Science 287, 443–449 (2000).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
22.Harding, K. C., Begon, M., Eriksson, A. & Wennberg, B. Increased migration in host–pathogen metapopulations can cause host extinction. J. Theor. Biol. 298, 1–7 (2012).MathSciNet 
CAS 
PubMed 
MATH 
Article 
PubMed Central 

Google Scholar 
23.Dowling, A. J., Hill, G. E. & Bonneaud, C. Multiple differences in pathogen-host cell interactions following a bacterial host shift. Sci. Rep. 10, 1–12 (2020).Article 
CAS 

Google Scholar 
24.Kuzmin, I. V. et al. Molecular inferences suggest multiple host shifts of rabies viruses from bats to mesocarnivores in Arizona during 2001–2009. PLoS Pathog 8, e1002786 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
25.Levine, R. S. et al. Supersuppression: Reservoir competency and timing of mosquito host shifts combine to reduce spillover of West Nile virus. Am. J. Trop. Med. Hyg. 95, 1174–1184 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
26.Langwig, K. E. et al. Context-dependent conservation responses to emerging wildlife diseases. Front. Ecol. Environ. 13, 195–202 (2015).Article 

Google Scholar 
27.Smith, K. F., Acevedo-Whitehouse, K. & Pedersen, A. B. The role of infectious diseases in biological conservation. Anim. Conserv. 12, 1–12 (2009).Article 

Google Scholar 
28.Xiao, Y., Tang, B., Wu, J., Cheke, R. A. & Tang, S. Linking key intervention timing to rapid decline of the COVID-19 effective reproductive number to quantify lessons from mainland China. Int. J. Infect. Dis. 97, 296–298 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
29.Cintrón-Arias, A., Castillo-Chávez, C., Betencourt, L., Lloyd, A. L. & Banks, H. T. The Estimation of the Effective Reproductive Number from Disease Outbreak Data. (North Carolina State University, Center for Research in Scientific Computation, 2008).30.Salpeter, E. E. & Salpeter, S. R. Mathematical model for the epidemiology of tuberculosis, with estimates of the reproductive number and infection-delay function. Am. J. Epidemiol. 147, 398–406 (1998).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
31.Grenfell, B. & Harwood, J. (Meta) population dynamics of infectious diseases. Trends Ecol. Evol. 12, 395–399 (1997).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
32.Millard, A. R., Roberts, C. A. & Hughes, S. S. Isotopic evidence for migration in Medieval England: The potential for tracking the introduction of disease. Soc. Biol. Human Affairs. 70, 9–13 (2005).
Google Scholar 
33.Chen, M. et al. The introduction of population migration to SEIAR for COVID-19 epidemic modeling with an efficient intervention strategy. Inf. Fusion 64, 252–258 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
34.Reed, K. D., Meece, J. K., Henkel, J. S. & Shukla, S. K. Birds, migration and emerging zoonoses: West Nile virus, Lyme disease, influenza A and enteropathogens. Clin. Med. Res. 1, 5–12 (2003).PubMed 
PubMed Central 
Article 

Google Scholar 
35.Roy, B. & Kirchner, J. Evolutionary dynamics of pathogen resistance and tolerance. Evolution 54, 51–63 (2000).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
36.Bliven, K. A. & Maurelli, A. T. Antivirulence genes: Insights into pathogen evolution through gene loss. Infect. Immun. 80, 4061–4070 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
37.Wild, G., Gardner, A. & West, S. A. Adaptation and the evolution of parasite virulence in a connected world. Nature 459, 983–986 (2009).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
38.Laine, A. L. Resistance variation within and among host populations in a plant–pathogen metapopulation: Implications for regional pathogen dynamics. J. Ecol. 92, 990–1000 (2004).Article 

Google Scholar 
39.Thrall, P. H. et al. Rapid genetic change underpins antagonistic coevolution in a natural host-pathogen metapopulation. Ecol. Lett. 15, 425–435 (2012).PubMed 
PubMed Central 
Article 

Google Scholar 
40.Juhas, M. Horizontal gene transfer in human pathogens. Crit. Rev. Microbiol. 41, 101–108 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
41.Soanes, D. & Richards, T. A. Horizontal gene transfer in eukaryotic plant pathogens. Annu. Rev. Phytopathol. 52, 583–614 (2014).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
42.Brunham, R. C., Plummer, F. A. & Stephens, R. S. Bacterial antigenic variation, host immune response, and pathogen-host coevolution. Infect. Immun. 61, 2273 (1993).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
43.Sasaki, A. Evolution of antigen drift/switching: Continuously evading pathogens. J. Theor. Biol. 168, 291–308 (1994).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
44.Lange, A. & Ferguson, N. M. Antigenic diversity, transmission mechanisms, and the evolution of pathogens. PLoS Comput. Biol. 5, 1000536 (2009).ADS 
MathSciNet 
Article 
CAS 

Google Scholar 
45.Alizon, S., Hurford, A., Mideo, N. & Van Baalen, M. Virulence evolution and the trade-off hypothesis: History, current state of affairs and the future. J. Evol. Biol. 22, 245–259 (2009).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
46.Messenger, S. L., Molineux, I. J. & Bull, J. Virulence evolution in a virus obeys a trade off. Proc. R. Soc. Lond. B 266, 397–404 (1999).CAS 
Article 

Google Scholar 
47.Alizon, S., de Roode, J. C. & Michalakis, Y. Multiple infections and the evolution of virulence. Ecol. Lett. 16, 556–567 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
48.Bull, J. J. & Lauring, A. S. Theory and empiricism in virulence evolution. PLoS Pathog 10, e1004387 (2014).49.Gray, M. J. & Chinchar, V. G. Ranaviruses: Lethal Pathogens of Ectothermic Vertebrates (Springer Science+ Business Media, New York, 2015).Book 

Google Scholar 
50.Dobbelaere, T., Muller, E. M., Gramer, L. J., Holstein, D. M. & Hanert, E. Coupled epidemio-hydrodynamic modeling to understand the spread of a deadly coral disease in Florida. Front. Mar. Sci. 7, 1016 (2020).Article 

Google Scholar 
51.Stoddard, S. T. et al. House-to-house human movement drives dengue virus transmission. Proc. Natl. Acad. Sci. USA 110, 994–999 (2013).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
52.Diekmann, O., Heesterbeek, J. & Roberts, M. G. The construction of next-generation matrices for compartmental epidemic models. J. R. Soc. Interface 7, 873–885 (2010).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
53.Smith, K. F., Sax, D. F. & Lafferty, K. D. Evidence for the role of infectious disease in species extinction and endangerment. Conserv. Biol. 20, 1349–1357 (2006).PubMed 
Article 
PubMed Central 

Google Scholar 
54.Anderson, R. M., Anderson, B. & May, R. M. Infectious Diseases of Humans: Dynamics and Control (Oxford University Press, Oxford, 1992).
Google Scholar 
55.O’Brien, S. J. et al. Genetic basis for species vulnerability in the cheetah. Science 227, 1428–1434 (1985).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
56.Ingvarsson, P. K. & Lundberg, S. The effect of a vector-borne disease on the dynamics of natural plant populations: A model for Ustilago violacea infection of Lychnis viscaria. J. Ecol. 81, 263–270 (1993).Article 

Google Scholar 
57.Carlson, S. M., Cunningham, C. J. & Westley, P. A. Evolutionary rescue in a changing world. Trends Ecol. Evol. 29, 521–530 (2014).PubMed 
Article 
PubMed Central 

Google Scholar 
58.Gonzalez, A., Ronce, O., Ferriere, R. & Hochberg, M. E. (The Royal Society, 2013).59.Fine, P. E. Herd immunity: History, theory, practice. Epidemiol. Rev. 15, 265–302 (1993).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
60.Fontanet, A. & Cauchemez, S. COVID-19 herd immunity: Where are we?. Nat. Rev. Immunol. 20, 583–584 (2020).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
61.Fine, P., Eames, K. & Heymann, D. L. “Herd immunity”: A rough guide. Clin. Infect. Dis. 52, 911–916 (2011).PubMed 
Article 
PubMed Central 

Google Scholar 
62.Barbarossa, M. V. & Röst, G. Immuno-epidemiology of a population structured by immune status: A mathematical study of waning immunity and immune system boosting. J. Math. Biol. 71, 1737–1770 (2015).MathSciNet 
CAS 
PubMed 
MATH 
Article 
PubMed Central 

Google Scholar 
63.Hamami, D., Cameron, R., Pollock, K. G. & Shankland, C. Waning immunity is associated with periodic large outbreaks of mumps: A mathematical modeling study of Scottish data. Front. Physiol. 8, 233 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
64.Klepac, P. & Caswell, H. The stage-structured epidemic: Linking disease and demography with a multi-state matrix approach model. Thyroid Res. 4, 301–319 (2011).
Google Scholar 
65.Anderson, R. M. & May, R. M. Population biology of infectious diseases: Part I. Nature 280, 361–367 (1979).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
66.Colizza, V. & Vespignani, A. Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: Theory and simulations. J. Theor. Biol. 251, 450–467 (2008).MathSciNet 
PubMed 
MATH 
Article 
PubMed Central 

Google Scholar 
67.Stone, C. M., Schwab, S. R., Fonseca, D. M. & Fefferman, N. H. Human movement, cooperation and the effectiveness of coordinated vector control strategies. J. R. Soc. Interface 14, 20170336 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
68.Jiao, J., Pilyugin, S. S., Riotte-Lambert, L. & Osenberg, C. W. Habitat-dependent movement rate can determine the efficacy of marine protected areas. Ecology 99, 2485–2495 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
69.Keeling, M. J., Rohani, P. & Grenfell, B. T. Seasonally forced disease dynamics explored as switching between attractors. Physica D 148, 317–335 (2001).ADS 
MATH 
Article 

Google Scholar  More

In March, Brazil’s Ministry of Agriculture took an alarming step to boost trade of artisanal animal products across states (see go.nature.com/3by9). It added reptiles and amphibians — already the most threatened vertebrates on Earth — to the list permitting the capture of fishes, crustaceans and molluscs for human consumption.Brazil has the fastest rate of decline of amphibian populations in South America, owing to habitat loss and infectious diseases (B. C. Scheele et al. Science 363, 1459–1463; 2019). If the policy takes effect in its current form, trade of amphibians will increase — compounding the spread of lethal pathogens such as Batrachochytrium species and ranavirus.We urge the government to align its policy with the Convention on Biological Diversity and other international commitments that are backed by substantial scientific evidence. More

Soil properties in different rice farming systemsFive treatments were designed in the three selected rice fields, including (1) rice monoculture field (RM); (2) planting area in the 1st year of rice–fish field (OP); (3) aquaculture area in the 1st year of rice–fish field (OA); (4) planting area in the 5th year of rice–fish field (FP); (5) aquaculture area in the 5th year of rice–fish field (FA). The soil properties of the five treatments were shown in Table 1. The highest soil available nitrogen (AN) content was observed in FP and was significantly higher than that in RM, OP and OA. The highest soil available phosphorus (AP) content was observed in RM and was significantly higher than that in the other 4 treatments. The highest soil available potassium (AK) content was measured in the 1st year of rice–fish field (OP and OA), followed by RM and the 5th year of rice–fish field (FP and FA), and significant differences were observed among different rice fields. The highest soil organic matter (OM) content appeared in the 5th year of rice–fish field (FP and FA), and was only significantly higher than that in OA. In addition, the soil pH in the 1st year of rice–fish field (OP and OA) was significantly lower than that in RM and the 5th year of rice–fish field (FP and FA). In summary, significant differences of soil properties were observed among the different rice farming systems.Table 1 Soil properties in different rice systems and areas.Full size tableSoil bacterial community compositionA total of 1,346,468 sequences were obtained by 16S rRNA MiSeq sequencing analysis after basal quality control (reads containing ambiguous bases were discarded; only overlapping sequences longer than 10 bp were assembled; Operational taxonomic units (OTUs) were clustered with 97% similarity). These sequences were classified as 46 phyla, 800 genera and 5335 OTUs. As shown in Fig. 1, the dominant bacterial phyla across different treatments were Proteobacteria (26.06–29.41%) and Chloroflexi (20.07–27.99%), followed by Actinobacteria (7.22–20.87%), Acidobacteria (11.36–14.46%) and Nitrospirae (3.11–8.50%). Since the implementation of rice–fish farming regime, the soil bacterial community composition has greatly changed. For example, Actinobacteria abundance decreased from 20.87% in RM to 7.22% in FA, while Nitrospirae abundance greatly increased from 3.11% in RM to 8.50% in FA. Between different areas in a same rice–fish field (i.e. OP vs OA or FP vs FA), the bacterial community composition were similar. The PCoA analysis on OTU level also showed that different areas within the same rice–fish field had high similarity in bacterial community composition. In contrast, the bacterial community composition differed distinctly among different rice farming systems (Fig. 2). Bacterial alpha diversity indices, as evaluated by Shannon, Simpson, ACE and Chao1, were shown in Table 2. Student’s t-test was adopted to evaluate the difference among treatments. The results showed that the alpha indices of FP were significantly lower than other treatments, except for Simpson index.Figure 1The average relative abundances on phylum level of soil bacterial communities in different rice systems and areas.Full size imageFigure 2PCoA analysis on OTU level based on bray_curtis distance algorithm (significance among treatments were conducted with ANOSIM test, R = 0.4294, P = 0.0010).Full size imageTable 2 Alpha diversity indices of soil bacterial in different rice systems and areas.Full size tableBased on the Kruskal–Wallis test, the statistical differences among treatments were evaluated in the abundances of the top 15 phyla. The results showed that 5 phyla, including Actinobacteria, Nitrospirae, Bacteroidetes, Unclassified_k_norank and SBR1093 were observed significant differences among treatments, and the most significant phylum was Nitrospirae (Fig. 3). In order to trace the source of the significant differences, the Wilcoxon tests were conducted between every two rice cultivation patterns separately (Fig. 4). The results indicated that the significant differences were mainly derived from the comparison between RM and F_group (FP & FA), as well as the comparison between the O_group (OP & OA) and F_group. In the comparison between the RM and O_group, only the phylum Gemmatimonadetes was observed to have a significant difference. Furthermore, we also compared the differences of the top 15 phyla between planting area (P_group) and aquaculture area (A_group) within rice–fish fields, and the results showed no phyla observed with significant differences in the abundances.Figure 3The differences with significance of the top 15 phyla in different rice systems and areas (* indicates 0.01  More

1.Cózar, A. et al. Plastic debris in the open ocean. Proc. Nat. Acad. Sci. USA 111, 10239–10244. https://doi.org/10.1073/pnas.1314705111 (2014).ADS 
CAS 
Article 
PubMed 

Google Scholar 
2.Lavers, J. L., Dicks, L., Dicks, M. R. & Finger, A. Significant plastic accumulation on the Cocos (Keeling) Islands, Australia. Sci. Rep. 9, 7102. https://doi.org/10.1038/s41598-019-43375-4 (2019).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
3.Cózar, A. et al. The arctic ocean as a dead end for floating plastics in the north atlantic branch of the thermohaline circulation. Sci. Adv. https://doi.org/10.1126/sciadv.1600582 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
4.Peeken, I. et al. Arctic sea ice is an important temporal sink and means of transport for microplastic. Nat. Commun. 9, 1505. https://doi.org/10.1038/s41467-018-03825-5 (2018).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
5.Woodall, L. C. et al. The deep sea is a major sink for microplastic debris. R. Soc. Open Sci. 1, 140317 (2014).ADS 
Article 

Google Scholar 
6.Chiba, S. et al. Human footprint in the abyss: 30 year records of deep-sea plastic debris. Mar. Policy 96, 204–212. https://doi.org/10.1016/j.marpol.2018.03.022 (2018).Article 

Google Scholar 
7.Bergmann, M., Tekman, M. & Gutow, L. Sea change for plastic pollution. Nature 544, 297 (2017).ADS 
CAS 
Article 

Google Scholar 
8.Jambeck, J. R. et al. Plastic waste inputs from land into the ocean. Science 347, 768–771. https://doi.org/10.1126/science.1260352 (2015).ADS 
CAS 
Article 
PubMed 

Google Scholar 
9.Gall, S. C. & Thompson, R. C. The impact of debris on marine life. Mar. Pollut. Bull. 92, 170–179. https://doi.org/10.1016/j.marpolbul.2014.12.041 (2015).CAS 
Article 
PubMed 

Google Scholar 
10.Camphuysen, C. J. Northern Gannets Morus bassanus found dead in the Netherlands, 1970–2000. Atlantic Seabirds 3, 15–30 (2001).
Google Scholar 
11.Gregory, M. R. Environmental implications of plastic debris in marine settings–entanglement, ingestion, smothering, hangers-on, hitch-hiking and alien invasions. Phil. Trans. R. Soc. B 364, 2013–2025 (2009).Article 

Google Scholar 
12.Ryan, P. G. The effects of ingested plastic on seabirds: Correlations between plastic load and body condition. Environ. Pollut. 46, 119–125 (1987).CAS 
Article 

Google Scholar 
13.Ryan, P. G. Effects of ingested plastic on seabird feeding: Evidence from chickens. Mar. Pollut. Bull. 19, 125–128 (1988).Article 

Google Scholar 
14.Pierce, K. E., Harris, R. J., Larned, L. S. & Pokras, M. A. Obstruction and starvation associated with plastic ingestion in a Northern Gannet Morus bassanus and a greater shearwater Puffinus gravis. Mar. Ornithol. 32, 187–189 (2004).
Google Scholar 
15.Ryan, P. G., Connell, A. D. & Gardner, B. D. Plastic ingestion and PCBs in seabirds: Is there a relationship?. Mar. Pollut. Bull. 19, 174–176 (1988).CAS 
Article 

Google Scholar 
16.Lavers, J. L., Bond, A. L. & Hutton, I. Plastic ingestion by Flesh-footed Shearwaters (Puffinus carneipes): Implications for chick body condition and the accumulation of plastic-derived chemicals. Environ. Pollut. 187, 124–129. https://doi.org/10.1016/j.envpol.2013.12.020 (2014).CAS 
Article 
PubMed 

Google Scholar 
17.Tanaka, K. et al. In vivo accumulation of plastic-derived chemicals into seabird tissues. Curr. Biol. 30, 723-728.e3. https://doi.org/10.1016/j.cub.2019.12.037 (2020).CAS 
Article 
PubMed 

Google Scholar 
18.Teuten, E. L. et al. Transport and release of chemicals from plastics to the environment and to wildlife. Phil. Trans. R. Soc. B 364, 2027–2045 (2009).CAS 
Article 

Google Scholar 
19.Tanaka, K., van Franeker, J. A., Deguchi, T. & Takada, H. Piece-by-piece analysis of additives and manufacturing byproducts in plastics ingested by seabirds: Implication for risk of exposure to seabirds. Mar. Pollut. Bull. 145, 36–41. https://doi.org/10.1016/j.marpolbul.2019.05.028 (2019).CAS 
Article 
PubMed 

Google Scholar 
20.Thiel, M. & Gutow, L. The ecology of rafting in the marine environment. I. The floating substrata. Oceanogr. Mar. Biol. Annu. Rev. 42, 181–264 (2005).
Google Scholar 
21.Kiessling, T., Gutow, L. & Thiel, M. Marine litter as habitat and dispersal vector. In: Bergmann M, Gutow L, Klages M, editors. Marine Anthropogenic Litter. p. 141–80 (2015).22.Day, R. H. & Shaw, D. G. Patterns of abundance of pelagic plastic and tar in the North Pacific Ocean, 1976–1985. Mar. Pollut. Bull. 18, 311–316 (1987).CAS 
Article 

Google Scholar 
23.Pichel, W. G. et al. Marine debris collects within the North Pacific Subtropical Convergence Zone. Mar. Pollut. Bull. 54, 1207–1211 (2007).CAS 
Article 

Google Scholar 
24.Yamashita, R. & Tanimura, A. Floating plastic in the Kuroshio Current area, western North Pacific Ocean. Mar. Pollut. Bull. 54, 485–488 (2007).CAS 
Article 

Google Scholar 
25.Titmus, A. J. & Hyrenbach, K. D. Habitat associations of floating debris and marine birds in the North East Pacific Ocean at coarse and meso spatial scales. Mar. Pollut. Bull. 62, 2496–2506 (2011).CAS 
Article 

Google Scholar 
26.Goldstein, M. C., Titmus, A. J. & Ford, M. Scales of spatial heterogeneity of plastic marine debris in the northeast pacific ocean. PLoS ONE 8, e80020 (2013).ADS 
Article 

Google Scholar 
27.Eriksen, M. et al. Plastic pollution in the world’s oceans: More than 5 trillion plastic pieces weighing over 250,000 tons afloat at sea. PLoS ONE 9, e111913 (2014).ADS 
Article 

Google Scholar 
28.IUCN. The IUCN Red List of Threatened Species. Version 2020–2. https://www.iucnredlist.org (2020).29.Lavers, J. L. & Bond, A. L. Ingested plastic as a route for trace metals in Laysan Albatross (Phoebastria immutabilis) and Bonin Petrel (Pterodroma hypoleuca) from Midway Atoll. Mar. Pollut. Bull. 110, 493–500. https://doi.org/10.1016/j.marpolbul.2016.06.001 (2016).CAS 
Article 
PubMed 

Google Scholar 
30.Roman, L., Hardesty, B. D., Hindell, M. A. & Wilcox, C. A quantitative analysis linking seabird mortality and marine debris ingestion. Sci. Rep. 9, 3202. https://doi.org/10.1038/s41598-018-36585-9 (2019).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
31.Jouventin, P. & Weimerskirch, H. Satellite tracking of wandering albatrosses. Nature 343, 746–748 (1990).ADS 
Article 

Google Scholar 
32.Kappes, M. A. et al. Hawaiian albatrosses track interannual variability of marine habitats in the North Pacific. Prog. Oceanogr. 86, 246–260 (2010).ADS 
Article 

Google Scholar 
33.Sakamoto, K. Q., Takahashi, A., Iwata, T. & Trathan, P. N. From the eye of the albatrosses: A bird-borne camera shows an association between albatrosses and a killer whale in the Southern Ocean. PLoS ONE 4, e7322 (2009).ADS 
Article 

Google Scholar 
34.Fukuoka, T. et al. The feeding habit of sea turtles influences their reaction to artificial marine debris. Sci. Rep. 6, 28015. https://doi.org/10.1038/srep28015 (2016).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
35.Nishizawa, B. et al. Albatross-borne loggers show feeding on deep-sea squids: Implications for the study of squid distributions. Mar. Ecol. Prog. Ser. 592, 257–265 (2018).ADS 
Article 

Google Scholar 
36.Hunt, G. L. Jr. & Schneider, D. Scale-dependent processes in the physical and biological environment of marine birds. In Seabirds: Feeding Ecology and Role in Marine Ecosystems (ed. Croxall, J. P.) 7–41 (Cambridge University Press, 1987).
Google Scholar 
37.Pinaud, D. & Weimerskirch, H. At-sea distribution and scale-dependent foraging behaviour of petrels and albatrosses: A comparative study. J. Anim. Ecol. 76, 9–19 (2007).Article 

Google Scholar 
38.Thiebot, J.-B., Nishizawa, B., Sato, F., Tomita, N. & Watanuki, Y. Albatross chicks reveal interactions of adults with artisanal longline fisheries within a short range. J. Ornithol. 159, 935–944 (2018).Article 

Google Scholar 
39.Froese, R. & Pauly, D. FishBase. World Wide Web electronic publication. www.fishbase.org, version (12/2019).40.Ryan, P. G. A simple technique for counting marine debris at sea reveals steep litter gradients between the Straits of Malacca and the Bay of Bengal. Mar. Pollut. Bull. 69, 128–136 (2013).CAS 
Article 

Google Scholar 
41.Mitani, Y. et al. Marine debris observed in the North Pacific during Oshoro-maru cruise in 2012. Bull. Fish. Sci. Hokkaido Univ. 64, 25–29 (2014).
Google Scholar 
42.Hyrenbach, K. D. et al. Plastic ingestion by Black-footed albatross from Kure Atoll, Hawai’i: linking chick loads and parental at-sea distributions. Mar. Ornithol. 45, 225–236 (2017).
Google Scholar 
43.Nevitt, G. A., Losekoot, M. & Weimerskirch, H. Evidence for olfactory search in wandering albatross, Diomedea Exulans. Proc. Nat. Acad. Sci. USA 105, 4576–4581 (2008).ADS 
CAS 
Article 

Google Scholar 
44.Savoca, M. S., Wohlfeil, M. E., Ebeler, S. E. & Nevitt, G. A. Marine plastic debris emits a keystone infochemical for olfactory foraging seabirds. Sci. Adv. 2, e1600395 (2016).ADS 
Article 

Google Scholar 
45.Santos, R. G., Andrades, R., Fardim, L. M. & Martins, A. S. Marine debris ingestion and Thayer’s law—The importance of plastic color. Environ. Pollut. 214, 585–588 (2016).CAS 
Article 

Google Scholar 
46.Castro, J. J., Santiago, J. A. & Santana-Ortega, A. T. A general theory on fish aggregation to floating objects: An alternative to the meeting point hypothesis. Rev. Fish Biol. Fish. 11, 255–277 (2002).Article 

Google Scholar 
47.Harrison, C. S., Hida, T. S. & Seki, M. P. Hawaiian seabird feeding ecology. Wildl. Monogr. 85, 1–71 (1983).
Google Scholar 
48.Hunte, W., Oxenford, H. A. & Mahon, R. Distribution and relative abundance of flyingfish (Exocoetidae) in the eastern Caribbean. II. Spawning substrata, eggs and larvae. Mar. Ecol. Prog. Ser. 117, 25–37 (1995).ADS 
Article 

Google Scholar 
49.Rapp, D. C., Youngren, S. M., Hartzell, P. & Hyrenbach, K. D. Community-wide patterns of plastic ingestion in seabirds breeding at French Frigate Shoals Northwestern Hawaiian Islands. Mar. Pollut. Bull. 123, 269–278 (2017).CAS 
Article 

Google Scholar 
50.Douglas, D. & Peucker, T. Algorithms for the reduction of the number of points required to represent a digitized line or its caricature. Cannadian Cartogr. 10, 112–122 (1973).Article 

Google Scholar 
51.Edelhoff, H., Signer, J. & Balkenhol, N. Path segmentation for beginners: an overview of current methods for detecting changes in animal movement patterns. Move. Ecol. 4, 21 (2016).Article 

Google Scholar 
52.R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. http://www.r-project.org/index.html (2020). More