Philippa Kaur

Soils and their microbial communities
The soil properties shown in Table 1 indicate variation in soil texture across the Broadbalk field, with less clay in the N0 and FYM plots, situated on the north side of Broadbalk field compared to N6 and woodland towards the south side. The soil pH ranged from 7.1 to 8.2, lowest in the mineral-nitrogen fertilized soil N6 and highest in the N0 soil that received no N fertilizer. The bulk density of woodland soil is much lower and the % SOC much higher compared to other soils; the FYM soil has lower bulk density and higher % SOC than the other arable soils. The ratio of SOC:total N was approximately 10:1 in the arable soils and 13:1 in the woodland soil.
The community structure of bacteria and archaea revealed by 16S rRNA amplicon sequencing of metagenomic DNA extracted from the soil samples, at collection from the field, shows significant differences, and distinct separation on a NMDS plot (Fig. 1). Of the 14 phyla (sub-phyla for the Proteobacteria) comprising  > 0.1% of the community in at least one of the soils, only the δ-Proteobacteria did not show significantly different (P ≤ 0.05) mean abundance in at least one soil, according to ANOVA (Fig. 1). For example, the woodland soil has more α-Proteobacteria and Verrucomicrobia but fewer Thaumarchaeota (archaea) and β-Proteobacteria than the other treatments. Both the FYM and woodland soil have more γ-Proteobacteria and fewer Gemmatimonadetes; the FYM soil has more Firmicutes than the other soils (Fig. 1).
Figure 1

Relative mean abundance of prokaryotic phyla/sub-phyla in soils of origin on collection from the field. Phyla with at least 0.1% of the total community present in at least one soil treatment are included. Proteobacteria sub-phyla: a = alpha, b = beta, d = delta, g = gamma; s.e.d. for each group is shown; letters indicate mean significantly different means within each group (P = 0.05, according to Tukey’s post-hoc test on ANOVA). Insert top right shows NMDS plot of OTU for prokaryotic communities – PERMANOVA F = 9.477, P (same) = 0.0001.

Full size image

