Philippa Kaur

1.Simões, R. A., Reis, L. G., Bento, J. M., Solter, L. F. & Delalibera, I. Jr. Biological and behavioral parameters of the parasitoid Cotesia flavipes (Hymenoptera: Braconidae) are altered by the pathogen Nosema sp. (Microsporidia: Nosematidae). Biol. Control 63, 164–171 (2012).Article 

Google Scholar 
2.Frago, E., Dicke, M. & Godfray, H. C. J. Insect symbionts as hidden players in insect–plant interactions. Trends Ecol. Evol. 27, 705–711 (2012).Article 

Google Scholar 
3.Himler, A. G. et al. Rapid spread of a bacterial symbiont in an invasive whitefly is driven by the fitness benefits and female bias. Science 332, 254–256 (2011).ADS 
CAS 
Article 

Google Scholar 
4.Lu, M., Wingfield, M. J., Gillette, N. & Sun, J. H. Do novel genotypes drive the success of an invasive bark beetle–fungus complex? Implications for potential reinvasion. Ecology 92, 2013–2019 (2011).Article 

Google Scholar 
5.Vilcinskas, A., Stoecker, K., Schmidtberg, H., Röhrich, C. R. & Vogel, H. Invasive harlequin ladybird carries biological weapons against native competitors. Science 340, 862–863 (2013).ADS 
CAS 
Article 

Google Scholar 
6.Zhao, L. et al. A native fungal symbiont facilitates the prevalence and development of an invasive pathogen–native vector symbiosis. Ecology 94, 2817–2826 (2013).Article 

Google Scholar 
7.Solter, L. F., Becnel, J. J. & Vávra, J. Research methods for entomopathogenic microsporidia and other protists. Manual of Techniques in Invertebrate Pathology 329–371 (2012).8.Maddox, J. V. Protozoan diseases. Epizootiol. Insect Dis. 1, 417–452 (1987).
Google Scholar 
9.Latchininsky, A. V. & VanDyke, K. A. Grasshopper and locust control with poisoned baits: a renaissance of the old strategy?. Outlooks Pest Manag. 17, 105–111 (2006).Article 

Google Scholar 
10.Sweeney, A. W. & Becnel, J. J. Potential of microsporidia for the control of mosquitoes. Parasitol. Today. 7, 217–220 (1991).CAS 
Article 

Google Scholar 
11.Capella-Gutierrez, S., Marcet-Houben, M. & Gabaldon, T. Phylogenomics supports microsporidia as the earliest diverging clade of sequenced fungi. BMC Biol. 10, 47–52 (2012).Article 

Google Scholar 
12.Tokarev, Y. S. et al. A formal redefinition of the genera Nosema and Vairimorpha (Microsporidia: Nosematidae) and reassignment of species based on molecular phylogenetics. J. Invertebr. Pathol. 169, 107279 (2020).CAS 
Article 

Google Scholar 
13.Tooke, F. G. C. The Eucalyptus Snout beetle, Gonipterus scutellatus Gyll. A study of its ecology and control by biological means Union of South Africa, Department of Agriculture. Entomol. Mem. 3, 1–184 (1955).
Google Scholar 
14.Mapondera, T. S., Burgess, T., Matsuki, M. & Oberprieler, R. G. Identification and molecular phylogenetics of the cryptic species of the Gonipterus scutellatus complex (Coleoptera: Curculionidae: Gonipterini). Aust. J. Entomol. 51, 175–188 (2012).Article 

Google Scholar 
15.Valente, C. et al. Economic outcome of classical biological control: a case study on the Eucalyptus snout beetle, Gonipterus platensis, and the parasitoid Anaphes nitens. Ecol Econ. 149, 40–47 (2018).Article 

Google Scholar 
16.Ansari, M. J., Al-Ghamdi, A., Nuru, A., Khan, K. A. & Alattal, Y. Geographical distribution and molecular detection of Nosema ceranae from indigenous honeybees of Saudi Arabia. Saudi J. Biol. Sci 24, 983–991 (2017).Article 

Google Scholar 
17.Ovcharenko, M., Świątek, P., Ironside, J. & Skalski, T. Orthosomella lipae sp. n. (Microsporidia) a parasite of the weevil, Liophloeus lentus Germar, 1824 (Coleoptera: Curculionidae). J. Invertebr. Pathol. 112, 33–40 (2013).Article 

Google Scholar 
18.Weiser, J. A new microsporidian from the bark beetle Pityokteines curvidens Germar (Coleoptera, Scolytidae) in Czechoslovakia. J. Invertebr. Pathol. 3, 324–329 (1961).
Google Scholar 
19.Malone, L. A. A new pathogen, Microsporidium itiiti n. sp. (Microsporida), from the Argentine Stem Weevil, Listronotus bonariensis (Coleoptera, Curculionidae). J. Protozool. 32, 535–541 (1985).Article 

Google Scholar 
20.Purrini, K. & Weiser, J. Ultrastructural study of the microsporidian Chytridiopsis typographi (Chytridiopsida: Microspora) infecting the bark beetle, Ips typographus (Scolytidae: Coleoptera), with new data on spore dimorphism. J. Invertebr. Pathol. 45, 66–74 (1985).Article 

