Philippa Kaur

The custom code used to clean occurrence records and construct SDMs is available at (github.com/RhettRautsaw/ VenomMaps). We used the following R16 packages for data cleaning, manipulation, species distribution modeling, and Shiny app creation: tidyverse17 readxl18, data.table19, sf20, sp21,22, rgdal23, raster24, smoothr25, ape26, phytools27, argparse28, parallel16, memuse29, dismo30, rJava31, concaveman32, spThin33, usdm34, ENMeval35, kuenm36, shiny37, leaflet38, leaflet.extras39, leaflet.extras240, RColorBrewer41, ggpubr42, ggtext43, and patchwork44.Updating occurrence record taxonomyOur goal was to update and reconstruct the distributions of New World pitvipers. We used the Reptile Database45 (May 2021) as our primary source for current taxonomy which included the following genera: Agkistrodon, Atropoides, Bothriechis, Bothrocophias, Bothrops, Cerrophidion, Crotalus, Lachesis, Metlapilcoatlus, Mixcoatlus, Ophryacus, Porthidium, and Sistrurus. However, to ensure we captured all New World pitvipers records, we incorporated all members of the family Viperidae (all vipers and pitvipers) into our pipeline for updating occurrence record taxonomy (i.e., to account for errors in the recorded latitude, longitude, or if subfamily was not recorded).First, we collected global occurrence records for “Viperidae” from GBIF (downloaded 2021-08-19)46, Bison (downloaded 2021-08-19)47, HerpMapper (only New World taxa; downloaded 2021-08-19)48, Brazilian Snake Atlas49, BioWeb (downloaded 2021-07-07)50, unpublished data/databases from RMR, GJV, EPH, LRVA, MM, and CLP, and georeferenced literature records totaling 373,673 species-level records, 292,425 of which are New World pitvipers. Given the fluidity of taxonomy, records were often associated with outdated names. For example, Crotalus mitchelli pyrrhus was elevated to Crotalus pyrrhus51, but may still be recorded as the former in a given repository (e.g., GBIF). To correct taxonomy in our database, we checked records against a list of synonyms found on the Reptile Database and compared them to current taxonomy. If species and subspecies columns matched the same taxon (or no subspecies was recorded), then species IDs were not altered. If species and subspecies IDs did not match the same taxon, we updated taxonomy by minimizing the number of changes required to a given character string. We then manually checked all changes.Constructing distribution mapsNext, we collected preliminary distribution maps from the International Union for Conservation of Nature (IUCN; downloaded 2018-11-27)52, Global Assessment of Reptile Distributions (GARD) v1.153, Heimes54, Campbell and Lamar55, and unpublished maps. We manually curated distribution maps for all New World pitvipers in QGIS using the occurrence records, previous distribution maps, and recent publications for each taxon (note that distributions for Old World Viperidae have not yet been updated). We used a digital relief map (maps-for-free.com) and The Nature Conservancy Terrestrial Ecoregions (TNG.org)56 to identify clear distribution boundaries (e.g., mountains). We then clipped the final distributions to a land boundary (GADM v3.6)57 and smoothed the distribution using the the “chaikin” method in the R package smoothr25.Occurrence-distribution overlapOur initial taxonomy check was only concerned with records for which a subspecies was recorded and had since been elevated to species status. Therefore, many records with no assigned subspecies likely remained associated with an incorrect or outdated generic and/or specific identification. Fortunately, taxonomic changes are typically associated with changes in the species’ expected distribution. For example, when Crotalus simus was resurrected from C. durissus, the distribution of C. durissus was split: the northern portion of its range in Central America now represented the resurrected species (C. simus) and the southern portion of its range remained C. durissus55. Yet, occurrence records in Central America often remain labelled as C. durissus in data repositories. Therefore, we spatially joined records with the newly reconstructed species distribution maps to determine if they overlapped with their expected distribution (Old World taxa were joined with the GARD 1.1 distributions53).Briefly, we developed a custom function (occ_cleaner.R) to perform the spatial join and update taxonomy. First, we calculated the distance for each record to the 20 nearest distributions within 50 km (full overlap resulted in a distance of 0 m). Next, we calculated the phylogenetic distance between the recorded species ID and each species with which that record overlapped using the tree from Zaher et al.58 and adding taxa based on recent clade-specific publications (bind.tip2.R; see github.com/RhettRautsaw/VenomMaps for full list of references and details). If records overlapped with their expected species, no changes were made. If records fell outside of their expected distribution, we filtered the potential overlapping and nearby species (within 50 km) to minimize phylogenetic distance. If multiple species were equally distant (i.e., share the same common ancestor), we attempted to minimize geographic distance. If multiple species remained equally distant in both phylogenetic and geographic distance, we flagged the record to be manually checked. We also flagged records if a species’ taxonomy had changed and records were additionally flagged as potentially dubious if the taxonomic change had a phylogenetic divergence greater than 5 million years. We manually checked all flagged records and returned records to their original species ID if species identity remained uncertain. We flagged these records as potentially dubious, along with records that fell outside of their expected distribution (within 50 km), and removed all flagged records for species distribution modeling. Our final cleaned database contained 344,998 global records, of which 275,087 were New World pitvipers.Species distribution modelingWe attempted to infer SDMs for the 158 species of New World pitvipers currently recognized by the Reptile Database (May 2021) and additionally modeled the three subspecies of Crotalus ravus separately based on recommendations for species status elevation by Blair et al.59 for a total of 160 species. We developed a unix-executable R script (autokuenm.R) designed to take occurrence records, distribution maps, and environmental data and prepare these data for species distribution modeling with kuenm36. We chose to use kuenm – and MaxEnt v3.4.460 – because it has been shown to have good predictive power61 and fine-tuning of this algorithm has performance comparable to more computationally intensive ensembles62,63. Additionally, MaxEnt allows for flexibility in parameter selection64 and can function entirely with presence data14.Prior to autokuenm, to account for sampling/spatial bias during SDM, we created a bias file by using the pooled New World pitviper occurrence records as representative background data65,66,67,68. Specifically, we converted occurrence records to a raster and performed two-dimensional kernel density estimation (kde2d) with the MASS package with default settings69 and rescaled the kernel density by a factor of 1000 and rounded to three decimal places. This was then used as input to factor out sampling bias by MaxEnt. We then ran autokuenm, which is designed to subset/partition the cleaned occurrence records for a given species and prepare additional files for SDM. We first defined M-areas – or areas accessible to a given species – using the World Wildlife Fund Terrestrial Ecoregions70. Biogeographic regions represent distributional limits for many species and are reasonable hypotheses for the areas accessible to a given species71,72. To do this, we created alpha hulls from the subset of occurrence records for a given species using concaveman32 with default settings. We then identified regions with at least 20% of the region covered by the alpha hull and merged these regions together to form our final M-area. All environmental layers and the bias file were cropped to this M-area which was used as the geographic extent for modeling. We then randomly selected 5% of records to function as an independent test set for final model evaluation. Next, we generated 2000 random background points across the cropped environmental layers and used ENMeval to partition occurrence records into four sets using the checkerboard2 pattern35. Note that the background points here were not used in MaxEnt. One of the four partitions was selected at random to be used as the testing set; the remaining three partitions were used for training the MaxEnt models. If the number of occurrence records in the independent test set was less than five, then we used the training partition for final model creation and used the testing partition for final model evaluation.We tested the top-contributing variables from three sets of environmental layers: (1) bioclimatic variables, (2) EarthEnv topographic variables73, and (3) a combination of these variables. To select the top-contributing variables in each set, we wrote a custom function (SelectVariables) which used a combination of MaxEnt permutation importance and Variable Inflation Factors (VIF) to remove collinearity while keeping the variables that contributed the most to the model. Compared with variable selection via principal component analysis loadings, the permutation importance and VIF methodology demonstrated significant improvement in MaxEnt model fit. First, we designed SelectVariables to run MaxEnt using dismo::maxent with default settings and then extracted the permutation importance. We removed variables if they had 0% permutation importance. Next, we calculated VIF with usdm::vif and then iteratively removed variables by selecting the variables with two highest VIF values and removing whichever variable had the lowest permutation importance. We then recalculated VIF and repeated the process until the maximum VIF value was less than 10. Finally, we recalculated permutation importance with the remaining variables using dismo::maxent with default settings and removed variables with less than 1% permutation importance to create the final variable sets. This process was done for each species independently.With the final environmental variable, testing, and training sets, we generated SDMs using kuenm. First, we created candidate calibration models with multiple combinations of regularization multipliers (0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 2, 3, 4, 5, 6, 8, 10), feature classes (l, q, h, lq, lp, lt, lh, qp, qt, qh, pt, ph, th, lqp, lqt, lqh, lpt, lph, lth, qpt, qph, qth, pth, lqpt, lqph, lqth, lpth, qpth, lqpth), and sets of environmental predictors (bioclimatic, topographic, combination) totaling 2,958 candidate models per species. We then ran each model in parallel using GNU Parallel74. Next, we evaluated the candidate models and selected the best models using statistical significance (partial ROC), prediction ability (omission rates; OR), and model complexity (AICc) with the “kuenm_ceval” function with default settings. Specifically, models were only considered if they were statistically significant and had an OR less than 5%. If no models passed the OR criteria, the models with the minimal OR were considered. Finally, any remaining models were filtered to those within 2 AICc of the top model (Supplementary Table 1). In addition to evaluating and comparing all models together, we evaluated bioclimatic-only and combination-only models separately since these two sets of environmental variables were expected to be the best performing models given the ubiquity of bioclimatic variables in species distribution modeling (Supplementary Table 1).We generated 10 bootstrap replicates for each of the “best” calibration models using the “kuenm_mod” function. We also performed jackknifing to assess variable importance and models were output in raw format. We evaluated the final models using “kuenm_feval” with default settings. To select the best model for each comparative set (i.e., all, bioclimatic-only, and combination-only sets), we filtered the final evaluation results to minimize the OR and maximize the AUC ratio (Supplementary Table 2). If multiple models remained and were considered equally competitive, we averaged these models together (Supplementary Table 3). Because we performed three different set of comparisons, there were three “best” models per species, so we again aimed to minimize the OR and maximize the AUC ratio to select a final model for each species (Supplementary Table 4). We then converted our final models into cloglog format for visualization and threshold the models using a 10th percentile training presence cutoff (Fig. S2). Both conversion and thresholding functions are provided as R functions (raw2log, raw2clog, raster_threshold in functions.R; github.com/RhettRautsaw/VenomMaps). More

