Philippa Kaur

Evolution of ACE2 in bats inhabiting urban or rural areas
We collected ACE2 orthologues from 46 bat species across the phylogeny (Fig. 1 and Supplementary Table 1). These species contained 28 species that roost or forage in urban areas near humans and 18 species more restricted to rural areas and hence likely to have minimal contact with humans (Supplementary Table 2). In total, we examined 46 species representing 11 bat families (Supplementary Table 3). After aligning the protein sequences of bat ACE2 orthologues, we examined 25 critical residues involved in the binding of the surface spike glycoprotein (S protein) of SARS-CoV-2 (ref. 9) (Extended Data Fig. 1). Genetic variations were observed in nearly all these 25 sites, which may have led to different abilities to support entry of SARS-CoV and SARS-CoV-2 (ref. 9). Furthermore, we detected at least 22 amino acid sites that are putatively under positive selection (Supplementary Table 4), which is indicative of heterogeneous selection pressure across sites. Notably, four of these positively selected sites are in the binding region of ACE2 to the SARS-CoV-2 S protein (Supplementary Table 4).
Fig. 1: Phylogenetic tree of 46 bat species in this study.

The labels of bat species in our experiments are indicated. Expression levels determined by western blot (Fig. 2a) are shown with asterisk symbols compared with human ACE2: the triple asterisk indicates high expression, the double asterisk indicates medium expression and the single asterisk indicates low but detectable expression. The ability of bat ACE2 to support SARS-CoV and SARS-CoV-2 pseudovirus entry is shown with different signs (Fig. 3a,b): infection data are presented as percentage mean values of bat ACE2 supporting infection compared with the infection supported by human ACE2. Infection efficiency 50% with a double plus sign. Bat phylogeny was taken from previous studies28,29,30.

Full size image

Interaction between bat ACE2 orthologues and SARS-CoV or SARS-CoV-2 receptor binding domain
Efficient binding between the S protein and the ACE2 receptor is essential for SARS-CoV and SARS-CoV-2 entry. This binding is mainly mediated by the interaction between the critical residues on the receptor-binding domain (RBD) and ACE2. To characterize the receptor function of ACE2 orthologues in a range of diverse bat species, we generated a stable cell library consisting of cell lines expressing the respective 46 bat ACE2 orthologues through lentiviral transduction of 293T cells lacking ACE2 expression10. All bat ACE2 orthologues were exogenously expressed at a comparable level after puromycin selection, as indicated by western blot and immunofluorescence assays detecting the C-terminal 3×FLAG-tag (Fig. 2a,b).
Fig. 2: Expression of bat ACE2 orthologues and their interaction with the SARS-CoV and SARS-CoV-2 RBD.

a, Western blot detected the expression levels of ACE2 orthologues on 293T stable cells by targeting the C-terminal 3×FLAG-tag. Glyceraldehyde 3-phosphate dehydrogenase was employed as a loading control. b, Visualization of the intracellular bat ACE2 expression level by immunofluorescence assay detecting the C-terminal 3×FLAG-tag. Scale bar, 100 μm. c,d, Assessment of the interaction between different ACE2 orthologues and SARS-CoV-RBD-hFc (c) or SARS-CoV-2-RBD-hFc (d) proteins. Species that do not support efficient binding are underlined. 293T cells stably expressing the different bat ACE2 orthologues were incubated with 5 μg ml−1 of the recombinant proteins at 37 °C for 1 h; binding efficiency was examined by Alexa Fluor 488 goat anti-human IgG via fluorescence assay. Scale bar, 200 μm.

Full size image

To analyse the interaction, we produced recombinant SARS-CoV or SARS-CoV-2 RBD human immunoglobulin G (IgG) Fc fusion proteins (RBD-hFc), previously reported to be sufficient to bind human ACE2 efficiently11,12. Protein binding efficiency was tested on the bat ACE2 cell library by means of immunofluorescence or flow cytometry targeting the human Fc. As expected, binding was almost undetectable on mock 293T cells but a strong binding signal was detected in the 293T cells expressing human ACE2 (Fig. 2c,d). Consistent with previous reports13,14, SARS-CoV-2 RBD showed higher binding to human ACE2 than SARS-CoV, which can also be observed on many bat ACE2 orthologues (Fig. 2c,d). Previous reports have shown that only a small fraction of ACE2 orthologues from tested mammalian species could not bind with SARS-CoV-2 S protein (n = 6 of 49 species7; n = 5 of 17 species15). However, our study revealed that many bat species (n = 32 and n = 28 of 46 species) do not support efficient binding with SARS-CoV-RBD and SARS-CoV-2-RBD, respectively (Fig. 2c,d). The overall profiles of bat ACE2 to bind to SARS-CoV and SARS-CoV-2 RBD are generally comparable; a few showed contrasting modes of binding preferences (Fig. 2c,d). For instance, Bat22 could bind to SARS-CoV but not SARS-CoV-2, whereas Bat14, 21 and 40 could bind to SARS-CoV-2 but not SARS-CoV (Fig. 2c,d). Flow cytometry analysis showed consistent results (Extended Data Fig. 2).
Overall, the RBD-hFc binding assays demonstrated that bat ACE2 orthologues showed different affinity and selectivity levels to SARS-CoV and SARS-CoV-2, indicating that the ACE2 receptors of many bat species may not support efficient SARS-CoV and SARS-CoV-2 infection.
Receptor function of bat ACE2 orthologues to support the entry of SARS-CoV and SARS-CoV-2 using pseudotyped and live viruses
To further evaluate the receptor function of different bat ACE2 orthologues, we employed a vesicular stomatitis virus (VSV)-based rhabdoviral pseudotyping system to mimic the coronavirus spike protein-mediated single-round entry15. SARS-CoV and SARS-CoV-2 pseudotypes were generated by assembling the coronavirus spike proteins and replication-deficient VSV with the VSV glycoprotein gene replaced with a fluorescence protein (VSV-dG-GFP) or a firefly luciferase (VSV-dG-Luc) reporter15. Both viruses showed minimal background infection on 293T cells, but efficient infection on 293T-human ACE2 cells (Extended Data Fig. 3). The susceptibility of the 293T cells expressing bat ACE2 orthologues was then examined with SARS-CoV and SARS-CoV-2 pseudotypes. The results showed that bat ACE2 orthologues have varying abilities to support coronavirus entry and different preferences for SARS-CoV and SARS-CoV-2. (Fig. 3a,b and Extended Data Fig. 4). Pseudotypes with green fluorescent protein (GFP) reporter showed similar results (Extended Data Fig. 5). Notably, we found that 24, 21 and 16 of the 46 bat species showed almost no entry for SARS-CoV, SARS-CoV-2 and both viruses, respectively (Figs. 1 and 3a,b and Supplementary Table 5), suggesting that these species are not likely to be potential hosts of either or both coronaviruses. The bat species showing no viral entry include those that occur in urban areas and those more restricted to rural areas (Fig. 1), suggesting that there is no correlation between proximity to humans and probability of being natural hosts of SARS-CoV or SARS-CoV-2. Although horseshoe bats were suggested as potential natural hosts of SARS-CoV and SARS-CoV-2 (refs. 1,2,3), only one of the three species examined (Rhinolophus sinicus) supported SARS-CoV entry; this species was suggested as the potential host of SARS-CoV3,16. None of these tested horseshoe bats showed entry for SARS-CoV-2 (Figs. 1 and 3). These results unambiguously indicate that ACE2 receptor usage is species-dependent.
Fig. 3: Characterization of bat ACE2 orthologues mediating entry of SARS-CoV and SARS-CoV-2 viruses.

a,b, Ability of bat ACE2 orthologues to support the entry of SARS-CoV and SARS-CoV-2 pseudovirus. 293T cells expressing bat ACE2 orthologues in a 96-well plate were infected with VSV-dG-Luc pseudotyped with SARS-CoV (a) and SARS-CoV-2 (b) spike proteins, respectively. Intracellular luciferase activity was determined at 20 h post-infection. RLU, relative light unit. c, 293T cells expressing bat ACE2 orthologues were inoculated with the SARS-CoV-2 live virus at an MOI = 0.01. N protein (red) in the infected cells was detected through immunofluorescence assay at 48 h post-infection. Scale bar, 200 μm. Samples expressing the indicated ACE2 orthologues that showed almost no entry for SARS-CoV-2 live virus are underlined. Data shown are representative results from 3 independent experiments and are presented as the mean ± s.d. (n = 3 for a and n = 2 for b).

Full size image

The SARS-CoV-2 S protein used in this study for pseudotyping contains a D614G mutation, which is currently a dominant variation17. The D614G mutation remarkably improved the in vitro infectivity of SARS-CoV-2 but may not significantly affect the receptor interaction since it is not in the RBD18. Indeed, we identified a very similar susceptibility profile using an original strain without D614G (Extended Data Fig. 4). We further demonstrated that the pseudotyped entry assay mimics the entry of live viruses through a SARS-CoV-2 infection assay (Fig. 3c). As expected, the profile of SARS-CoV-2 N protein expression is highly consistent with the results from the VSV-dG-based pseudotyped virus entry assay, except for some ACE2 that showed relatively higher infection efficiency (for example, Bat43–46) compared with the pseudovirus infection assay, which may be attributed to the different virus strains used (Fig. 3c). In addition, the live virus infection resulted in the phenotype of plaque formation, while the pseudotypes showed evenly distributed, single-round infection (Extended Data Fig. 5), which also partially explains why some bat ACE2 showed higher infection in the live virus infection assay.
When comparing the RBD-hFc binding and pseudotyped entry profiles, we found that binding and susceptibility are not always consistent, although the phenotypes were reproducible. For instance, some species (Bat12, 13, 14) were able to bind to SARS-CoV-2 RBD-hFc efficiently but could not support infection of the same virus, indicating that high binding affinity does not guarantee efficient viral entry (Figs. 2 and 3). In contrast, some species (Bat3–8) were defective or less efficient in SARS-CoV RBD-hFc binding but supported the entry of the same virus to some degree (Figs. 2 and 3). We hypothesize that such minimal binding may be sufficient for viral entry mediated by those ACE2 orthologues; alternatively, additional residues outside the traditional RBD region might be required for efficient interaction. These hypotheses should be tested in the future. Together, our results demonstrated dramatic variation of susceptibility to SARS-CoV and SARS-CoV-2 infection among bat species, suggesting that SARS-CoV and SARS-CoV-2 can selectively use some bat ACE2 as functional receptors for viral entry and many—if not most—bat ACE2 are not favoured by one or both viruses.
Evaluation of critical residues in bat ACE2 orthologues affecting viral binding and entry efficiency or specificity
We comprehensively analysed the relationship between critical RBD binding sites in bat ACE2 sequences and their ability to support SARS-CoV and SARS-CoV-2 RBD binding and viral entry. Several critical residues were identified that may play critical roles in the determination of species specificity (Extended Data Fig. 1). According to the sequence alignment, two species pairs (Bat33 and Bat34 and Bat38 and Bat40) were selected to demonstrate the role of critical residues in RBD binding and viral entry because they were phylogenetically close but showed contrasting phenotypes for supporting RBD binding and viral entry. Specifically, Bat34 and Bat38 do not support SARS-CoV and SARS-CoV-2 RBD binding and infection, while Bat33 supports efficient binding and infection of both viruses and Bat40 supports infection of both viruses and to a lesser degree SARS-RBD binding (Figs. 2 and 3). We compared their protein sequences and highlighted the residues that may affect RBD interaction. For example, substitutions I27K, N31G and K42E were observed when comparing Bat33 with Bat34, while Q24L, E30K, K35Q and G354N were present between Bat38 and Bat40 (Fig. 4a). We hypothesized that the discrepancy in binding and infection phenotype is determined by their differences in critical residues for RBD interaction. To test this hypothesis, we designed a residue swap mutagenesis assay to investigate the role of critical residues on RBD binding and virus entry (Fig. 4a). We generated four swap mutations and corresponding 293T stable cell lines to test whether these substitutions could achieve gain-of-function and loss-of-function. All bat ACE2 orthologues and related mutants were expressed at a comparable level after lentiviral transduction, as indicated by the immunofluorescence of the C-terminal 3×FLAG-tag (Fig. 4b). Recombinant SARS-CoV and SARS-CoV-2 RBD-hFc proteins were applied to the cells expressing different ACE2 and binding efficiency was evaluated by immunofluorescence (Fig. 4c) and flow cytometry assays (Fig. 4d). As expected, the swap of critical residues on the selected four bat ACE2 changed their receptor function to the opposite, except for Bat38 mutant, which remained unable to bind SARS-CoV RBD-hFc (Fig. 4c,d). GFP (Fig. 4e) and luciferase levels (Fig. 4f) from the pseudotyped virus entry assay and the N protein staining from the live SARS-CoV-2 infection assay (Fig. 4g) further confirmed our hypothesis at the viral entry level. Structure modelling of bat ACE2/SARS-CoV-2-RBD complexes showed that the substitutions of I27K and N31G between Bat33 and Bat34 lead to a reduced packing interaction and the substitution of K42E disrupts the hydrogen bond with Y449, which may be related to the difference of susceptibility between Bat33 and Bat34 (Fig. 4h,i and Extended Data Fig. 6). In comparison, the substitutions of Q24L and E30K between Bat38 and Bat40 destroyed the favourable hydrophilic interactions with N487 and K417, respectively (Extended Data Fig. 6).
Fig. 4: Evaluation of the critical binding sites determining the species-specific restriction of SARS-CoV and SARS-CoV-2 binding and entry.