Google Scholar 
21.Yaman, M., Radek, R., Aslan, I. & Erturk, O. Characteristic features of Nosema phyllotretae Weiser 1961, a microsporidian parasite of Phyllotreta atra (Coleoptera: Chrysomelidae) in Turkey. Zool. Stud. Taipei. 44, 368 (2005).
Google Scholar 
22.Zhu, F. et al. A new isolate of Nosema sp. (Microsporidia, Nosematidae) from Phyllobrotica armata Baly (Coleoptera, Chrysomelidae) from China. Jour J. Invertebr. Pathol. 106, 339–342 (2011).CAS 
Article 

Google Scholar 
23.Andreadis, T. G., Takaoka, H., Otsuka, Y. & Vossbrinck, C. R. Morphological and molecular characterization of a microsporidian parasite, Takaokaspora nipponicus n. gen. n. sp. from the invasive rock pool mosquito, Ochlerotatus japonicus japonicus. J. Invertebr. Pathol. 114, 161–172 (2013).CAS 
Article 

Google Scholar 
24.Sapir, A. et al. Microsporidia-nematode associations in methane seeps reveal basal fungal parasitism in the deep sea. Front. Microbiol. 5, 43–52 (2014).Article 

Google Scholar 
25.Solter, L. F., Maddox, J. V. & McManus, M. L. Host specificity of microsporidia (Protista: Microspora) from European populations of Lymantria dispar (Lepidoptera: Lymantriidae) to indigenous North American Lepidoptera. J. Invertebr. Pathol. 69, 135–150 (1997).CAS 
Article 

Google Scholar 
26.Knell, J. D., Allen, G. E. & Hazard, E. I. Light and electron microscope study of Thelohania solenopsae n. sp. (Microsporida: Protozoa) in the red imported fire ant Solenopsis invict. J. Invertebr. Pathol. 29, 192–200 (1977).CAS 
Article 

Google Scholar 
27.Henry, J. E., & Oma, E. A. Pest control by Nosema locustae, a pathogen of grasshoppers and crickets. Microbial Control of Pests and Plant Diseases 1970–1980 (1981).28.Vávra, J. & Maddox, J. V. Methods in microsporidiology. In Biology of the Microsporidia 281–319 (Springer, Boston, 1976).
Google Scholar 
29.Simões, R. A., Feliciano, J. R., Solter, L. F. & Delalibera, I. Jr. Impacts of Nosema sp. (Microsporidia: Nosematidae) on the sugarcane borer, Diatraea saccharalis (Lepidoptera: Crambidae). J. Invertebr. Pathol. 129, 7–12 (2015).Article 

Google Scholar 
30.Inglis, G. D., Lawrence, A. M. & Davis, F. M. Impact of a novel species of Nosema on the southwestern corn borer (Lepidoptera: Crambidae). J. Econ. Entomol. 96, 12–20 (2003).CAS 
Article 

Google Scholar 
31.Zheng, H. Q. et al. Spore loads may not be used alone as a direct indicator of the severity of Nosema ceranae infection in honey bees Apis mellifera (Hymenoptera: Apidae). J. Econ. Entomol. 107, 2037–2044 (2014).Article 

Google Scholar 
32.Goettel, M. S., Inglis, G. D. & Lacey, L. A. Manual of Techniques in Invertebrate Pathology (Academic Press, 2012).
Google Scholar 
33.Canning, E. U., Curry, A., Cheney, S., Lafranchi-Tristem, N. J., Haque, M. A. Vairimorpha imperfecta n. sp., a microsporidian exhibiting an abortive octosporous sporogony in Plutella xylostella L. (Lepidoptera: Yponomeutidae). Parasitology 119, 273–286 (1999).34.Tsai, S. J., Lo, C. F., Soichi, Y. & Wang, C. H. The characterization of microsporidian isolates (Nosematidae: Nosema) from five important lepidopteran pests in Taiwan. J. Invertebr. Pathol. 83, 51–59 (2003).CAS 
Article 

Google Scholar 
35.Cai, S. F., Lu, X. M., Qiu, H. H., Li, M. Q. & Feng, Z. Z. Phagocytic uptake of Nosema bombycis (Microsporidia) spores by insect cell lines. J. Integr. Agric. 11, 1321–1326 (2012).Article 

Google Scholar 
36.Dong, S., Shen, Z., Xu, L. & Zhu, F. Sequence and phylogenetic analysis of SSU rRNA gene of five microsporidia. Curr. Microbiol. 60, 30 (2010).CAS 
Article 

Google Scholar 
37.Becnel, J. J. & Andreadis, T. G. Microsporidia in insects. The microsporidia and microsporidiosis 447-501 (1999).38.Knell, R. J. & Webberley, K. M. Sexually transmitted diseases of insects: Distribution, evolution, ecology and host behaviour. Biol. Rev. 79, 557–581 (2004). (PERMANECE)39.Bell, H. A., Down, R. E., Kirkbride‐Smith, A. E. & Edwards, J. P. Effect of microsporidian infection in Lacanobia oleracea (Lep., Noctuidae) on prey selection and consumption by the spined soldier bug Podisus maculiventris (Het., Pentatomidae). J. Appl. Entomol. 128(8), 548–553 (2004).40.Dakhel, W. H., Latchininsky, A. V. & Jaronski, S. T. Efficacy of two entomopathogenic fungi, Metarhizium brunneum, strain F52 alone and combined with Paranosema locustae against the migratory grasshopper, Melanoplus sanguinipes, under laboratory and greenhouse conditions. Insects 10(4), 94–102 (2019).Article 