16S rRNA and denitrification gene abundance
At the end of the experiment, DNA was extracted and amplified from all samples but sufficient RNA for further analysis was obtained only from the FYM and woodland soils which contained more organic matter and larger microbial communities. ANOVA comparing the abundance for each set of genes and transcripts measured using qPCR showed that the soil of origin had a significant influence in all cases (Table 2). However, other factors (presence/absence of wheat and addition or not of N-fertilizer) and interactions between them were not significant, except for nosZI which was significantly influenced by the plant. Bacterial abundance indicated by 16S rRNA gene copy number was 2 × 109 g−1 soil in the N0 and N6 soils and significantly higher in the FYM and woodland soil, 5 × 109 and 7 × 109 copies g−1, respectively (Fig. 2). This pattern of relative abundance was seen for nirK (7 × 108–4 × 109 copies), nosZI (5 × 107–2 × 108 copies) and nosZII (4 × 106–1 × 107 copies g−1 soil). The exception was nirS where N0, N6 and woodland soil had similar gene abundance (1 × 107 copies g−1 soil) and FYM significantly more with 4 × 107 copies g−1soil (Fig. 2). The ratio nirK:nirS gene copies in the woodland soil was 300:1, significantly more than the mean of 55:1 in the arable soils (F3,32 = 102.63, P  More

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript. More

1.
Krishnan, G., Holshouser, D. L. & Nissen, S. J. Weed control in soybean (Glycine max) with green manure crops. Weed Technol. 12, 97–102 (1998).
Article  Google Scholar 
2.
Parsons, W. T., Parsons, W. T. & Cuthbertson, E. G. Noxious weeds of Australia (CSIRO Publishing, Canberra, 2001).
Google Scholar 

3.
Whish, J. P. M., Sindel, B. M., Jessop, R. S. & Felton, W. L. The effect of row spacing and weed density on yield loss of chickpea. Aust. J. Agric. Res. 53, 1335–1340 (2002).
Article  Google Scholar 

4.
McDonald, G. K., Hollaway, K. L. & McMurray, L. Increasing plant density improves weed competition in lentil (Lens culinaris). Aust. J. Exp. Agric. 47, 48–56 (2007).
Article  Google Scholar 

5.
Manalil, S., Werth, J., Jackson, R., Chauhan, B. S. & Preston, C. An assessment of weed flora 14 years after the introduction of glyphosate-tolerant cotton in Australia. Crop Pasture Sci. 68, 773–780 (2017).
Article  Google Scholar 

6.
Adkins, S. W. et al. Weeds resistant to chlorsulfuron and atrazine from the north-east grain region of Australia. Weed Res. 37, 343–349 (1997).
Article  Google Scholar 

7.
Manalil, S., Ali, H. H. & Chauhan, B. S. Germination ecology of turnip weed (Rapistrum rugosum (L.) All.) in the northern regions of Australia. PLoS ONE 13, e0201023 (2018).
PubMed  Article  CAS  PubMed Central  Google Scholar 

8.
Llewellyn, R. et al. Impact of Weeds in Australian Grain Production (Grains Research and Development Corporation, Canberra, 2016).
Google Scholar 

9.
Chauhan, B. S., Gill, G. & Preston, C. African mustard (Brassica tournefortii) germination in southern Australia. Weed Sci. 54, 891–897 (2006).
CAS  Article  Google Scholar 

10.
Chauhan, B. S., Gill, G. & Preston, C. Factors affecting turnipweed (Rapistrum rugosum) seed germination in southern Australia. Weed Sci. 54, 1032–1036 (2006).
CAS  Article  Google Scholar 

11.
Wilson, B. J. & Wilson, J. T. Effect of time of seedling emergence on seed production and time to flowering of eight weeds, in 6th Australian Weeds Conference.’ Gold Coast, Queensland, Australia, Queensland Weed Society. https://www.cawsorgau/awc/1981/awc198110351.pdf (1981).

12.
Minnich, R. A. & Sanders, A. C. Brassica tournefortii Gouan. In Invasive Plants of California’s Wildlands (eds Bossard, C. C. et al.) 68–72 (University of California Press, Berkeley, 2000).
Google Scholar 

13.
Patterson, D. T. Research on exotic weeds. In Exotic Plant Pests and North American Agriculture (eds Wilson, C. L. & Graham, C. L.) 381–393 (Academic Press, New York, 1983).
Google Scholar 

14.
Meekins, J. F. & McCarthy, B. C. Competitive ability of Alliaria petiolata (garlic mustard, Brassicaceae), an invasive, nonindigenous forest herb. Int J Plant Sci. 160, 743–752 (1999).
Article  Google Scholar 

15.
Heap, I. International Survey of Herbicide Resistant Weeds. https://weedscience.org (NOAA, 2019).

16.
Manalil, S. & Chauhan, B. S. Interference of turnipweed (Rapistrum rugosum) and Mexican pricklepoppy (Argemone mexicana) in wheat. Weed Sci. 67, 666–672 (2019).
Article  Google Scholar 

17.
Gill, G. S. & Davidson, R. M. Weed interference. In Australian Weed Management Systems (ed. Sindel, B. M.) 61–80 (R.G. and F.J. Richardson, Meredith, 2000).
Google Scholar 

18.
Mahajan, G., Mutti, N. K., Jha, P., Walsh, M. & Chauhan, B. S. Evaluation of dormancy breaking methods for enhanced germination in four biotypes of Brassica tournefortii. Sci. Rep. 8, 17103 (2018).
ADS  PubMed  Article  CAS  PubMed Central  Google Scholar 

19.
Trader, M. R., Brooks, M. L. & Draper, J. V. Seed production by the non-native Brassica tournefortii (Sahara mustard) along desert roadsides. Madrono 53, 313–321 (2006).
Article  Google Scholar 

20.
Lindström, J. & Kokko, H. Cohort effects and population dynamics. Ecol. Lett. 5, 338–344 (2002).
Google Scholar 

21.
Hussain, S., Khaliq, A., Matloob, A., Fahad, S. & Tanveer, A. Interference and economic threshold level of little seed canary grass in wheat under different sowing times. Environ. Sci. Pollut. Res. 22, 441–449 (2015).
Google Scholar 

22.
Estorninos, L. E., Gealy, D. R., Gbur, E. E., Talbert, R. E. & McClelland, M. R. Rice and red rice interference. II. Rice response to population densities of three red rice (Oryza sativa) ecotypes. Weed Sci. 53, 683–689 (2005).
CAS  Google Scholar 

23.
Knezevic, S. Z., Weise, S. F. & Swanton, C. J. Interference of redroot pigweed (Amaranthus retroflexus) in corn (Zea mays). Weed Sci. 42, 568–573 (1994).
CAS  Article  Google Scholar 

24.
Bosnic, A. C. & Swanton, C. J. Influence of barnyardgrass (Echinochloa crus-galli) time of emergence and density on corn (Zea mays). Weed Sci. 45, 276–282 (1997).
CAS  Article  Google Scholar 

25.
Willenborg, C. J., May, W. E., Gulden, R. H., Lafond, G. P. & Shirtliffe, S. J. Influence of wild oat (Avena fatua) relative time of emergence and density on cultivated oat yield, wild oat seed production, and wild oat contamination. Weed Sci. 53, 342–352 (2005).
CAS  Article  Google Scholar 

26.
Borchert, R. et al. Photoperiodic induction of synchronous flowering near the Equator. Nature 433, 627 (2005).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

27.
Morellato, L. P. C. et al. Linking plant phenology to conservation biology. Biol. Conserv. 195, 60–72 (2016).
Article  Google Scholar 

28.
Ekeleme, F., Forcella, F., Archer, D. W., Akobundu, I. O. & Chikoye, D. Seedling emergence model for tropic ageratum (Ageratum conyzoides). Weed Sci. 53, 55–61 (2005).
CAS  Article  Google Scholar 

29.
Hock, S. M., Knezevic, S. Z., Martin, A. R. & Lindquist, J. L. Soybean row spacing and weed emergence time influence weed competitiveness and competitive indices. Weed Sci. 54, 38–46 (2006).
CAS  Article  Google Scholar 

30.
Spaunhorst, D. J. et al. Phenology of five Palmer amaranth (Amaranthus palmeri) populations grown in northern Indiana and Arkansas. Weed Sci. 66, 457–469 (2018).
Article  Google Scholar 

31.
Mobli, A., Mijani, S., Ghanbari, A. & Rastgoo, M. Seed germination and emergence of two flax-leaf alyssum (Alyssum linifolium Steph. ex. Willd.) populations in response to environmental factors. Crop Pasture Sci. 70, 807–813 (2019).
Article  Google Scholar 

32.
Faraji, A. Determination of phenological response of spring canola (Brassica napus L.) genotypes to sowing date, temperature and photoperiod. Seed Plant Prod. J. 1, 25–41 (2010).
Google Scholar 

33.
Salehian, H., Jamshidi, M., Karamzade, H. & Modaresi, H. An investigation of the phenology some weed species. Int. J. Agron. Plant Prod. 4, 1511–1520 (2013).
Google Scholar 

34.
Baskin, J. M. & Baskin, C. C. A classification system for seed dormancy. Seed Sci Res. 14, 1–16 (2004).
ADS  Article  Google Scholar 

35.
Baskin, J. M., Baskin, C. C. & Li, X. Taxonomy, anatomy and evolution of physical dormancy in seeds. Plant Species Biol. 15, 139–152 (2000).
Article  Google Scholar 

36.
Huang, J. Z. et al. Effect of temperature and photoperiod on the phenological development of wild mustard (Sinapis arvensis L.). Field Crops Res. 70, 75–86 (2001).
Article  Google Scholar 

37.
Hatfield, J. L. et al. Climate impacts on agriculture: implications for crop production. Agron. J. 103, 351–370 (2011).
Article  Google Scholar 

38.
Hatfield, J. L. & Prueger, J. H. Temperature extremes: effect on plant growth and development. Weather Clim. Extremes 10, 4–10 (2015).
Article  Google Scholar 

39.
Ghersa, C. M. & Holt, J. S. Using phenology prediction in weed management: a review. Weed Res. 35, 461–470 (1995).
Article  Google Scholar 

40.
Deen, W., Hunt, T. & Swanton, C. J. Influence of temperature, photoperiod, and irradiance on the phenological development of common ragweed (Ambrosia artemisiifolia). Weed Sci. 46, 555–560 (1998).
CAS  Google Scholar 

41.
O’Donovan, J. T., Remy, E. A. D. S., O’Sullivan, P. A., Dew, D. A. & Sharma, A. K. Influence of the relative time of emergence of wild oat (Avena fatua) on yield loss of barley (Hordeum vulgare) and wheat (Triticum aestivum). Weed Sci. 33, 498–503 (1985).
Google Scholar 

42.
Chikoye, D., Weise, S. F. & Swanton, C. J. Influence of common ragweed (Ambrosia artemisiifolia) time of emergence and density on white bean (Phaseolus vulgaris). Weed Sci. 43, 375–380 (1995).
CAS  Google Scholar 

43.
Hartzler, R. G., Battles, B. A. & Nordby, D. Effect of common waterhemp (Amaranthus rudis) emergence date on growth and fecundity in soybean. Weed Sci. 52, 242–245 (2004).
CAS  Google Scholar 

44.
Gioria, M. & Pyšek, P. Early bird catches the worm: germination as a critical step in plant invasion. Biol. Invasions. 19(4), 1055–1080 (2017).
Google Scholar 

45.
Bajwa, A. A., Chauhan, B. S. & Adkins, S. W. Germination ecology of two Australian biotypes of ragweed parthenium (Parthenium hysterophorus) relates to their invasiveness. Weed Sci. 66, 62–70 (2018).
Article  Google Scholar 

46.
Harrison, S. K. Interference and seed production by common lambsquarters (Chenopodium album) in soybeans (Glycine max). Weed Sci. 38, 113–118 (1990).
ADS  Article  Google Scholar 

47.
Chauhan, B. S. & Johnson, D. E. Implications of narrow crop row spacing and delayed Echinochloa colona and Echinochloa crus-galli emergence for weed growth and crop yield loss in aerobic rice. Field Crops Res. 117, 177–182 (2010).
Article  Google Scholar 

48.
Simard, M. J. & Benoit, D. L. Potential pollen and seed production from early-and late-emerging common ragweed in corn and soybean. Weed Technol. 26, 510–516 (2012).
Article  Google Scholar 

49.
Fahad, S. et al. Consequences of narrow crop row spacing and delayed Echinochloa colona and Trianthema portulacastrum emergence for weed growth and crop yield loss in maize. Weed Res. 54, 475–483 (2014).
Article  Google Scholar 

50.
Thullen, R. J. & Keeley, P. E. Germination, growth, and seed production of Ipomoea hederacea when planted at monthly intervals. Weed Sci. 31, 837–840 (1983).
Article  Google Scholar  More

Site description and sampling
This experimental site was established in July 2009 at the Kessler Atmospheric and Ecological Field Station (KAEFS) in the US Great Plains in McClain County, Oklahoma (34̊ 59ʹN, 97̊ 31ʹW)14,48. Experimental design and site description were described in detail previously25. Briefly, Ambrosia trifida, Solanum carolinense and Euphorbia dentate belonging to C3 forbs, and Tridens flavus, Sporobolus compositus and Sorghum halapense belonging to C4 grasses are dominant in the site25,48. Annual mean temperature is 16.3 °C and annual precipitation is 914 mm, based on Oklahoma Climatological Survey data from 1948 to 1999. The soil type of this site is Port–Pulaski–Keokuk complex with 51% of sand, 35% of silt and 13% of clay, which is a well-drained soil that is formed in loamy sediment on flood plains. The soil has a high available water holding capacity (37%), neutral pH and 1.2 g cm−3 bulk density with 1.9% total organic matter and 0.1% total nitrogen (N)25,48. Four blocks were used in the field site experiment, in which warming is a primary factor. Two levels of warming (ambient and +3 °C) were set for four pairs of 2.5 m × 1.75 m plots by utilizing a real or dummy infrared radiator (Kalglo Electronics, Bethlehem, PA, USA). In the warmed plots, a real infrared radiator was suspended 1.5 m above the ground, and the dummy infrared radiator was suspended to simulate a shading effect of the device in the control plots.
In this study, eight surface (0–15 cm) soil samples, four from the warmed and four from the control plots, were collected annually at approximately the date of peak plant biomass (September or October) from 2010 to 2016. Three soil cores (2.5 cm diameter × 15 cm depth) were taken by using a soil sampler tube in each plot and composited to have enough samples for soil chemistry, microbiology and molecular biology analyses. A total of 56 soil samples were analyzed in this study.
Environmental and soil chemical measurements
Precipitation data were obtained from the Oklahoma Mesonet Station (Washington Station)48 located 200 m away from our experiment site, and 12-month version of the standardized precipitation-evapotranspiration index (SPEI-12) was used as annual drought index49. Air temperature, soil temperature and volumetric soil water content were described in detail previously25. Specifically, air temperature and soil temperature at the depth of 7.5 cm in the center of each field plot were measured by using Constantan-copper thermocouples wired to a Campbell Scientific CR10x data logger (Campbell Scientific, UT, USA). A portable time domain reflectometer (Soil Moisture Equipment Corp.) was used to measure soil moisture from the soil surface to a 15-cm depth once or twice a month. Three measurements of soil moisture were performed in each plot and the average of three technical replicates were used in further analyses.
All soil samples were analyzed to determine soil total organic carbon (TOC), total nitrogen (TN), soil nitrate (NO3−) and ammonia (NH4+) by the Soil, Water, and Forage Analytical Laboratory at Oklahoma State University (Stillwater, OK, USA). Soil pH was measured using a pH meter with a calibrated combined glass electrode50.
Aboveground plant communities
Aboveground plant community investigations were annually conducted at peak biomass (usually September)48,51. Aboveground plant biomass, separated into C3 and C4 species, was indirectly estimated by a modified pin-touch method48,51. Detailed description of biomass estimation is provided by Sherry et al.52. A pin frame used in this study is 1 m long and have 10 pins 10 cm apart at 30° from vertical. Pins with a 0.75 m length were raised within the frame to count hits up to 1 m high (hits over 1 m are negligible at this site). The pin frame was placed in the center of each plot to record the contact numbers of the pins separately with C3 and C4 plants (e.g., leaves and stems). The contact numbers of C3 and C4 plants were then used to estimate plant biomass using calibration equations derived from calibration plots, which were located near the experimental plots. Biomass in the calibration plots was clipped at a height of 10 cm above the ground at approximately the date of peak plant biomass (September or October). All of the species in plant community within each plot were identified to estimate species richness. Clipped plant materials were oven-dried and then correlated with the total contact number. C3 and C4 plant biomasses were estimated by using the calibration equation of contact number and plant biomass. All of the species within each plot were identified to estimate species richness of plants.
Ecosystem C fluxes and soil respiration
Ecosystem C fluxes and soil respirations were measured once or twice a month between 10:00 and 15:00 (local time) from January 2010 to December 2016 by following previous methods14,48. One square aluminum frame (0.5 m × 0.5 m) was inserted in the soil at 2 cm depth in each plot to provide a flat base between the soil surface and the CO2 sampling chamber. NEE and ecosystem respiration (ER) were measured using LI-6400 portable photosynthesis system (LI-COR). Gross primary productivity (GPP) was estimated as the difference between NEE and ER. Meanwhile, soil surface respiration was monthly measured using a LI-8100A soil flux system attached to a soil CO2 flux chamber (LI-COR). Measurements were taken above a PVC collar (80 cm2 in area and 5 cm in depth) and a PVC tube (80 cm2 in area and 70 cm in depth) in each plot. The PVC tube was permanently fixed on the ground to cut off old plant roots and prevent new roots from growing inside the tube. Any aboveground parts of living plants were removed from the PVC tubes and collars before each measurement. The CO2 efflux measured above the PVC tubes represented heterotrophic respiration (Rh) from soil microbes, while that measured above the PVC collars represented soil total respiration (Rt) including heterotrophic and autotrophic respiration (Rh and Ra) from soil microbes and plant root, respectively.
Soil decomposition rate
Weighted cellulose filter paper (Whatman CAT No. 1442-090) was placed into fiberglass mesh bags and placed vertically at 0–10 cm soil depth in each plot in March 2016. All of decomposition bags were collected back in September 2016, rinsed and dried at 60 °C for weighing. The percentage of mass loss was calculated to represent soil decomposition rate.
Molecular analyses of soil samples
The C substrate utilization patterns of soil microbial communities in 2016 were analyzed by BIOLOG EcoPlateTM (BIOLOG). The BIOLOG EcoPlateTM contains 31 of the most useful labile carbon sources for soil community analysis, which are repeated three times in each plate. In this study, the plates with diluted soil supernatant (0.5 g soil with 45 mL 0.85% NaCl) were incubated in a BIOLOG OmniLog PM System at 25 °C for 4.5 days. The color change of each well was shown as absorbance curve. The net area under the absorbance versus time curve was calculated to represent physiological activity of various C sources53. The average value from three replicates was used for analyses in this study.
Soil total DNA was extracted from 1.5 g soil by freeze-grinding and SDS-based lysis54, and purified with a MoBio PowerSoil DNA isolation kit (MoBio Laboratories)25. Then, 10 ng DNA per sample were used for library construction and amplicon sequencing. Amplicons sequencing was performed with cautions in terms of experimental preparations and data analyses to ensure sequence representativeness and semi-quantitative nature55. The V4 region of bacterial and archaeal 16S rRNA genes were amplified with the primer set 515F (5ʹ-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5ʹ-GGACTACHVGGGTWTCTAAT-3ʹ), and fungal ITSs between 5.8S and 28S rRNA genes were amplified with the primer set ITS7F (5ʹ-GTGARTCATCGARTCTTTG-3ʹ) and ITS4R (5ʹ-TCCTCCGCTTATTGATATGC-3ʹ). PCR products from different samples were sequenced on a MiSeq platform (Illumina, Inc.) using 2 × 250 pair-end sequencing kit. Raw sequences were submitted to our Galaxy sequence analysis pipeline (http://zhoulab5.rccc.ou.edu:8080) to further analyze according to the protocol in the pipeline25. Finally, OTUs were clustered by UPARSE56 at 97% identity for both 16S rRNA gene and ITS. All sequences were randomly resampled to 30,000 sequences for 16S rRNA gene and 10,000 sequences for ITS per sample. Representative sequences of OTUs were annotated taxonomically by the Ribosomal Database Project (RDP) Classifier with 50% confidence estimates.
GeoChip 5.0 M, a functional gene array57, was used for all 56 samples from 2010 to 2016. GeoChip hybridization, scanning and data processing were performed in the Institute for Environmental Genomics, University of Oklahoma57,58. Specifically, 800 ng of purified soil DNA of each sample was mixed with 5.5 µl random primers (Life Technologies, random hexamers, 3 µg/µl), diluted with nuclease-free water to 35 µl, heated to 99 °C for 5 min, and placed on ice immediately. The labeling master mix (15 µl), including 0.5 µl of Cy-3 dUTP (25 nM; GE Healthcare), 2.5 µl of dNTP (2.5 mM dTTP, 5 mM dAGC-TP), 1 µl of Klenow (imer; San Diego, CA; 40 U ml−1), 5 µl Klenow buffer, and 2.5 µl of water, was added in the sample mixed solution. The samples were incubated at 37 °C for 6 h in a thermocycler, and then incubated at 95 °C for 3 min to inactivate the enzyme. Subsequently, samples were protected from the light as much as possible. Labeled DNA was cleaned using a QIAquick purification kit (Qiagen) according the manufacturer’s instructions and then dried thoroughly in a SpeedVac (45 °C, 45 min; ThermoSavant).
Labeled DNA was resuspended into 27.5 µl of DNase-free water, and then mixed completely with 99.4 µl of hybridization solution, containing 63.5 µl of formamide (10% final concentration), 2 × HI-RPM hybridization buffer, 12.7 µl of 10 × aCGH blocking agent, 0.05 μg/µl Cot-1 DNA, and 10 pM CORS58. The mixed solution was denatured at 95 °C for 3 min, and then incubated at 37 °C for 30 min. The DNA solution was centrifuged at 6000 × g for 1 min to collect liquid at the bottom of the tube. 110 µl of the solution was pipetted into the center of the well of the gasket slide. The array slide was placed on the gasket slide, sealed using a SureHyb chamber, hybridized at 67 °C for 24 h at 20 rpm in a hybridization oven. After hybridization, slides were washed in room temperature with Wash Buffer 1 (Agilent) and Wash Buffer 2 (Agilent).
The slides were imaged as a Multi-TIFF with a NimbleGen MS200 Microarray Scanner (Roche NimbleGen, Inc., Madison, WI, United States). The raw signals from NimbleGen were submitted to the Microarray Data Manager on our website (http://ieg.ou.edu/microarray), cleaned, normalized and analyzed using the data-analysis pipeline. Briefly, probe quality was assessed, and poor or low signal probes were removed. Probe spots with coefficient of variance (CV; probe signal SD/signal) >0.8 were removed. Then, the signal-to-noise ratio (SNR) was calculated. As suggested by Agilent, the average signal of Agilent’s negative control probes within each subarray was used as the background signal for the probes in that subarray instead of the local background typically used. The signal intensity for each spot was corrected by subtracting the background signal intensity. If the net difference was30%, aligned length >20 a.a., and e-value More

1.
Thuiller, W. et al. Consequences of climate change on the tree of life in Europe. Nature 470, 531–534 (2011).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 
2.
IPCC. Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (Cambridge University Press, Cambridge, 2013).
Google Scholar 

3.
Sohngen, B. & Tian, X. Global climate change impacts on forests and markets. For. Policy Econ. 72, 18–26 (2016).
Article  Google Scholar 

4.
Fawzy, S., Osman, A. I., Doran, J. & Rooney, D. W. Strategies for mitigation of climate change: a review. Environ. Chem. Lett. https://doi.org/10.1007/s10311-020-01059-w (2020).
Article  Google Scholar 

5.
Luyssaert, S. et al. Trade-offs in using European forests to meet climate objectives. Nature 562, 259 (2018).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

6.
Seidl, R., Schelhaas, M.-J., Rammer, W. & Verkerk, P. J. Increasing forest disturbances in Europe and their impact on carbon storage. Nat. Clim. Change 4, 806–810 (2014).
ADS  CAS  Article  Google Scholar 

7.
Osman, A. I. Mass spectrometry study of lignocellulosic biomass combustion and pyrolysis with NOx removal. Renew. Energy 146, 484–496 (2020).
CAS  Article  Google Scholar 

8.
Woolf, D., Amonette, J. E., Street-Perrott, F. A., Lehmann, J. & Joseph, S. Sustainable biochar to mitigate global climate change. Nat. Commun. 1, 56 (2010).
ADS  PubMed  Article  CAS  PubMed Central  Google Scholar 

9.
Qambrani, N. A., Rahman, Md. M., Won, S., Shim, S. & Ra, C. Biochar properties and eco-friendly applications for climate change mitigation, waste management, and wastewater treatment: a review. Renew. Sustain. Energy Rev. 79, 255–273 (2017).
CAS  Article  Google Scholar 

10.
Choi, S. H. & Manousiouthakis, V. I. On the carbon cycle impact of combustion of harvested plant biomass vs. fossil carbon resources. Comput. Chem. Eng. 140, 106942 (2020).
CAS  Article  Google Scholar 

11.
Pan, Y. et al. A large and persistent carbon sink in the world’s forests. Science 333, 988–993 (2011).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

12.
Laiho, R. & Laine, J. Tree stand biomass and carbon content in an age sequence of drained pine mires in southern Finland. For. Ecol. Manag. 93, 161–169 (1997).
Article  Google Scholar 

13.
Martin, A. R. & Thomas, S. C. A reassessment of carbon content in tropical trees. PLoS ONE 6, e23533 (2011).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

14.
Jagodziński, A. M., Dyderski, M. K., Gęsikiewicz, K. & Horodecki, P. Tree- and stand-level biomass estimation in a Larix decidua Mill. Chronosequence. Forests 9, 587 (2018).
Article  Google Scholar 

15.
Teobaldelli, M., Somogyi, Z., Migliavacca, M. & Usoltsev, V. A. Generalized functions of biomass expansion factors for conifers and broadleaved by stand age, growing stock and site index. For. Ecol. Manag. 257, 1004–1013 (2009).
Article  Google Scholar 

16.
Forrester, D. I. et al. Generalized biomass and leaf area allometric equations for European tree species incorporating stand structure, tree age and climate. For. Ecol. Manag. 396, 160–175 (2017).
Article  Google Scholar 

17.
Jagodziński, A. M. et al. How do tree stand parameters affect young Scots pine biomass? Allometric equations and biomass conversion and expansion factors. For. Ecol. Manag. 409, 74–83 (2018).
Article  Google Scholar 

18.
Picard, N., Saint-Andre, L. & Henry, M. Manual for building tree volume and biomass allometric equations: from field measurement to prediction. (Food and Agricultural Organization of the United Nations and Centre de Coopération Internationale en Recherche Agronomique pour le Développement, 2012).

19.
Grote, R., Schuck, J., Block, J. & Pretzsch, H. Oberirdische holzige Biomasse in Kiefern-/Buchen- und Eichen-/Buchen-Mischbeständen. Forstwiss. Cent. Ver. Mit Tharandter Forstl. Jahrb. 122, 287–301 (2003).
Google Scholar 

20.
Zhang, L., Zhang, Y., Wang, H., Zou, J. & Siemann, E. Chinese Tallow trees (Triadica sebifera) from the invasive range outperform those from the native range with an active soil community or phosphorus fertilization. PLoS ONE 8, e74233 (2013).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

21.
Frouz, J. et al. Development of canopy cover and woody vegetation biomass on reclaimed and unreclaimed post-mining sites. Ecol. Eng. 84, 233–239 (2015).
Article  Google Scholar 

22.
Mangla, S., Sheley, R. L., James, J. J. & Radosevich, S. R. Intra and interspecific competition among invasive and native species during early stages of plant growth. Plant Ecol. 212, 531–542 (2011).
Article  Google Scholar 

23.
Gómez-García, E., Crecente-Campo, F. & Diéguez-Aranda, U. Tarifas de biomasa aérea para abedul (Betula pubescens Ehrh.) y roble (Quercus robur L.) en el noroeste de España. Madera Bosques 19, 71–91 (2013).
Google Scholar 

24.
Albert, K., Annighöfer, P., Schumacher, J. & Ammer, C. Biomass equations for seven different tree species growing in coppice-with-standards forests in Central Germany. Scand. J. For. Res. 29, 210–221 (2014).
Article  Google Scholar 

25.
Repola, J. Biomass equations for birch in Finland. Silva Fenn. 42, 605–624 (2008).
Article  Google Scholar 

26.
Uri, V. et al. Biomass production and carbon sequestration in a fertile silver birch (Betula pendula Roth) forest chronosequence. For. Ecol. Manag. 267, 117–126 (2012).
Article  Google Scholar 

27.
Sellin, A. et al. Elevated air humidity affects hydraulic traits and tree size but not biomass allocation in young silver birches (Betula pendula). Front. Plant Sci. 6, 1–10 (2015).
Article  Google Scholar 

28.
Xiao, C.-W. et al. Above- and belowground biomass and net primary production in a 73-year-old Scots pine forest. Tree Physiol. 23, 505–516 (2003).
PubMed  Article  PubMed Central  Google Scholar 

29.
Yuste, J. C. et al. Contrasting net primary productivity and carbon distribution between neighboring stands of Quercus robur and Pinus sylvestris. Tree Physiol. 25, 701–712 (2005).
CAS  Article  Google Scholar 

30.
Petersson, H. & Ståhl, G. Functions for below-ground biomass of Pinus sylvestris, Picea abies, Betula pendula and Betula pubescens in Sweden. Scand. J. For. Res. 21, 84–93 (2006).
Article  Google Scholar 

31.
Repola, J. & AhnlundUlvcrona, K. Modelling biomass of young and dense Scots pine (Pinus sylvestris L.) dominated mixed forests in northern Sweden. Silva Fenn. 48, 1190 (2014).
Article  Google Scholar 

32.
Ozolinčius, R., Mikšys, V. & Stakénas, V. Above-ground phytomass and light regime in Norway spruce stands planted with different initial density. Biomass Bioenergy 11, 201–206 (1996).
Article  Google Scholar 

33.
Johansson, T. Sprouting ability and biomass production of downy and silver birch stumps of different diameters. Biomass Bioenergy 32, 944–951 (2008).
Article  Google Scholar 

34.
Pajtík, J., Konôpka, B. & Lukac, M. Individual biomass factors for beech, oak and pine in Slovakia: a comparative study in young naturally regenerated stands. Trees 25, 277–288 (2011).
Article  Google Scholar 

35.
Gezici-Koç, Ö, Erich, S. J. F., Huinink, H. P., van der Ven, L. G. J. & Adan, O. C. G. Bound and free water distribution in wood during water uptake and drying as measured by 1D magnetic resonance imaging. Cellulose 24, 535–553 (2017).
Article  CAS  Google Scholar 

36.
Samuelsson, R., Nilsson, C. & Burvall, J. Sampling and GC-MS as a method for analysis of volatile organic compounds (VOC) emitted during oven drying of biomass materials. Biomass Bioenergy 30, 923–928 (2006).
CAS  Article  Google Scholar 

37.
Samuelsson, R., Burvall, J. & Jirjis, R. Comparison of different methods for the determination of moisture content in biomass. Biomass Bioenergy 30, 929–934 (2006).
CAS  Article  Google Scholar 

38.
R Core Team. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, Vienna, 2020).

39.
Jabłoński, M. & Budniak, P. Estimating above-ground woody biomass of forests in Poland for UNECE/FAO and UNFCCC reporting. For. Res. Pap. 75, 277–289 (2014).
Google Scholar 

40.
Claessens, H., Oosterbaan, A., Savill, P. & Rondeux, J. A review of the characteristics of black alder (Alnus glutinosa (L.) Gaertn.) and their implications for silvicultural practices. Forestry 83, 163–175 (2010).
Article  Google Scholar 

41.
Horodecki, P. & Jagodziński, A. M. Tree species effects on litter decomposition in pure stands on afforested post-mining sites. For. Ecol. Manag. 406, 1–11 (2017).
Article  Google Scholar 

42.
Horodecki, P., Nowiński, M. & Jagodziński, A. M. Advantage of mixed tree stands in restoration of upper soil layers on post-mining sites: a five-year leaf litter decomposition experiment. Land Degrad. Dev. 30, 3–13 (2019).
Article  Google Scholar 

43.
Repola, J. Models for vertical wood density of Scots pine, Norway spruce and birch stems, and their application to determine average wood density. Silva Fenn. 40, 4 (2006).
Google Scholar 

44.
Poorter, H. et al. How does biomass distribution change with size and differ among species? An analysis for 1200 plant species from five continents. New Phytol. 208, 736–749 (2015).
PubMed  PubMed Central  Article  Google Scholar 

45.
Annighöfer, P. et al. Species-specific and generic biomass equations for seedlings and saplings of European tree species. Eur. J. For. Res. 135, 313–329 (2016).
Article  Google Scholar 

46.
Schepaschenko, D. et al. A dataset of forest biomass structure for Eurasia. Sci. Data 4, 201770 (2017).
Article  Google Scholar 

47.
Muukkonen, P. Generalized allometric volume and biomass equations for some tree species in Europe. Eur. J. For. Res. 126, 157–166 (2007).
Article  Google Scholar 

48.
Chave, J. et al. Towards a worldwide wood economics spectrum. Ecol. Lett. 12, 351–366 (2009).
PubMed  Article  PubMed Central  Google Scholar 

49.
Flores, O. et al. An evolutionary perspective on leaf economics: phylogenetics of leaf mass per area in vascular plants. Ecol. Evol. 4, 2799–2811 (2014).
PubMed  PubMed Central  Article  Google Scholar 

50.
Paź-Dyderska, S., Dyderski, M. K., Nowak, K. & Jagodziński, A. M. On the sunny side of the crown: quantification of intra-canopy SLA variation among 179 taxa. For. Ecol. Manag. 472, 118254 (2020).
Article  Google Scholar 

51.
Zanne, A. E. et al. Global Wood Density Database. (2009) https://datadryad.org/handle/10255/dryad.235.

52.
Kleyer, M. et al. The LEDA Traitbase: a database of life-history traits of the Northwest European flora. J. Ecol. 96, 1266–1274 (2008).
Article  Google Scholar 

53.
FAO. Global Forest Resources Assessment. (UN Food and Agriculture Organization, 2015).

54.
Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).
Article  Google Scholar 

55.
Kuznetsova, A., Brockhoff, P. B. & Christensen, R. H. B. lmerTest: tests in linear mixed effects models. J. Stat. Softw. 82, 13 (2017).
Article  Google Scholar 

56.
Bartoń, K. MuMIn: Multi-Model Inference. (2017). https://cran.r-project.org/package=MuMIn (Accessed 12 Aug 2020).

57.
Jagodziński, A. M., Dyderski, M. K., Gęsikiewicz, K. & Horodecki, P. Effects of stand features on aboveground biomass and biomass conversion and expansion factors based on a Pinus sylvestris L. chronosequence in Western Poland. Eur. J. For. Res. 138, 673–683 (2019).
Article  Google Scholar 

58.
Jagodziński, A. M., Dyderski, M. K. & Horodecki, P. Differences in biomass production and carbon sequestration between highland and lowland stands of Picea abies (L.) H. Karst. and Fagus sylvatica L. For. Ecol. Manag. 474, 118329 (2020).
Article  Google Scholar  More

Website sampling
We used five different search terms, all translations of ‘reptiles for sale’ (reptile à vendre, reptilien zu verkaufen, 爬虫類の販売, reptil para la venta), on the appropriately localized version of two search engines (Google: https://www.google.com/, https://www.google.fr/, https://www.google.de/, https://www.google.jp/, https://www.google.es/ and Bing: https://www.bing.com/?cc=en, https://www.bing.com/?cc=fr, https://www.bing.com/?cc=de, https://www.bing.com/?cc=jp, https://www.bing.com/?cc=es) to retrieve a list of reptile selling websites, extracting URLs using the XML v.3.99.0.339, assertthat v.0.2.140 and stringr v.1.4.0 packages41 in R v.3.5.342 and R studio v.1.2.133543 (Supplementary Code 1). We completed searches in Firefox44, while signed out of search engine accounts and in a private window to minimise browsing history’s impact on searches.
Once we had generated a list of search result URLs, we manually reviewed each URL’s content—679 websites led to 151 searchable reptile selling websites. Our review had three goals: ensure the website was selling reptiles, check whether the website terms and conditions did not explicitly forbid automated data collection, and identify the most appropriate method of searching the content of the website (see Supplementary Data 1 for example of review datasheet).
We employed a hierarchy of five search methods depending on the structure of the website and the display of stocklists (Supplementary Code 3). The hierarchical approach minimised the analysis of irrelevant pages and minimised server load.
(1) Searched only a single html or pdf page. Where the seller has supplied a full list of stock, we could review all animals sold with only a single page using the downloadr v.0.4 package45. Occasionally animals were listed on large pages that displayed a store’s full non-animal stock as well. For stocklists supplied as a pdf, we manually downloaded them and accessed the text using the pdftools v.2.2 package46.
(2) Systematic cycling through search results. Forty-nine websites with adequate search functions allowed us to request all reptiles for sale, then examine the pages of search results one-by-one. We employed this method when website search results contained the complete details of the reptiles for sale on the search page. We ceased cycling through search pages when a URL returned a 404 error, or when 100 pages had been cycled through. Hundred pages were surveyed to prevent endless cycling back onto initial pages, or deriving errors from misinterpreting the number of search pages returned, while still exceeding the number of pages on most sites. We performed a post hoc review of ten sites searched using a cycle search method to check whether species ordering could have led systematic biases for names near the beginning of the alphabet or price. For four websites we could not determine how species were ordered, for six websites species listings were ordered by date, and for one website, species were ordered by popularity. Thus even for sites with more pages, we feel the results will not be impacted by biases given the inconsistency of approach for ordering entries on different sites. The 100-page limit may have led to missing species on large websites, but undercounting likely only affected a small portion of the websites searched via cycling methods and overlap between websites species lists mitigate suboptimal sampling on any particular website (see species-accumulation curves, Supplementary Figs. 1 and 2).
(3) Systematic cycling followed by level 1 crawl. We employed this method when sites had adequate search functionality but the details or full names of species for sale were buried one level deeper into the search results. In these instances, we ran a level 1 crawl on every search results page (Rcrawler v.0.1.9.1 package47). We followed the same stop criteria as the single page systematic cycling, 404 error or 100 search result pages.
(4) Basic level 1 crawl. Some sites had a full species list but split between clades or categories. In this case we passed the page containing links to all the different clade lists and completed a level 1 crawl40.
(5) Basic level 2 crawl. We required a level 2 crawl47 when the subsection of reptiles for sale was more specific. For example, to detect a Boa constrictor on a site that divided its stock multiple times would require moving from the ‘snake’ section, through to the ‘Boa’ section where the details of the stock are listed.
We employed these five search methods hierarchically, 1–5, and included 20 s delays between crawled requests to minimise server load on reptile selling websites. For search method 3, there was a significant chance of duplicated pages being returned; we removed duplicated pages prior to keyword searching. A few sites required multiple methods to extract complete stocklists. For search methods 3, 4 and 5, we limited the crawl further by selecting (where possible) keywords within a website’s URL that must be included for a page to be searched. For example, a website may list animals on pages that all include the pattern ‘/category=reptiles/’, therefore limiting the search of irrelevant non-stock pages.
We augmented our 2019 snapshot sampling by exploring the archived web pages stored on the Internet Archive48. For the most species-rich site (from the 2019 snapshot) we retrieved all archived web pages using the Internet Archive’s Wayback machine API49, adapting code from the wayback v.0.4.0 package50, with functions from httr v.1.4.151, jsonlite v.1.6.152, downloader v.0.445, lubridate v.1.7.453 and tibble v.2.1.3 packages54. We limited the search to pages directly pertaining to sales.
Though our online search analysis provided the number of mentions per species per page, we do not detail these numbers because sellers may list multiple individuals at once, sellers may post the same advertisement numerous times, or that advertisements can be repeated on different pages within the same website. Therefore, numbers derived from online analysis did not provide a reliable estimate of numbers per species for sale, and we elected to restrict analysis to binary species appearances.
Keyword generation
We use the complete list of 11,050 reptiles created by Reptile Database16, updated 14 August 2019, as our naming standard. We downloaded the complete list from http://www.reptile-database.org/data/, then fed the list of species into code designed to query and extract all common names, historic scientific names and locality information for each species from Reptile Database. The extraction code made use of functions from stringr v.1.4.041, XML v.3.99.0.339, xml2 v.1.2.255 and rvest v.0.3.5 packages56. We combined the resulting list with names, both common and scientific, supplied by CITES (http://checklist.cites.org/#/en [accessed 6 September 2019]) using the dplyr v.0.8.4 package57 (Supplementary Code 2). Five CITES listed species had no matching counterpart in the Reptile Database; we determined that this was caused by minor spelling mistakes. We included both spellings in our complete list of species keywords. Overall, our species keyword list comprised all scientific and common names from both Reptile Database and CITES (Supplementary Data 2), with an average of 5.82 ± 0.06 s.e. (standard error) per species and grand total of 64,342 terms (s.e. calculated using the pracma v.2.2.5 package58). Common names were predominantly English, French, German or Spanish, but occasionally included local names. We compared the number of species detected via scientific names, to the number of species detected via a combination of scientific and common names because previous work highlighted that, while correlated, they can produce different search results59.
Keyword searching and species comparison
On a site by site basis, we cleaned each page’s html code of extraneous punctuation, numbers and spacing, replacing them with single spaces. That way, two-word keywords split by line breaks, punctuation or double-spacing appeared the same as those split only by spaces.
After cleaning the html code, we searched each page for keyword matches using the stringr v.1.4.0 package41 (Supplementary Code 4, 5). Because of the large quantity of keywords and the high computational cost of collation string matching, we used a fixed string matching set to be case insensitive. Fixed string matching has the disadvantage of being sensitive to variation in how diacritics or ligatures are displayed as single or multiple characters. Our keyword searches returned the website, page number (as an index relating to the total number of pages retrieved from a given website), the keywords detected, and the corresponding Reptile Database name (Supplementary Data 3).
We searched the pages obtained from The Internet Archive48 using the same list of species (Supplementary Code 4, 7). Because there is likely a connection between the number of pages available and the number of species detected, we regressed the number of species detected in a year against the number of pages searched (n = 15, intercept = 483.72, gradient = 1.65). We excluded 2002, 2003 and 2019 for this regression because they had considerably fewer pages than all other years (mean of 3.7 ± 1.2 pages, compared to mean of 296.6 ± 48.8 pages). We plotted the residuals from the regression alongside counts of unique species per year and the number of species included in CITES appendices. To show the sensitivity to the keywords used, we counted the number of unique species in two ways: (1) counting all species detected using either scientific or common name keywords, (2) counting species only detected using scientific name keywords. The two keyword groups produce slightly different yearly species lists; therefore, changing the number of unique species per year and yearly residuals.
We compared the list of species names generated by our keyword searches to those listed in CITES. Because names of species have changed, we first converted the CITES scientific names to the most recent Reptile Database used name. However, due to species synonymisations, splits and name changes, comparisons between the list of traded species and CITES species contain some ambiguity. The ambiguity can be seen in the variation between the number of traded species covered by CITES when comparing to only the top Reptile Database name, versus when comparing the CITES list to any historically used name of traded species. For general reporting we used the more generous matching using any historic name, boosting the estimations of CITES coverage. For examination of counts of CITES covered species traded over time (Fig. 2d) we used the more stringent single name matching because of the added complexity of a changing list of CITES species and the assumption that new CITES listings would use the most recently accepted name.
Data exploration and display
We used forcats v.0.4.060 and dplyr v.0.8.457 to manipulate data, and ggplot2 v.3.2.161, scico v.1.1.062, ggpubr v.0.263 and ggforce v.0.3.164 to generate the plots. We undertook keyword searching in R v.3.5.342 and R studio v.1.2.133543. Silhouette images were obtained from http://phylopic.org/, in cases where the images were not public domain for free from attribution they were produced by Aline M. Ghilardi (CC BY-NC 3.0) and Roberto Díaz Sibaja (CC BY 3.0).
We explored the completeness of our samples—2019 snapshot and temporal—in two ways. The first was only applied to the snapshot data. We built an accumulation curve illustrating the relationship between the number of sites sampled and the number of species detected. We accomplished this by randomly resampling a subset of websites, increasing the subsample by one website until all were included. We repeated the resampling process 100 times, and plotted the results with a loess smoothed curve. The second method we applied to both snapshot and temporal data. Using the iNEXT v.2.0.19 package65,66, treating our data as raw incidences, we calculated both sample-size and coverage-based rarefaction and extrapolation metrics providing us with estimates of total species richness and sample completeness. For snapshot data we used ‘website’ (n = 151) as the resampling method, for the temporal data we used ‘year’ (2002–2019, n = 18).
We compared our data to two international trade databases (compiled species list is available in Supplementary Data 5, code for review of data sources is available in Supplementary Code 9): CITES and United States Fish and Wildlife Service’s Law Enforcement Management Information System (LEMIS). Following the online web scraping, the same types of analysis and cleaning were applied to all three databases. CITES data was retrieved from https://trade.cites.org/# on 13 May 2020) using the comparative tabulations for all ‘reptilia’ and the appropriate years (the snapshot of 2019, and 2004–2019) to download all reptile species traded over this time. We retrieved LEMIS data67,68 (v.1.1.0) via R using the lemis package69 (Supplementary Code 6). LEMIS data includes records of imports to the USA, alongside information pertaining to purpose, quantity, origin, date among other metadata, therefore quantitative data on imports for each species, or based on origins and source could be calculated in R using dplyr v.0.8.457. As for the CITES species lists, the unstandardised LEMIS names were matched to those present in Reptile Database (operating as our backbone nomenclature), leading to both synonymisations and splits. A LEMIS name was converted to a Reptile Database name if it matched any current, common, or historically used name. Names would fail to match if misspelled. By LEMIS naming, there were 639 instances of genus level listing, that were matched to 510 Reptile Database names. Of the 510 converted names, 442 appeared in other sources, suggesting genus level listings in LEMIS did not inflate species counts. Outside of generic level listings, 83 full names could not be matched. We compared the 83 names to the traded list from other sources, looking for names with fewer than 5 different characters (using the similiars v.0.1.0 package70); 56 species were found to be present in other sources by this metric. Those that failed to be converted were not included in total species counts; as final counts were entirely based on Reptile Database explicit species naming. Final species counts from all data sources are based on unique Reptile Database names and do not include any remaining generic identifiers after this synonymization/split process.
Though the research focused on the percentages of species vulnerable to trade based on various forms of IUCN and CITES categorisation, we made some efforts to quantify the proportions of items with different statuses within CITES and LEMIS. Quantifications were made using a number of different approaches. Online assessments were not directly quantified due to the possibility of listing the same individuals multiple times, or having mixed batches of specimens with variable numbers. For CITES we used the summary statistics tool in ArcMap 10.3 to quantify the means and totals for the numbers exported and imported (and listings of both are provided throughout where the numbers differ), and the range for each species or endangerment status is provided in text (or a single number if they were the same). RedList status was associated with the data by joining the scientific name field between the two databases. Sums were made for various sources, purposes and endangerment statuses for CITES data using this same approach, based on the 2004–2019 data from the CITES trade portal. ‘Terms’ (i.e., skins) were also explored, recategorising the standard terms (57 were used for reptiles) into nine (i.e., fashion, live, food, decorative, medicinal, specimen, egg, body, other uses), then summing the total item number imported and exported and determining the percentage. In addition to this we tried to quantify the trade in wild captured individuals within CITES. To try to represent individuals, terms from the CITES trade database were filtered to only include bodies, carapaces, eggs, live, shells, skeletons, skins, specimens, trophies, as most of these are mutually exclusive, though the huge quantity of reptile leather and meat could not be converted to representative individuals, skins or bodies listed as weights were also removed. Following from this the individuals from each source imported and exported could be calculated to percentages of individuals from the wild or captive bred within CITES, though these percentage values were very similar to summed total values showing the results are consistent. To investigate the extent of wild capture in LEMIS data, we restricted our summaries to items that represent individuals (whole dead bodies, live eggs, dead specimens, live individuals, full specimens, substantially whole skins, and full animal trophies), filtering out 75.6% other reptile items (79,812,310/105,536,941) leaving us 25,724,631 items to review source and purpose. The filter terms are close to those used in other recent publications which also quantified elements of trade (“live”, “bodies”, “skins”, “gall bladder”, “skulls”, “heads”, “tails”, “trophies” and “skeletons”)71, but we also excluded body-parts that may have come from the same individual (i.e., skin and skull) which may otherwise inflate numbers (79,812,310 items including skulls and skeletons; 79,796,472 excluding skulls and skeletons). The filtering to individuals made negligible difference in summaries of origin (wild or captive): 58.17% wild-sourced without the filter, 58.05% wild-sourced with the filter (61,390,757/105,536,941 without; 14,933,888/25,724,631 with); 41.32% captive sourced without the filter, 41.23% captive sourced with the filter (43,611,039/10,5536,941 without; 10,605,330/25,724,631 with). Our quantification of non-commercial trade was calculated by the number of individual animal items listed as Scientific, Reintroduction, or Biomedical research; our quantification of captive sourced trade was calculated by the number of individual animal items listed as being bred/born in captivity, commercially bred, or from ranching operations. We excluded all instances of NA in either purpose or source filters (127,881 reptile items had a missing source, purpose, or description). We additionally include clade-based analysis of source, as some taxa (i.e., crocodylia) may be more impacted for fashion trade and are imported in greater numbers. For clade-based summaries of wild capture, we summarised the quantity of traded items by genus, and further simplified the genus-summary to clade using Reptile Database genera and family information. For genera missing from Reptile Database (e.g., where genus information was family such as Varanidae), we manually assigned the clade.
Maps were created using the Global Assessment of Reptile Distributions (GARD) database72 combined with each species list as appropriate using join field and then connecting by scientific names in ArcMap 10.3 based on the corrected lists. Join by field was also used to connect species to their RedList status (downloaded from https://www.iucnredlist.org/search) and CITES appendix (from http://checklist.cites.org/#/en). To create hotspot vulnerability maps we extracted each group with different IUCN classes then used ‘count overlapping polygons’ to count the number of species with each status in any given area. This was then repeated separately for the species listed within each of the three data sources, to map the species listed as traded within each separately in addition to the total number in trade.
To obtain overall number of species and percentage of species we separated each species polygon for species in trade, and all species using QGIS, then converted them to rasters with a resolution of ~1 km using ArcCatalog. Mosaic to new raster was then used in groups of 200 species, then all mosaics added to determine overall richness for reptiles, and richness for reptiles in trade, and the percentage of reptiles in trade determined using the raster calculator ((traded_species/all_species)*100). Other trends, i.e., the percentage of species coming from different sources or with different statuses was calculated in Excel using basic approaches to quantify listings with different qualities (i.e., seized, wild, commercial and personal use) and the percentage with that status within CITES based on the number of exports and imports. For more extensive analysis of multiple factors, summary statistics were used in ArcMap after joining fields to connect species data from traded specimens of the three data sources with RedList assessments. This provided some simple statistics to further understand patterns as detailed in text, as CITES data lacks the detail of some other data sources; it was largely used to understand what species were in trade relative to existing regulations and threat.
To determine trends on a country basis we joined the CITES appendix and RedList status to the GARD layer. We used QGIS to separate a global country layer (http://thematicmapping.org/downloads/world_borders.php) into constituent countries, then clipped the GARD layer into each country with trade status noted. ISO2 codes were added to each of the country layers, then each country merged again to list each species and country and thus provide a species list for all countries within the GARD database. The number of species in and out of trade, and with and without an appendix was then calculated for each species using summary statistics, and this was repeated for each RedList status (Supplementary Data 6). For all species listed in Reptile Database but with no GARD layer, countries were listed separately and the process repeated based on the listing on the website, then the total combined with that from the GARD layer to map country richness, and the number of species with each trade status and endangerment to provide an understanding of the level of potential threat to the reptile faunas of different regions based on the trade, threat and CITES appendix of species listed in those countries.
For exploration of the time lag between species descriptions and their detection in the trade, we relied on the date of description from the author details supplied by Reptile Database. We extracted the earliest description date for each species using the stringr v.1.4.0 package41 (Supplementary Code 8), and compared this to the year reported alongside the archived pages or trade date in LEMIS. For species detected only in the snapshot data we used 2019 as their date of appearance in the trade and did not include them in the calculation of mean lag time. We only included species that had been detected directly with the scientific name of the new descriptions, subsequent name changes or common names were ignored for this analysis. We also excluded species listed as only being traded for LEMIS non-commercial purposes in this part of analysis.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

1.
Freeman, K. H. & Hayes, J. M. Fractionation of carbon isotopes by phytoplankton and estimates of ancient CO2 levels. Glob. Biogeochem. Cycles 6, 185–198 (1992).
Google Scholar 
2.
Pancost, R. D. et al. Further evidence for the development of photic-zone euxinic conditions during Mesozoic oceanic anoxic events. J. Geol. Soc. London 161, 353–364 (2004).
Google Scholar 

3.
Monteiro, F. M., Pancost, R. D., Ridgwell, A. & Donnadieu, Y. Nutrients as the dominant control on the spread of anoxia and euxinia across the Cenomanian-Turonian oceanic anoxic event (OAE2): model-data comparison. Paleoceanography 27, PA4209 (2012).
Google Scholar 

4.
Schlanger, S. O. & Jenkyns, H. C. Cretaceous oceanic anoxic events: causes and consequences. Geol. Mijnbouw 55, 179–184 (1976).
Google Scholar 

5.
Jones, M. M. et al. Astronomical pacing of relative sea level during Oceanic Anoxic Event 2: preliminary studies of the expanded SH#1 Core, Utah. Geol. Soc. Am. Bull. 131, 1702–1722 (2019).
Google Scholar 

6.
Gale, A. S. & Christenson, W. K. Occurrence of the belemnite Actinocamax plenus in the Cenomanian of SE France and its significance. Bull. Geol. Soc. Den. 43, 68–77 (1996).
Google Scholar 

7.
O’Connor, L. K. et al. A re-evaluation of the Plenus Cold Event, and the links between CO2, temperature, and seawater chemistry during OAE 2. Paleoceanogr. Paleoclimatol. 35, e2019PA003631 (2019).
Google Scholar 

8.
Kuhnt, W. et al. Unravelling the onset of Cretaceous Oceanic Anoxic Event 2 in an extended sediment archive from the Tarfaya-Laayoune Basin, Morocco. Paleoceanogr. Paleoclimatol. 32, 923–946 (2017).
Google Scholar 

9.
Kuroda, J. & Ohkouchi, N. Implications of spatiotemporal distribution of black shales deposited during the Cretaceous oceanic anoxic event-2. Paleontol. Res. 10, 345–358 (2006).
Google Scholar 

10.
Owens, J. D., Lyons, T. W. & Lowery, C. M. Quantifying the missing sink for global organic carbon burial during a Cretaceous oceanic anoxic event. Earth Planet. Sci. Lett. 499, 83–94 (2018).
Google Scholar 

11.
Berner, R. A. Phanerozoic atmospheric oxygen: new results using the GEOCARBSULF model. Am. J. Sci. 309, 603–606 (2009).
Google Scholar 

12.
Baker, S. J., Hesselbo, S. P., Lenton, T. M., Duarte, L. V. & Belcher, C. M. Charcoal evidence that rising atmospheric oxygen terminated Early Jurassic ocean anoxia. Nat. Commun. 8, 15018 (2017).
Google Scholar 

13.
Kump, L. R. Terrestrial feedback in atmosphere oxygen regulation by fire and phosphorus. Nature 335, 152–154 (1988).
Google Scholar 

14.
Watson, A., Lovelock, J. E. & Margulis, L. Methanogenesis, fires and the regulation of atmospheric oxygen. Biosystems 10, 293–298 (1978).
Google Scholar 

15.
Bond, W. J. & Scott, A. C. Fire and the spread of flowering plants in the Cretaceous. New Phytol. 188, 1137–1150 (2010).
Google Scholar 

16.
Brown, S. A. E., Scott, A. C., Glasspool, I. J. & Collinson, W. E. Cretaceous wildfires and their impact on the Earth system. Cretac. Res. 36, 162–190 (2012).
Google Scholar 

17.
Glasspool, I. J. & Scott, A. C. Phanerozoic concentrations of atmospheric oxygen reconstructed from sedimentary charcoal. Nat. Geosci. 3, 627–630 (2010).
Google Scholar 

18.
Baker, S. J. et al. CO2-induced climate forcing on the fire record during the initiation of Cretaceous oceanic anoxic event 2. Geol. Soc. Am. Bull. 132, 321–333 (2019).
Google Scholar 

19.
Zhang, M., Dai, S., Du, B., Ji, L. & Hu, S. Mid-Cretaceous hothouse climate and the expansion of early angiosperms. Acta Geol. Sin. Engl. 92, 2004–2025 (2018).
Google Scholar 

20.
Blumer, M. Polycyclic aromatic compounds in nature. Sci. Am. 234, 35–45 (1976).
Google Scholar 

21.
Lima, A. L. C., Farrington, J. W. & Reddy, C. M. Combustion-derived polycyclic aromatic hydrocarbons in the environment—a review. Environ. Forensics 6, 109–113 (2005).
Google Scholar 

22.
Youngblood, W. W. & Blumer, M. Polycyclic aromatic hydrocarbons in the environment: homologous series in soils and recent marine sediments. Geochim. Cosmochim. Acta 39, 1303–1314 (1975).
Google Scholar 

23.
Killops, S. D. & Massoud, M. S. Polycyclic aromatic hydrocarbons of pyrolytic origin in ancient sediments: evidence for Jurassic vegetation fires. Org. Geochem. 18, 1–7 (1992).
Google Scholar 

24.
Finkelstein, D. B., Pratt, L. M., Curtin, T. M. & Brassell, S. C. Wildfires and seasonal aridity recorded in Late Cretaceous strata from south-eastern Arizona, USA. Sedimentology 52, 587–599 (2005).
Google Scholar 

25.
Belcher, C. M., Finch, P., Collinson, M. E., Scott, A. C. & Grassineau, N. V. Geochemical evidence for combustion of hydrocarbons during the K-T impact event. Proc. Natl Acad. Sci. USA 106, 4112–4117 (2009).
Google Scholar 

26.
Tsikos, H. et al. Carbon-isotope stratigraphy recorded by the Cenomanian-Turonian Oceanic Anoxic Event: correlation and implications based on three key localities. J. Geol. Soc. London 161, 711–719 (2004).
Google Scholar 

27.
Jarvis, I., Lignum, J. S., Grocke, D. R., Jenkyns, H. C. & Pearce, M. A. Black shale deposition, atmospheric CO2 drawdown, and cooling during the Cenomanian-Turonian Oceanic Anoxic Event. Paleoceanography 26, PA3201 (2011).
Google Scholar 

28.
Joo, Y. J. & Sageman, B. B. Cenomanian to Campanian carbon isotope chemostratigraphy from the western interior basin, USA. J. Sediment. Res. 84, 529–542 (2014).
Google Scholar 

29.
Jenkyns, H. C., Dickson, A. J., Ruhl, M. & van den Boorn, S. H. J. M. Basalt-seawater interaction, the Plenus Cold Event, enhanced weathering and geochemical change: deconstructing Oceanic Anoxic Event 2 (Cenomanian-Turonian, Late Cretaceous). Sedimentology 64, 16–43 (2017).
Google Scholar 

30.
Heimhofer, U. et al. Vegetation response to exceptional global warmth during Oceanic Anoxic Event 2. Nat. Commun. 9, 3832 (2018).
Google Scholar 

31.
Elder, W. P. Geometry of Upper Cretaceous bentonite beds: implications about volcanic source areas and paleowind patterns, western interior, United States. Geology 16, 835–838 (1988).
Google Scholar 

32.
He, T., Pausas, J. G., Belcher, C. M., Schwilk, D. W. & Lamont, B. B. Fire-adapted traits of Pinus arose in the fiery Cretaceous. New Phytol. 194, 751–759 (2012).
Google Scholar 

33.
Belcher, C. M. & Hudspith, V. A. Changes to Cretaceous surface fire behavior influenced the spread of the early angiosperms. New Phytol. 213, 1521–1532 (2016).
Google Scholar 

34.
Chumakov, N. M. et al. Climate belts of the mid-Cretaceous time. Stratigr. Geol. Correl. 3, 241–260 (1995).
Google Scholar 

35.
Hasegawa, H. et al. Drastic shrinking of the Hadley circulation during the mid-Cretaceous Supergreenhouse. Clim. Past 8, 1323–1337 (2012).
Google Scholar 

36.
Hay, W. W. Possible solutions to several enigmas of Cretaceous climate. Int. J. Earth Sci. 108, 587–620 (2018).
Google Scholar 

37.
Hay, W. W. & Floegel, S. New thoughts about the Cretaceous climate and oceans. Earth Sci. Rev. 115, 262–272 (2012).
Google Scholar 

38.
Scopelliti, G. et al. High-resolution geochemical and biotic records of the Tethyan ‘Bonarelli Level’ (OAE2, latest Cenomanian) from the Calabianca-Guidaloca composite section, northwestern Sicily, Italy. Palaeogeogr. Palaeoclimatol. Palaeoecol. 208, 293–317 (2004).
Google Scholar 

39.
Charbonnier, G. et al. Obliquity pacing of the hydrological cycle during the Oceanic Anoxic Event 2. Earth Planet. Sci. Lett. 499, 266–277 (2018).
Google Scholar 

40.
Van Helmond, N. A. G. M. et al. A perturbed hydrological cycle during Oceanic Anoxic Event 2. Geology 42, 123–126 (2014).
Google Scholar 

41.
Carr, A. S. et al. Leaf wax n-alkane distributions in arid zone South African flora: environmental controls, chemotaxonomy and palaeoecological implications. Org. Geochem. 67, 72–84 (2014).
Google Scholar 

42.
Denis, E. H., Pedentchouk, N., Schouten, S., Pagani, M. & Freeman, K. H. Fire and ecosystem change in the Arctic across the Paleocene-Eocene Thermal Maximum. Earth Planet. Sci. Lett. 467, 149–156 (2017).
Google Scholar 

43.
Mills, B. J. E., Belcher, C. M., Lenton, T. M. & Newton, R. J. A modeling case for high atmospheric oxygen concentrations during the Mesozoic and Cenozoic. Geology 22, 1023–1026 (2016).
Google Scholar 

44.
Bergman, N. M., Lenton, T. M. & Watson, A. J. COPSE: a new model of biogeochemical cycling over Phanerozoic time. Am. J. Sci. 304, 397–437 (2004).
Google Scholar 

45.
Kump, L. Chemical stability of the atmosphere and ocean. Palaeogeogr. Palaeoclimatol. Palaeoecol. 75, 123–136 (1989).
Google Scholar 

46.
Saltzman, M. R. et al. Pulse of atmospheric oxygen during the late Cambrian. Proc. Natl Acad. Sci. USA 108, 3876–3881 (2011).
Google Scholar 

47.
Huang, J. et al. The global oxygen budget and its future projection. Sci. Bull. 63, 1180–1186 (2018).
Google Scholar 

48.
Klages, J. P. et al. Temperature rainforests near the South Pole during peak Cretaceous warmth. Nature 580, 81–86 (2020).
Google Scholar 

49.
Turgeon, S. C. & Creaser, R. A. Cretaceous oceanic anoxic event 2 triggered by a massive magmatic episode. Nature 454, 323–326 (2008).
Google Scholar 

50.
Jones, M. M., Sageman, B. B. & Selby, D. Stratigraphic record of OAE2 from the Western Interior Basin (N. America): new insights from osmium isotopes (OSi) and the expanded Big Water, UT site. In Society for Sedimentary Geology (SEPM) Research Conference on Oceanic Anoxic Events (Oral Presentation) (2016).

51.
Arinobu, T., Ishiwatari, R., Kaiho, K. & Lamolda, M. A. Spike of pyrosynthetic polycyclic aromatic hydrocarbons associated with an abrupt decrease in δ13C of a terrestrial biomarker at the Cretaceous-Tertiary boundary at Caravaca, Spain. Geology 27, 723–726 (1999).
Google Scholar 

52.
Finkelstein, D. B., Pratt, L. M. & Brassell, S. C. Can biomass burning produce a globally significant carbon-isotope excursion in the sedimentary record? Earth Planet. Sci. Lett. 250, 501–510 (2006).
Google Scholar 

53.
Barclay, R. S., McElwain, J. C. & Sageman, B. B. Carbon sequestration activated by a volcanic CO2 pulse during Ocean Anoxic Event 2. Nat. Geosci. 3, 205–208 (2010).
Google Scholar 

54.
van Bentum, E. C., Reichart, G.-J., Forster, A. & Sinninghe Damsté, J. S. Latitudinal differences in the amplitude of the OAE-2 carbon isotopic excursion: pCO2 and paleo productivity. Biogeosciences 9, 717–731 (2012).
Google Scholar 

55.
Abatzoglou, J. T. & Williams, A. P. Impact of anthropogenic climate change on wildfire across western US forests. Proc. Natl Acad. Sci. USA 113, 11770–11775 (2016).
Google Scholar 

56.
Shen, W., Sun, Y., Lin, Y., Liu, D. & Chai, P. Evidence for wildfire in the Meishan section and implications for Permian-Triassic events. Geochim. Cosmochim. Acta 75, 1992–2006 (2011).
Google Scholar 

57.
Raison, R. J. Modification of the soil environment by vegetation fires, with particular reference to nitrogen transformations: a review. Plant Soil 51, 73–108 (1979).
Google Scholar 

58.
Spencer, C. N. & Hauer, F. R. Phosphorus and nitrogen dynamics in streams during a wildfire. J. North Am. Benthol. Soc. 10, 24–30 (1991).
Google Scholar 

59.
Moody, J. A. & Martin, D. A. Initial hydrologic and geomorphic response following a wildfire in the Colorado Front Range. Earth Surf. Process. Landf. 26, 1049–1070 (2001).
Google Scholar 

60.
Guieu, C., Bonnet, S., Wagener, T. & Loye-Piot, M.-D. Biomass burning as a source of dissolved iron to the open ocean? Geophys. Res. Lett. 32, L19608 (2005).
Google Scholar 

61.
Shakesby, R. A. & Doerr, S. H. Wildfire as a hydrological and geomorphological agent. Earth Sci. Rev. 74, 269–307 (2006).
Google Scholar 

62.
Kaiho, K. et al. A forest fire and soil erosion event during the Late Devonian mass extinction. Palaeogeogr. Palaeoclimatol. Palaeoecol. 392, 272–280 (2013).
Google Scholar 

63.
Barkley, A. E. et al. African biomass burning is a substantial source of phosphorus deposition to the Amazon, tropical Atlantic Ocean, and Southern Ocean. Proc. Natl Acad. Sci. USA 116, 16216–16221 (2019).
Google Scholar 

64.
Leckie, R. M., Yuretich, R. F., West, O. L. O., Finkelstein, D. & Schmidt, M. in Stratigraphy and Paleoenvironments of the Cretaceous Western Interior Seaway, USA Vol. 6 (eds Dean, W. E. & Arthur, M. A.) 101–126 (Society for Sedimentary Geology, 1998).

65.
Pogge von Strandmann, P. A. E., Jenkyns, H. C. & Woodfine, R. G. Lithium isotope evidence for enhanced weathering during Oceanic Anoxic Event 2. Nat. Geosci. 6, 668–672 (2013).
Google Scholar 

66.
Blättler, C. L., Jenkyns, H. C., Reynard, L. M. & Henderson, G. H. Significant increases in global weathering during Oceanic Anoxic Events 1a and 2 indicated by calcium isotopes. Earth Planet. Sci. Lett. 309, 77–88 (2011).
Google Scholar 

67.
Knoll, M. A. & James, W. C. Effect of the advent and diversification of vascular land plants on mineral weathering through geologic time. Geology 15, 1099–1102 (1987).
Google Scholar 

68.
Lenton, T. M. & Watson, A. J. Redfield revisited: what regulates the oxygen content of the atmosphere? Glob. Biogeochem. Cycles 14, 149–168 (2000).
Google Scholar 

69.
Likens, G. E., Bormann, F. H. & Johnson, N. M. in Some Perspectives of the Major Biogeochemical Cycles (ed. Likens, G. E.) 93–112 (John Wiley & Sons, 1981).

70.
Boudinot, F. G. et al. Neritic ecosystem response to Oceanic Anoxic Event 2 in the Cretaceous Western Interior Seaway, USA. Palaeogeogr. Palaeoclimaol. Palaeoecol. 546, 109673 (2020).
Google Scholar 

71.
Sinninghe Damsté, J. S., van Bentum, E. C., Reichart, G.-J., Pross, J. & Schouten, S. A CO2 decrease-driven cooling and increased latitudinal temperature gradient during the mid-Cretaceous Oceanic Anoxic Event 2. Earth Planet. Sci. Lett. 293, 97–103 (2010).
Google Scholar 

72.
Van Helmond, N. A. G. M. et al. Equatorward phytoplankton migration during a cold spell within the Late Cretaceous super-greenhouse. Biogeosciences 13, 2856–2872 (2016).
Google Scholar 

73.
Forster, A., Schouten, S., Moriya, K., Wilson, P. A. & Sinninghe Damsté, J. S. Tropical warming and intermittent cooling during the Cenomanian/Turonian oceanic anoxic event 2: sea surface temperature records from the equatorial Atlantic. Paleoceanogr. 22, PA1219 (2007).
Google Scholar 

74.
Boudinot, F. G. and Sepúlveda, J. Organic geochemistry of SH#1 core: fires. PANGAEA https://doi.pangaea.de/10.1594/PANGAEA.921198 (2020). More

1.
Gienapp, P., Teplitsky, C., Alho, J. S., Mills, J. A. & Merilä, J. Climate change and evolution: disentangling environmental and genetic responses. Mol. Ecol. 17, 167–178 (2008).
CAS  Article  Google Scholar 
2.
Anderson, J. T., Inouye, D. W., McKinney, A. M., Colautti, R. I. & Mitchell-Olds, T. Phenotypic plasticity and adaptive evolution contribute to advancing flowering phenology in response to climate change. Proc. R. Soc. B Biol. Sci. 279, 3843–3852 (2012).
Article  Google Scholar 

3.
Siemann, E. & Rogers, W. E. Genetic differences in growth of an invasive tree species. Ecol. Lett. 4, 514–518 (2001).
Article  Google Scholar 

4.
Bossdorf, O., Prati, D., Auge, H. & Schmid, B. Reduced competitive ability in an invasive plant. Ecol. Lett. 7, 346–353 (2004).
Article  Google Scholar 

5.
Maron, J. L., Vilà, M., Bommarco, R., Elmendorf, S. & Beardsley, P. Rapid evolution of an invasive plant. Ecol. Monogr. 74, 261–280 (2004).
Article  Google Scholar 

6.
Byrne, K. & Nichols, R. A. Culex pipiens in London underground tunnels: differentiation between surface and subterranean populations. Heredity 82, 7–15 (1999).
Article  Google Scholar 

7.
Lee, C. E. Rapid and repeated invasions of fresh water by the copepod Eurytemora affinis. Evolution 53, 1423–1434 (1999).
Article  Google Scholar 

8.
Linnen, C. R. et al. Adaptive evolution of multiple traits through multiple mutations at a single gene. Science 339, 1312–1316 (2013).
ADS  CAS  Article  Google Scholar 

9.
Yeh, P. J. & Price, T. D. Adaptive phenotypic plasticity and the successful colonization of a novel environment. Am. Nat. 164, 531–542 (2004).
Article  Google Scholar 

10.
Price, T. D., Yeh, P. J. & Harr, B. Phenotypic plasticity and the evolution of a socially selected trait following colonization of a novel environment. Am. Nat. 172, S49–S62 (2008).
Article  Google Scholar 

11.
Lande, R. Evolution of phenotypic plasticity in colonizing species. Mol. Ecol. 24, 2038–2045 (2015).
Article  Google Scholar 

12.
Chevin, L. M. & Lande, R. Adaptation to marginal habitats by evolution of increased phenotypic plasticity. J. Evol. Biol. 24, 1462–1476 (2011).
Article  Google Scholar 

13.
Orizaola, G. & Laurila, A. Developmental plasticity increases at the northern range margin in a warm-dependent amphibian. Evol. Appl. 9, 471–478 (2016).
Article  Google Scholar 

14.
Nyamukondiwa, C., Kleynhans, E. & Terblanche, J. S. Phenotypic plasticity of thermal tolerance contributes to the invasion potential of Mediterranean fruit flies (Ceratitis capitata). Ecol. Entomol. 35, 565–575 (2010).
Article  Google Scholar 

15.
Richards, C. L., Bossdorf, O., Muth, N. Z., Gurevitch, J. & Pigliucci, M. Jack of all trades, master of some? On the role of phenotypic plasticity in plant invasions. Ecol. Lett. 9, 981–993 (2006).
Article  Google Scholar 

16.
Crispo, E. Modifying effects of phenotypic plasticity on interactions among natural selection, adaptation and gene flow. J. Evol. Biol. 21, 1460–1469 (2008).
CAS  Article  Google Scholar 

17.
Baldwin, J. M. A new factor in evolution. Am. Nat. 30(441–451), 536–553 (1896).
Article  Google Scholar 

18.
Waddington, C. H. Genetic assimilation. Adv. Genet. 10, 257–293 (1961).
CAS  Article  Google Scholar 

19.
Price, T. D., Qvarnström, A. & Irwin, D. E. The role of phenotypic plasticity in driving genetic evolution. Proc. R. Soc. B Biol. Sci. 270, 1433–1440 (2003).
Article  Google Scholar 

20.
Lande, R. Adaptation to an extraordinary environment by evolution of phenotypic plasticity and genetic assimilation. J. Evol. Biol. 22, 1435–1446 (2009).
Article  Google Scholar 

21.
Levis, N. A. & Pfennig, D. W. Evaluating ‘Plasticity-First’ evolution in nature: key criteria and empirical approaches. Trends Ecol. Evol. 31, 563–574 (2016).
Article  Google Scholar 

22.
Charmantier, G. Ontogeny of osmoregulation in crustaceans: a review. Invertebr. Reprod. Dev. 33, 177–190 (1998).
CAS  Article  Google Scholar 

23.
Cervetto, G., Gaudy, R. & Pagano, M. Influence of salinity on the distribution of Acartia tonsa (Copepoda, Calanoida). J. Exp. Mar. Bio. Ecol. 239, 33–45 (1999).
Article  Google Scholar 

24.
Ho, P.-T. et al. Impacts of salt stress on locomotor and transcriptomic responses in the intertidal gastropod Batillaria attramentaria. Biol. Bull. 236, 224–241 (2019).
Article  Google Scholar 

25.
Yang, S. et al. The salinity tolerance of the invasive golden apple snail (Pomacea canaliculata). Molluscan Res. 38, 90–98 (2018).
Article  Google Scholar 

26.
Deaton, L. E., Derby, J. G. S., Subhedar, N. & Greenberg, M. J. Osmoregulation and salinity tolerance in two species of bivalve mollusc: Limnoperna fortunei and Mytilopsis leucophaeta. J. Exp. Mar. Bio. Ecol. 133, 67–79 (1989).
Article  Google Scholar 

27.
Jordan, P. J. & Deaton, L. E. Osmotic regulation and salinity tolerance in the freshwater snail Pomacea bridgesi and the freshwater clam Lampsilis teres. Comp. Biochem. Physiol. A Mol. Integr. Physiol. 122, 199–205 (1999).
Article  Google Scholar 

28.
Bouétard, A., Côte, J., Besnard, A. L., Collinet, M. & Coutellec, M. A. Environmental versus anthropogenic effects on population adaptive divergence in the freshwater snail Lymnaea stagnalis. PLoS ONE https://doi.org/10.1371/journal.pone.0106670 (2014).
Article  PubMed  PubMed Central  Google Scholar 

29.
Sinclair, C. S. Surfing snails: population genetics of the land snail Ventridens ligera (Stylommatophora: Zonitidae) in the Potomac Gorge. Am. Malacol. Bull. 28, 105–112 (2010).
Article  Google Scholar 

30.
Hartl, D. L. & Clark, A. G. Principles of Population Genetics 4th edn. (Sinauer, Sunderland, 2007).
Google Scholar 

31.
Dmitrieva, N. I. & Burg, M. B. Elevation of extracellular NaCl increases secretion of von Willebrand Factor from endothelial cells. FASEB J. 27, 686.3 (2013).
Google Scholar 

32.
Mansour, M. M. F. Nitrogen containing compounds and adaptation of plants to salinity stress. Biol. Plant. 43, 491–500 (2000).
CAS  Article  Google Scholar 

33.
Somero, G. N. & Bowlus, R. D. Osmolytes and metabolic end products of molluscs: the design of compatible solute systems. in Mollusca, Vol. 2. Environ. Biochem. Physiol. 77–100 (1983).

34.
Lv, J. et al. Transcriptome analysis of Portunus trituberculatus in response to salinity stress provides insights into the molecular basis of osmoregulation. PLoS ONE https://doi.org/10.1371/journal.pone.0082155 (2013).
Article  PubMed  PubMed Central  Google Scholar 

35.
Wiesenthal, A. A., Müller, C., Harder, K. & Hildebrandt, J. P. Alanine, proline and urea are major organic osmolytes in the snail Theodoxus fluviatilis under hyperosmotic stress. J. Exp. Biol. https://doi.org/10.1242/jeb.193557 (2019).
Article  PubMed  Google Scholar 

36.
Yamanaka, O. & Takeuchi, R. UMATracker: an intuitive image-based tracking platform. J. Exp. Biol. https://doi.org/10.1242/jeb.182469 (2018).
Article  PubMed  Google Scholar 

37.
Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for illumina sequence data. Bioinformatics 30, 2114–2120 (2014).
CAS  Article  Google Scholar 

38.
Grabherr, M. G. et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644–652 (2011).
CAS  Article  Google Scholar 

39.
Li, B. & Dewey, C. N. Assembly of non-unique insertion content using next-generation sequencing. BMC Bioinform. https://doi.org/10.1186/1471-2105-12-S6-S3 (2011).
Article  Google Scholar 

40.
Sun, J., Nishiyama, T., Shimizu, K. & Kadota, K. TCC: An R package for comparing tag count data with robust normalization strategies. BMC Bioinform. https://doi.org/10.1186/1471-2105-14-219 (2013).
Article  Google Scholar 

41.
Tang, M., Sun, J., Shimizu, K. & Kadota, K. Evaluation of methods for differential expression analysis on multi-group RNA-seq count data. BMC Bioinform. https://doi.org/10.1186/s12859-015-0794-7 (2015).
Article  Google Scholar 

42.
Conesa, A. et al. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676 (2005).
CAS  Article  Google Scholar 

43.
Li, W. & Godzik, A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 22, 1658–1659 (2006).
CAS  Article  Google Scholar 

44.
McKenna, A. et al. The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303 (2010).
CAS  Article  Google Scholar 

45.
Foll, M. & Gaggiotti, O. A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: a Bayesian perspective. Genetics 180, 977–993 (2008).
Article  Google Scholar 

46.
Excoffier, L. & Lischer, H. E. L. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564–567 (2010).
Article  Google Scholar  More