Literature reviewWe conducted a comprehensive review of relevant literature published since 1995. Studies were extracted from the China National Knowledge Infrastructure and Web of Science using the following keywords: “N (nitrogen) loss OR NO (nitric oxide) emission OR N2O (nitrous oxide) emission OR NH3 (ammonia volatilization) emission OR NO3− (nitric leaching) OR N (nitrogen) runoff AND wheat AND China”. We excluded the following types of experiment: experiments not covering the entire wheat growing season, experiments conducted in greenhouses or laboratories, experiments without zero-N control, and experiments including manure, controlled release fertilizer, or inhibitors. In total, we extracted 941 observations from 138 articles, consisting of 121 observations of NO emission, 383 of N2O emission, 185 of NH3 emission, 188 of NO3− leaching, and 64 of Nr runoff. We also extracted data on N application rates, and climate and soil variables (Fig. 1). Missing climate data were obtained from China Meteorological Data Network (https://data.cma.cn/), miss values of soil organic carbon (SOC) and total N content were obtained from the National Scientific Fertilizer Network (http://kxsf.soilbd.com/), and missing soil silt, clay, sand content, bulk density, cation exchange capacity (CEC), and pH data were obtained from the Harmonized World Soil Database (HWSD) v. 1.2 (http://www.fao.org/soils-portal/soil-survey/soilmaps-and-databases/harmonized-world-soildatabase-v12/en). Based on this dataset, the EFs of Nr loss pathways were calculated by the following equation:$$E{F}_{i}=left({E}_{treatment}{rm{-}}{E}_{control}right){rm{/}}N;applied$$
(1)
where i = 1–5, represented NO, N2O, NH3, NO3− leaching and Nr runoff, respectively. Etreatment is the loss rate of experimental treatments with applied N fertilizer, Econtrol is the loss rate of experimental control without applied N fertilizer, and N applied is the N application rate corresponding to Etreatment. The resulting data was used to develop RF models to predict EFs of the five Nr loss pathways.Fig. 1The generate framework of the Nr loss from Chinese wheat system (Nr-Wheat) 1.0 database.Full size imageRF modelsRF models outperformed empirical models in previous studies15,18,19. We employed RF models to predict the EFs of NO, N2O, NH3, NO3− leaching, and Nr runoff. Environmental factors were selected via redundancy analysis20. Redundancy analysis, a basic ordination technique for gradients analysis, produces an ordination summarizing the variation in several response variables that can be best explained by a matrix of explanatory variables based on multiple linear regression. We conducted redundancy analysis using Canoco 5 to further analyze the effects of 10 environmental factors, including 4 soil physical factors (bulk density, silt, clay, and sand content), 4 soil chemical factors (pH, SOC, CEC and total N content), and 2 weather factors (total rainfall and mean temperature during the wheat growing period) of different EFs. Ultimately, the dataset of each pathway contained an ensemble of different environmental factors (Table 1).Table 1 Environmental factors were employed to build RF model for each pathway and total explanatory rates.Full size tableWhen establishing the RF model, the first step was to select k features from a total of m (k  More

Ambitious new targets are needed to conserve nature by protecting parks and species.Credit: Tang Dehong/VCG/Getty

It will take ample time and money to slow the world’s catastrophic loss of plant and animal species — and right now, both are running dangerously low. This year, nations are due to agree to an action plan to protect global biodiversity at the 15th Conference of the Parties (COP15) to the United Nations Convention on Biological Diversity. But the meeting is already two years late because of the pandemic, and China, which will host the conference in Kunming, has yet to set a new date.Now, conflicts over financing are adding to the tension. Conservation groups and advocates suggest that rich nations must donate at least US$60 billion annually to help less-affluent ones to fund projects such as protecting areas where wildlife can thrive and tackling the illegal wildlife trade that is driving hundreds of species to extinction. This is much more than the $4 billion to $10 billion that they are estimated to be spending today, and well below the amount they are giving low- and middle-income countries (LMICs) to fight climate change, which reached around $50 billion in 2019 according to one estimate. Yet limited overseas development funds are spread ever thinner as donors deal with the pandemic and now the fallout from Russia’s invasion of Ukraine. This is where COP15 is meant to deliver: as well as agreeing to the action plan, called the Global Biodiversity Framework, nations will be encouraged to pledge more money.A mix of public and private money has started to trickle in. Currently, biodiversity funding on the table ahead of COP15 amounts to roughly $5.2 billion per year, according to estimates by a group of five leading conservation organizations. Most comes from six governments, including France, the United Kingdom and Japan, and the European Union. In April, the Global Environment Facility (GEF) — a multilateral fund to support international environmental agreements — announced that, over the next four years, around $1.9 billion will go to projects dedicated to biodiversity. However, it’s unclear how much of this will come from the coffers that donor countries have already pledged.Some cash for conservation is coming from private philanthropic donors — such as $2 billion committed by entrepreneur Jeff Bezos last year. And starting in 2020, a group of financial institutions (now 89 of them) promised to annually report their financing activities and investments that affect biodiversity, and to move away from those that do harm — a form of ecological accounting that could help to shrink the budget needed to protect biodiversity. Donors will need to reach much deeper into their pockets to meet the demands of LMICs, the custodians of much of the world’s biodiversity. In March, a group of LMICs, led by Gabon, asked for $100 billion per year in new funding when officials met in Geneva, Switzerland, to discuss progress on the Global Biodiversity Framework. The LMICs want the money placed in a new multilateral fund for biodiversity, separate from, but complementary to, the GEF.Aside from cash, the fund will need to find a new home and structure — and there are a few options. A proposal from Brazil, circulated at the Geneva meeting, suggests the fund be governed by a board of 24 members, with an equal number from rich and lower-income nations. The board would be responsible for funding decisions and would prioritize projects that help to achieve the biodiversity convention’s goals. The pitch generated interest among some countries, but also concerns that it’s an attempt by Brazil to divert attention from its failure over the past few years to protect the Amazon rainforest and prevent other environmental harm.Another option is the Kunming Biodiversity Fund, which China announced in October last year to help LMICs to safeguard their ecosystems. It allocated 1.5 billion yuan (US$223 million) to seed the fund and invited other countries to contribute, but so far none has. Sources knowledgeable about the fund say that donor countries are reluctant to pitch in because China is holding on too tightly to the reins and is not involving others in its deliberations. Details of how the fund will operate are scarce, but Nature has learnt that China is floating the idea of housing it at the Asian Infrastructure Investment Bank (AIIB), based in Beijing. Set up in 2016, the AIIB has $100 billion in total capital and 105 members, including Germany, France and the United Kingdom. The AIIB has big green plans. By 2025, it wants half of all infrastructure projects it finances to focus on climate issues. With rigorous oversight and transparency, the AIIB would make a good home for the Kunming fund.As countries prepare to meet in Nairobi on 20–26 June in a last-ditch attempt to push the biodiversity framework forwards before COP15, China, as the host, must urgently provide stronger leadership on financing, including more transparency and engagement. Progress will require quick, generous contributions from donor nations — which should prioritize grants, not loans, for biodiversity projects.Holding the COP15 meeting must be a priority, too. As China tightens restrictions in the face of a COVID-19 surge, some researchers fear that delays will stretch on, stalling conservation work and leaving less time to meet biodiversity targets. China must either commit to holding the meeting this year or let it proceed elsewhere. One option being quietly discussed is moving the meeting to Canada — home of the United Nations biodiversity convention’s secretariat — and this deserves consideration. The world needs an ambitious biodiversity plan now — nature cannot wait. More

Khamraev, Sh. R. & Bezborodov, Yu. G. Results of research on the reduction of physical evaporation of moisture from the cotton fields. Sci. World 2(33), 86–93 (2016).
Google Scholar 
Khan, A. U. et al. Production of organic fertilizers from rocket seed (Eruca sativa L.), chicken peat and Moringa oleifera leaves for growing linseed under water deficit stress. Sustainability 13(1), 1–19 (2021).CAS 

Google Scholar 
Patil Shirish, S., Kelkar Tushar, S. & Bhalerao Satish, A. Mulching: A soil and water conservation practice. Res. J. Agric For. Sci. 1(3), 26–29 (2013).
Google Scholar 
Matkovic, A. et al. Mulching as a physical weed control method applicable in medicinal plants cultivations. J. Lekovite Sirovine 35, 37–51 (2015).Article 

Google Scholar 
Nawaz, A., Lal, R., Shrestha, R. K. & Farooq, M. Mulching affects soil properties and greenhouse gas emissions under long-term no-till and plough-till systems in alfisol of Central Ohio. Land Degrad. Dev. 28(2), 673–681 (2016).Article 

Google Scholar 
Brant, V. et al. Splash erosion in maize crops under conservation management in combination with Shallow Strip-tillage before Sowing. Soil Water Res. 12(2), 106–116 (2017).CAS 
Article 

Google Scholar 
Kumar, R. et al. Effect of plant spacing and organic mulch on growth, yield and quality of natural sweetener plant Stevia and soil fertility in western Himalayas. Int. J. Plant Prod. 8(3), 311–334 (2014).ADS 

Google Scholar 
Seleiman, M. F. & Kheir, A. M. S. Maize productivity, heavy metals uptake and their availability in contaminated clay and sandy alkaline soils as affected by inorganic and organic amendments. Chemosphere 204, 514–522 (2018).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Seleiman, M. F. & Kheir, A. M. S. Saline soil properties, quality and productivity of wheat grown with bagasse ash and thiourea in different climatic zones. Chemosphere 193, 538–546 (2018).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Chakraborty, D. et al. Effect of mulching on soil and plant water status, and the growth and yield of wheat (Triticum aestivum L.) in a semi-arid environment. Agric. Water Manag. 95(12), 1323–1334 (2008).Article 

Google Scholar 
Ahmad, Z. I., Ansar, M., Iqbal, M. & Minhas, N. M. Effect of planting geometry and mulching on moisture conservation, weed control and wheat growth under rainfed conditions. Pak. J. Bot. 39(4), 1189–1195 (2007).
Google Scholar 
Teame, G. Effect of organic mulches and land preparation methods on soil moisture and sesame productivity. Afr. J. Agric. Res. 12(38), 2836–2843 (2017).Article 

Google Scholar 
Lehar, L., Wardiyati, T., Moch Dawam, M. & Suryanto, A. Influence of mulch and plant spacing on yield of Solanum tuberosum L. cv. Nadiya at medium altitude. Int. Food Res. J. 24(3), 1338–1344 (2017).CAS 

Google Scholar 
Arash, K. The evaluation of water use efficiency in common bean (Phaseolus vulgaris L.) in irrigation condition and mulch. Sci. Agric. 2(3), 60–64 (2013).
Google Scholar 
Artyszak, A., Gozdowski, D. & Kucińska, K. The yield and technological quality of sugar beet roots cultivated in mulches. Plant Soil Environ. 60(10), 464–469 (2014).Article 

Google Scholar 
Brittaine, R. & Lutaladio, N. Jatropha: A Smallholder Bioenergy Crop. The Potential for Pro-poor Development Integrated Crop Management, Vol. 8 (IFAD/FAO, 2010). http://www.fao.orgElbehri, A., Segerstedt, A. & Liu, P. Biofuels and the sustainability challenge: A global assessment of sustainability issues, trends and policies for biofuels and related feedstocks. Food and Agric. Organ. United Nations (FAO) xvi-174 (2013).King, A. J. et al. Potential of Jatropha curcas as a source of renewable oil and animal feed. J. Exp. Bot. 60(10), 2897–2905 (2009).CAS 
PubMed 
Article 

Google Scholar 
Raheman, H. 14 Jatropha. Handbook of Bioenergy Crop Plants, 315–345 (2012).Ullah, F., Bano, A. & Nosheen, A. Sustainable measures for biodiesel. Effects 36(23), 2621–2628 (2014).CAS 

Google Scholar 
Irshad, M. et al. Evaluation of Jatropha curcas L. leaves mulching on wheat growth and biochemical attributes under water stress. BMC Plant Biol. 21(1), 1–12 (2021).Article 
CAS 

Google Scholar 
Dieye, T. et al. The effect of Jatropha curcas L. leaf litter decomposition on soil carbon and nitrogen status and bacterial community structure (Senegal). J. Soil Sci. Environ Manag. 7(3), 32–44 (2016).CAS 
Article 

Google Scholar 
Kafi, M. & Salehi, M. Kochia scoparia as a model plant to explore the impact of water deficit on halophytic communities. Pak. J. Bot. 44, 257–262 (2012).
Google Scholar 
Yang, Y. M., Liu, X. J., Li, W. Q. & Li, C. Z. Effect of different mulch materials on winter wheat production in desalinized soil in Heilonggang region of North China. J. Zhejiang Univ. Sci. B 7(11), 858–867 (2006).PubMed 
PubMed Central 
Article 

Google Scholar 
Xie, Z. K., Wang, Y. J. & Li, F. M. Effect of plastic mulching on soil water use and spring wheat yield in arid region of northwest China. Agric. Water Manag. 75(1), 71–83 (2005).Article 

Google Scholar 
Khan, R. H., Anwar-ul-Haq, K. & Sajjad, M. R. Effect of different types of mulches on grain yield and yield components of wheat (Triticum aestivum) under rainfed condition. J. Biol. Agric. Healthc. 4(12), 85–91 (2014).
Google Scholar 
Weidhuner, A., Afshar, R. K., Luo, Y., Battaglia, M. & Sadeghpour, A. Particle size affects nitrogen and carbon estimate of a wheat cover crop. Agron. J. 111(6), 3398–3402 (2019).CAS 
Article 

Google Scholar 
Ding, Z. et al. The integrated effect of salinity, organic amendments, phosphorus fertilizers, and deficit irrigation on soil properties, phosphorus fractionation and wheat productivity. Sci. Rep. 10(1), 1–13 (2020).Article 
CAS 

Google Scholar 
Rummana, S., Amin, A. K. M. R., Islam, M. S. & Faruk, G. M. Effect of irrigation and mulch materials on growth and yield of wheat. Bang. Agron. J. 21(1), 71–76 (2018).Article 

Google Scholar 
Richard, L. A. Diagnosis and improvement of saline and alkaline soils. Handbook No. 60 (US Depart. Agric., 1954).McLean, E. O. Soil pH and lime requirement. Methods of Soil Analysis: Part 2 Chemical and Microbiological Properties, Vol. 9, 199–224 (1983).Walkley, A. A critical examination of a rapid method for determining organic carbon in soils—Effect of variations in digestion conditions and of inorganic soil constituents. Soil Sci. 63, 251–264 (1947).ADS 
CAS 
Article 

Google Scholar 
Singleton, V. L., Orthofer, R. & Lamuela-Raventos, R. M. Analysis of total phenols and other oxidation substrates and antioxidants by means of Folin–Ciocalteu reagent. Methods Enzymol. 299, 152–178 (1999).CAS 
Article 

Google Scholar 
Vance, E. D., Brookes, P. C. & Jenkinson, D. S. An extraction method for measuring soil microbial biomass C. Soil Biol. Biochem. 19, 703–707 (1987).CAS 
Article 

Google Scholar 
Bremner, J. M. & Mulvaney, C. S. Nitrogen-total. In Methods of Soil Analysis. Part 2. Chemical and Microbiological Properties (eds Page, A. L. et al.) 595–624 (Soil Sci. Society America, 1982).
Google Scholar 
Steel, R. G. D., Torrie, J. H. & Dickey, D. A. Principles and Procedures of Statistics: A Biometrical Approach 3rd edn, 246 (McGraw-Hill, 1997).
Google Scholar 
Brady, N. C. & Weil, R. R. Soil colloids: Seat of soil chemical and physical acidity. Nat. Prop. Soils 5(13), 311–358 (2008).
Google Scholar 
Scharenbroch, B. C. & Lloyd, J. E. Particulate organic matter and soil nitrogen availability in urban landscapes. Arboricul. Urb. For. 32(4), 180–191 (2006).Article 

Google Scholar 
Bhadha, J. H., Capasso, J. M., Khatiwada, R., Swanson, S. & LaBorde, C. Raising soil organic matter content to improve water holding capacity. UF/IFAS 1–5 (2017).Chalker-Scott, L. Impact of mulches on landscape plants and the environment—A review. J. Environ. Hortic. 25(4), 239–249 (2007).Article 

Google Scholar 
Liu, Z., Fu, B., Zheng, X. & Liu, G. Plant biomass, soil water content and soil N:P ratio regulating soil microbial functional diversity in a temperate steppe: A regional scale study. Soil Biol. Biochem. 42(3), 445–450 (2010).CAS 
Article 

Google Scholar 
Bai, S. H., Blumfield, T. J. & Reverchon, F. The impact of mulch type on soil organic carbon and nitrogen pools in a sloping site. Biol. Fertil. Soils 50(1), 37–44 (2014).Article 

Google Scholar 
Yang, H. et al. The combined effects of maize straw mulch and no-tillage on grain yield and water and nitrogen use efficiency of dry-land winter wheat (Triticum aestivum L.). Soil Tillage Res. 197, 104485 (2020).Article 

Google Scholar 
Li, X. J. et al. Abscisic acid pretreatment enhances salt tolerance of rice seedlings: Proteomic evidence. Biochim. Biophys. Acta (BBA) Proteins Proteomics 1804(4), 929–940 (2010).CAS 
Article 

Google Scholar 
Fang, S., Xie, B., Liu, D. & Liu, J. Effects of mulching materials on nitrogen mineralization, nitrogen availability and poplar growth on degraded agricultural soil. New For. 41(2), 147–162 (2011).Article 

Google Scholar 
Houghton, J. T. Climate Change 2001: The Scientific Basis 419–470 (2001).Johnson, D. et al. Plant community composition affects the biomass, activity and diversity of microorganisms in limestone grassland soil. Eur. J. Soil Sci. 54(4), 671–678 (2003).Article 

Google Scholar 
Johnson, M. J., Lee, K. Y. & Scow, K. M. DNA finger printing reveals links among agricultural crops, soil properties, and the composition of soil microbial communities. Geoderma 114, 279–303 (2003).ADS 
Article 

Google Scholar 
Nielsen, N. M., Winding, A. & Binnerup, S. Microorganisms as Indicators of Soil Health 15–16 (Ministry of the Environment, National Environ. Res. Inst., 2002).
Google Scholar 
Wilkinson, S. C. et al. PLFA profiles of microbial communities in decomposing conifer litters subject to moisture stress. Soil Biol. Biochem. 34(2), 189–200 (2002).CAS 
Article 

Google Scholar 
Drenovsky, R. E., Vo, D., Graham, K. J. & Scow, K. M. Soil water content and organic carbon availability are major determinants of soil microbial community composition. Microb. Ecol. 48(3), 424–430 (2004).CAS 
PubMed 
Article 

Google Scholar 
Liu, Y. Y., Yao, H. Y. & Huang, C. Y. Influence of soil moisture regime on microbial community diversity and activity in a paddy soil. Acta Pedol. Sin. 43, 828–834 (2006).
Google Scholar 
Jensen, K. D., Beier, C., Michelsen, A. & Emmett, B. A. Effects of experimental drought on microbial processes in two temperate heathlands at contrasting water conditions. Appl. Soil Ecol. 24(2), 165–176 (2003).Article 

Google Scholar 
Stoklosa, A., Hura, T., Stupnicka-Rodzynkiewicz, E., Dabkowska, T. & Lepiarczyk, A. The influence of plant mulches on the content of phenolic compounds in soil and primary weed infestation of maize. Acta. Agron. Bot. 61(2), 205–219 (2008).
Google Scholar 
Ohno, T. Oxidation of phenolic acid derivatives by soil and its relevance to allelopathic activity. J. Environ. Qual. 30(5), 1631–1635 (2001).CAS 
PubMed 
Article 

Google Scholar 
Farooq, S., Shahid, M., Khan, M. B., Hussain, M. & Farooq, M. Improving the productivity of bread wheat by good management practices under terminal drought. J. Agric. Crop Sci. 201(3), 173–188 (2015).Article 

Google Scholar 
Madani, A., Rad, A. S., Pazoki, A., Nourmohammadi, G. & Zarghami, R. Wheat (Triticum aestivum L.) grain filling and dry matter partitioning responses to source: Sink modifications under postanthesis water and nitrogen deficiency. Acta Sci. Agron. 32, 145–151 (2010).CAS 
Article 

Google Scholar 
Deng, X. P., Shan, L., Zhang, H. & Turner, N. C. Improving agricultural water use efficiency in arid and semiarid areas of China. Agric. Water Manag. 80(1–3), 23–40 (2006).Article 

Google Scholar 
Athar, H. R., Khan, A. & Ashraf, M. Inducing salt tolerance in wheat by exogenously applied ascorbic acid through different modes. J. Plant Nutr. 32, 1799–1817 (2009).CAS 
Article 

Google Scholar 
Luo, et al. Dual plastic film and straw mulching boosts wheat productivity and soil quality under the El Nino in semiarid Kenya. Sci. Total Environ. 738, 139808 (2020).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Duan, et al. Improvement of wheat productivity and soil quality by arbuscular mycorrhizal fungi is density-and moisture-dependent. Agron. Sustain. Dev. 41(1), 1–12 (2021).Article 
CAS 

Google Scholar  More

Mesocosm experimentsSi:Nexport measurementsBetween 2010 and 2014, we conducted five in situ mesocosm experiments to assess impacts of OA on natural plankton communities. Study sites covered a large latitudinal gradient (28 °N–79 °N) and diverse oceanic environments/ecosystems (Extended Data Fig. 1 and Extended Data Table 1). Sample collection and processing was conducted every 1 or 2 days throughout the experiments. Sinking particulate matter was obtained from sediment traps attached to the bottom of each mesocosm, thereby collecting the entire material sinking down in the enclosed water column36. Processing of sediment trap samples followed a previous protocol37. Samples for particulate matter suspended in the water column were collected with depth-integrating water samplers (HYDRO-BIOS) and filtered following standard procedures. Biogenic silica was leached from the sediment trap samples and filters by alkaline pulping (0.1 M NaOH at 85 °C). After 135 min the leaching process was terminated with 0.05 M H2SO4 and dissolved silica was measured spectrophotometrically38. Carbon and nitrogen content were determined using an elemental CN analyser (EuroEA)39.Analysis of OA impactsTo test for a systemic influence of OA on Si:Nexport, we synthesized the datasets from the different experiments and (i) conducted a meta-analysis to quantify effect sizes, and (ii) computed probability density estimates. Because the experimental design, the range of CO2 treatments, and the time periods for our analysis of Si:Nexport varied to some extent among experiments (Extended Data Table 1), we pooled mesocosms for ambient conditions and in the ({{p}}_{{{rm{CO}}}_{2}}) range of ~700–1,000 μatm (‘OA treatment’), corresponding to end-of-century values according to RCP 6.0 and 8.5 emission scenarios15. Effect sizes were calculated as log-transformed response ratios lnRR, an approach commonly used in meta-analysis40:$${rm{l}}{rm{n}}{rm{R}}{rm{R}}={rm{l}}{rm{n}}{X}_{{rm{O}}{rm{A}}}-{rm{l}}{rm{n}}{X}_{{rm{c}}{rm{o}}{rm{n}}{rm{t}}{rm{r}}{rm{o}}{rm{l}}},$$where X is the arithmetic mean of Si:Nexport ratios under OA and ambient conditions (Extended Data Table 1). Effect sizes 0 indicate that the effect was positive. Effects are considered statistically significant when 95% confidence intervals (calculated from pooled standard deviations) do not overlap with zero. The overall effect size across all studies was computed by weighing individual effect sizes according to their variance, following the common methodology for meta-analyses40. In addition, we computed probability densities of Si:Nexport based on kernel density estimation, which better accounts for data with skewed or multimodal distributions41. Another advantage of this approach is that it does not require the calculation of temporal means. Instead, the entire data timeseries can be incorporated into the analysis, thus retaining information about temporal variability. Confidence intervals of the density estimates were calculated with a bootstrapping approach using data resampling (1,000 permutations)41. The resulting probability density plots can be interpreted analogously to histograms. Differences among ambient and OA conditions are considered statistically significant when confidence intervals of the probability density distributions do not overlap. Numbers for suspended and sinking Si, C and N (and their respective ratios) for the analysis period are given in Extended Data Table 2.Analysis of pH effects on Si:N in global sediment trap dataWe analysed a recent compilation of global sediment trap data (674 locations collected between 1976 and 2012)35. The aim of this analysis was to assess the influence of pH on opal dissolution in the world ocean. In contrast to the mesocosm experiments, where export fluxes were measured only at one depth, the global dataset provides depth-resolved information, enabling us to examine the vertical change in the Si:N ratio of sinking particulate matter and how this correlates with pH. It has long been known that the silica content of sinking particles increases with depth, as opal dissolution is less efficient than organic matter remineralization25,42. The resulting accumulation of Si relative to N can be quantified as the change in Si:N with increasing depth, that is, the slope of the relationship of depth versus Si:N (ΔSi:N, in units of m−1). Our approach is analogous to previous studies, which used vertical profiles of Si:C as a proxy for differential dissolution/remineralization of opal and organic matter, and its regional variability in the ocean24,42. We extracted all data that (I) included simultaneous measurements of Si and N, and (II) contained vertical profiles with at least three depth levels (so that ΔSi:N [m−1] can be calculated). We then calculated linear regressions for individual Si:N profiles and subsequently extracted those for which Si:N displayed a statistically significant relationship with depth (p  More

Storey, M. et al. Timing of hot spot—Related volcanism and the breakup of Madagascar and India. Science (80-) 267, 852–855 (1995).CAS 
Article 
ADS 

Google Scholar 
Wilmé, L., Goodman, S. M. & Ganzhorn, J. U. Biogeographic evolution of Madagascar’s microendemic biota. Science (80-) 312, 1063–1065 (2006).Article 
ADS 
CAS 

Google Scholar 
Myers, N., Mittermeler, R. A., Mittermeler, C. G., Da Fonseca, G. A. B. & Kent, J. Biodiversity hotspots for conservation priorities. Nature 403, 853–858 (2000).CAS 
PubMed 
Article 
ADS 

Google Scholar 
Vences, M., Wollenberg, K. C., Vieites, D. R. & Lees, D. C. Madagascar as a model region of species diversification. Trends Ecol. Evol. 24, 456–465 (2009).PubMed 
Article 

Google Scholar 
Rakotomanana, H., Jenkins, R. K. B. & Ratsimbazafy, J. Conservation challenges for Madagascar in the next decade. In Conservation Biology: Voices from the Tropics (eds Raven, P. H., Sodhi, N. S. & Gibson, L.) 33–39 (Wiley-Blackwell, 2013). https://doi.org/10.1002/9781118679838.ch5.Jenkins, R. K. B. et al. Extinction risks and the conservation of Madagascar’s reptiles. PLoS ONE 9, 1 – 14 (2014). https://doi.org/10.1371/journal.pone.0100173Velosoa, J. et al. An integrated research, management, and community conservation program for the Rere (Madagascar Big-headed turtle), Erymnochelys madagascariensis. In Chelonian Research Monographs, Contributions in Turtle and Tortoise Research (eds Rhodin, A. G. J.) 171–177 (Chelonian Research Foundation, 2014). https://doi.org/10.3854/crm.6.a27p171.Leuteritz, T., Kuchling, G., Garcia, G. & Velosoa, J. Erymnochelys madagascariensis. In Chelonian Research Monographs, Contributions in Turtle and Tortoise Research (eds Rhodin, A. G. J.) 56–58 (Chelonian Research Foundation, 2014). https://doi.org/10.3854/crm.6.a11p56.Rafeliarisoa, T., Shore, G., Engberg, S., Louis, E. & Brenneman, R. Characterization of 11 microsatellite marker loci in the Malagasy big-headed turtle (Erymnochelys madagascariensis). Mol. Ecol. Notes 6, 1228–1230 (2006).CAS 
Article 

Google Scholar 
Roca, V., García, G. & Montesinos, A. Gastrointestinal helminths found in the three freshwater turtles (Erymnochelys madagascariensis, Pelomedusa subrufa and Pelusios castanoides) from Ankarafantsika National Park, Madagascar. Helminthologia 44, 177–182 (2007).Article 

Google Scholar 
Kuchling, G. & Garcia, G. Pelomedusidae, freshwater turtles. In The Natural History of Madagascar (eds Goodman, S. M. & Benstead, J. P.) 956–960 (University of Chicago Press, 2003).
Google Scholar 
Pedrono, M. & Smith, L. Overview of the natural history of Madagascar’s endemic tortoises and freshwater turtles: Essential components for effective conservation. In Chelonian Research Monographs, Contributions in Turtle and Tortoise Research (eds Rhodin, A. G. J.) 59–66 (Chelonian Research Foundation, 2014). https://doi.org/10.3854/crm.6.a12p59.Kuchling, G. Population structure, reproductive potential and increasing exploitation of the freshwater turtle Erymnochelys madagascariensis. Biol. Conserv. 43, 107–113 (1988).Article 

Google Scholar 
Allnutt, T. F. et al. A method for quantifying biodiversity loss and its application to a 50-year record of deforestation across Madagascar. Conserv. Lett. 1, 173–181 (2008).Article 

Google Scholar 
Leuteritz, T., Kuchling, G., Garcia, G. & Velosoa, J. Erymnochelys madagascariensis (errata version published in 2016). The IUCN Red List of Threatened Species. 2008, 1–3 (2008).Kuchling, G. Concept and design of the Madagascar side-necked turtle Erymnochelys madagascariensis breeding facility at Ampijoroa, Madagascar. Dodo 36, 62–74 (2000).
Google Scholar 
Witzenberger, K. A. & Hochkirch, A. Ex situ conservation genetics: A review of molecular studies on the genetic consequences of captive breeding programmes for endangered animal species. Biodivers. Conserv. 20, 1843–1861 (2011).Article 

Google Scholar 
Stanton, D. W. G. et al. Genetic structure of captive and free-ranging okapi (Okapia johnstoni) with implications for management. Conserv. Genet. 16, 1115–1126 (2015).Article 

Google Scholar 
Boumans, L., Vieites, D. R., Glaw, F. & Vences, M. Geographical patterns of deep mitochondrial differentiation in widespread Malagasy reptiles. Mol. Phylogenet. Evol. 45, 822–839 (2007).CAS 
PubMed 
Article 

Google Scholar 
Orozco-Terwengel, P., Andreone, F., Louis, E. & Vences, M. Mitochondrial introgressive hybridization following a demographic expansion in the tomato frogs of Madagascar, genus Dyscophus. Mol. Ecol. 22, 6074–6090 (2013).CAS 
PubMed 
Article 

Google Scholar 
Pearson, R. G. & Raxworthy, C. J. The evolution of local endemism in Madagascar: Watershed versus climatic gradient hypotheses evaluated by null biogeographic models. Evolution (New York) 63, 959–967 (2009).
Google Scholar 
Sunde, J., Yıldırım, Y., Tibblin, P. & Forsman, A. Comparing the performance of microsatellites and RADseq in population genetic studies: Analysis of data for pike (Esox lucius) and a synthesis of previous studies. Front. Genet. 11, 218 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Hulce, D., Li, X., Snyder-Leiby, T. & Liu, J. GeneMarker® genotyping software: Tools to increase the statistical power of DNA fragment analysis. J. Biomol. Tech. https://doi.org/10.1002/wps.20394 (2011).Article 
PubMed Central 

Google Scholar 
Van Oosterhout, C., Hutchinson, W. F., Wills, D. P. M. & Shipley, P. MICRO-CHECKER: Software for identifying and correcting genotyping errors in microsatellite data. Mol. Ecol. Notes 4, 535–538 (2004).Article 
CAS 

Google Scholar 
Carlsson, J. Effects of microsatellite null alleles on assignment testing. J. Hered. 99, 616–623 (2008).CAS 
PubMed 
Article 

Google Scholar 
Bossuyt, F. & Milinkovitch, M. C. Convergent adaptive radiations in Madagascan and Asian ranid frogs reveal covariation between larval and adult traits. Proc. Natl. Acad. Sci. U. S. A. 97, 6585–6590 (2000).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
Rousset, F. GENEPOP’007: A complete re-implementation of the GENEPOP software for Windows and Linux. Mol. Ecol. Resour. 8, 103–106 (2008).PubMed 
Article 

Google Scholar 
Beaumont, M. A. Detecting population expansion and decline using microsatellites. Genetics 153, 2013–2029 (1999).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Bulut, Z. et al. Microsatellite mutation rates in the eastern tiger salamander (Ambystoma tigrinum tigrinum) differ 10-fold across loci. Genetica 136, 501–504 (2009).CAS 
PubMed 
Article 

Google Scholar 
Brooks, S. P. & Gelman, A. General methods for monitoring convergence of iterative simulations. J. Comput. Graph. Stat. 7, 434–455 (1998).MathSciNet 

Google Scholar 
Plummer, M. & Murrell, P. CODA: Convergence Diagnosis and Output Analysis for MCMC. R News. 6, 7–11 (2006).R Development Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2008). https://doi.org/10.1017/CBO9781107415324.004.Book 

Google Scholar 
Pritchard, J. K., Stephens, M. & Donnelly, P. Inference of Population Structure Using Multilocus Genotype Data. Genetics 155, 945–959 (2000). https://doi.org/10.1111/j.1471-8286.2007.01758.x.Puechmaille, S. J. The program structure does not reliably recover the correct population structure when sampling is uneven: Subsampling and new estimators alleviate the problem. Mol. Ecol. Resour. 16, 608–627 (2016).PubMed 
Article 

Google Scholar 
Hale, M. L., Burg, T. M. & Steeves, T. E. Sampling for microsatellite-based population genetic studies: 25 to 30 individuals per population is enough to accurately estimate allele frequencies. PLoS ONE 7, e45170 (2012).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
Francis, R. M. Pophelper: An R package and web app to analyse and visualize population structure. Mol. Ecol. Res. 17, 27–32 (2017).Evanno, G., Regnaut, S. & Goudet, J. Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. Mol. Ecol. 14, 2611–2620 (2005).CAS 
PubMed 
Article 

Google Scholar 
Kearse, M. et al. Geneious basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28, 1647–1649 (2012).PubMed 
PubMed Central 
Article 

Google Scholar 
Dieringer, D. & Schlötterer, C. Microsatellite analyser (MSA): A platform independent analysis tool for large microsatellite data sets. Mol. Ecol. Notes 3, 167–169 (2003).CAS 
Article 

Google Scholar 
Narum, S. R. Beyond Bonferroni: Less conservative analyses for conservation genetics. Conserv. Genet. 7, 783–787 (2006).CAS 
Article 

Google Scholar 
Goudet, J. FSTAT (version 1.2): A computer program to calculate F-statistics. J. Hered. 86, 485–486 (1995).Article 

Google Scholar 
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).PubMed 
Article 

Google Scholar 
Librado, P. & Rozas, J. DnaSP v5: A software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452 (2009).CAS 
PubMed 
Article 

Google Scholar 
Prost, S. & Anderson, C. N. K. TempNet: A method to display statistical parsimony networks for heterochronous DNA sequence data. Methods Ecol. Evol. 2, 663–667 (2011).Article 

Google Scholar 
Paquette, S. R. et al. Riverbeds demarcate distinct conservation units of the radiated tortoise (Geochelone radiata) in southern Madagascar. Conserv. Genet. 8, 797–807 (2007).CAS 
Article 

Google Scholar 
Bouchard, C., Tessier, N. & Lapointe, F. J. Watersheds influence the wood turtle’s (Glyptemys insculpta) genetic structure. Conserv. Genet. 20, 653–664 (2019).Article 

Google Scholar 
Perlman, S. J., Hodson, C. N., Hamilton, P. T., Opit, G. P. & Gowen, B. E. Maternal transmission, sex ratio distortion, and mitochondria. Proc. Natl. Acad. Sci. U. S. A. 112, 10162–10168 (2015).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
Pearse, D. E. et al. Estimating population structure under nonequilibrium conditions in a conservation context: Continent-wide population genetics of the giant Amazon river turtle, Podocnemis expansa (Chelonia; Podocnemididae). Mol. Ecol. 15, 985–1006 (2006).CAS 
PubMed 
Article 

Google Scholar 
Pearse, D. E. & Avise, J. C. Turtle mating systems: Behavior, sperm storage, and genetic paternity. J. Hered. 92, 206–211 (2001).CAS 
PubMed 
Article 

Google Scholar 
Claussen, M. et al. Simulation of an abrupt change in Saharan vegetation in the mid-Holocene. Geophys. Res. Lett. 26, 2037–2040 (1999).Article 
ADS 

Google Scholar 
Virah-Sawmy, M., Willis, K. J. & Gillson, L. Threshold response of Madagascar’s littoral forest to sea-level rise. Glob. Ecol. Biogeogr. 18, 98–110 (2009).Article 

Google Scholar 
Wahlund, S. Zusammensetzung von populationen und korrelationserscheinungen vom standpunkt der vererbungslehre aus betrachtet. Hereditas 11, 65–106 (1928).Article 

Google Scholar 
Hurst, G. D. D. & Jiggins, F. M. Problems with mitochondrial DNA as a marker in population, phylogeographic and phylogenetic studies: The effects of inherited symbionts. Proc. R. Soc. B Biol. Sci. 272, 1525–1534 (2005).CAS 
Article 

Google Scholar 
Hill, W. G. & Robertson, A. The effect of linkage on limits to artificial selection. Genet. Res. (Camb.) 89, 311–336 (2008).Article 

Google Scholar 
Valenzuela, N. Multiple paternity in side-neck turtles Podocnemis expansa: Evidence from microsatellite DNA data. Mol. Ecol. 9, 99–105 (2000).CAS 
PubMed 
Article 

Google Scholar 
Moritz, C. Defining ‘Evolutionarily Significant Units’ for conservation. Trends Ecol. Evol. 9, 373–375 (1994).CAS 
PubMed 
Article 

Google Scholar 
Volkmann, L., Martyn, I., Moulton, V., Spillner, A. & Mooers, A. O. Prioritizing populations for conservation using phylogenetic networks. PLoS ONE 9, 1–10 (2014). https://doi.org/10.1371/journal.pone.0088945Article 
CAS 

Google Scholar 
García-Dorado, A. & Caballero, A. Neutral genetic diversity as a useful tool for conservation biology. Conserv. Genet. 22, 541–545 (2021).Article 

Google Scholar 
Frankham, R. Genetic rescue of small inbred populations: Meta-analysis reveals large and consistent benefits of gene flow. Mol. Ecol. 24, 2610–2618 (2015).PubMed 
Article 

Google Scholar 
Teixeira, J. C. & Huber, C. D. The inflated significance of neutral genetic diversity in conservation genetics. Proc. Natl. Acad. Sci. U. S. A. 118, 1–10 (2021). https://doi.org/10.1073/pnas.2015096118Araki, H., Cooper, B. & Blouin, M. S. Genetic effects of captive breeding cause a rapid, cumulative fitness decline in the wild. Science (80-) 318, 100–103 (2007).CAS 
Article 
ADS 

Google Scholar  More

In this study, we aimed to determine if tropical cyclones in the northwestern GOM differentially affected mesozooplankton community structure. We found that multivariate community structure varied between storm and non-storm years. However, among the three hurricanes, only the post-Harvey mesozooplankton communities were distinct from years where no storms occurred. Multivariate dispersion, a measure of variance within community structure, did not differ between storm and non-storm years. This result refutes our hypothesis that variability in zooplankton community structure would be higher during storm years opposed to non-storm years. Our prediction that post-storm mesozooplankton communities would differ from non-storm communities was supported, as was our expectation that mesozooplankton community structure varied among storms. We hypothesized that due to the major flooding and rainfall of Harvey, reduced salinity would likely be the main driver of mesozooplankton community differences relative to non-storm years. This expectation was partially met; differences in mesozooplankton community structure between Harvey and years with no storms were driven by salinity, stratification, and ‘distance-from-shore’ rather than solely salinity. This indicates that NWGOM zooplankton community structure varies holistically with biophysical conditions rather than being primarily driven by one or two dominating factors24. Moreover, we found that the presence of Hurricane Harvey, rather than temperature, explained the second greatest amount of variance in zooplankton community structure reflecting the importance of considering how complex disturbance mechanisms might compound and result in a unique ecological responses following a tropical cyclone. The overall PSEM showed that high salinity was directly associated with reduced fluorescence and depressed zooplankton abundance. Higher abundances were found to result in higher community dominance (i.e., lower evenness). Conversely, high salinities indirectly reduced zooplankton evenness via greater water column stratification. Lower fluorescence at higher salinities supports the spatial patterns in NWGOM phytoplankton identified by Kurtay et al.20 following Hurricane Harvey. Those authors observed declines in overall phytoplankton abundance and communities increasingly dominated by cells  More

Data preparationHistorical dataThe primary source of data is historical parish registers, which have been transcribed under the supervision of many of the study authors over a number of decades, primarily for evolutionary demographic research. Our dataset (Supplementary Data 1) includes nine European populations, including some for which the positive relationship between maternal lifetime twinning status and maternal total births has been described previously5,6,8,10. Details for the populations used in this study are given in Table 1 and in Supplementary Table 14. The sourcing of each dataset and the socio-ecological background of each population have already been described in previous studies (see Table 1 for references). Overall, there is no reason to suspect a high level of consanguinity in these populations62, so our analyses do not account for the variable level of relatedness between individuals. The datasets cover pre-industrial periods in which the lifetime reproductive success was high and the majority of people were living and working in agrarian communities, except for the Samis (from northern Finland and Sweden) who made their living fully or in part from a combination of herding, fishing and hunting. The smooth decline of the probability of parity progression with parity (Fig. 4a) suggests that mothers did not effectively limit their reproductive success with the aim of achieving a small family size, as found in populations that have undergone a demographic transition.Data selectionWe use the term family to describe a mother and all individuals to whom she gave birth over her life. For our analyses, all families considered met the following criteria: the mother’s age was known at a monthly resolution and her life course traced until at least age 45 (approximating full reproductive life), the birth year and month of all offspring must have been recorded and consecutive births were all at least nine months apart from one another. In the case of one population (Norway) and of a few observations in the other populations, the month of birth was not available. These data were thus not considered in the results presented in main text because some of our analyses require an accurate estimation of the interbirth interval. Most analyses are thus based on data from eight populations. Nevertheless, the slope of the negative relationship between twinning and total birth remained very similar irrespectively of whether or not such data were included (Supplementary Fig. 5), which suggests that the exclusion of Norwegian data does not alter our main conclusions. Information on the populations considering also the data for which the birth months were missing is provided in Supplementary Table 14.Twin identificationIn our data, the maximum number of offspring to constitute a multiple birth was three. We use the term twin(s) to refer to offspring who were the result of the same multiple birth (including 1745 sets of twins sensu stricto and 19 sets of triplets in the filtered dataset and, respectively, 1915 and 20 in the dataset including observations that lack birth months). Although twins are sometimes explicitly indicated in the data sources, this is not always the case. Thus, for the sake of consistency across our populations, twin births were identified when at least two individuals born to the same mother appeared with similar birth dates, according to strict criteria: if the exact birth dates were available, then offspring were identified as twins if their birth dates were no more than one day apart. If the exact birth dates were not available, then an identical birth year and month were considered sufficient for positive twin identification.Analyses and simulationsCharacterisation of the relationship between twinning and fertilityWe began by characterising the relationship between lifetime twinning status and maternal total births by fitting two models. For the first, we used a Generalised Linear Mixed-effects Model (GLMM) to investigate whether the mothers of twins (twinners) had experienced a larger or smaller number of births than mothers who only had singletons (non-twinners). We fitted this GLMM on the mother-level data with the R package spaMM63 using the call:$$begin{array}{c}{{{rm{fitme}}}}left({{{{rm{births}}}}}_{{{{rm{total}}}}}sim 1+{{{rm{twinner}}}}+(1|{{{rm{pop}}}}),right.\ {{{rm{data}}}}={{{rm{mother}}}}_{{{rm{level}}}}_{{{rm{data}}}},\ left.{{{rm{family}}}}={{{rm{Tnegbin}}}}({{{rm{link}}}}={hbox{“}}{{{mathrm{log}}}}{hbox{”}})right)end{array}$$
(1)
The response variable births_total refers to all births recorded over a mother’s observed lifetime (count data). The term 1 informs the function to fit an intercept (which happens by default, but is indicated here for clarity). The predictor twinner refers to maternal lifetime twinning status (binary: twinner vs non-twinner) and is modelled as a fixed effect. The term pop refers to the population identity (qualitative variable with eight levels) and is modelled by a Gaussian random effect acting on the intercept, which allows for the modelling of the heterogeneity between populations that is not captured by the fixed effects. The argument family is used to define the error structure and the link function of the GLMM (more on this below).In a second model, we reversed which variable is used as a response and which is used as the fixed-effect predictor. This allowed us to analyse how maternal total births predicted the probability of a mother producing twins during her lifetime using the call:$$begin{array}{c}{{{{{rm{fitme}}}}}}left({{{{{rm{twinner}}}}}} sim 1+{{{{{rm{births}}}}}}_{{{{{rm{total}}}}}}+(1|{{{{{rm{pop}}}}}}),right.\ {{{{{rm{data}}}}}}={{{{{rm{mother}}}}}}_{{{{{rm{level}}}}}}_{{{{{rm{data}}}}}},\ left.{{{{{rm{family}}}}}}={{{{{rm{binomial}}}}}}({{{{{rm{link}}}}}}={hbox{“}}{{{{{rm{logit}}}}}}{hbox{”}})right)end{array}$$
(2)
The fitted models 1 and 2 are depicted in Fig. 1a, b and Supplementary Tables 1,2. While models 1 and 2 represent two sides of the same coin, the fit of both models is justified because each model formulation provides complementary information: expressing the effect of twinning on fertility relates to previous work5,6,7,8,9,10,57 and expressing the effect of fertility on twinning is a first step toward identifying what shapes twinning propensity, the focus of this paper.For the model predicting total births (model 1), we chose to use a negative binomial error structure. Using this error structure produced a fit of the data that was better than a (truncated) Poisson—the usual alternative for count data—as evidenced by much smaller marginal and conditional AIC values64. Here we specifically used a truncated negative binomial distribution because the data do not possess zeros by construction (only mothers are present in the dataset, i.e. there are no nulliparous women). For the model predicting lifetime twinning status (model 2), we chose a binomial error structure which is appropriate for binary data.Modelling the proportion of twin births among all births per mother is an effective way to avoid biases caused by differences in exposure to the risk of having twins affecting the relationship between twinning and fertility. For this, we fitted the following third model:$$begin{array}{l}{{{{{rm{fitme}}}}}}left({{{{{rm{cbind}}}}}};({{{{{rm{twin}}}}}}_{{{{{rm{total}}}}}},{{{{{rm{singleton}}}}}}_{{{{{rm{total}}}}}})right.\ sim ,1+{{{{{rm{births}}}}}}_{{{{{rm{total}}}}}}+(1|{{{{{rm{pop}}}}}}),\ ,{{{{{rm{data}}}}}}={{{{{rm{mother}}}}}}_{{{{{rm{level}}}}}}_{{{{{rm{data}}}}}},\ ,left.{{{{{rm{family}}}}}}={{{{{rm{binomial}}}}}};({{{{{rm{link}}}}}}={hbox{“}}{{{{{rm{logit}}}}}}{hbox{”}})right)end{array}$$
(3)
In this model, the variable twin_total refers to the mother’s total number of twin births (i.e. one for each twinning event), singleton_total refers to the lifetime number of singleton births, and the cbind() function serves to indicate the fitting function to model the frequency of twinning events based on these two variables, which is interpreted as number of successes and failures of a binomial experience. The fitted model 3 is depicted in Fig. 2 and Supplementary Table 3.We modified model 3 so as to test whether the effect of total births differed significantly between populations. To do so, we considered that the effect of populations on total births could either be modelled as an interaction between fixed effects or as a random slope. For the former representation, we thus compared the fit of a model with linear predictor structure defined in spaMM as 1+births_total*pop to that of a model with the structure 1+births_total+pop. For the latter representation, we compared the fit with linear predictor structure 1+births_total + (1|pop) (i.e. model 3 as introduced above) to that of 1+births_total + (1+births_total|pop). We performed this testing procedure by comparing the likelihood ratio between each pair of alternative fits to the expectation of such a ratio under the null hypothesis. The distribution of the statistics used for the test was computed using 1000 parametric bootstrap replicates, which we generated using the function anova() provided by spaMM63. The test revealed a small non-significant variation in slopes between populations (see Results). For the sake of simplicity, we thus considered the effect of births_total the same across populations in all other analyses.Modelling life-history events using GLMMsTo reveal the biological mechanisms responsible for the relationship between twinning and fertility, we first fitted statistical models describing how age, parity and twin/singleton status, as well as individual and population differences influenced three key life-history events: parity progression (PP), the duration of interbirth intervals (IBI) and the twinning outcome of births (T). These models were fitted on birth-level data by the following calls:$$begin{array}{l}{{{{{rm{fitme}}}}}}left({{{{{rm{PP}}}}}} sim 1+{{{{{rm{twin}}}}}}+{{{{{rm{poly}}}}}}({{{{{rm{cbind}}}}}}({{{{{rm{age}}}}}},,{{{{{rm{parity}}}}}}),,{{{{{rm{best}}}}}}_{{{{{rm{order}}}}}})right.\ ,+,(1|{{{{{rm{maternal}}}}}}_{{{{{rm{id}}}}}})+(1|{{{{{rm{pop}}}}}}),\ ,{{{{{rm{data}}}}}}={{{{{rm{birth}}}}}}_{{{{{rm{level}}}}}}_{{{{{rm{data}}}}}},\ ,{{{{{rm{family}}}}}}={{{{{rm{binomial}}}}}};({{{{{rm{link}}}}}}={hbox{“}}{{{{{rm{logit}}}}}}{hbox{”}})end{array}$$
(4)
$$begin{array}{l}{{{{{rm{fitme}}}}}}left({{{{{rm{IBI}}}}}} sim 1+{{{{{rm{twin}}}}}}+{{{{{rm{poly}}}}}};({{{{{rm{cbind}}}}}}({{{{{rm{age}}}}}},,{{{{{rm{parity}}}}}}),{{{{{rm{best}}}}}}_{{{{{rm{order}}}}}})right.\ ,+,(1|{{{{{rm{maternal}}}}}}_{{{{{rm{id}}}}}})+(1|{{{{{rm{pop}}}}}}),\ ,{{{{{rm{data}}}}}}={{{{{rm{birth}}}}}}_{{{{{rm{level}}}}}}_{{{{{rm{data}}}}}},\ ,left.{{{{{rm{family}}}}}}={{{{{rm{negbin}}}}}}({{{{{rm{link}}}}}}={hbox{“}}{rm log} {hbox{”}})right)end{array}$$
(5)
$$begin{array}{l}{{{{{rm{fitme}}}}}}left({{{{{rm{T}}}}}} sim 1+{{{{{rm{poly}}}}}}({{{{{rm{cbind}}}}}}({{{{{rm{age}}}}}},,{{{{{rm{parity}}}}}}),,{{{{{rm{best}}}}}}_{{{{{rm{order}}}}}})right.\ ,+,(1|{{{{{rm{maternal}}}}}}_{{{{{rm{id}}}}}})+(1|{{{{{rm{pop}}}}}}),\ ,{{{{{rm{data}}}}}}={{{{{rm{birth}}}}}}_{{{{{rm{level}}}}}}_{{{{{rm{data}}}}}},\ ,left.{{{{{rm{family}}}}}}={{{{{rm{binomial}}}}}};({{{{{rm{link}}}}}}={hbox{“}}{{{{{rm{logit}}}}}}{hbox{”}})right).end{array}$$
(6)
The response variables of models 4, 5 and 6 are thus PP, IBI and T, which refer to whether the mother went on to reproduce again or not (a boolean), the duration of the interbirth interval between the focal birth and the next (a discrete number of months) and whether the birth resulted in twins or not (a boolean), respectively. In addition to the terms that have already been defined, we now have the term poly(cbind(age, parity), best_order) to code for a polynomial describing the effect of maternal age, parity and their possible interaction. The two-variable polynomial function was applied on maternal age (with a monthly resolution) and parity (i.e. the current birth rank). Such a polynomial term allowed us to explore the influences of maternal age and parity on each response variable while encompassing the non-linearity of these predictors. We also have the predictor variable twin, which is a boolean that indicates if the previous birth event of a given mother resulted in twins or not (the variable twin and T are the same, but we used two different names to clarify when it is used as a response or as a predictor). Finally, we have the random effect “maternal identity” (maternal_id), which is used to represent intrinsic variation among mothers, that is, heterogeneity of expected response among individuals, beyond that due to the fixed effects and the population random effect. This random effect therefore measures maternal intrinsic fertility (in models 4 & 5) and twinning propensity (in model 6).To determine the best polynomial order (best_order) for the polynomial term we attempted orders from 0 to 6 and selected, for each model, the order leading to the model fit associated with the smallest marginal AIC. A polynomial of order 6 is sufficient to fit very complex shapes. Polynomial orders obtained by this procedure are given in the summary tables of the model fits given in Supplementary Tables. Importantly, maternal age and parity are highly correlated together (Spearman’s rho = 0.69), unequally correlated to response variables and exert non-linear effects. These are precisely the conditions in which collinearity issues are the most severe65. This justifies why we considered them jointly in all statistical models, as well as why we did not attempt to partition their respective biological effects in our analyses (except for the visual representation in Fig. 4).In order to study how the lifetime twinning status influenced maternal age at first birth, we also fitted the following model:$$begin{array}{l}{{{{{rm{fitme}}}}}}left({{{{{rm{AFB}}}}}} sim 1+{{{{{rm{twinner}}}}}}* {{{{{rm{births}}}}}}_{{{{{rm{total}}}}}}_{{{{{rm{fac}}}}}}+(1|{{{{{rm{pop}}}}}}),right.\ ,{{{{{rm{data}}}}}}={{{{{rm{mother}}}}}}_{{{{{rm{level}}}}}}_{{{{{rm{data}}}}}},\ ,left.{{{{{rm{family}}}}}}={{{{{rm{negbin}}}}}}({{{{{rm{link}}}}}}={hbox{“}}{rm log} {hbox{”}})right)end{array}$$
(7)
In this model, the response variable AFB corresponds to the age at first birth expressed as a number of months (discrete data) and the predictor births_total_fac corresponds to a qualitative variable referring to maternal total births (10 levels: 1, 2, …, 9, 10 + ). We here considered a possible interaction between twinner and births_total_fac. We used the negative binomial family as in model 1 but as for model 5 there is no need to consider here the truncated form of the distribution. All other terms have already been defined. The fitted model is depicted in Supplementary Fig. 1 and Supplementary Table 7.Marginal predictions for GLMMsAll predictions shown in plots or given in text represent marginal predictions. This means that the predictions for the quantities of interest (maternal lifetime births, twinning probabilities and age at first birth) are a function of coefficients of the fixed effects, and of the variance of the random effects. To be more precise, we averaged, over the fitted distribution of random effects, the predictions expressed on the scale of the response (i.e. back-transformed from the scale of the linear predictor) and conditional on the fixed and random effects. Unlike the traditional conditional predictions computed for a specific value of the random effects (often 0), such computation provides unbiased predictions and should be favoured in the context of GLMMs where random effects act non-additively on the expected response (which is the case when the link function of the model is not identity66). We estimated 95% intervals for these marginal predictions (CI95%) using parametric bootstraps with the help of the function spaMM_boot() from the R package spaMM and boot.ci() from the R package boot. More details can be found by looking at the code of the functions compute_predictions() and compare_prediction() in our supporting R package twinR (see Code availability).Simulating the life history of mothersWe produced an individual-based simulation model of human female life history to investigate the contribution of four mechanisms to the relationship, shown in Fig. 2, between per-birth twinning probability and maternal total births—an approach generally known as pattern-oriented modelling67. Each simulation proceeds in the following way: first, we initialised the simulation with representations of the exact same mothers present in the observed dataset, setting their population and maternal identities as the real ones, their starting ages at the observed values for age at first birth and their parity to one. Following this initialisation, the virtual lives of mothers proceeded as multiple iterations of a sequence of three life-history events, informed by statistical models (see below) and subject to the hypothesis being tested (Supplementary Figs. 3, 4). Specifically, for each mother, the twin/singleton status (T) of the current birth was first determined using a GLMM predicting T. Then, whether or not she will go on to reproduce at least once more was determined by simulating her parity progression status (PP) using a GLMM predicting the parity progression probability. For mothers who do continue reproducing, we finally used a third GLMM to determine the length of the interval between the current birth and the next one (IBI). At each iteration, a mother’s parity is increased by one, and age is increased by the simulated length of the interbirth interval. All predictions were performed conditionally on the value for the predictor characterised by both fixed and random effects. The process of simulating PP, IBI and T was then reiterated until all women had ceased reproduction, which happens necessarily since the probability of parity progression is lower than one. We also set this probability to zero once women reached 60 years old to save computation time in particular simulation scenarios leading to unrealistic life histories (and bad goodness of fit). Note that the maximum recorded age at which a mother gave birth was 55.1 years in our data. For the same reason, we also capped the maximum simulated duration for the interbirth interval to 30 years.Drawing life-history events from the fit of the model formulas shown above for models 4, 5 and 6 corresponds to simulating the scenario PISH (i.e. all four hypothetical mechanisms are activated). For simulating other simulation scenarios, we had to fit additional GLMMs derived from models presented above. Specifically, the term twin was dropped from model 4 to deactivate mechanism P (model 8; Supplementary Table 8); the term twin was dropped from model 5 to deactivate mechanism I (model 9; Supplementary Table 9); the term poly(cbind(age, parity), best_order) was dropped from model 6 to deactivate the mechanism S (model 10 and 11; Supplementary Tables 10, 11); and the term (1|maternal_id) was dropped from model 6 to deactivate mechanism H (model 11 and 12; Supplementary Tables 11, 12).Testing candidate mechanisms using simulationsTo test how each mechanism or association of mechanisms influenced the relationship between twinning and fertility, we ran simulations under each possible set of activated or inactivated mechanisms. We tested all possible sets and we thus built a total of 42 = 16 simulation scenarios (Supplementary Figs. 3, 4).For each simulation scenario, we ran simulation replicates (see Supplementary Notes for details and information on the numbers of replicates), then fitted model 3 on the dataset produced by each replicate and extracted the estimate for the slope associated with the term births_total in that model (β*). We then consider, in turn, that each simulation scenario may have generated the data. Each simulation scenario is thus considered as a null hypothesis which we aim at testing. Such a test is traditionally referred to as a goodness-of-fit test68. The result of such a test, a p-value, answers the question: what is the probability of obtaining a value equal to, or more extreme than, the statistic of interest, if the null hypothesis were true? The rejection of the null hypothesis by the test (i.e. a p-value ≤ 0.05) signifies a rejection of the null hypothesis, and thus, here, the rejection of a simulation scenario which represents a particular mechanism, or combination of mechanisms. In contrast, a large (i.e. non-significant) p-value would here denote support for the simulation scenario under consideration.A first candidate, as a statistic of interest to build our goodness-of-fit test, is the slope β*. However, when viewed as a goodness-of-fit test, the direct comparison of the observed and simulated slopes may be conservative when other life-history parameters are fitted to the data. This is because the data tend to be more likely given parameter values fitted to the data than given the actual (unknown) parameter values that generated the data. The goodness-of-fit test is however only guaranteed to provide uniformly distributed p-values (a feature necessary for the correctness of any null hypothesis testing) when samples are drawn under the latter parameter values. This is a general issue in statistics which has also been discussed long ago, for example, when the data-generating process is the normal distribution and a Kolmogorov–Smirnov test of goodness-of-fit is applied69. We thus designed and validated a specific procedure to correct for such bias while testing each simulation scenario (Supplementary Notes). In the text, we only report outcomes from this unbiased goodness-of-fit test (for details, see Supplementary Notes and Supplementary Table 13).Studying the effect of twinning propensity on the number of offspring using simulationsTo study how twinning propensity influences the total number of offspring that mothers produced during their lifetime, we ran two sets of simulations, each with 100 replicates. In the first set, we ran the simulation as described in the section “Simulating the life history of mothers” using the fits of the models associated with the simulation scenario PIS (i.e. fits of models 4, 5 and 12). In the second set, we did the same, except that we modified the intercept of the model predicting twinning events (fit of the model 12) by adding 2.5 to its intercept. We also tried other values, some smaller (e.g. 0.25), some larger (e.g. 5), to make sure that the magnitudes of the change of the intercept did not impact our qualitative statements. For each set of replicates, we extracted the twinning rate, the twinner rate, the mean number of offspring produced and the mean total number of births. We report the means of these metrics, as well as the 95% Central Range from simulation replicates (CR95%), which we directly computed by extracting the corresponding quantiles from the distribution generated by the replicates.Realism of the simulationsWe checked that the simulated life history closely matched that of the real mothers represented in our dataset beyond what is captured by the relationship between twinning and fertility. To do so, we compared different metrics related to fertility and twinning between the real and simulated data. We chose to perform this comparison under the simulation scenario PIS since it produces the best goodness-of-fit. The results of this quality check confirm that our simulations represent the reproductive lives of the mothers appropriately (Supplementary Fig. 6).Studying the effect of mortality on the number of offspring using simulationsTo account for the fact that not all offspring have the same expected survival, we also applied a survival weight to each simulated offspring before averaging the numbers for a given simulation set (baseline twinning propensity or enhanced, see Results). We used as weights the estimates for the probability of offspring survival between birth and adulthood provided by two publications associated with some of the data we used there. Specifically, following Helle et al.8, we used a weight of 0.603 for twins, 0.838 for singletons from twinners and 0.815 for singletons from non-twinners. Alternatively, following Haukioja et al.3, we used a weight of 0.337 for twins and 0.706 for singletons from all mothers.Implementation detailsAll statistical analyses were performed in R version 4.170. The main R packages we used were spaMM63 version 3.9.40 for the fit of all the statistical models, boot71,72 version 1.3-28 for the computation of confidence intervals based on parametric bootstraps, and R673 version 2.5.1 for defining the object used to run the simulations. The DESCRIPTION file from our package twinR (see Code availability) lists the additional R packages required for this project (e.g. those used for plotting and data manipulation).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More