Google Scholar 
41.Guo, Y., An, Z. & Shi, W. Control of grasshoppers by combined application of Paranosema locustae and an insect growth regulator (IGR) (cascade) in rangelands in China. J. Econ. Entomol. 105(6), 1915–1920 (2012).Article 

Google Scholar 
42.Lockwood, J. A., Bomar, C. R. & Ewen, A. B. The history of biological control with Nosema locustae: Lessons for locust management. Int. J. Trop. Insect Sci. 19(4), 333–350 (1999).Article 

Google Scholar 
43.Larem, A., Fritsch, E., Undorf-Spahn, K., Kleespies, E. G. & Jehle, J. A. Interaction of Phthorimaea operculella granulovirus with a Nosema sp. microsporidium in larvae of Phthorimaea operculella. J. Invertebr. Pathol. 160, 76–86 (2019).Article 

Google Scholar 
44.Tokarev, Y. S., Grizanova, E. V., Ignatieva, A. N. & Dubovskiy, I. M. Greater wax moth Galleria mellonella (Lepidoptera: Pyralidae) as a resistant model host for Nosema pyrausta (Microsporidia: Nosematidae). J. Invertebr. Pathol. 157, 1–3 (2018).Article 

Google Scholar 
45.Coombs, N. J., Gough, A. C. & Primrose, J. N. Optimisation of DNA and RNA extraction from archival formalin-fixed tissue. Nucleic Acids Res. 27, e12-I (1999).46.Huang, W. F., Tsai, S. J., Lo, C. F., Soichi, Y. & Wang, C. H. The novel organization and complete sequence of the ribosomal RNA gene of Nosema bombycis. Fungal Genet. Biol 41, 473–481 (2004).CAS 
Article 

Google Scholar 
47.Karnovsky, M. J. A formaldehyde glutaraldehyde fixative of high osmolality for use in electron microscopy. J. Cell. Biol. 27, 1A-149A (1965).Article 

Google Scholar  More

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