a, Swap mutagenesis assay to investigate the role of critical residues on bat ACE2 orthologues for tropism determination. Residues involved in RBD (according to the structure between SARS2-RBD and human ACE2, Protein Data Bank 6M0J) interaction are shown in the table. Residues that changed in the mutagenesis assay are marked in red. b, The expression level of the bat ACE2 orthologues and related mutants in transduced 293T cells was determined by an immunofluorescence assay recognizing the 3×FLAG-tag. Scale bar, 200 μm. c,d, Binding efficiency of SARS2-RBD-hFc and SARS2-RBD-hFc on 293T cells expressing bat ACE2 and related mutants. Cells were incubated with 5 μg ml−1 of recombinant proteins at 37 °C for 1 h and then washed and incubated with a secondary antibody recognizing human Fc. Immunostaining (c) and flow cytometry (d) were conducted to show binding efficiency. Scale bar, 200 μm. e,f, Ability of the indicated ACE2 and related mutants to support the entry of coronavirus pseudotypes. The 293T cells expressing the indicated ACE2 and their mutants were infected with SARS-CoV and SARS-CoV-2 pseudotypes expressing GFP (e) and luciferase (f). Infection was analysed at 20 h post-infection. Scale bar, 200 μm. Data are presented as the mean with s.d. (n = 2). g, 293T cells infected by the SARS-CoV-2 live virus at an MOI = 0.01; the infection was examined at 48 h post-infection through N protein (red) immunostaining. Nuclei were stained with Hoechst 33342 (blue). Scale bar, 200 μm. h,i, Comparison of the interface between Bat33/SARS-CoV-2-RBD and Bat34/SARS-CoV-2-RBD. Bat33 and its complexed RBD are coloured cyan and gold, respectively (h); Bat34 and its complexed RBD are coloured wheat and green, respectively (i). The mutated residues in ACE2 and the corresponding residues in SARS-CoV-2-RBD are shown and labelled. The red dotted lines between residues indicate hydrogen or ionic bonds.

Full size image

In addition, two bat cell lines, the lung epithelial cell line Tb 1 Lu of Tadarida brasiliensis (Bat31) and the kidney epithelial cell line of Pteropus alecto (Bat2), were used to validate our findings derived from human HEK293T cells. Endogenous ACE2 expression was almost undetectable in these two cell lines, accounting for at least 1,000 folds lower than the susceptible Vero-E6 cells (Extended Data Fig. 7a). Therefore, these cells cannot support the entry of SARS-CoV and SARS-CoV-2. We successfully generated Tb 1 Lu stable cell lines expressing human ACE2 and bat ACE2 (Bat2, 3, 31, 32) since the transduction efficiency of Tb 1 Lu is much higher than that of PakiT03 cells (Extended Data Fig. 7b). As expected, Tu 1 Lu were susceptible to both SARS-CoV and SARS-CoV-2 when human ACE2 or some bat ACE2 orthologues (Bat2, 3 and 31) were expressed, yet remained non-susceptible when an ACE2 of a closely related species (Bat32) was expressed (Extended Data Fig. 7c–e). Furthermore, we conducted SARS-CoV and SARS-CoV-2 pseudovirus entry assays on the two bat cell lines transiently transfected with various bat ACE2 (Bat2, 3, 31, 32, 33, 34, 38, 40) and their mutants (mutant Bat33, 34, 38 and 40m). The results were consistent with those derived from human cells, further confirming that ACE2 is the main receptor for the species-specific entry of SARS-CoV and SARS-CoV-2 in these bat cells (Extended Data Fig. 7f,g). More

Genome sequencing and assembly
We performed a de novo genome sequencing and assembly of CBB using a hybrid approach by combining 454-FLX and Illumina reads from female and male individuals. A total of 3.02 Gb of high-quality 454-FLX sequences and 26 Gb of Illumina sequences were obtained in this study (Table S1), which represent approximate 19 × and 160 × genome coverage respectively based on a previously estimated CBB genome size of 163Mb21. The genome hybrid assembly approach we used involved an initial pre-assembly of the 454FLX data with Newbler and the Illumina data with ABySS22, followed by merging of these two pre-assemblies into a single genome consensus with Metassembler23. Our final hybrid H. hampei CENICAFE_Hham1.1 (Hham1.1) genome assembly had a size of 162.57 Mb, comprising 8198 genome scaffolds (Table 1). This assembly represents an improvement in sequence contiguity, containing a 36.3-Kb contig-N50; 340.2-Kb scaffold-N50 and 4.9 Mb for the largest genome scaffold, compared with a previously published CBB genome assembly21, which resulted in contig and scaffold N50 of 10.5-Kb and 44.7-Kb respectively and largest genome scaffold of 440-Kb. The Hham1.1 genome assembly completeness was assessed using Benchmarking Universal Single-Copy Orthologs (BUSCO)24. BUSCO recovered 98.22% of the 1066 Arthropoda core gene set, from which 96.25% were complete genes and 2% were fragmented genes (Fig. S1). BUSCO results indicate that almost the entire genome of H. hampei was sequenced and de novo assembled in this study.
Table 1 Hypothenemus hampei genome assembly (CENICAFE_Hham1.1) statistics.
Full size table