MDE hypothesisGradient features of species diversity of plant communities refer to regular changes in species diversity along a gradient of environmental factors at the community level12,43. The altitudinal gradient includes gradient effects of multiple environmental factors. It is therefore important to study altitudinal patterns of species diversity to reveal changes in biodiversity along environmental gradients. The width and range of species distribution along geographical gradients reflect the ecological adaptability, diffusivity, and evolutionary history of species44. To some extent, geographical distribution patterns of species diversity can be interpreted as outcomes of synthetic actions across altitudinal gradients resulting from eurychoric species with a greater distribution width and stenochoric species with a smaller distribution range along geographical gradients44. Hence, the MDE, environmental gradient, distribution area, human disturbance, and habitat heterogeneity all have effects on the vertical distribution patterns of species diversity45,46. According to the MDE hypothesis, there is overlap in the distribution range of species along altitudinal gradients; the highest overlap intensity occurs at intermediate elevations. There is relatively low overlap intensity at low and high elevation areas, and the peak values of species diversity occur at intermediate elevations.In this study, forest ecosystems on the Loess Plateau were separated into three communities: tree, shrub, and herb communities. The altitudinal patterns and factors that influence species diversity of mountainous vegetation were then determined at the plant community level related to the form and family of the plant species. We discovered that the family numbers of the herb, shrub, and tree communities reached their greatest values at intermediate, intermediate, and lower elevations, respectively. We also discovered that correlations of species diversity indices with elevations conformed to unimodal change patterns for herb, shrub, and tree communities, which presented their greatest values at higher, lower, and intermediate elevations, respectively. This showed that MDE is another important factor that affects the distribution patterns of species diversity along regional altitudinal gradients, in addition to temperature, precipitation, and the terrain.A large number of studies have already verified that MDE is a significant mechanism that influences the gradient patterns of species diversity. MDE not only functions along altitudinal gradients but also acts along latitudinal and temporal gradients28,29,30,47. However, the effects of MDE on species diversity patterns are highly controversial. Some studies considered MDE to be the main factor that results in unimodal patterns of species diversity47,48, whereas other studies affirmed that the effects of MDE are smaller in contrast to the functions of the distribution area, environmental gradient, and other factors29,30. Besides MDE, other factors may also lead to unimodal vertical gradient patterns of species diversity, such as plant productivity, human disturbance, and the regional climate45,49. Our research indicated that the relationships of species diversity conformed to unimodal change patterns along various elevations for mountainous herb, shrub, and tree communities in a semi-arid region of the temperate zone. It can be concluded that vertical patterns of species diversity with a unimodal type may be a more universal phenomenon, relative to monotonic decreasing or increasing patterns of species diversity with increasing elevations.Factors influencing MDE at the plant community levelIn this study, more forbs and grasses were found at higher elevations, whereas more sedges occurred at lower elevations. The responses of the importance values of tree species to the altitudinal gradient demonstrated the following variation patterns: evergreen coniferous trees had higher importance values than deciduous coniferous trees, followed by deciduous broad-leaved trees. This showed that MDE was influenced by species life form. The species diversity of different life forms responded differently to the environment, and plant species with different life forms presented different diversity patterns along altitudinal gradients50. In New Zealand, the number of mountainous plants species decreased with increasing elevations and the total species number of all plants also decreased significantly, while species diversity had no significant distribution trend in response to elevation when plants with different life forms were considered under different layers of plant communities51.MDE is a common pattern of species diversity of mountainous plants with changing elevations. Our study area, located on the Loess Plateau of China, belongs to a semi-arid mountainous region in the temperate zone where the maximum species diversity of the tree community occurred at intermediate elevations. This finding was in accordance with the MDE hypothesis. Research from the Kinabalu Mountains in Sabah, Malaysia, indicated an obvious MDE pattern of species number linked to elevations52. On the Haleakala Mountains, Hawaii, USA, the highest species diversity occurred at intermediate elevations, which was also in accordance with the MDE hypothesis53. The MDE hypothesis was also proved by studies conducted in the Yu Mountains of Taiwan and the Emei Mountains of Szechwan in China54. However, the MDE pattern of species diversity in tree communities is caused by precipitation, which is the highest at intermediate elevations52. This situation also occurred in the herb community.There are many factors that affect the distribution of herbaceous plants, so the variations in species diversity in the herb community with elevations are complex55. In this study, we found that the herb community exhibited higher species diversity at higher elevations; more forbs and grasses were distributed at higher elevations, whereas more sedges were distributed at lower elevations. In the Siskiyou Mountains in Oregon, USA, the species diversity of herbaceous plants had a significantly positive correlation with elevation. This correlation occurred mainly due to an increase in the number of grass species, which was the primary reason that radiation was enhanced by a drastic reduction in community coverage as a result of increasing elevations. Consequently, there was an increase in the species diversity of herbaceous plants24. A decrease in species diversity with increasing elevations is a more familiar pattern for herbaceous plants in temperate56 and tropical22 forests.We also discovered that the family numbers of herb and shrub species all showed unimodal change patterns with high values at their central elevations in this semi-arid region. This conformed to the MDE hypothesis as well. In arid temperate grasslands, species diversity indicated an MDE distribution pattern. For example, the species diversity of herbaceous plants presented an MDE pattern in drought areas of the Siskiyou Mountains57. However, in semi-humid mountains in the temperate zone, the species diversity of herbaceous plants was principally in control of the community structure, and community coverage did not respond uniformly to elevation. Studies in New Zealand showed that there were no evident distribution trends for species diversity of herbaceous plants along elevations51. In low bush communities of Chile, the species diversity of herbaceous plants declined with increasing elevations after longstanding succession, but it increased during the early stage of succession58. Therefore, relationships between the species diversity of herb species and elevations were not completely clear in semi-humid regions.The major factors that control the distribution areas of species differ among different families and genera, and thus the vertical distribution patterns of species diversity differ22. We concluded that the family number of the tree species had a maximum distribution at lower elevations, unlike herb and shrub species; meanwhile, the responses of the importance values to the altitudinal gradient in the tree community were also different among evergreen coniferous trees, deciduous coniferous trees, and deciduous broad-leaved trees. These differences may have been related to environmental factors. Due to various distribution patterns of environmental parameters with elevations, the distribution patterns of species diversity showed large changes along elevations59. For example, the distribution of fern and Melastomataceae species is principally related to humidity, that of Acanthaceae and Bromeliaceae species is correlated with temperature, and that of Araceae species is related to transpiration59. Research conducted in the Gongga Mountains, China, showed that the species diversity with different floral components exhibited different distribution patterns along elevations due to differences in the environment and species origin60. We also discovered that the importance values of dominant families in the shrub (Rosaceae) and tree (Pinaceae) communities exhibited changing patterns in contrast to MDE. In our study, only the family numbers in the herb and shrub communities, as well as species diversity in the tree community, conformed to the MDE hypothesis. Therefore, we concluded that the MDE hypothesis of species diversity of mountainous vegetation is influenced by the species life form and family of different plant communities in the temperate semi-arid region of China.Factors influencing plant species diversity at the environment levelThe altitudinal distribution patterns of the plant community diversity had greater discrepancies in mountainous regions and between different community types, which might be connected to regional environmental conditions, relative heights of mountains, and the geological landscape35. Concerning the altitudes of mountains, serious human disturbances (e.g., deforestation, grazing, and land-use conversion) had negative effects on biodiversity in low-altitude regions61. In high-altitude regions, the cold climate slowed down plant growth and soil development, while other harsh environments exceeded the tolerance limits for growth of the majority of species, such as by intense solar radiation or large temperature differences between day and night62. In the middle-altitude regions, the species diversity was relatively higher due to less human disturbances and the formation of transition zones of plant species differentiation between the low- and high-altitude regions62. Hence, the plant community diversity and its altitudinal gradient patterns in mountainous regions were largely influenced by regional climate and human disturbances.Comparisons of the diversity at different levels indicated that the responses of the plant community diversity to the environment were not the same for diverse gradations, and different species exhibited different gradient patterns owing to restrictions from environmental factors63,64. The primary factors leading to the altitudinal differentiation of diversity included the temperature, moisture, soil nutrients, and succession process65. In our study area, compared with Guancen Mountain and Guandi Mountain, Wulu Mountain at the lower latitude of the Lvliang Mountains had a lower altitude and was located in the continental monsoon subhumid climate region of the warm temperate zone, making it suitable for the growth of secondary forests and shrub vegetation. However, the vegetation growth in the herb layer was restricted in Wulu Mountain, making diversity in the herb layer the greatest on Guandi Mountain at the middle latitude of the Lvliang Mountains35. This showed that the plant species diversity in the east of the Loess Plateau changed with the altitude, while being affected by complicated habitat conditions such as latitude and human disturbances. This characterization of the study area was correlated with the unimodal patterns observed. In this study, we observed that the family numbers of the herb and shrub communities presented unimodal patterns across an altitudinal gradient; the importance values of dominant families also presented unimodal patterns in the shrub and tree communities; and the species diversity indices of the herb, shrub, and tree communities conformed to unimodal change patterns following an altitudinal gradient as well.In our recent studies on the species diversity of herbaceous communities in the Lvliang Mountains66, we found that the results calculated for β-diversity using different indices revealed the highest value for the Cody index and the lowest value for the Bray–Curtis index at altitudes between 1900 and 2000 m, indicating that areas located between 1900 and 2000 m form a transition zone in which the herbaceous community undergoes a rapid process of species renewal and changes in its composition. The results for γ-diversity indicated a pattern of unimodal variation in relation to altitude. Changes in altitude gradient had highly significant impacts on changes in temperature and humidity, indicating that various environmental factors (notably humidity and temperature) and human disturbances had combined effects on changes in the values of the α-diversity indices.At present, it is widely believed that the formation of herbs in different life forms was principally impacted by precipitation, whereas in areas with similar rainfall, water, heat, and light conditions need to be considered. These conditions chiefly included average annual precipitation, accumulated temperature, and illumination time67. In this research, we observed that herb and tree species in different life forms showed different trends with altitudinal gradients in the Lvliang Mountains. At higher elevations, forbs and grasses grew well, whereas sedges grew well at lower elevations. The responses of different tree life forms to the altitudinal gradient were greater for evergreen coniferous tree species than for deciduous coniferous tree species and deciduous broad-leaved tree species. From the north to the south in the Lvliang Mountains, increases in the average annual precipitation increased the number of species and components of the annual herbs, while the hydrothermal matching requirements of Guandi Mountain at the middle latitude were preferred for annual herb growth in comparison with Guancen Mountain at the higher latitude67. However, considering whole mountains, the Lvliang Mountains located in the continental monsoon climate region of the warm temperate zone had four distinctive seasons with drought and wind in the spring, and a quick rise in air temperatures, and larger diurnal temperature differences35. These conditions conformed to the habitat features of herbs and trees. Hence, the hydrothermal distribution status affected by latitude and human disturbances determined the altitudinal distribution patterns of plant community diversity in the Lvliang Mountains.The MDE at different elevationsIn this study, we discussed the MDE of mountainous vegetation on the Loess Plateau with an elevational range from 1324 to 2745 m, including tree, shrub and herb community. This range was a very large elevation range for a case study, but a very short range in comparison to global elevation ranges, which extended from the sea level to well over 8000 m (though the highest locations did not have any vegetation). Therefore, owing to this limitation, the results we obtained in this research were suitable for lower elevation mountains in semi-arid areas.The MDE was changed with different elevations and vegetation layers. In studies on the Daiyun Mountains with an elevation from 900 to 1600 m, the phylogenetic diversity and species diversity of tree community indicated an intermediate high expansion pattern along elevations and their peak values all appeared at the elevation of 1200 m68. This conclusion conformed to the MDE pattern. In our studies on the Lvliang Mountains with an elevation from 1459 to 2610 m, higher species diversity of tree community was observed at intermediate elevations with a peak value being at the 2000 m, which conformed to the MDE pattern either. Therefore, at smaller elevations less than 2600 m, species diversity of tree conformed to the MDE pattern.When an elevation reached 2700 m on the Lvliang Mountains, the vegetation types changed to shrub and grass, and only their family numbers followed the MDE pattern across an altitudinal gradient. Slimily, the species richness of shrub and herb layer showed an obvious “lateral pattern” on an elevational gradient from 2950 to 3750 m on the Three River Headwater, both reaching the maximum value at the 3150 m; while with the rise of altitude, α diversity of shrub layer and herb layer showed a “wave”-shaped changed trend, reaching the lowest value at the 3550 m69. It illustrated that species diversity of shrub and grass did not conform to the MDE pattern completely at medium elevations from 2700 to 3700 m.As for an elevation extending from 3000 to 4400 m on the Gongga Mountains, the vegetation type was alpine meadow, and the species richness index presented an obvious unimodal pattern with a peak value at the 3850 m, which accorded with the MDE pattern70. Similarly, studies on an alpine meadow on the Gannan revealed that the number of richness, Shannon-Weiner index and phylogenetic diversity of plant community all showed a “humped-back” relationship with the increase of altitude from 3000 to 4000 m, and reached the maximum value at the 3800 m71,72. Thereby, at greater elevations more than 3800 m, species diversity of alpine meadow conformed to the MDE pattern.However, when an elevation exceeded 5000 m, research object on species diversity were not vegetation but animals along an altitudinal gradient. For example, on the Himalaya Mountains with an elevation from 3755 to 5016 m, the ant species richness illustrated a “unimodal curve” pattern along the rise of altitude, and the Shannon–Wiener index and Fisher α index of ant community commonly expressed the “Multi-Domain Effect” phenomenon73. Another research on mammalian richness was also conducted on the Himalaya Mountains. It concluded that most of elevational species richness patterns were hump-shaped from 100 to 6000 m on the Himalayas Mountains74. As a result, the MDE pattern was also extremely common in animal communities. More

The catenary sequence of soilsThe catena of Andosols down a slope near a hot spring in the Valley of Geysers was subdivided into four thermal zones (Fig. 1a–e), which are described below.Figure 1Location of study area (a–c), soil pits along the catena (d, e) and photos of soil pits (f–i). (a b) Study area location. Soil map is from133 (open access) with additions and corrections by I.N. Semenkov based on the map of soil temperature at a depth of 15 cm in the Valley of Geysers131 and the soil names from44 using CorelDraw X7 software (https://www.coreldraw.com/). (c) Top view of the left side of the Geysernaya River with the body of a catastrophic landslide, a visitor center (in the left lower part) and the location of the transect studied (1–15). (d) The location of transect studied. (e) A schematic profile of the catena studied with numbered soil pits. The main soil pits selected for comprehensive analyses (see section ‘Soil analyses’) are in red. (f) Non-heated Eutrosilic Silandic Andosols (Arenic, Cutanic) on pyroclastic material (pit no 16, Zone I), within levelled parts of the interfluve, under tall-herb meadow communities with local patches of Erman’s birch woods. (g) Slightly heated Eutrosilic Aluandic Andosols (Cutanic, Loamic, Natric) in the upper part of the catena, on hydrothermally altered sandy-loamy pyroclastic material (pit no 12, Zone II), on slightly heated slopes, under tall-herb meadows. (h) Moderately heated Eutrosilic Gleyic Aluandic Andosols (Loamic, Reductic, Protosalic, Hyperthionic) in the middle part of the catena, on hydrothermally altered clayey pyroclastic material (pit no 9.1, Zone III), under different moss and ‘microzonal’ communities. (i) Hot Gleyic Aluandic Andosols (Clayic, Reductic, Salic, Hyperthionic) in the lower part of the catena, on hydrothermal clays (pit no 4, Zone IV), on most heated bare slopes.Full size imageZone I. Non-heated Eutrosilic Silandic Andosols (Arenic, Cutanic) under Kamchatka’s tall herb communities and fragmented Erman’s birch woodsNon-heated Andosols with temperatures of 50 with very small amounts of PM5–50 ( More

1.Moller, H. Lessons for invasion theory from social insects. Biol. Conserv. 78, 125–142 (1996).Article 

Google Scholar 
2.Pimentel, D. et al. Economic and environmental threats of alien plant, animal, and microbe invasions. Agric. Ecosyst. Environ. 84, 1–20 (2001).Article 

Google Scholar 
3.Kenis, M. et al. Ecological effects of invasive alien insects. Biol. Invasion 11, 21–45 (2009).Article 

Google Scholar 
4.McGlynn, T. P. The worldwide transfer of ants: Geographical distribution and ecological invasions. J. Biogeog. 26, 535–548 (1999).Article 

Google Scholar 
5.Morrison, L. W., Porter, S. D., Daniels, E. & Korzukhin, M. D. Potential global range expansion of the invasive fire ant, Solenopsis invicta. Biol. Invasions 6, 183–191 (2004).Article 

Google Scholar 
6.Fitzpatrick, M. C., Weltzin, J. F., Sanders, N. J. & Dunn, R. R. The biogeography of prediction error: Why does the introduced range of the fire ant over-predict its native range?. Global Ecol. Biogeog. 16, 24–33 (2007).Article 

Google Scholar 
7.Bertelsmeier, C., Luque, G. M., Hoffmann, B. D. & Courchamp, F. Worldwide ant invasions under climate change. Biodivers. Conserv. 24, 117–128 (2015).Article 

Google Scholar 
8.Rowles, A. D. & O’Dowd, D. J. Impacts of the invasive Argentine ant on native ants and other invertebrates in coastal scrub in south-eastern Australia. Austral. Ecol. 34, 239–248 (2009).Article 

Google Scholar 
9.Naumann, K. & Higgins, R. J. The European fire ant (Hymenoptera: Formicidae) as an invasive species: Impact on local ant species and other epigaeic arthropods. Can. Entomol. 147, 592–601 (2015).Article 

Google Scholar 
10.Uchida, S. et al. Effects of an invasive ant on land snails in the Ogasawara Islands. Conserv. Biol. 30, 1330–1337 (2016).PubMed 
Article 

Google Scholar 
11.Helanterä, H., Strassmann, J. E., Carrillo, J. & Queller, D. C. Unicolonial ants: Where do they come from, what are they and where are they going?. Trends Ecol. Evol. 24, 341–349 (2009).PubMed 
Article 

Google Scholar 
12.Holway, D. A., Lach, L., Suarez, A. V., Tsutsui, N. D. & Case, T. J. The causes and consequences of ant invasions. Ann. Rev. Ecol. Syst. 33, 181–233 (2002).Article 

Google Scholar 
13.Stafford, C. T. Hypersensitivity to fire ant venom. Ann. Allerg. Asthma Im. 77, 87–99 (1996).CAS 
Article 

Google Scholar 
14.Kemp, S. F., deShazo, R. D., Moffitt, J. E., Williams, D. F. & Buhner, W. A. II. Expanding habitat of the imported fire ant (Solenopsis invicta): A public health concern. J. Allergy Clin. Immun. 105, 683–691 (2000).CAS 
PubMed 
Article 

Google Scholar 
15.Morrison, J. E. Jr., Williams, D. F., Oi, D. H. & Potter, K. N. Damage to dry crop seed by red imported fire ant (Hymenoptera: Formicidae). J. Econ. Entomol. 90, 218–222 (1997).Article 

Google Scholar 
16.Mikissa, J. B., Jeffery, K., Fresneau, D. & Mercier, J. L. Impact of an invasive alien ant, Wasmannia auropunctata Roger., on a specialised plant-ant mutualism, Barteria fistulosa Mast. and Tetraponera aethiops F. Smith., in a Gabon forest. Ecol. Entomol. 38, 580–584 (2013).Article 

Google Scholar 
17.Keller, R. P., Lodge, D. M. & Finnoff, D. C. Risk assessment for invasive species produces net bioeconomic benefits. Proc. Natl. Acad. Sci. 104, 203–207 (2007).ADS 
CAS 
PubMed 
Article 

Google Scholar 
18.Hulme, P. E., Nentwig, W., Pyšek, P. & Vilà, M. Common market, shared problems: Time for a coordinated response to biological invasions in Europe. Neobiota 8, 3–19 (2009).
Google Scholar 
19.Pluess, T. et al. Which factors affect the success or failure of eradication campaigns against Alien species?. PLoS One 7, 48157 (2012).ADS 
Article 
CAS 

Google Scholar 
20.Sakamoto, Y., Kumagai, N. H. & Goka, K. Declaration of local chemical eradication of the Argentine ant: Bayesian estimation with a multinomial-mixture model. Sci. Rep. 7, 1–8 (2017).Article 
CAS 

Google Scholar 
21.Jiménez-Carmona, F., Carpintero, S. & Reyes-López, J. L. The digging-in effect on ant studies with pitfall traps: Influence of type of habitat and sampling time. Entomol. Exp. Appl. 167, 906–914 (2019).Article 

Google Scholar 
22.Paknia, O., Bergmann, T. & Hadrys, H. Some ‘ant’swers: Application of a layered barcode approach to problems in ant taxonomy. Mol. Ecol. Resour. 15, 1262–1274 (2015).CAS 
PubMed 
Article 

Google Scholar 
23.Chen, Y. & Zhou, S. Phylogenetic relationships based on DNA barcoding among 16 species of the Ant Genus Formica (Hymenoptera: Formicidae) from China. J. Insect Sci. 17, 1–7 (2017).CAS 
Article 

Google Scholar 
24.Giraud, T., Pedersen, J. S. & Keller, L. Evolution of supercolonies: Argentine ants of southern Europe. Proc. Natl. Acad. Sci. 99, 6075–6079 (2002).ADS 
CAS 
PubMed 
Article 

Google Scholar 
25.Touyama, Y., Ogata, K. & Sugiyama, T. The Argentine ant, Linepithema humile, in Japan: Assessment of impact on species diversity of ant communities in urban environments. Entomol. Sci. 6, 57–62 (2003).Article 

Google Scholar 
26.Okaue, M. et al. Distribution of the Argentine ant, Linepithema humile, along the Seto Inland Sea, western Japan: Result of surveys in 2003–2005. Entomol. Sci. 10, 337–342 (2007).Article 

Google Scholar 
27.Sunamura, E., Nishisue, K., Terayama, M. & Tatsuki, S. Invasion of four Argentine ant supercolonies into Kobe Port, Japan: Their distributions and effects on indigenous ants (Hymenoptera: Formicidae). Sociobiol. 50, 659–674 (2007).
Google Scholar 
28.Inoue, M. N. et al. Recent range expansion of Argentine ant in Japan. Divers. Distrib. 19, 2937 (2013).Article 

Google Scholar 
29.Park, S. H., Hosoishi, S. & Ogata, K. Long-term impacts of Argentine ant invasion of urban parks in Hiroshima, Japan. J. Ecol. Environ. 37, 123–129 (2014).Article 

Google Scholar 
30.Sugiyama, T. & Onishi, O. Invasion of Argentinean ants in Kyoto City Ari. J. Myrmecol. Soci. Jpn. 32, 127–129 (2009).
Google Scholar 
31.Murakami, K. Exotic ants in PortIsland, Kobe City Ari. J. Myrmecol. Soci. Jpn. 26, 45–46 (2002).
Google Scholar 
32.Ficetola, G. F., Miaud, C., Pompanon, F. & Taberlet, P. Species detection using environmental DNA from water samples. Biol. Lett. 4, 423–425 (2008).PubMed 
PubMed Central 
Article 

Google Scholar 
33.Jerde, C. L., Mahon, A. R., Chadderton, W. L. & Lodge, D. M. “Sight-unseen” detection of rare aquatic species using environmental DNA. Conserv. Lett. 4, 150–157 (2011).Article 

Google Scholar 
34.Minamoto, T., Yamanaka, H., Takahara, T., Honjo, M. N. & Kawabata, Z. Surveillance of fish species composition using environmental DNA. Limnology 13, 193–197 (2012).CAS 
Article 

Google Scholar 
35.Bienert, F. et al. Tracking earthworm communities from soil DNA. Mol. Ecol. 21, 2017–2030 (2012).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
36.Turner, C. R., Uy, K. L. & Everhart, R. C. Fish environmental DNA is more concentrated in aquatic sediments than surface water. Biol. Conserv. 183, 93–102 (2015).Article 

Google Scholar 
37.Goldberg, C. S. et al. Critical considerations for the application of environmental DNA methods to detect aquatic species. Methods Ecol. Evol. 7, 1299–1307 (2016).Article 

Google Scholar 
38.Dejean, T. et al. Improved detection of an alien invasive species through environmental DNA barcoding: The example of the American bullfrog Lithobates catesbeianus. J. Appl. Ecol. 49, 953–959 (2012).Article 

Google Scholar 
39.Takahara, T., Minamoto, T. & Doi, H. Using environmental DNA to estimate the distribution of an invasive fish species in Ponds. PLoS One 8, e56584 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
40.Geerts, A. N., Boets, P., Van den Heede, S., Goethals, P. & Van der Heyden, C. A search for standardized protocols to detect alien invasive crayfish based on environmental DNA (eDNA): A lab and field evaluation. Ecol. Indic. 84, 564–572 (2018).CAS 
Article 

Google Scholar 
41.Valentin, R. E., Fonseca, D. M., Nielsen, A. L., Leskey, T. C. & Lockwood, J. L. Early detection of invasive exotic insect infestations using eDNA from crop surfaces. Front. Ecol. Environ. 16, 265–270 (2018).Article 

Google Scholar 
42.Shogren, A. J. et al. Controls on eDNA movement in streams: Transport, retention, and resuspension. Sci. Rep. 7, 1–11 (2017).CAS 
Article 

Google Scholar 
43.Jo, T., Murakami, H., Yamamoto, S., Masuda, R. & Minamoto, T. Effect of water temperature and fish biomass on environmental DNA shedding, degradation, and size distribution. Ecol. Evol. 9, 1135–1146 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
44.Takahara, T., Minamoto, T., Yamanaka, H., Doi, H. & Kawabata, Z. Estimation of fish biomass using environmental DNA. PLoS One 7, e35868 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
45.Tillotson, M. D. et al. Concentrations of environmental DNA (eDNA) reflect spawning salmon abundance at fine spatial and temporal scales. Biol. Conserv. 220, 1–11 (2018).Article 

Google Scholar 
46.Sakata, M. K. et al. Sedimentary eDNA provides different information on timescale and fish species composition compared with aqueous eDNA. Environ. DNA 2, 505–518 (2020).Article 

Google Scholar 
47.Pietramellara, G. et al. Extracellular DNA in soil and sediment: Fate and ecological relevance. Biol. Fert. Soils 45, 219–235 (2009).CAS 
Article 

Google Scholar 
48.Matisoo-Smith, E. et al. Recovery of DNA and pollen from New Zealand lake sediments. Quat. Int. 184, 139–149 (2008).Article 

Google Scholar 
49.Foucher, A. et al. Persistence of environmental DNA in cultivated soils: Implication of this memory effect for reconstructing the dynamics of land use and cover changes. Sci. Rep. 10, 10502 (2020).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
50.Sirois, S. H. & Buckley, D. H. Factors governing extracellular DNA degradation dynamics in soil. Environ. Microbiol. Rep. 11, 173–184 (2019).CAS 
PubMed 
Article 

Google Scholar 
51.Ogram, A. V., Mathot, M. L., Harsh, J. B., Boyle, J. & Pettigrew, C. A. Effects of DNA polymer length on its adsorption to soils. Appl. Environ. Microbiol. 60, 393–396 (1994).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
52.Braid, M. D., Daniels, L. M. & Kitts, C. L. Removal of PCR inhibitors from soil DNA by chemical flocculation. J. Microbiol. Methods 52, 389–393 (2003).CAS 
PubMed 
Article 

Google Scholar 
53.Shogren, A. J. et al. Modelling the transport of environmental DNA through a porous substrate using continuous flow-through column experiments. J. R. Soc. Interface 13, 20160290 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
54.Oliverio, A. M., Gan, H., Wickings, K. & Fierer, N. A DNA metabarcoding approach to characterize soil arthropod communities. Soil Biol. Biochem. 125, 37–43 (2018).CAS 
Article 

Google Scholar 
55.Hayami, K. et al. Effects of sampling seasons and locations on fish environmental DNA metabarcoding in dam reservoirs. Ecol. Evol. 10, 5354–5367 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
56.Sato, K. et al. Relationship among establishment durations, kin relatedness, aggressiveness, and distance between populations of eight invasive Argentine Ant (Hymenoptera: Formicidae) supercolonies in Japan. J. Econ. Entomol. 110, 1676–1684 (2017).CAS 
PubMed 
Article 

Google Scholar 
57.Nakajima, S. et al. A list of ant species on the Yard of Kyoto Prefectural Institute of Public Health and Environment, Fushimi, Kyoto Annual report of Kyoto Pref. Inst. Hyg. Environ. Sci. 58, 47–50 (2013).
Google Scholar 
58.Nakajima, S. et al. Time saving improvement to the time-unit sampling method for quantitative surveys of ant faunas. Jpn. J. Environ. Entomol. Zool. 24, 39–50 (2013).
Google Scholar 
59.Kouduka, M. et al. A new DNA extraction method by controlled alkaline treatments from consolidated subsurface sediments. FEMS Microbiol. Lett. 326, 47–54 (2012).CAS 
PubMed 
Article 

Google Scholar 
60.The eDNA Society. Environmental DNA Sampling and Experiment Manual (ver. 2.1). https://ednasociety.org/eDNA_manual_Eng_v2_1_3b.pdf (2019).61.R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. Retrieved from https://www.R-project.org/ (2018). More