Transcriptome assembly
Illumina RNA-seq data obtained from whole-body female and male adults were de novo assembled using rnaSPades25 and sequence redundancy reduced by CD-HIT26. The resulting transcript assembly was composed of 64,244 contigs (available at NCBI TSA accession: GIPB00000000.1). The average transcript length was 1103-bp, transcript N50 of 2145-bp and largest transcript of 26,019-bp. The transcript assembly completeness with BUSCO recovered 99.6% (98.97% completed and 0.65% fragmented genes) of the 1066 Arthropoda core gene set. (Fig. S1). Using TransDecoder27, we extracted 35,558 protein-encoding transcripts with full Open Reading Frames (ORFs), from which 33,378 (95%) were annotated against InterPro and NCBI NR proteins. As expected, top BLAST hits were against the Coleoptera species, including D. ponderosae (61%) Sitophilus orizae (22%), Anoplophora glabripennis (3%) and Tribolium castaneum (5.7%); whereas the remaining hits were against other insect species (14%).
Gene prediction and functional assignations
We identified 18,765 gene models encoding 20,801 proteins on the Hham1.1 genome assembly using BRAKER2 gene predictor and all available RNA-seq evidence for H. hampei at NCBI. The number of gene models found here for our Hham1.1 assembly is slightly smaller than the previous gene prediction (19,222) performed on the first published H. hampei genome draft21. Completeness of the Hham1.1 gene set using BUSCO recovered 97.2% (94.1% completed and 3.1% fragmented genes) of the Arthropoda core gene set (Fig. S1). BLASTP found 18,364 (88.3%) Hham1.1 predicted proteins similar (e-value  More

Study site
Small passerines were caught using a 13 m mist net placed between trees and shrubs on the edge of the village of Tarbet, Argyll & Bute (56.21 N 4.71 W), adjacent to a large forestry plantation. Siskins can forage up to at least 5 km from their nest during the breeding season5,18, and the area within 5 km was therefore considered likely to include breeding siskins that would move through the catching site. Plantation forestry species, age, and area were determined from maps from the Forestry Commission compartment data base. There were 1152 ha of plantation forestry within 5 km of the catching site, comprising 79% Sitka spruce, 7% Norway spruce Picea abies, 13% larch and  More

National yield stability
We used the FAOSTAT database (http://www.fao.org/faostat, visited in September 2019) to obtain data on annual crop production (in tons) and area harvested (in hectares) from 1961 to 2010 for 138 crops in 91 populous nations. Following Renard and Tilman16, we accounted for differences among nations in data quality and excluded five nations, namely North Korea, Guinea, Kenya, Mozambique and Zambia, for which at least 20% of the data on area harvested or production were extrapolated by the FAO (see details in16). We calculated for each nation and each year the total annual caloric yield (millions of kcal ha-1). To do so, we first calculated the kcal production of each crop by multiplying the production of each crop by its commodity-specific kilocalorie conversion factor from the USDA Nutrient Database32. In doing so, we were able to compare the production of different crops. Then, we summed these kcal harvests across all crops and divided this value by the sum of harvested area for all crops. We calculated national yield stability (S) as the ratio of mean total annual caloric yield (µT) over its time-detrended standard deviation (σT) for fifty consecutive years (1961–2010). We accounted for a temporal trend of increasing total annual crop yield by implementing a loess regression between annual crop yield and years. σT corresponds to the standard deviation of the residuals of this regression. Finally, we compared this stability index (largely used in the biodiversity-ecological functioning research, e.g.14,16,17,18) with the resilience index used by Zampieri et al.22. Both indices were strongly correlated (r = 0.992), strengthening our findings.
Individual crop yield stability and yield asynchrony
For each country, we quantified the average stability of yields of individual crops as the mean of the inverse of the coefficient of variation of yield of each crop:

$$ {{left( {mathop sum limits_{i = 1}^{N} frac{{mu_{i} }}{{sigma_{i} }}} right)} mathord{left/ {vphantom {{left( {mathop sum limits_{i = 1}^{N} frac{{mu_{i} }}{{sigma_{i} }}} right)} N}} right. kern-nulldelimiterspace} N} $$
(1)

where (mu_{i}) is the temporal mean of crop’s annual kcal yield and (sigma_{i}) its time-detrended standard deviation. Time-detrended crop yield was computed through a loess regression between individual, annual crop yield and years.
We computed the asynchrony between crop yield fluctuations following the index developed by Loreau and De Mazancourt11:

$$ Phi = 1 – frac{{sigma^{2}_{T} }}{{left( {mathop sum nolimits_{i = 1}^{N} sigma_{i} } right)^{2} }} $$
(2)

where Φ is the asynchrony of crop species based on annual caloric yield (millions of kcal ha−1) with (sigma_{T}^{2}) the temporal variance of the time-detrended national yield and (sigma_{i}) the time-detrended standard deviation of each crop’s annual kcal yield. The value of asynchrony varies between zero (perfect synchrony) and one (perfect asynchronous temporal fluctuations).
To test whether yield fluctuations of the most abundant crops have a greater impact on the stability of national food production, we weighted the annual yield of each crop by the proportion of total harvested area occupied by that crop. Average stability of yields of individual crops and yield asynchrony were computed on both the non-weighted and abundance-weighted yields.
Crop diversity
For each country and year, we used both the total number of crop commodities (i.e. crop richness) and the Shannon information index (H′) to quantify crop diversity. H′ weights each crop in a nation by the proportion of total cropland it occupies (pi):

$$ H^{prime } = – mathop sum limits_{i = 1}^{N} left( {p_{i} lnleft[ {p_{i} } right]} right) $$
(3)

with N being the total number of crops grown in a country each year.
The exponential form of the Shannon diversity index gives the effective crop diversity that is the number of crops representing an equal share of harvested area24. In other words, the exponential of the Shannon diversity index weighs all species by their frequency, without favouring either common or rare species24. We averaged the annual effective diversity of crop across the fifty years studied to test the effect of crop diversity on national yield stability.
Agricultural inputs
We extracted the annual national application of nitrogen and the annual cropland area equipped for irrigation from the FAOSTAT database. Because Ireland, New Zealand and Netherlands use much of their fertilizers on pastures rather than croplands, we excluded these nations from our analysis. Similarly, we excluded Egypt because it has 100% of cropland equipped for irrigation. We calculated the annual rates of nitrogen application and irrigation per hectare by dividing their use by the total annual cropland area.
Climate variability
We used global gridded climatic data from the Climate Research Unit of the University of East Anglia33 to compute the year-to-year variability of growing season precipitation and temperature for each country, both strongly affecting the stability of national food production16. From these data, we derived annual precipitation and temperature for each grid cell in a country by taking the sum of monthly precipitation and the mean of monthly temperature values weighted by the proportion of cropland in each grid cell34. We then computed the year-to-year coefficient of variation of cropland-based temperature and precipitation for each country.
Statistical analysis
We used structural equation models (SEMs) to evaluate how irrigation, intensity of use of nitrogen fertilizers and crop diversity affected national yield stability through changes in the average stability of yields of individual crops and asynchrony of yields. SEMs represent a powerful way to disentangle complex mechanisms controlling crop diversity-stability relationships, as previously done in natural ecosystems (e.g.14,15,35,36). We set up two different structural equation models, one based on non-weighted indices of stability of individual crops and asynchrony, the other based on the same indices weighted by the proportion of total harvested area accounted for by each crop. We firstly considered the effects of agricultural inputs and crop diversity on the stability of national food production via the path of average yield stability. The second path quantified the indirect effects of agricultural inputs and crop diversity on national stability via their impacts on crop yield asynchrony. We also accounted for the direct effects of agricultural inputs and crop diversity on national yield stability. Finally, we controlled for the effects of climate variability on total, national yield stability, individual crop yield stability and yield asynchrony. SEMs were run with the lavaan R library37. We used the standardized estimates to compare the relative importance of the different paths. The model fit was evaluated using the Fisher C’score and its associated p values. Because the structural equation model assumes linear relationships between predictors and the dependent variable, we also plotted the relationships between total national yield stability and both asynchrony and average stability of individual crop yield to control for linearity (Fig. 2). Similarly, we investigated the relationships between crop diversity and asynchrony (Fig. 3), as well as between irrigation rate and the average stability of individual crop yield (Fig. 4). More

MDA-MB-231 generally out-compete MCF-7 cells under all pH, glucose and glutamine levels, and starting frequencies
Fluorescently labeled MCF-7 (GFP) and MDA-MB-231 (RFP) cells allowed us to non-destructively measure population dynamics within each micro-spheroid. We used this spheroid system to: (1) observe the outcome of competition (mixed culture), (2) compare population growth models (mono-cultures), (3) measure intrinsic growth rates (DI), (4) measure carrying capacities (DD), and (5) determine competition coefficients (FD). In experiments 1 and 2, spheroids started with a plating density of 20,000 and 10,000 cells, respectively, with initial conditions of 0%, 20%, 40%, 60%, 80% and 100% MDA-MB-231 cells relative to MCF-7 cells. In experiment 1, additional treatments included neutral (7.4) versus low (6.5) pH and varying levels of glucose (0, 1, 2, 4.5 g/L) (Fig. 1b). Experiment 2 had treatments of pH (7.4, 6.5), glucose (0, 4.5 g/L), and glutamine positive and starved media (Fig. 1c). The ranges of these values encompass both physiologic and extreme tumor conditions. For example, pH ~ 6.5 and very low levels of glucose have been measured near the necrotic core of the tumor. Conversely, high glucose and physiologic pH are present at the tumor edge and near vasculature, and also reflect standard culture conditions for our cell lines. In experiment 1, microspheres were grown for 670 h with fluorescence measured every 24–72 h. In experiment 2, microspheres were grown for 550 h with fluorescence measured every 6–12 h. Growth medium was replaced every 4 days without disturbing the integrity of the microsphere.
Counter to our hypothesis, MDA-MB-231 outcompete MCF-7 cells in physiologic (2 g/L glucose, 7.4 pH) conditions (Fig. 2). However, in more acidic conditions (2 g/L glucose, 6.4 pH) the MCF-7 cells appear to persist, suggesting coexistence (Fig. 3, Supplementary Fig. 1). Furthermore MCF-7 cells persist in either pH condition with high 4 g/L glucose (Fig. 4). This raises two questions: (1) Is the success of MDA-MB-231 in most conditions a consequence of higher initial growth rates, increased carrying capacity, or consistently favorable competition coefficients? and, (2) what characteristics of MCF-7 cells allow them to competitively persist at high glucose conditions?
Figure 2

Photographs of spheroids illustrative of the progression and outcome of competition in normal pH. This series shows the progression of spheroids grown physiological conditions (2 g/L glucose, 7.4 pH). Note that, in this case, the MDA-MB-231 cells appear to take over the spheroids by day 20.

Full size image

Figure 3

Photographs of spheroids illustrative of the progression and outcome of competition in acidic pH. This series shows the progression of spheroids grown in acidic conditions (6.5 pH) with 2 g/L glucose. Note that, in this case, the MCF-7 cells appear to persist to the end of the experiment.

Full size image

Figure 4

Three time points spanning from the middle to the end for all culture conditions.

Full size image

Analysis of growth rate
The growth rate of MCF-7 and MDA-MB-231 cells was quantified by measuring the instantaneous exponential growth rate (in units of per day) for each monoculture during the first 72 h. In general, MCF-7 cells had equal or higher growth rates under all growth conditions when compared to MDA-MB-231 cells. MCF-7 cell growth rate with lower (10,000 initial cells) plating density was reduced when grown in low pH 6.5 (Fig. 5b,c) environments compared to physiologic pH 7.4 (Fig. 5e,f). The higher (20,000 initial cells) plating density resulted in reduced growth rates in normal pH compared to low 6.5 pH (Fig. 5a,d). Not surprisingly, in the harshest environment of low pH 6.5, no glucose, and no glutamine, MCF-7 cells had the lowest, even negative, growth rate (Fig. 5c). Interestingly, the absence of glutamine (Fig. 5c,f) reduced the difference in growth rates between the two cell types, regardless of pH or glucose concentration, consistent with glutamine’s role as a substrate for respiration, which is more pronounced in MCF-7 cells.
Figure 5

Growth rates (per day) of both MCF-7 and MDA-MB-231 cells as measured by the instantaneous exponential growth rate for each monoculture in the first 72 h under all experimental conditions. ANOVA analysis is provided in Supplementary Tables 1–3.

Full size image

While the growth rates of MDA-MB-231 cells during these early time points were generally unaffected by changes in environmental conditions, the initial plating density had a significant effect. The higher (20,000 initial cells) plating density (Fig. 5a,d) resulted in slower, even negative, growth rates when compared to the lower (10,000 initial cells) plating density. This suggests that these plating densities may be approaching MDA-MB-231’s carrying capacities.
Analysis of carrying capacity
Due to small differences in the RFU intensity between the GFP and RFP used for the MCF-7 and MDA-MB-231 cells respectively, a direct comparison of carrying capacity using these RFU measurements cannot be made. In this way we split the results in to two figures (Figs. 6 and 7). We see RFU’s higher than 60 in the MCF-7 data while the maximum RFU for MDA-MB-231 data is never greater than 20 (please note the varying y-axis limits in the corresponding figures). Even with this difference in RFU intensity, it is not unreasonable to assume, due to the large difference, the MCF-7 cells do indeed have a higher carrying capacity across all experimental conditions.
Figure 6

Carrying capacities of MCF-7 cells under all experimental conditions, as measured by the mean fluorescence of the final ten time points. ANOVA analysis is provided in Supplementary Tables 4, 5.

Full size image

Figure 7

Carrying capacities of MDA-MB-231 cells under all experimental conditions, as measured by the mean fluorescence of the final ten time points. ANOVA analysis is provided in Supplementary Tables 4, 5.

Full size image

For MCF-7 cells, an increase in glucose resulted in an increase in carrying capacity across experimental conditions, as would be expected (Fig. 6a–f). Surprisingly, the MCF-7 cells experience a decrease in carrying capacity under normal pH conditions when compared to acidic pH 6.5 conditions suggesting an advantage in acidic conditions.
For the MDA-MB-231 cells, after a high seeding density of 20,000 cells, an increase in glucose resulted in an increase in carrying capacity, just like the MCF-7 cells (Fig. 7a,d). For the low seeding density of 10,000 cells, the opposite occurred (Fig. 7b,c,e, and f), though not as dramatically. Furthermore, the MDA-MB-231 cells experience an increase in carrying capacity in normal pH when compared to acidic pH 6.5 conditions suggesting an advantage in normal pH conditions.
Analysis of initial seeding frequencies
For the analysis of frequency dependent interactions of the two cells lines we analyzed the instantaneous growth rates (in units of per day) during the first 72 h in the mixed spheroids (Figs. 8 and 9). The figures are split between low glucose conditions of 0 g/L (Fig. 8) and high glucose conditions of 4.5 g/L (Fig. 9). Interestingly, under both glucose environments, MDA-MB-231 cells have their maximum growth rates when there are high initial frequencies of MCF-7 cells in the spheroid. This suggests some benefit to the MDA-MB-231 cells from the presence of MCF-7 cells, or less inhibition of growth rates from MCF-7 cells. Furthermore, we see that the MCF-7 cells generally show a maximum growth rate when the seeding frequencies are 40%/60% or 60%/40% (Figs. 8 and 9b,c,e, and f), again suggesting possible benefits to having both cell types in the spheroid.
Figure 8

Based on the first 72 hours, estimates of growth rates (per day) for mixed population tumor spheroids under glucose starved conditions of 0 g/L. ANOVA analysis is provided in Supplementary Tables 1–3.

Full size image

Figure 9

Based on the first 72 hours, estimates of growth rates (per day) for mixed population tumor spheroids under glucose rich conditions of 4.5 g/L. ANOVA analysis is provided in Supplementary Tables 1–3.

Full size image

In low glucose and at the high seeding density of 20,000 cells, the MCF-7 cells experience a significant decrease in growth rate, even negative growth rates, as the seeding frequency of MDA-MB-231 cells increase. This points to the MDA-MB-231 cells potentially approaching the carrying capacity provided within these nutrient starved conditions, essentially fully depleting all the resources for other cell types. When the glucose is increased, these MCF-7 growth rates increase, suggesting there are now sufficient nutrients for both cell types (Figs. 8 and 9a,d). Upon seeding it appears that both cell types generally have higher initial growth rates when the starting frequencies of MDA-MB-231 are low.
The logistic model best describes individual cell growth
Frequency-dependent effects (here measured as competition coefficients) can only be quantified by modeling the growth of multiple interacting types or species. In order to obtain an accurate estimate of these, the mode of growth must first be determined. In this way, we use the monoculture spheroids of either MCF-7 or MDA-MB-231 alone as a data set to first understand the growth dynamics without the added complexity of the mixed cultured spheroids. Although several growth models have been proposed for spheroids, tumors, and organoids, no specific consensus has been reached34,35. Because of this, we agnostically compared four population growth models: exponential, logistic, Gompertz, and Monod-like (Table 1).
Table 1 The four models considered for estimating the mode of growth of the cancer cells based on data from monoculture spheroids seeded at a low density of 10,000 cells.
Full size table

The exponential growth model is commonly used in bacterial studies and assumes unlimited resources. Logistic population growth is the simplest constrained growth model and assumes that per capita population growth rate declines linearly with population size or density36. Although it gives a similar shape, Gompertz growth is exponentially constrained by increasing population. Several authors and works have suggested that the Gompertz equation provides a better fit to the growth of tumors in vivo37. The Monod-like equation imagines resource matching where per capita growth rates are proportional to the amount of resources supplied per individual. It provides a good fit to the population growth curves of bacteria and other single cell organisms grown in chemostats38. Here, we consider all four models when estimating the DI parameter of intrinsic growth rate (r) and the latter three models for estimating the DD parameter of carrying capacity (K).
Using constrained parameter optimization, all four models were fit to the mono-culture spheroid growth of MCF-7 or MDA-MB-231 cells under the eight environmental conditions (Supplementary Fig. 1). Due to the lower number of time points in the experiments with a high seeding density of 20,000 cells, these experiments were not used. Quality of fit for each model was evaluated using adjusted R2, Root Mean Square Error (RMSE), and Akaike Information Criterion (AIC) (Table 2, Supplementary Table 6).
Table 2 Average adjusted R2, RMSE, and AIC when fitting the four growth models to monoculture spheroids.
Full size table

Two-way ANOVAs (one for each cell type with pH and nutrient levels as independent variables) showed that adjusted R2’s were similar for the three density-dependent models and significantly lower for the exponential model (Supplementary Table 6). Conversely, the RMSE indicated the best fit for exponential growth (lowest value), lowest quality fit for the Monod-like equation, and similar, intermediate values for the logistic and Gompertz growth models. AIC analysis suggests logistic fits to be the best fit in the majority of scenarios.
While any of these four models could be used to model this system, this analysis suggests either an exponential or logistic model provides the best fit to the data. The principal difference between these two models is that logistic growth includes density-dependence and the growth curve approaches a carrying capacity, whereas exponential growth ultimately leads to extinction (if negative growth rates) or infinitely large population sizes (if positive growth rates). In this way, we evaluated whether the spheroids showed evidence of approaching or reaching a carrying capacity. To do this, we examined the slope of the final ten time points for both the red and green fluorescent values (RFU’s) for all experimental conditions (Supplementary Table 7). As most mono-culture spheroids showed declining slopes, or slopes near to zero, we favor the logistic model for this study.
Parameter estimation for Lotka–Volterra competition model
With the logistic equation showing the best fit for the mono-culture experimental data, we expanded to the Lotka–Volterra competition equations to model the mixed population tumor spheroids. The Lotka–Volterra equations represent a two species extension of logistic growth with the addition of competition coefficients.

$$begin{aligned} frac{{dN_{1} }}{dt} & = N_{1} r_{1} left( {1 – frac{{N_{1} + alpha N_{2} }}{{K_{1} }}} right) \ frac{{dN_{2} }}{dt} & = N_{2} r_{2} left( {1 – frac{{N_{2} + beta N_{1} }}{{K_{2} }}} right) \ end{aligned}$$

where Ni are population sizes, ri are intrinsic growth rates and Ki are carrying capacities. We let species 1 and 2 be MCF-7 and MDA-MB-231, respectively. The competition coefficients (α and β) scale the effects of inter-cell-type competition where α is the effect of MDA-MB-231 on the growth rate of MCF-7 in units of MCF-7, and vice-versa for β. If α (or β) is close to zero, then there is no inter-cell-type suppression; if α (or β) is close to 1 then intra-cell-type competition and inter-cell-type competition are comparable; if α (or β) is > 1 then inter-cell-type competition is severe, and if α (or β) is  (K_{MDA})/β and (K_{MCF})/α  > (K_{MDA}). Similarly, MDA-MB-231 cells are expected to outcompete MCF-7 cells if (K_{MDA})/β  > (K_{MCF}) and (K_{MDA})  > (K_{MCF})/α. Either MCF-7 or MDA-MB-231 will out-compete the other depending upon initial conditions (alternate stable states) when (K_{MCF})  > (K_{MDA})/β and (K_{MDA})  > (K_{MCF})/α. The stable coexistence of MCF-7 and MDA-MB-231 cells is expected when (K_{MDA})/β  > (K_{MCF}) and (K_{MCF})/α  > (K_{MDA}).
To estimate the values for (K_{MCF}), (K_{MDA} ,{ }) α, and β, we fit the co-culture data with low seeding density to the two-species Lotka–Volterra equation creating different fits for each combination of experimental conditions. It is important to note here that the Lotka–Volterra competition equations are not spatially explicit. From Figs. 2, 3 and 4, it can be seen that the cells are segregated within the tumor spheroid. In our application, the Lotka–Volterra competition equations simply provide a first-order linear approximation of between cell-type competition, not a mechanistic model of the consumer-resource dynamics or the spatial dynamics that likely occur in the spheroids.
Using nonlinear constrained optimization, we varied the values for instantaneous growth rates, carrying capacities, and competition coefficients α and β to minimize the RMSE between the Lotka–Volterra model and the experimental data (Fig. 10). In addition to estimating the competition coefficients, this uses the mixed culture data to re-estimate the growth rates and carrying capacities previously estimated from the monoculture spheroids as growth rate and carrying capacity depend on initial seeding density (Figs. 8 and 9).
Figure 10

(a) Estimates for growth rates, carrying capacities, and competition coefficients using nonlinear constrained optimization fitting to the Lotka–Volterra competition growth curves using the low (10,000 cell) seeding density experiments. The average of the four replicates for each experimental condition are shown along with the optimized fit to the Lotka–Volterra model. (b) Optimized parameters for growth rates (per day), carrying capacities, and competition coefficients are shown for each experimental condition. The predicted outcome of competition is given in the last column.

Full size image

As in the monoculture parameter estimation, MCF-7 had higher intrinsic growth rates and carrying capacities than MDA-MB-231 under all culture conditions. For all culture conditions, α  > 1 (as high as 12.6) and β ≈ 0 (the magnitude of six of the eight β’s was less than 0.1). Under two conditions (normal Ph and glutamine), the β’s are less than zero allowing for the possibility that MCF-7 cells are actually facilitating the per capita growth rates of the MDA-MB-231. Thus, across all treatments MDA-MB-231 cells had large competitive effects on MCF-7 growth, while the MCF-7 cells had virtually no inhibitory effects on the growth of MDA-MB-231 cells. It appears that this frequency dependent (FD) effect provides a competitive advantage for the MDA-MB-231 cells that allows them to remain in the tumor spheroids despite MCF-7’s higher intrinsic growth rates (DI) and carrying capacities (DD).
The greatest frequency dependent effect of MDA-MB-231 on MCF-7 (highest α) occurred under conditions of 4.5 g/L glucose, 6.5pH, and no glutamine. The greatest effect of MCF-7 cells on MDA-MB-231 (highest β) occurred in the absence of glucose and glutamine at neutral pH. This fits our original predictions; however, the effect of β is not strong enough to rescue MCF-7 cells from competitive exclusion. In accordance with the hypothesis that glucose permits MDA-MB-231 cells to be highly competitive, α’s (the effect of MDA-MB-231 on MCF-7) were generally lower under glucose-starved conditions. Further, both cell types appear to experience decreased growth efficiency and competition when glucose/glutamine starved, and under acidic conditions. This is consistent with more general ecological studies showing that competition between species is generally less under harsh physical conditions39,40.
Using the optimized parameters values calculated by the constrained optimization analysis of the Lotka–Volterra equations can predict the outcome of competition. This analysis suggests that under physiologic pH, regardless of glucose or glutamine concentration, MDA-MB-231 cells will outcompete the MCF-7 cells. Furthermore, under acidic pH, regardless of glucose or glutamine concentration, MDA-MB-231 cells and MCF-7 cells will actually coexist. This is observed experimentally in Figs. 2 and 3 where acidic pH results in coexistence of the two cell types and normal pH results in no MCF-7 cells remaining at the end of the experiment. While a robust result, this does not accord with our original expectations.
In vivo exploration of cell competition
MDA-MB-231 and MCF-7 cells were grown in the mammary fat pads of 8–10 week old nu/nu female mice. Three tumor models were evaluated: only MCF-7, only MDA-MB-231, and 1:1 mix of both cell types with three replicates for each (total of nine mice). Tumor volumes were measured weekly. At 5 weeks, the tumors were harvested and evaluated for ER expression which marks MCF-7 (ER positive) cells providing estimates of the ratio of MDA-MB-231 to MCF-7 cells.
H&E staining indicates that MCF-7 tumors had the greatest viability and lowest necrosis while MDA-MB-231 tumors displayed decreased viability and increased necrosis (Fig. 11a). The small differences in necrotic and viable tissue between the MDA-MB-231 and mixed tumors are indicative of the expansive effect that MDA-MB-231 phenotype has on tumor progression (Fig. 11c).
Figure 11

Results of in vivo mono- and co-culture tumors. (a) Histology slides of MCF-7, MDA-MB-231, and 1:1 co-cultured tumors. Slides were stained with H&E and ER immune stain. (b) Changes in the standardized tumor volume with time (3 mice per treatment). (c) Tumor viability for all three tumor types. Notice a decrease in viability and increase in necrosis in the MDA-MB-231 and mixed tumors. (d) Measures of CAIX, Glut1, and ER biomarker staining for each tumor type. (e) The presence of small clumps of ER staining in mixed tumors reflect the clumping of MCF-7 cells in the progression of the mixed tumor spheroids (these are enlarged photos of the indicated portions of the co-cultured tumors shown in a), suggesting the coexistence of these two phenotypes in certain tumor microenvironments. All error bars show standard error of the mean.

Full size image

Changes in tumor size over 5 weeks indicated that MCF-7 cell tumors grew more slowly than MDA-MB-231 tumors and mixed tumors (Fig. 11b; ANCOVA of natural logarithm of standardized size with day as a covariate, F1,49 = 248.0 p  More

1.
Liste, H.-H. & Felgentreu, D. Crop growth, culturable bacteria, and degradation of petrol hydrocarbons (PHCs) in a long-term contaminated field soil. Appl. Soil. Ecol. 31, 43–52 (2006).
Article  Google Scholar 
2.
Smith, M. J., Flowers, T. H., Duncan, H. J. & Alder, J. Effects of polycyclic aromatic hydrocarbons on germination and subsequent growth of grasses and legumes in freshly contaminated soil and soil with aged PAHs residues. Environ. Pollut. 141, 519–525 (2006).
CAS  PubMed  Article  Google Scholar 

3.
Meudec, A., Poupart, N., Dussauze, J. & Deslandes, E. Relationship between heavy fuel oil phytotoxicity and polycyclic aromatic hydrocarbon contamination in Salicornia fragilis. Sci. Total Environ. 381, 146–156 (2007).
ADS  CAS  PubMed  Article  Google Scholar 

4.
Euliss, K., Ho, C.-H., Schwab, A. P., Rock, S. & Banks, M. K. Greenhouse and field assessment of phytoremediation for petroleum contaminants in a riparian zone. Bioresour. Technol. 99, 1961–1971 (2008).
CAS  PubMed  Article  Google Scholar 

5.
Hutchinson, T. C. & Freedman, W. Effects of experimental crude oil spills on subarctic boreal forest vegetation near Norman Wells, N.W.T., Canada. Can. J. Bot. 56, 2424–2433 (1978).
Article  Google Scholar 

6.
Lin, Q. & Mendelssohn, I. A. The combined effects of phytoremediation and biostimulation in enhancing habitat restoration and oil degradation of petroleum contaminated wetlands. Ecol. Eng. 10, 263–274 (1998).
Article  Google Scholar 

7.
Racine, C. H. Long-term recovery of vegetation on two experimental crude oil spills in interior Alaska black spruce taiga. Can. J. Bot. 72, 1171–1177 (1994).
Article  Google Scholar 

8.
Fatima, K., Afzal, M., Imran, A. & Khan, Q. M. Bacterial rhizosphere and endosphere populations associated with grasses and trees to be used for phytoremediation of crude oil contaminated soil. Bull. Environ. Contam. Toxicol. 94, 314–320 (2015).
CAS  PubMed  Article  Google Scholar 

9.
Hashmat, A. J. et al. Characterization of hydrocarbon-degrading bacteria in constructed wetland microcosms used to treat crude oil polluted water. Bull. Environ. Contam. Toxicol. 102, 358–364 (2019).
CAS  PubMed  Article  Google Scholar 

10.
Khan, F. I., Husain, T. & Hejazi, R. An overview and analysis of site remediation technologies. J. Environ. Manag. 71, 95–122 (2004).
Article  Google Scholar 

11.
Sarkar, D., Ferguson, M., Datta, R. & Birnbaum, S. Bioremediation of petroleum hydrocarbons in contaminated soils: comparison of biosolids addition, carbon supplementation, and monitored natural attenuation. Environ. Pollut. 136, 187–195 (2005).
CAS  PubMed  Article  Google Scholar 

12.
Gan, S., Lau, E. V. & Ng, H. K. Remediation of soils contaminated with polycyclic aromatic hydrocarbons (PAHs). J. Hazard. Mater. 172, 532–549 (2009).
CAS  PubMed  Article  Google Scholar 

13.
Anchugova, E. M., Melekhina, E. N., Markarova, MYu. & Shchemelinina, T. N. Approaches to the assessment of the efficiency of remediation of oil-polluted soils. Eurasian Soil Sci. 49, 234–237 (2016).
ADS  CAS  Article  Google Scholar 

14.
Erkenova, M. I., Tolpeshta, I. I., Trofimov, S. Y., Aptikaev, R. S. & Lazarev, A. S. Changes of the content of oil products in the oil-polluted peat soil of a high-moor bog in a field experiment with application of lime and fertilizers. Eurasian Soil Sci. 49, 1310–1318 (2016).
ADS  CAS  Article  Google Scholar 

15.
Sorkhoh, N. A. et al. Bioremediation of volatile oil hydrocarbons by epiphytic bacteria associated with American grass (Cynodon sp.) and broad bean (Vicia faba) leaves. Int. Biodeterior. Biodegrad. 65, 797–802 (2011).
CAS  Article  Google Scholar 

16.
Roy, A. S. et al. Bioremediation potential of native hydrocarbon degrading bacterial strains in crude oil contaminated soil under microcosm study. Int. Biodeterior. Biodegrad. 94, 79–89 (2014).
CAS  Article  Google Scholar 

17.
Cai, B. et al. Comparison of phytoremediation, bioaugmentation and natural attenuation for remediating saline soil contaminated by heavy crude oil. Biochem. Eng. J. 112, 170–177 (2016).
CAS  Article  Google Scholar 

18.
Murygina, V., Gaydamaka, S., Gladchenko, M. & Zubaydullin, A. Method of aerobic-anaerobic bioremediation of a raised bog in Western Siberia affected by old oil pollution. A pilot test. Int. Biodeterior. Biodegrad. 114, 150–156 (2016).
CAS  Article  Google Scholar 

19.
Tahseen, R. et al. Rhamnolipids and nutrients boost remediation of crude oil-contaminated soil by enhancing bacterial colonization and metabolic activities. Int. Biodeterior. Biodegrad. 115, 192–198 (2016).
CAS  Article  Google Scholar 

20.
Baoune, H. et al. Effectiveness of the Zea mays-Streptomyces association for the phytoremediation of petroleum hydrocarbons impacted soils. Ecotoxicol. Environ. Saf. 184, 109591 (2019).
CAS  PubMed  Article  Google Scholar 

21.
Ra, T., Zhao, Y. & Zheng, M. Comparative study on the petroleum crude oil degradation potential of microbes from petroleum-contaminated soil and non-contaminated soil. Int. J. Environ. Sci. Technol. 16, 7127–7136 (2019).
CAS  Article  Google Scholar 

22.
Rajkumari, J., Bhuyan, B., Das, N. & Pandey, P. Environmental applications of microbial extremophiles in the degradation of petroleum hydrocarbons in extreme environments. Environ. Sustain. 2, 311–328 (2019).
CAS  Article  Google Scholar 

23.
Newman, L. A. & Reynolds, C. M. Phytodegradation of organic compounds. Curr. Opin. Biotechnol. 15, 225–230 (2004).
CAS  PubMed  Article  Google Scholar 

24.
Unterbrunner, R. et al. Plant and fertiliser effects on rhizodegradation of crude oil in two soils with different nutrient status. Plant Soil 300, 117–126 (2007).
CAS  Article  Google Scholar 

25.
Muratova, A. Y., Dmitrieva, T. V., Panchenko, L. V. & Turkovskaya, O. V. Phytoremediation of oil-sludge-contaminated soil. Int. J. Phytorem. 10, 486–502 (2008).
CAS  Article  Google Scholar 

26.
Shirdam, R., Zand, A., Bidhendi, G. & Mehrdadi, N. Phytoremediation of hydrocarbon-contaminated soils with emphasis on the effect of petroleum hydrocarbons on the growth of plant species. Phyto 89, 21–29 (2008).
CAS  Article  Google Scholar 

27.
Farias, V. et al. Phytodegradation Potential of Erythrina crista-galli L., Fabaceae petroleum-contaminated soil. Appl. Biochem. Biotechnol. 157, 10–22 (2009).
PubMed  Article  CAS  Google Scholar 

28.
Peng, S., Zhou, Q., Cai, Z. & Zhang, Z. Phytoremediation of petroleum contaminated soils by Mirabilis jalapa L. in a greenhouse plot experiment. J. Hazard. Mater. 168, 1490–1496 (2009).
CAS  PubMed  Article  Google Scholar 

29.
Basumatary, B., Saikia, R., Bordoloi, S., Das, H. C. & Sarma, H. P. Assessment of potential plant species for phytoremediation of hydrocarbon-contaminated areas of upper Assam, India. J. Chem. Technol. Biotechnol. 87, 1329–1334 (2012).
CAS  Article  Google Scholar 

30.
Moubasher, H. A. et al. Phytoremediation of soils polluted with crude petroleum oil using Bassia scoparia and its associated rhizosphere microorganisms. Int. Biodeterior. Biodegrad. 98, 113–120 (2015).
CAS  Article  Google Scholar 

31.
Fatima, K., Imran, A., Amin, I., Khan, Q. M. & Afzal, M. Plant species affect colonization patterns and metabolic activity of associated endophytes during phytoremediation of crude oil-contaminated soil. Environ. Sci. Pollut. Res. Int. 23, 6188–6196 (2016).
CAS  PubMed  Article  Google Scholar 

32.
Khan, S., Afzal, M., Iqbal, S. & Khan, Q. M. Plant–bacteria partnerships for the remediation of hydrocarbon contaminated soils. Chemosphere 90, 1317–1332 (2013).
ADS  CAS  PubMed  Article  Google Scholar 

33.
Yavari, S., Malakahmad, A. & Sapari, N. B. A review on phytoremediation of crude oil spills. Water Air Soil Pollut. 226, 279 (2015).
ADS  Article  CAS  Google Scholar 

34.
Okoh, E., Yelebe, Z. R., Oruabena, B., Nelson, E. S. & Indiamaowei, O. P. Clean-up of crude oil-contaminated soils: bioremediation option. Int. J. Environ. Sci. Technol. 17, 1185–1198 (2020).
CAS  Article  Google Scholar 

35.
Naeem, U. & Qazi, M. A. Leading edges in bioremediation technologies for removal of petroleum hydrocarbons. Environ. Sci. Pollut. Res. 27, 27370–27382 (2019).

36.
Tyagi, M., da Fonseca, M. M. R. & de Carvalho, C. C. C. R. Bioaugmentation and biostimulation strategies to improve the effectiveness of bioremediation processes. Biodegradation 22, 231–241 (2011).
CAS  PubMed  Article  Google Scholar 

37.
Fatima, K., Imran, A., Amin, I., Khan, Q. M. & Afzal, M. Successful phytoremediation of crude-oil contaminated soil at an oil exploration and production company by plants-bacterial synergism. Int. J. Phytoremediat. 20, 675–681 (2018).
CAS  Article  Google Scholar 

38.
Walker, D. A. et al. Cumulative impacts of oil fields on Northern Alaskan Landscapes. Science 238, 757–761 (1987).
ADS  CAS  PubMed  Article  Google Scholar 

39.
Maganov, R. U., Markarova, M. Y., Mulyak, V. V., Zagvozdkin, V. E. & Zaikin, I. A. Nature conservation management at the oil and gas companies. Part 1. Reclamation of oil-polluted lands in the Usinsky district of the Komi Republic (Komi Science Center Ural Branch of RAS, Syktyvkar, 2006) ((in Russian)).
Google Scholar 

40.
Vervaeke, P. et al. Phytoremediation prospects of willow stands on contaminated sediment: a field trial. Environ. Pollut. 126, 275–282 (2003).
CAS  PubMed  Article  Google Scholar 

41.
Huang, X.-D., El-Alawi, Y., Gurska, J., Glick, B. R. & Greenberg, B. M. A multi-process phytoremediation system for decontamination of persistent total petroleum hydrocarbons (TPHs) from soils. Microchem. J. 81, 139–147 (2005).
CAS  Article  Google Scholar 

42.
Robson, D. B., Knight, J. D., Farrell, R. E. & Germida, J. J. Natural revegetation of hydrocarbon-contaminated soil in semi-arid grasslands. Can. J. Bot. 82, 22–30 (2004).
Article  Google Scholar 

43.
Grime, J. P. & Pierce, S. The evolutionary strategies that shape ecosystems (Wiley-Blackwell, Chichester, 2012).
Google Scholar 

44.
Ramenskiy, L. G. On principal rules, basic concepts, and terms of land typology, geobotany, and ecology. Sov. Bot. 4, 25–42 (1935).
Google Scholar 

45.
Grime, J. P., Hodgson, J. G. & Hunt, R. Comparative plant ecology: a functional approach to common British species (Unwin Hyman, London, 1988).
Google Scholar 

46.
Thompson, K. Predicting the fate of temperate species in response to human disturbance and global change. In Biodiversity, temperate ecosystems, and global change (eds Boyle, T. J. B. & Boyle, C. E. B.) 61–76 (springer, Berlin, 1994).
Google Scholar 

47.
Massant, W., Godefroid, S. & Koedam, N. Clustering of plant life strategies on meso-scale. Plant Ecol. 205, 47–56 (2009).
Article  Google Scholar 

48.
Novakovsky, A. B. & Panyukov, A. N. Analysis of successional dynamics of a sown meadow using Ramenskii–Grime’s System of ecological strategies. Russ. J. Ecol. 49, 119–127 (2018).
Article  Google Scholar 

49.
Novakovskiy, A. B. & Elsakov, V. V. Hydrometeorological database (HMDB) for practical research in ecology. Data Sci. J. 13, 57–63 (2014).
Article  Google Scholar 

50.
PND F 16.1: 2.21-98. Quantitative chemical analysis of soils. The method of measuring the mass fraction of oil products in soil and soil samples by the fluorimetric method on the Fluorat-02 fluid analyzer. https://www.russiangost.com/p-275219-fr131201213170.aspx

51.
Archegova, I. B., Markarova, M. Y. & Gromova, O. V. Method for reclaiming posttechnogenic lands and lands in remote districts of extreme North. US Patent RU2093974C1 (1997).

52.
Archegova, I. B., Markarova, M. Y. & Gromova, O. V. Method for producing granular fertilizing-seeding material. US Patent RU2099917C1 (1997).

53.
Murygina, V. P., Vojshvillo, N. E. & Kaljuzhnyj, S. V. Biological preparation ‘Roder’ for cleaning soils, soil grounds, sweet and mineralized waters to remove crude oil and petroleum products. US Patent RU2174496C2 (2001).

54.
Murygina, V. P., Markarova, M. Y. & Kalyuzhnyi, S. V. Application of biopreparation “Rhoder” for remediation of oil polluted polar marshy wetlands in Komi Republic. Environ. Int. 31, 163–166 (2005).
CAS  PubMed  Article  Google Scholar 

55.
Lavorel, S. et al. Assessing functional diversity in the field—methodology matters!. Funct. Ecol. 22, 134–147 (2008).
Google Scholar 

56.
Kindt, R. & Coe, R. Tree diversity analysis: a manual and software for common statistical methods for ecological and biodiversity studies (World Agroforestry Centre, Nairobi, 2006).

57.
Hodgson, J. G., Wilson, P. J., Hunt, R., Grime, J. P. & Thompson, K. Allocating C-S-R plant functional types: a soft approach to a hard problem. Oikos 85, 282–294 (1999).
Article  Google Scholar 

58.
Pierce, S., Brusa, G., Vagge, I. & Cerabolini, B. E. L. Allocating CSR plant functional types: the use of leaf economics and size traits to classify woody and herbaceous vascular plants. Funct. Ecol. 27, 1002–1010 (2013).
Article  Google Scholar 

59.
Magguran, A. E. Measuring biological diversity (Blackwell Publishing, Oxford, 2004).
Google Scholar 

60.
Peet, R. K. The measurement of species diversity. Annu. Rev. Ecol. Syst. 5, 285–307 (1974).
Article  Google Scholar 

61.
Kruskal, J. B. Nonmetric multidimensional scaling: a numerical method. Psychometrika 29, 115–129 (1964).
MathSciNet  Article  MATH  Google Scholar 

62.
McCune, B. & Grace, J. B. Analysis of ecological communities (MjM Software Design, Gleneden Beach, 2002).
Google Scholar 

63.
Melekhina, E. N., Markarova, MYu., Shchemelinina, T. N., Anchugova, E. M. & Kanev, V. A. Secondary successions of biota in oil-polluted peat soil upon different biological remediation methods. Eurasian Soil Sci. 48, 643–653 (2015).
ADS  Article  Google Scholar 

64.
Borowik, A., Wyszkowska, J., Gałązka, A. & Kucharski, J. Role of Festuca rubra and Festuca arundinacea in determinig the functional and genetic diversity of microorganisms and of the enzymatic activity in the soil polluted with diesel oil. Environ. Sci. Pollut. Res. Int. 26, 27738–27751 (2019).
CAS  PubMed  PubMed Central  Article  Google Scholar 

65.
Wyszkowska, J., Borowik, A. & Kucharski, J. The resistance of Lolium perenne L. × hybridum, Poa pratensis, Festuca rubra, F. arundinacea, Phleum pratense and Dactylis glomerata to soil pollution by diesel oil and petroleum. Plant Soil Environ. 65, 307–312 (2019).
CAS  Article  Google Scholar 

66.
Freedman, W. & Hutchinson, T. Physical and biological effects of experimental crude-oil spills on Low Arctic Tundra in Vicinity of Tuktoyaktuk, Nwt, Canada. Can. J. Bot.-Rev. Can. Bot. 54, 2219–2230 (1976).
Article  Google Scholar 

67.
Kazantseva, M. N. The effect of oil extraction on ground cover of West Siberian taiga forests. Contemp. Probl. Ecol. 4, 582–587 (2011).
Article  Google Scholar 

68.
Lapshina, E. D. & Bleuten, W. Types of deterioration and self-restoration of vegetation of olygotrophic bogs in oil-production areas of Tomsk province, Krylovia. Siberian Bot. J. 1, 129–140 (1999).
Google Scholar 

69.
Cook, R. L., Landmeyer, J. E., Atkinson, B., Messier, J.-P. & Nichols, E. G. Field note: successful establishment of a phytoremediation system at a petroleum hydrocarbon contaminated shallow aquifer: trends, trials, and tribulations. Int. J. Phytorem. 12, 716–732 (2010).
Article  Google Scholar 

70.
Nichols, E. G. et al. Phytoremediation of a petroleum-hydrocarbon contaminated shallow aquifer in Elizabeth City, North Carolina, USA. Remediat. J. 24, 29–46 (2014).
Article  Google Scholar 

71.
Seburn, D. C., Kershaw, G. P. & Kershaw, L. J. Vegetation response to a subsurface crude oil spill on a subarctic right-of-way, Tulita (Fort Norman), Northwest Territories, Canada. Arctic 49, 321–327 (1996).
Article  Google Scholar 

72.
Melekhina, E. N., Markarova, M. Y. U., Anchugova, E. M., Shchemelinina, T. N. & Kanev, V. A. The efficiency assessment of oil polluted soil recultivation methods. Bull. Komi Sci. Center Ural Branch of RAS 27, 61–70 (2016) ((in Russian)).
Google Scholar 

73.
Ma, X. & Burken, J. G. VOCs fate and partitioning in vegetation: use of tree cores in groundwater analysis. Environ. Sci. Technol. 36, 4663–4668 (2002).
ADS  PubMed  Article  Google Scholar 

74.
Haroni, N. N., Badehian, Z., Zarafshar, M. & Bazot, S. The effect of oil sludge contamination on morphological and physiological characteristics of some tree species. Ecotoxicology 28, 507–519 (2019).
CAS  PubMed  Article  Google Scholar 

75.
Banks, M. K., Kulakow, P., Schwab, A. P., Chen, Z. & Rathbone, K. Degradation of crude oil in the rhizosphere of sorghum bicolor. Int. J. Phytorem. 5, 225–234 (2003).
CAS  Article  Google Scholar  More

The following museum acronyms are used (mostly following55):
BMNH, NHM, BM: Natural History Museum, London, formerly British Museum of Natural History (UK);
CAS: California Academy of Sciences, San Francisco (USA);
MHNC: Musée d’histoire naturelle, La Chaux-de-Fonds (CH);
MHNG: Muséum d’Histoire Naturelle, Genève (CH);
MNHN: Muséum national d’Histoire naturelle, Paris (F);
SMNS: Staatliches Museum für Naturkunde, Stuttgart (D);
UWBM The Washington State Museum of Natural History and Culture/Burke Museum, University of Washington (USA);
ZMB: Museum für Naturkunde Berlin, formerly Zoologisches Museum Berlin (D).
Study organism and taxonomic background
Acanthodactylus was first described as a subgenus of Lacerta Cuvier (sic)56, the type species is A. boskianus (Daudin, 1802). The few meristic characters described to be diagnostic for Acanthodactylus56 include: collar connected in the center but free on the sides, tempora squamosal, i.e., temporal region covered by small scales rather than larger shields, ventral scales rectangular and arranged in longitudinal rows, digits acutely fimbriate-denticulate forming “toe fringes”. Fringed toes have evolved in various shapes multiple times in lizards and in Lacertidae, and when triangular, projectional or conical as in Acanthodactylus they are commonly seen as adaptation to windblown sand substrate57.
The enigmatic Acanthodactylus guineensis is among the lesser-known species of the Lacertidae. Only limited information is available regarding the species’ morphology, habitat and distribution (see Meinig & Böhme44 for a review and references therein; and 34, 38, 45, 58). Based on one young specimen, A. guineensis was described as member of the genus Eremias Fitzinger, 1834 by Boulenger40. Boulenger does not mention the slightly projecting third row of scales around the toes and fingers (“fringed toes”) of the type specimen in his quite comprehensive description40, probably because he did not notice it in such a small (SVL 24 mm) specimen and because the projection is indeed rather minor compared to other members of the genus. This is probably also the reason why he did not place guineensis in the genus Acanthodactylus Wiegmann, 1834 but in Eremias Fitzinger, 1834 (in56).
About 30 years later, Boulenger revised the genus Eremias and divided it into five sections he assumed to be natural associations59: (1) Taenieremias Boulenger, 1918—monotypic, type species Eremias guineensis Boulenger, 1887 (currently Acanthodactylus guineensis); (2) Lampreremias Boulenger, 1918—type species Eremias nitida Günther, 1872 (currently Heliobolus nitidus); (3) Pseuderemias Boettger, 1883—type species Eremias mucronata (Blandford, 1870) (currently Pseuderemias mucronata); (4) Mesalina Gray, 1838—type species Eremias rubropunctata (Lichtenstein, 1823) (currently Mesalina rubropunctata); (5) Eremias s. str. Fitzinger, 1834—type Eremias velox (Pallas, 1771) (currently Eremias velox).
Monard47 described Eremias (Taenieremias) benuensis from Cameroon based on a few minor morphological differences compared to E. guineensis, but in 1969 E. benuensis was synonymized with E. guineensis60.
Salvador17 and one year later also Arnold18 in their respective major revisions of the genus Acanthodactylus both found that Eremias (Taenieremias) guineensis agrees with all the characteristic morphological features of Acanthodactylus with the exception of the arrangement of scales around the nostril. However, it was suggested that the E. guineensis condition with an extra suture across the area occupied by the first upper labial scale to produce a smaller first upper labial is easily derived from that found in Acanthodactylus, evidenced by BMNH 1966.430, a juvenile A. erythrurus lineomaculatus (sic) with a similar scale arrangement18. Within Acanthodactylus, both authors placed guineensis in the Western clade and in the Acanthodactylus erythrurus group which was assumed to consist of A. erythrurus, A. savignyi, A. boueti, A. guineensis, and A. blanci which was considered a subspecies of A. savignyi at the time14, 20.
Molecular data and phylogenetic analyses
We aimed at reconstructing a well-supported phylogeny of the genus Acanthodactylus using whole mitochondrial DNA sequences, assembled by means of a shotgun next generation sequencing strategy. Analyses of whole mitogenomes have been shown to resolve many nodes of the lacertid tree with high statistical support (e.g.31). We retrieved muscle tissue from a museum voucher of A. guineensis (ZFMK 59511 from Daroha, near Bobo Dioulasso, Burkina Faso43) that likely was never in contact with formalin for preservation and therefore offered good chances to obtain sufficient amounts of DNA of decent quality for sequencing. We extracted genomic DNA using the Qiagen DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) following the protocol provided by the manufacturer. We additionally sequenced a tissue sample of a freshly caught individual of Acanthodactylus schmidti (SB 642) from Abu Dhabi (UAE) to increase sampling of Acanthodactylus spp. in the mitogenomic tree. The lizard was euthanized by injection of an aqueous solution of benzocaine (20%) into the body cavity. Subsequently, a sample of muscle tissue was taken from the thigh, preserved in 98% Ethanol and stored in − 80 °C. Handling, euthanizing and collection of tissue samples of A. schmidti individuals was approved by the NYUAD Institutional Animal Care and Use Committee (IACUC 19-0002) and UAE No Objection Certificate (NOC 8416), and all applied methods were performed in accordance with the relevant guidelines and regulations.
Genomic DNA of SB 642 was extracted from ethanol-preserved muscle tissue using the Qiagen MagAttract HMW DNA Kit (Qiagen, Hilden, Germany) for high molecular weight DNA. We determined DNA yields on a Qubit fluorometer (Qubit, London, UK) with a dsDNA high sensitivity kit. Totals of 52 ng and 80 ng of DNA (diluted in 26 µl 10 mM Tris·Cl, 0.5 mM EDTA, pH 9.0 (AE buffer)) from A. guineensis and A. schmidti, respectively, were used for library preparation.
ZFMK 59511 libraries were prepared with NEB Ultra II FS DNA (New England Biolabs, Ipswich, MA, USA) kit as per protocol instructions with input below 100 ng. For sample SB642, linked reads were generated on a 10X Genomics Chromium following Genome reagent kits v2 instructions. Resulting libraries’ concentration, size distribution and quality were assessed on a Qubit fluorometer (Qubit, London, UK) with a dsDNA high sensitivity kit and on an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA) using a High Sensitivity DNA kit. Subsequently, libraries were normalized and pooled, and pools quantified with a KAPA Library quantification kit for Illumina platforms (Roche Sequencing, Pleasanton, CA, USA) on a ABI StepOnePlus qPCR machine (Thermo Fisher Scientific Inc., Waltham, MA, USA), then loaded on a SP flowcell and paired-end sequenced (2 × 150 bp) on an Illumina NovaSeq 6000 next generation sequencer (Illumina, San Diego, CA, USA), and a S2 flowcell for linked read library. Raw reads were deposited in the Sequence Read Archive (SRA) under BioProject ID PRJNA700414 . All mitogenome assemblies and original alignments were deposited in Figshare under  https://doi.org/10.6084/m9.figshare.13754083.v1 .
Raw FASTQ sequenced reads were first assessed for quality using FastQC v0.11.561. The reads where then passed through Trimmomatic v0.3662 (parameters ILLUMINACLIP: trimmomatic_adapter.fa:2:30:10 TRAILING:3 LEADING:3 SLIDINGWINDOW:4:15 MINLEN:36) for quality trimming and adapter sequence removal. Following the quality trimming, the reads were assessed again using FastQC. The executed workflows were performed using BioSAILs63.
For SB 642 (A. schmidti), we received 436,370,000 single-end reads (average read length 139.5 bases (b)) with a raw coverage of 29.29X and a scaffold N50 of 21.12 kilobases (kb). For the MITObim analysis, read number for SB 642 was randomly reduced to 1,000,000 reads (360,463,035 b) using the awk command. For ZFMK 59511 (A. guineensis), we received 57,936,345 paired-end reads (average read length 151 b, raw coverage ~ 10x). Subsequently, the two sets of reads were converted to interleaved format. We used the quality filtered reads to assemble the mitogenomes of A. guineensis and A. schmidti using an iterative mapping strategy in MITObim v. 1.9.164. We used the Acanthodactylus aureus mitogenome (GB accession number xxxxx; assembly method is provided in the next paragraph) as seed for both samples; this rendered an initial mapping of a conserved region from a more distantly related individual unnecessary64. We therefore applied the –quick option in MITObim and iterations were run until no additional reads could be incorporated into the assembly (14 in A. guineensis, eight in A. schmidti).
We also assembled complete or nearly complete mitogenomes for additional twelve species of Gallotiinae and Eremiadini using anchored hybrid enrichment sequence data from Garcia-Porta et al.6 (see Supplementary methods in6 for extraction protocol and sequencing methods and Supplementary Table S1 online in6) with the Podarcis muralis complete mitogenome (GB accession number NC_011607) as reference sequence. The raw data of each individual was quality filtered using Trimmomatic v0.3662 (parameter MINLEN:45) and assembled using MITObim v. 1.9.1. All multiple alignment files generated in the final MITObim iteration were imported to Geneious R11 (https://www.geneious.com) to check for assembly quality and coverage.
The resulting assemblies were annotated with MITOS65 using defaults settings with the vertebrata database as a reference. The annotated assemblies were imported into PhyloSuite66 together with existing mitogenomic sequences of Lacertidae, and Blanus cinereus (Amphisbaenia; outgroup) available in GenBank (as of July 2020). Details of all mitogenomic sequences included in this study and the corresponding GenBank accession numbers are provided in Supplementary Table S1 online. Using PhyloSuite, we exported all protein coding sequences plus the two rRNAs. The resulting files were aligned with MAFFT67 using “auto” settings. Alignments of coding sequences were refined with MACSE68 to account for open reading frame structure during the alignments. The protein coding and rRNAs alignments were finally concatenated into a single alignment delimiting partitions by marker and, for the protein-coding genes, by codons within markers. Mitogenome phylogenies were inferred using a maximum likelihood approach with IQ-TREE v. 1.5.4 software69 and Bayesian inference (BI) analysis with MrBayes 3.270.
For the maximum likelihood approach, best-fitting partitioning schemes and substitution models were selected based on the Akaike Information Criterion (AIC) using the heuristic algorithms implemented in the MFP + MERGE option (Supplementary Table S3 online). After ML inference, branch support was assessed with 1000 ultrafast bootstrap replicates. As the alignments used for phylogenetic inference contained more than one sequence for some species, these terminals in the resulting tree were collapsed a-posteriori for aesthetic reasons using FigTree v1.4 (as depicted in Fig. 1). An uncollapsed version of this tree is presented in Supplementary Fig. S1 online. For the BI analysis, the best-fitting partition/substitution model scheme, as selected by the ModelFinder algorithm in iQTree (Supplementary Table S4 online), was implemented with MrBayes 3.270. Results of two independent runs of 10 million generations, each comprising four Markov Chains (three heated and one cold), were sampled every 1000 generation. Chain mixing and stationery was assessed by examining the standard deviation of split frequencies and by plotting the –lnL per generation using Tracer 1.5 software71. Samples corresponding to the initial phase of the Markov chain (25%) were discarded as burn-in and the remaining results were combined to obtain a majority rule consensus tree and the respective posterior probabilities of nodes.
An additional, more comprehensive, alignment was produced containing all currently known Acanthodactylus species to obtain a higher resolution perspective on the phylogenetic placement of the target species. Since no nuclear data was available for A. guineensis, we extracted from the newly obtained sequences the four best represented mitochondrial gene fragments across Lacertidae, as compiled by Garcia-Porta et al.6, despite the known shortcomings of partial mitochondrial datasets to resolve lacertid relationships9, 72, 73. The corresponding sequences of the ribosomal RNAs (12S, 16S), COB, and NADH-dehydrogenase subunit (ND4) were later added to the alignment of6 using –add and –keeplength options in MAFFT (Supplementary Table S5 online).
Species distribution modeling
As basis for characterizing the species’ bioclimatic envelopes, we compiled an updated distribution map for A. guineensis and A. boueti including new discoveries in museum collections from the current study (see Supplementary Material online), all known records from the literature16, 34, 37, 38, 40, 42, 44, 45, 47, 58, 60, 74,75,76,77 and from the Global Biodiversity Information Facility (GBIF78). We provide all coordinates as latitude (decimal degrees), longitude (decimal degrees). Due to the general scarcity of records for A. guineensis and A. boueti, we added two localities for A. guineensis (12.5°, − 2.5°; 7.5°, 13.5°) and one locality for A. boueti (− 2.5°, 7.5°) published in Trape et al.45 that were not found on GBIF or in other publications, by extracting the center coordinates of the respective 1 × 1 degree grid. These additional localities have coordinates with very low accuracy and the corresponding environmental variables extracted from a 1 × 1 km resolution grid (30 arc-seconds; see below) therefore have high degrees of uncertainty. When examined, the majority of values fell within the total range for the respective environmental variable (the two exceptions are mentioned in the Results section), we consequently decided to keep them for our analysis. We reduced the final locality dataset to one record per km2 (the resolution of the environmental input data, see below) to avoid pseudoreplication using the R package spThin79. We then tested whether the remaining presence points are randomly dispersed or clustered using the Nearest Neighbor Index NNI in the R packages sp80 and spatialEco81. NNI was 1.4 for A. guineensis and 1.2 for A. boueti, indicating that the filtered point dataset consists of randomly distributed points which are not spatially autocorrelated.
As environmental parameters we downloaded spatial layers of the Terrestrial Ecoregions of the World (TEOW)82, global bioclimate and elevation layers39, data on potential evapotranspiration83 and global aridity index84 with a spatial resolution of 30 arc-seconds (~ 1 km2 near the equator). Although the climate datasets are interpolated from data from weather stations much farther apart from one another and are consequently not measurements of actual environmental conditions, they were rigorously cross-validated with observed data (including satellite data) during the development of the dataset and generally showed high correspondence with observations (especially temperature variables)39. Since for many applications such as our species distribution models data at high spatial resolution (i.e., 1 × 1 km) are preferable over lower resolution to capture variation for example across steep climate gradients in mountains39 we chose the 30 arc-seconds dataset over the 5 arc-minutes dataset (~ 9 × 9 km near the equator) for the ecoregions and climate datasets. In addition, we downloaded spatial layers with forest and grassland/scrub/woodland cover from the harmonized soil database (only available in 5 arc-minutes spatial resolution85). We downscaled the forest and grassland/scrub/woodland dataset to 30 arc-seconds resolution despite the resulting slight inaccuracy, which we kept in mind during interpretation but did not consider relevant for the vegetation datasets.
In order to compare the prevailing climate within the distribution range of A. guineensis and A. boueti to all other species of the genus Acanthodactylus we downloaded all available records of the other currently known Acanthodactylus species from GBIF. We cleaned the downloaded datasheet by deleting all entries with missing data for either latitude, longitude, species epithet, or all of those, and removed duplicates or spelling mistakes. We further removed all GBIF entries with imprecise coordinates, i.e. records with coordinates that when plotted landed just outside of coastal areas in the ocean instead of on land (usually coordinates with only two decimals). In addition, we added several curated locality data from Garcia-Porta et al.6 as available from Figshare (https://doi.org/10.6084/m9.figshare.8866271.v1/). Following Tamar et al.14 we treated A. lineomaculatus as a junior synonym of A. erythrurus and changed the respective records in our database accordingly. For species that did not have records published on GBIF we added further locality data from the literature (Supplementary Table S6 online). The final database contained 4286 records from all 44 currently recognized species of Acanthodactylus. We acknowledge that our database is by no means complete, however, we wanted to follow the most conservative approach and use only records that are supported by a voucher specimen and can be traced back using published databases. We trust that our record list covers a nearly complete representation of the ecological conditions inhabited by all currently accepted Acanthodactylus species.
Using the bioclim dataset39 we plotted annual mean temperature (bio1) and annual precipitation (bio12) for all locality records in our Acanthodactylus species database and compared the respective data for A. guineensis and A. boueti with its congeners. We developed species distribution models using Maxent 3.4.146 for A. guineensis and A. boueti. Maxent applies the maximum entropy principle86 for model fitting under the basic premise that the estimated species distribution deviates from a uniform distribution as minimally as required to explain the observations46. We used the cleaned datasets of A. guineensis and A. boueti with only one record per km2 as input presence locality data. We extracted environmental data for all presence sites from the 19 bioclim variables (bio1-19), elevation (alt), aridity index (AI), percentage cover of forest (forest) and grassland, scrub and woodland (grass) per grid cell, as well as four parameters comprising potential evapotranspiration (PET): PET of the wettest, driest, warmest and coldest quarter of the year (PETwet, PETdry, PETwarm, PETcold). To avoid issues resulting from correlated parameters we reduced the number of environmental predictors by performing multicollinearity analyses using the Variance Inflation Factor (VIF) with a threshold  More

1.
Bakker, E. S. & Hilt, S. Impact of water-level fluctuations on cyanobacterial blooms: options for management. Aquat. Ecol. 50, 485–498 (2016).
CAS  Article  Google Scholar 
2.
Li, Q. et al. Impact of water level fluctuations on the development of phytoplankton in a large subtropical reservoir: implications for the management of cyanobacteria. Environ. Sci. Pollut. R. 25, 1306–1318 (2017).
Article  CAS  Google Scholar 

3.
Sellheim, K. L., Vaghti, M. & Merz, J. E. Vegetation recruitment in an enhanced floodplain: ancillary benefits of salmonid habitat enhancement. Limnol. Ecol. Manag. Inland Waters 58, 94–102 (2016).
Article  Google Scholar 

4.
Yuan, S., Yang, Z., Liu, X. & Wang, H. Key parameters of water level fluctuations determining the distribution of Carex in shallow lakes. Wetlands 37, 1005–1014 (2017).
Article  Google Scholar 

5.
Tharme, R. E. A global perspective on environmental flow assessment: emerging trends in the development and application of environmental flow methodologies for rivers. River Res. Appl. 19, 397–441 (2003).
Article  Google Scholar 

6.
Wang, M., Liu, Z., Luo, F., Lei, G. & Li, H. Do amplitudes of water level fluctuations affect the growth and community structure of submerged macrophytes. PLoS ONE 11, e0146528 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

7.
Wang, M., Qi, S. & Zhang, X. Wetland loss and degradation in the Yellow River Delta, Shandong Province of China. Environ. Earth Sci. 67, 185–188 (2011).
Article  Google Scholar 

8.
Zhou, D., Gong, H., Wang, Y., Khan, S. & Zhao, K. Driving Forces for the Marsh Wetland Degradation in the Honghe National Nature Reserve in Sanjiang Plain Northeast China. Environ. Model. Assess. 14, 101–111 (2008).
Article  Google Scholar 

9.
Zhang, J., Ma, K. & Fu, B. Wetland loss under the impact of agricultural development in the Sanjiang Plain NE China. Environ. Monit. Assess. 166, 139–148 (2009).
PubMed  Article  Google Scholar 

10.
Coops, H., Beklioglu, M. & Crisman, T. L. The role of water-level fluctuations in shallow lake ecosystems – workshop conclusions. Hydrobiologia 506–509, 23–27 (2003).
Article  Google Scholar 

11.
Leira, M. & Cantonati, M. Effects of water-level fluctuations on lakes: an annotated bibliography. Ecol. Effect Water-Level Fluctuat. Lakes 613, 171–184 (2008).
Article  Google Scholar 

12.
Cooke, G. D., Welch, E. B. & Peterson, S. Restoration and management of lakes and reservoirs 2nd edn. (Lewis Publishers, Boca Raton, FL, 1994).
Google Scholar 

13.
Blindow, I. Long- and short-term dynamics of submerged macrophytes in two shallow eutrophic lakes. Freshwater Biol. 28, 15–27 (2010).
Article  Google Scholar 

14.
Gafny, S. Spatially and temporally sporadic appearance of macrophytes in the littoral zone of Lake Kinneret, Israel: taking advantage of a window of opportunity. Aquat. Bot. 62, 249–267 (1999).
Article  Google Scholar 

15.
Lécrivain, N. et al. Water-level fluctuation enhances sediment and trace metal mobility in lake littoral. Chemosphere 264, 128451 (2021).
ADS  PubMed  Article  CAS  Google Scholar 

16.
Yin, D., Peng, F., He, T., Xu, Y. & Wang, Y. Ecological risks of heavy metals as influenced by water-level fluctuations in a polluted plateau wetland, southwest China. Sci. Total Environ. 742, 140319 (2020).
ADS  CAS  PubMed  Article  Google Scholar 

17.
Yin, D. et al. Production and migration of methylmercury in water-level-fluctuating zone of the Three Gorges Reservoir, China: Dual roles of flooding-tolerant perennial herb. J. Hazard. Mater. 381, 120962 (2020).
CAS  PubMed  Article  Google Scholar 

18.
Gownaris, N. J. et al. Water level fluctuations and the ecosystem functioning of lakes. J. Great Lakes Res. 44, 1154–1163 (2018).
Article  Google Scholar 

19.
Liu, J. F. et al. Water-level fluctuations are key for phytoplankton taxonomic communities and functional groups in Poyang Lake. Ecol. Indic. 104, 470–478 (2019).
Article  Google Scholar 

20.
Liu, Q. et al. Vegetation dynamics under water-level fluctuations: Implications for wetland restoration. J. Hydrol. 581, 124418 (2019).
Article  Google Scholar 

21.
Wang, P., Zhang, Q., Xu, Y. & Yu, F. Effects of water level fluctuation on the growth of submerged macrophyte communities. Flora Morphol. Distrib. Funct. Ecol. Plants. 223, 83–89 (2016).
Article  Google Scholar 

22.
Yang, Z. et al. Discharge and water level fluctuations in response to flow regulation in impounded rivers: an analytical study. J. Hydrol. 590, 125519 (2020).
Article  Google Scholar 

23.
Yuan, S., Yang, Z., Liu, X. & Wang, H. Water level requirements of a Carex hygrophyte in Yangtze floodplain lakes. Ecol. Eng. 129, 29–37 (2019).
Article  Google Scholar 

24.
Zweig, C. L. & Kitchens, W. M. Multi-state succession in wetlands: a novel use of state and transition models. Ecology 90, 1900–1909 (2009).
CAS  PubMed  Article  Google Scholar 

25.
Bovo-Scomparin, V. M. & Train, S. Long-term variability of the phytoplankton community in an isolated floodplain lake of the Ivinhema River State Park Brazil. Hydrobiologia 610, 331–344 (2008).
Article  Google Scholar 

26.
Yang, Y., Cao, Y. & Zhang, S. Effects of soil moisture regime on rhizomatic germination and young shoot growth of Carex cinerascens. J. Ecol. Rural Environ. 31, 180–187 (2015).
CAS  Google Scholar 

27.
Liu, L., Liu, D., Johnson, D. M., Yi, Z. Q. & Huang, Y. L. Effects of vertical mixing on phytoplankton blooms in Xiangxi Bay of Three Gorges Reservoir: Implications for management. Water Res. 46, 2121–2130 (2012).
CAS  PubMed  Article  Google Scholar 

28.
Souza, L. Morphology-based functional groups as the best tool to characterize shallow lake-dwelling phytoplankton on an Amazonian floodplain. Ecol. Indic. 95, 579–588 (2018).
Article  CAS  Google Scholar 

29.
Stević, F., Mihaljević, M. & Špoljarić, D. Changes of phytoplankton functional groups in a floodplain lake associated with hydrological perturbations. Hydrobiologia 709, 143–158 (2013).
Article  CAS  Google Scholar 

30.
Lytle, D. A. & Poff, L. R. Adaptation to natural flow regimes. Trends Ecol. Evol. 19, 94–100 (2004).
Article  Google Scholar 

31.
Visser, E. J. W., Bogemann, G. M., Steeg, H. M. V. D., Pierik, R. & Blom, C. W. P. M. Flooding tolerance of Carex species in relation to field distribution and aerenchyma formation. New Phytol. 148, 93–103 (2000).
Article  Google Scholar 

32.
Huber, H. et al. Plasticity as a plastic response: How submergence-induced leaf elongation in rumex palustris depends on light and nutrient availability in its early life stage. New Phytol. 194, 572–582 (2012).
PubMed  Article  Google Scholar 

33.
Fraser, L. H. & Karnezis, J. P. A comparative assessment of seedling survival and biomass accumulation for fourteen wetland plant species grown under minor water-depth differences. Wetlands. 25, 520–530 (2005).
Article  Google Scholar 

34.
Deng, H. et al. Morphology and physiological characteristics of Stachys lanata seedling under water stress. Acta Bot. Boreal 38, 1099–1108 (2018).
Google Scholar 

35.
Tan, Z., Zhang, Q., Li, Y., Xu, X. & Jiang, J. Distribution of typical vegetation communities along elevation in Poyang Lake wetlands. Wetland Sci. 14, 506–515 (2016).
Google Scholar 

36.
Scholte, P. Maximum flood depth characterizes above-Ground biomass in african seasonally shallowly flooded grasslands. J. Trop. Ecol. 23, 63–72 (2007).
Article  Google Scholar 

37.
Yan, B., Dai, Q., Liu, X., Huang, S. & Wang, Z. Flooding-induced membrane damage, lipid oxidation and activated oxygen generation in corn leaves. Plant Soil. 179, 261–268 (2015).
Article  Google Scholar 

38.
Shelford, V. E. Physiological animal geography. J. Morphol. 22, 551–618 (1911).
Article  Google Scholar 

39.
Luan, Z., Wang, Z., Yan, D., Liu, G., Xu, Y. The ecological response of Carex lasiocarpa community in the Riparian Wetlands to the environmental gradient of water depth in Sanjiang Plain, Northeast China. Sci. World J. 1–7 (2013).

40.
Zhang J.T. Quantitative Ecology. Beijing: Science press. 592 (2004).

41.
Zhang, M. et al. Characteristics of the soil seed bank of planted and natural restored draw-down zones in the Three Gorges Reservoir Region. Ecol. Eng. 103, 127–133 (2017).
Article  Google Scholar 

42.
Ma, Y.-R. et al. Effects of flooding on seed viability and nutrient composition in three riparian shrubs and implications for restoration. J. Freshwater Ecol. 33, 449–460 (2018).
CAS  Article  Google Scholar 

43.
Chen, F. et al. Impact of regulated water level fluctuations on the sexual reproduction of remnant Myricaria laxiflora populations. Glob. Ecol. Conserv. 18, e00628 (2019).
Article  Google Scholar 

44.
Webb, J. A., Wallis, E. M. & Stewardson, M. J. A. systematic review of published evidence linking wetland plants to water regime components. Aquat. Bot. 103, 1–14 (2012).
Article  Google Scholar 

45.
Deegan, B. M., White, S. D. & Ganf, G. G. The influence of water level fluctuations on the growth of four emergent macrophyte species. Aquat. Bot. 86, 309–315 (2007).
Article  Google Scholar 

46.
Sarneel, J. M., Janssen, R. H., Rip, W. J., Bender, I. & Bakker, E. S. Windows of opportunity for germination of riparian species after restoring water level fluctuations: a field experiment with controlled seed banks. J. Appl. Ecol. 51, 1006–1014 (2014).
Article  Google Scholar 

47.
Wei, G. et al. Growth responses of eight wetland species to water level fluctuation with different ranges and frequencies. PLoS ONE 14, e0220231 (2019).
CAS  PubMed  PubMed Central  Article  Google Scholar 

48.
Crosby, S. C. et al. Flowering and biomass allocation in US Atlantic coast Spartina alterniflora. Am. J. Bot. 102, 669–676 (2015).
PubMed  Article  Google Scholar 

49.
Mony, C., Mercier, E., Bonis, A. & Bouzillé, J. B. Reproductive strategies may explain plant tolerance to inundation: a mesocosm experiment using six marsh species. Aquat. Bot. 92, 99–104 (2010).
Article  Google Scholar 

50.
Arias, M. E., Wittmann, F., Parolin, P., Murray-Hudson, M. & Cochrane, T. A. Interactions between flooding and upland disturbance drives species diversity in large river floodplains. Hydrobiologia 814, 5–17 (2016).
Article  Google Scholar 

51.
Gattringer, J. P., Ludewig, K., Harvolk-Schöning, S., Donath, T. W. & Otte, A. Interaction between depth and duration matters: flooding tolerance of 12 floodplain meadow species. Plant Ecol. 219, 973–984 (2018).
Article  Google Scholar 

52.
Feng, W., Xu, L., Wang, X., Li, H. & Jiang, J. Response of Carex cinerascens populations to groundwater level gradients in the Poyang Lake wetland. Acta Ecol. Sin. 36, 5109–5115 (2016).
Google Scholar 

53.
Chen, D. et al. A multi-species comparison of selective placement patterns of ramets in invasive alien and native clonal plants to light, soil nutrient and water heterogeneity. Sci. Total Environ. 657, 1568–1577 (2019).
ADS  CAS  PubMed  Article  Google Scholar 

54.
Liu, B. et al. Differential flooding impacts on echinochloa caudata and scirpus planiculmis: implications for weed control in wetlands. Wetlands 36, 1–6 (2016).
Google Scholar 

55.
Hao, S., Cao, H., Wang, H. & Pan, X. The physiological responses of tomato to water stress and re-water in different growth periods. Sci. Hortic. Amsterdam. 249, 143–154 (2019).
CAS  Article  Google Scholar 

56.
Liang, F. et al. Growth and physiological response of Barringtonia acutangular to freshwater flooding stress. J. Southw. For. Univ. (Nat. Sci.). 39, 24–31 (2019).
Google Scholar 

57.
Plazas, M. et al. Comparative analysis of the responses to water stress in eggplant (Solanum melongena) cultivars. Plant Physiol. Biochem. 143, 72–78 (2019).
CAS  PubMed  Article  Google Scholar 

58.
Juan-Ovejero, R., Benito, E., Barreal, M. E., Rodeiro, J. & Briones, M. J. I. Tolerance to fluctuating water regimes drives changes in mesofauna community structure and vertical stratification in peatlands. Pedobiologia 76, 150571 (2019).
Article  Google Scholar 

59.
Wang, H., Liu, X. & Wang, H. The Yangtze River floodplain: threats and rehabilitation fishery resources, environment, and conservation in the Mississippi and Yangtze (Changjiang). River Basins 84, 263–291 (2016).
Google Scholar 

60.
Casanova, M. T. & Brock, M. A. How do depth, duration and frequency of flooding influence the establishment of wetland plant communities?. Plant Ecol. 147, 237–250 (2000).
Article  Google Scholar 

61.
Wang, Q. et al. The growth responses of two emergent plants to the water depth. Acta Hydrobiol. Sin. 3, 583–587 (2012).
Google Scholar 

62.
Zhang, Q. An investigation of enhanced recessions in Poyang Lake: comparison of Yangtze River and local catchment impacts. J. Hydrol. 517, 425–434 (2014).
ADS  Article  Google Scholar 

63.
Fang, C., Cao, W., Mao, J. & Li, H. Relationship between Poyang Lake and Yangtze River and influence of Three Georges Reservoir. J. Hydraul. Eng. 39, 175–181 (2018).
Google Scholar 

64.
Guo, H., Hu, Q. & Wang, Y. Annual variations in climatic and hydrological processes and related flood and drought occurrences in the Poyang Lake Basin. Acta Geogr. Sin. 67, 699–709 (2012).
Google Scholar 

65.
Hu, Z. & Fu, J. Quantitative study on hydrology relationship between the Yangtze River and Poyang Lake and its changes. J. Hydraul. Eng. 49, 570–579 (2018).
Google Scholar 

66.
Dai, X., Wan, R., Yang, G. & Wang, X. Temporal variation of hydrological rhythm in Poyang Lake and the associated water exchange with the Yangtze River. Sci. Geogr. Sin. 34, 1488–1496 (2014).
Google Scholar 

67.
Feng, L. et al. MODIS observations of the bottom topography and its inter-annual variability of Poyang Lake. Remote Sens. Environ. 115, 2729–2741 (2011).
ADS  Article  Google Scholar 

68.
Xie, F. & Wang, S. The meadow of poyang lake. Chin. J. Grassland 1, 8–14 (1984).
Google Scholar 

69.
Wu, Z. & Chen, X. Flora reipublicae popularis sinicae (Science Press, Beijing, 2004).
Google Scholar 

70.
Cao, Y., Wang, J., Guo, Z., Wu, H. & Luo, S. Study on natural distribution and influential factors of Phalaris Arundinacea in Poyang Lake National Wetland Park. Yangtze River 50, 59–64 (2018).
Google Scholar 

71.
Cao, Y., Guo, Z., Yang, Y., Wang, G. & Xie, Z. The ecological amplitude of acorus calamus young shoots under water level gradient. Pol. J. Ecol. 63, 585–592 (2015).
Article  Google Scholar 

72.
Liang, Q., Wang, Z., Lin, B. & Li, C. Effects of flooding stress on seedling physiological indexes and leaf fluorescence characteristics of Sorbus pohuashanensis. J. Jilin For. Sci. Technol. 48, 15–17 (2019).
Google Scholar 

73.
Mommer, L., Lenssen, J. P. M., Huber, H., Visser, E. J. W. & Kroon, H. D. Ecophysiological determinants of plant performance under flooding: a comparative study of seven plant families. J. Ecol. 94, 1117–1129 (2006).
Article  Google Scholar 

74.
Wright, S. D. & Mcconnaughay, K. D. M. Interpreting phenotypic plasticity: the importance of ontogeny. Plant Spec. Biol. 17, 119–131 (2002).
Article  Google Scholar 

75.
Coleman, M. C. S. Biomass allocation in plants: ontogeny or optimality? A test along three resource gradients. Ecology 80, 2581–2593 (1999).
Article  Google Scholar 

76.
Parolin, P. Submergence tolerance vs. escape from submergence: two strategies of seedling establishment in Amazonian floodplains. Environ. Exp. Bot. 48, 177–186 (2002).
Article  Google Scholar 

77.
Zhang, D. et al. Effects of drought and re-flooding on growth and photosynthesis of Carex schmidtii Meinsh: Implication for tussock restoration. Ecol. Indic. 103, 134–144 (2019).
CAS  Article  Google Scholar 

78.
Deegan, B. M., White, S. D. & Ganf, G. G. Nutrients and water level fluctuations: a study of three aquatic plants. River Res. Appl. 28, 359–368 (2012).
Article  Google Scholar 

79.
Warwick, N. W. M. & Brock, M. A. Plant reproduction in temporary wetlands: the effects of seasonal timing, depth, and duration of flooding. Aquat Bot. 77, 153–167 (2003).
Article  Google Scholar 

80.
André, M., Blanch, S. & Grillas, P. Effects of submergence on the growth of Phragmites australis seedlings. Aquat Bot. 69, 147–164 (2001).
Article  Google Scholar 

81.
Venterink, H. O., Wassen, M. J., Belgers, J. D. M. & Verhoeven, J. T. A. Control of environmental variables on species density in fens and meadows: importance of direct effects and effects through community biomass. J. Ecol. 89, 1033–1040 (2001).
Article  Google Scholar 

82.
Hu, Z., Ge, G., Liu, C., Chen, F. & Li, S. Structure of Poyang Lake wetland plants ecosystem and influence of lake water level for the structure. Resour. Environ. Yangtze Basin 19, 597–605 (2010).
Google Scholar 

83.
Gause, G. F. The influence of ecological factors on the size of population. Am. Nat. 65, 70–76 (1931).
Article  Google Scholar 

84.
Zhang, L., Yin, J., Jiang, Y. & Wang, H. Relationship between hydrological conditions and vegetation communities in Poyang Lake national nature reserve of China. Adv. Water Sci. 23, 755–768 (2012).
Google Scholar 

85.
Wu, J. et al. Structure analysis of beach vegetation in Poyang Lake in autumn. Jiangxi Sci. 28, 549–554 (2010).
Google Scholar 

86.
Qi, Q. et al. The driving mechanisms for community expansion in a restored Carex tussock wetland. Ecol. Ind. 121, 107040 (2021).
Article  Google Scholar  More