Ecology

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

Physicochemical properties of water
Physicochemical properties of water and associated environmental factors according to the flow directions of surface flow wetlands are shown in Table 1. Analysis of variance for indexes with p values of less than 0.05 showed that all indexes exhibited large variations during different months.
Table 1 Physicochemical characteristics of the surface-flow constructed wetlands in each unit.
Full size table

Some indexes exhibited large variations during different months according to the flow directions of surface flow wetlands. For example, DO was first reduced and then increased in May, but increased in August and showed differences compared with that in May and October. The salinity was first reduced and then increased in May but then remained stable from August to October. Some indexes showed similar changes according to the flow directions of surface flow wetlands. For example, ORP showed an initial decrease followed by an increase. At the same time, pH, SpCond, TDSs, and TN showed reduced variability over time. The changes in temperature were minimal, although the temperature was higher in the summer and autumn than in the spring.
Surface flow wetlands are in direct contact with the environment and are greatly influenced by outside environmental factors17. Thus, most physicochemical factors of the water samples showed large variability.
Denitrifying bacteria diversity and abundance
For 10 samples from different seasons showing 97% similarity in clustering analysis, the numbers of OTUs differed in May, August, and October (575, 869, and 741, respectively), and the Fig. 2 showed that the denitrifying bacterial abundance indexes (ACEs) were 686.8, 996.2, and 887.3 in May, August, and October, respectively. Additionally, the Shannon-Weiner indexes (H′) were 3.718, 4.303, and 4.432, respectively, indicating that the abundance tended to increase initially, followed by a decrease, and diversity tended to increase. The different seasons affected both the denitrifying bacteria abundance and diversity. Abundance was the largest in August, but its diversity was lower than that in October. These data suggested that the main species became dominant during August, affecting the structure of the denitrifying bacteria.
Figure 2

Biodiversity and abundance of the surface-flow constructed wetland in each unit.

Full size image

For different processing units, the abundance and diversity of denitrifying bacteria varied slightly; both the ACE and H′ index showed low variability. The units F, H, and J showed greater declines than the initial values. In May, the ACE index peaked, with a value of 841 at location I. In August, the ACE index peaked at location E (1251), and that in October peaked at location G (1042). In different months, denitrifying bacteria abundances showed similar changes. Because bacterial diversity in the flowing water and static water were affected by different factors, the surface flow wetlands will be susceptible to various factors, and the bacterial community interactions with internal and external environmental factors will be important for bacterial survival18. Additionally, the number of constructed wetland bacteria decreases as the depth and distance increases19,20, suggesting that denitrifying bacteria may be affected by physical and chemical indicators of changes in water.
Community structure of denitrifying bacteria
Similar OTUs (97% similarity) were identified by sequencing. Database analysis of sequence alignment results revealed that there were many bacteria in the environmental samples that could not be cultivated but that showed high similarity; thus, the denitrifying bacteria were mostly present in the surface flow wetlands and were not cultured. Figure 3 shows statistical analysis of the denitrifying bacterial categories in a histogram format. During the different months, OTUs mainly belonged to seven genera: unclassified bacteria (37.12%), unclassified Proteobacteria (18.16%), Dechloromonas (16.21%), unranked environmental samples (12.51%), unclassified Betaproteobacteria (9.73%), unclassified Rhodocyclaceae (2.14%), Rhodanobacter (1.51%), and other genera (2.62%, representing less than 1% each). Several genera have also been found in surface flow wetlands21 and other types of constructed wetlands22,23,24, albeit with different proportions.
Figure 3

Community structure of the surface-flow constructed wetland in each unit.

Full size image

The same processing units showed different denitrifying bacterial community structures during different seasons and were always changing. Unclassified bacteria showed a greater weight during May for the A processing unit, although its weight was lower than that of Dechloromonas in August. In October, unclassified bacteria had become the most dominant group, and the proportion of Dechloromonas was extremely low. For the B processing unit, from May to October, the proportion of Dechloromonas was decreased, and the proportion of unclassified Proteobacteria was increased, overtaking Dechloromonas. For the C and D processing units, unclassified Proteobacteria were dominant in May, and unclassified bacteria were dominant in August and October. For the E, G, H, I, and J processing units, unclassified bacteria were dominant at all time points. For the F processing unit, the bacterial groups were similar to those of the B processing unit, with proportion of Dechloromonas decreasing and the proportion of unclassified bacteria increasing in October.
Figure 4 shows the means and variances of denitrifying bacteria genus proportions among different processing units and seasons. The means and variances of the dominant genus were large at the same time during different seasons. Thus, the dominant genus often determined the changes in denitrifying bacteria community structures during different seasons in the same unit. However, the greatest variance was observed in the genus Dechloromonas, which was the second most dominant genus in May. This suggested that this genus showed greater changes in different processing units than others. In August, the largest variances were observed in the genus Dechloromonas and in unclassified Proteobacteria, which had lower means than unclassified bacteria. Similar results were observed in October. The largest variance was observed in unclassified Proteobacteria, indicating that the denitrifying bacterial community structures were affected by the second dominant genus over time in the different processing units.
Figure 4

Mean and variance of different denitrifying bacterial taxa in the surface-flow constructed wetland (May, August, and October).

Full size image

For the denitrifying bacteria community structures in different months, we used nonmetric multidimensional scaling to determine the similarities between different processing units during the same months. As shown in Fig. 5, in May, A, B, C, and E showed high similarity, whereas other samples were more dispersed. The distances between D and H and between G and I were shorter than the other distances. F and J were alone in a group. Sample distributions were concentrated in August; C, D, E, F, G, H, and I were relatively similar, and D and E showed maximum similarity. A and B showed some similarity. In contrast, J was distinct. In October, distributions were more dispersed, and the distances between two points were not relatively similar, whereas the differences between the various processing units were higher.
Figure 5

Nonmetric multidimensional scaling map (May, August, and October).

Full size image

Relationships between denitrifying bacteria and environmental factors
Next, we carried out PCA analysis to determine the main factors affecting denitrifying bacteria. After maximum variance orthogonal rotating (p = 0.05), there were two principal component eigenvalues that were greater than the average. The two top principal components contributed to 53.9% and 28.7% of the variance. The first principal component mainly reflected SpCond, TDSs, ORP, and salinity (factor loading was 0.409, 0.403, 0.398, and 0.403, respectively), and the second principal component reflected DO, TN, pH, and temperature (factor loading was 0.449, 0.449, 0.465, and 0.446, respectively). The load distribution characteristics of different environmental factors showed that the surface flow wetlands were affected by the main environmental factors, including temperature, SpCond, DO, pH, ORP, and TN (Fig. 6).
Figure 6

PCA of various environmental factors in the surface-flow constructed wetlands.

Full size image

Table 2 shows bacterial abundance indexes and H′ index for different external environmental factors, as analyzed by Pearson correlation analysis. The results showed that the bacterial abundance was strongly correlated with temperature, DO, and pH, and H′ was strongly correlated with all parameters except TN.
Table 2 Relationships between biodiversity and environment factors.
Full size table

RDA was performed (Fig. 7) for analysis of community distributions and the relationships among environmental factors. For screening of the physicochemical factors of water and the proportions of denitrifying bacterial genera, standardization to center (Monte Carlo permutation) tests were used, and refinement of the information extracted from the first and second axes showed that the total explained variance rate was 80.94%. The results showed that all denitrifying bacterial genera were greatly affected by environmental factors, including temperature and pH, and that the effects of SpCond and ORP were similar. The predominance of unclassified bacteria and unclassified Proteobacter could be explained by positive correlations with temperature, pH, ORP, and SpCond and negative correlations with TN and DO. Dechloromonas showed the opposite trends. In contrast, unranked environmental samples were similar to unclassified Betaproteobacteria, with positive correlations for temperature and pH but negative correlations for TN and DO.
Figure 7

Relationships between denitrifying bacterial community structures and environment factors.

Full size image

Denitrifying bacterial diversity is affected by water nutrient elements and other environmental factors. Most denitrifying bacteria were heterotrophic bacteria. In this study, the autotrophic denitrifying bacteria Dechloromonas accounted for a large proportion in each processing unit25,26,these bacteria can accumulate phosphate and exhibit denitrification activity, partly explaining the lack of TOC removal in association with the observed TN removal. The SpCond of the water reflected its salinity and could be explained by positive correlations with a high proportion of unclassified Proteobacteria. However, SpCond was not generally correlated with denitrifying bacterial abundance. Our results showed that the water SpCond in surface-flow constructed wetlands affected salinity-related denitrifying bacteria but did not affect other denitrifying bacteria. The ORP was positively correlated with denitrifying bacterial genera that were suitable for survival in a strong oxidizing environment, such as unclassified Proteobacteria.
Different physical and chemical properties can influence the structure of the bacterial community owing to the influence of different species on the living environment27,28. In this study, we assessed environmental factors that differed according to season and showed that denitrifying bacteria varied according to some environmental parameters. A comprehensive analysis of the trend of physical and chemical properties of water showed that all parameters except DO and salinity were not highly affected by season and that the trend of the abundances of denitrifying bacteria communities did not change with season along the flow direction of different processing units. However, environmental indicators have a more significant impact on different denitrifying bacteria, which also changes the diversity of denitrifying bacteria community. Accordingly, these results, combined with prediction models of the effects of environmental factors on nitrogen and phosphorus29,30,31, could be used to predict changes in the denitrifying bacterial community structure. More

1.
Sanchez-Gonzalez, A. & Lopez-Mata, L. Plant species richness and diversity along an altitudinal gradient in the Sierra Nevada, Mexico. Divers. Distrib. 11, 567–575 (2005).
Article  Google Scholar 
2.
Dai, L., Feng, Y., Luo, G., Li, Y. & Xu, W. The relationship between soil, climate and forest development in the mid-mountain zone of the Sangong River watershed in the northern Tianshan Mountains, China. J. Arid Land 7, 63–72 (2014).
Article  Google Scholar 

3.
Baiping, Z., Ya, T. & Senguo, M. O. Digital spectrum and analysis of altitudinal belts in the Tianshan Mountains. J. Mt. Res. 1, 18–28 (2004).
Google Scholar 

4.
Pretzsch, H., Biber, P., Schutze, G., Uhl, E. & Rotzer, T. Forest stand growth dynamics in Central Europe have accelerated since 1870. Nat. Commun. 5, 4967 (2014).
ADS  CAS  Article  Google Scholar 

5.
Li, W. J. et al. Effects of climate change on potential habitats of the cold temperate coniferous forest in Yunnan province, southwestern China. J. Mt. Sci. Engl. 13, 1411–1422 (2016).
Article  Google Scholar 

6.
Donat, M. G., Lowry, A. L., Alexander, L. V., O’Gorman, P. A. & Maher, N. Addendum: More extreme precipitation in the world’s dry and wet regions. Nat. Clim. Change 7, 154–158 (2017).
ADS  Article  Google Scholar 

7.
Sun, J., Qin, X. J. & Yang, J. The response of vegetation dynamics of the different alpine grassland types to temperature and precipitation on the Tibetan Plateau. Environ. Monit. Assess. 188, 20 (2016).
Article  Google Scholar 

8.
Windmaisser, T. & Reisch, C. Long-term study of an alpine grassland: Local constancy in times of global change. Alpine Bot. 123, 1–6 (2013).
Article  Google Scholar 

9.
Mahdavi, P., Akhani, H. & Van der Maarel, E. Species diversity and life-form patterns in steppe vegetation along a 3000 m altitudinal gradient in the Alborz Mountains, Iran. Folia Geobot. 48, 7–22 (2013).
Article  Google Scholar 

10.
Rumpf, S. B. et al. Extinction debts and colonization credits of non-forest plants in the European Alps. Nat. Commun. 10, 4293 (2019).
ADS  Article  Google Scholar 

11.
Lamprecht, A., Semenchuk, P. R., Steinbauer, K., Winkler, M. & Pauli, H. Climate change leads to accelerated transformation of high-elevation vegetation in the central Alps. New Phytol. 220, 447–459 (2018).
Article  Google Scholar 

12.
Lenoir, J., Gegout, J. C., Marquet, P. A., de Ruffray, P. & Brisse, H. A significant upward shift in plant species optimum elevation during the 20th century. Science 320, 1768–1771 (2008).
ADS  CAS  Article  Google Scholar 

13.
Kueppers, L. M. et al. Warming and provenance limit tree recruitment across and beyond the elevation range of subalpine forest. Glob. Change Biol. 23, 2383–2395 (2017).
ADS  Article  Google Scholar 

14.
Sedmakova, D. et al. Growth-climate responses indicate shifts in the competitive ability of European beech and Norway spruce under recent climate warming in East-Central Europe. Dendrochronologia 54, 37–48 (2019).
Article  Google Scholar 

15.
Fadrique, B. & Feeley, K. J. Commentary: Novel competitors shape species’ responses to climate change. Front. Ecol. Evol. 4, 33 (2016).
Article  Google Scholar 

16.
Cavieres, L. A. et al. Facilitative plant interactions and climate simultaneously drive alpine plant diversity. Ecol. Lett. 17, 193–202 (2014).
ADS  Article  Google Scholar 

17.
Li, B. F., Chen, Y. N., Chen, Z. S., Xiong, H. G. & Lian, L. S. Why does precipitation in northwest China show a significant increasing trend from 1960 to 2010?. Atmos. Res. 167, 275–284 (2016).
Article  Google Scholar 

18.
Peng, D. D. & Zhou, T. J. Why was the arid and semiarid northwest China getting wetter in the recent decades?. J. Geophys. Res. Atmos. 122, 9060–9075 (2017).
ADS  Article  Google Scholar 

19.
Hong, C. P. et al. Impacts of climate change on future air quality and human health in China. Proc. Natl. Acad. Sci. U.S.A. 116, 17193–17200 (2019).
ADS  CAS  Article  Google Scholar 

20.
Xu, C. C., Chen, Y. N., Chen, Y. P., Zhao, R. F. & Ding, H. Responses of surface runoff to climate change and human activities in the arid region of Central Asia: A case study in the Tarim River Basin, China. Environ Manag. 51, 926–938 (2013).
ADS  Article  Google Scholar 

21.
Deng, H. J., Chen, Y. N., Wang, H. J. & Zhang, S. H. Climate change with elevation and its potential impact on water resources in the Tianshan Mountains, Central Asia. Glob. Planet. Change 135, 28–37 (2015).
ADS  Article  Google Scholar 

22.
Luo, M. et al. Identifying climate change impacts on water resources in Xinjiang, China. Sci. Total Environ. 676, 613–626 (2019).
ADS  CAS  Article  Google Scholar 

23.
Yue, X. Y., Liu, G., Chen, J. M. & Zhou, C. Y. Synergistic regulation of the interdecadal variability in summer precipitation over the Tianshan mountains by sea surface temperature anomalies in the high-latitude Northwest Atlantic Ocean and the Mediterranean Sea. Atmos. Res. 233, UNSP 104717 (2020).
Article  Google Scholar 

24.
Zhang, H. K. K. & Roy, D. P. Using the 500 m MODIS land cover product to derive a consistent continental scale 30 m Landsat land cover classification. Remote Sens. Environ. 197, 15–34 (2017).
ADS  Article  Google Scholar 

25.
Hu, Z. Y., Dietz, A. & Kuenzer, C. The potential of retrieving snow line dynamics from Landsat during the end of the ablation seasons between 1982 and 2017 in European mountains. Int. J. Appl. Earth Obs. 78, 138–148 (2019).
Article  Google Scholar 

26.
Geng, L. Y., Che, T., Wang, X. F. & Wang, H. B. Detecting spatiotemporal changes in vegetation with the BFAST model in the Qilian Mountain region during 2000–2017. Remote Sens. Basel 11, 103 (2019).
ADS  Article  Google Scholar 

27.
Pham, H. T., Marshall, L., Johnson, F. & Sharma, A. A method for combining SRTM DEM and ASTER GDEM2 to improve topography estimation in regions without reference data. Remote Sens. Environ. 210, 229–241 (2018).
ADS  Article  Google Scholar 

28.
Piloyan, A. & Milan, K. Semi-automated classification of landform elements in Armenia based on SRTM DEM using K-means unsupervised classification. Quaest. Geogr. 36, 93–103 (2017).
Article  Google Scholar 

29.
Gonzalez-Moradas, M. D. R. & Viveen, W. Evaluation of ASTER GDEM2, SRTMv3.0, ALOS AW3D30 and TanDEM-X DEMs for the Peruvian Andes against highly accurate GNSS ground control points and geomorphological-hydrological metrics. Remote Sens. Environ. 237, 111509 (2020).
ADS  Article  Google Scholar 

30.
Florinsky, I., Skrypitsyna, T. & Luschikova, O. Comparative accuracy of the AW3D30DSM, ASTER GDEM, and SRTM1 DEM: A case study on the Zaoksky testing ground, Central European Russia. Remote Sens. Lett. 9, 706–714 (2018).
Article  Google Scholar 

31.
Xu, M., Kang, S. C., Wu, H. & Yuan, X. Detection of spatio-temporal variability of air temperature and precipitation based on long-term meteorological station observations over Tianshan Mountains, Central Asia. Atmos. Res. 203, 141–163 (2018).
Article  Google Scholar 

32.
Wu, P., Ding, Y. H., Liu, Y. J. & Li, X. C. The characteristics of moisture recycling and its impact on regional precipitation against the background of climate warming over Northwest China. Int. J. Climatol. 39, 5241–5255 (2019).
Article  Google Scholar 

33.
Lutz, A. F., Immerzeel, W. W., Shrestha, A. B. & Bierkens, M. F. P. Consistent increase in High Asia’s runoff due to increasing glacier melt and precipitation. Nat. Clim. Change. 4, 587–592 (2014).
ADS  Article  Google Scholar 

34.
Chen, Y. N., Li, W. H., Deng, H. J., Fang, G. H. & Li, Z. Changes in Central Asia’s water tower: Past, present and future. Sci. Rep. 6, 35458 (2016).
ADS  CAS  Article  Google Scholar  More

1.
Hutchinson, G. E. Cold spring harbor symposium on quantitative biology. Cold Spring Harb. Symp. Quant. Biol. 22, 415–427 (1957).
Article  Google Scholar 
2.
Chown, S., Gaston, K. & Robinson, D. Macrophysiology: large‐scale patterns in physiological traits and their ecological implications. Funct. Ecol. 18, 159–167 (2004).
Article  Google Scholar 

3.
Spicer, J. I., Morley, S. A. & Bozinovic, F. Physiological diversity, biodiversity patterns and global climate change: testing key hypotheses involving temperature and oxygen. Philos. Trans. R. Soc. B 374 (2019).

4.
Chown, S. L. & Gaston, K. J. Macrophysiology—progress and prospects. Funct. Ecol. 30, 330–344 (2016).
Article  Google Scholar 

5.
Sunday, J. M., Bates, A. E. & Dulvy, N. K. Thermal tolerance and the global redistribution of animals. Nat. Clim. Change 2, 686–690 (2012).
Article  Google Scholar 

6.
Deutsch, C. A. et al. Impacts of climate warming on terrestrial ectotherms across latitude. PNAS 105, 6668–6672 (2008).
CAS  Article  Google Scholar 

7.
Waldock, C., Stuart‐Smith, R. D., Edgar, G. J., Bird, T. J. & Bates, A. E. The shape of abundance distributions across temperature gradients in reef fishes. Ecol. Lett. 22, 685–696 (2019).
Article  Google Scholar 

8.
Payne, N. L. et al. Temperature dependence of fish performance in the wild: links with species biogeography and physiological thermal tolerance. Funct. Ecol. 30, 903–912 (2016).
Article  Google Scholar 

9.
Sunday, J. M., Bates, A. E. & Dulvy, N. K. Global analysis of thermal tolerance and latitude in ectotherms. Proc. R. Soc. B. 278, 1823–1830 (2011).
Article  Google Scholar 

10.
Sunday, J. M. et al. Thermal-safety margins and the necessity of thermoregulatory behavior across latitude and elevation. PNAS 111, 5610–5615 (2014).
CAS  Article  Google Scholar 

11.
Araujo, M. B. et al. Heat freezes niche evolution. Ecol. Lett. 16, 1206–1219 (2013).
Article  Google Scholar 

12.
Payne, N. L. & Smith, J. A. An alternative explanation for global trends in thermal tolerance. Ecol. Lett. 20, 70–77 (2017).
Article  Google Scholar 

13.
Sinclair, B. J. et al. Can we predict ectotherm responses to climate change using thermal performance curves and body temperatures? Ecol. Lett. 19, 1372–1385 (2016).
Article  Google Scholar 

14.
Rezende, E. L. & Bozinovic, F. Thermal performance across levels of biological organization. Philos. Trans. R. Soc. B 374, 20180549 (2019).
CAS  Article  Google Scholar 

15.
Barnes, D. K., Peck, L. S. & Morley, S. A. Ecological relevance of laboratory determined temperature limits: colonization potential, biogeography and resilience of Antarctic invertebrates to environmental change. Glob. Change Biol. 16, 3164–3169 (2010).
Article  Google Scholar 

16.
Clark, T. D., Sandblom, E. & Jutfelt, F. Aerobic scope measurements of fishes in an era of climate change: respirometry, relevance and recommendations. J. Exp. Biol. 216, 2771–2782 (2013).
Article  Google Scholar 

17.
Norin, T., Malte, H. & Clark, T. D. Aerobic scope does not predict the performance of a tropical eurythermal fish at elevated temperatures. J. Exp. Biol. 217, 244–251 (2014).
Article  Google Scholar 

18.
Kearney, M. & Porter, W. Mechanistic niche modelling: combining physiological and spatial data to predict species’ ranges. Ecol. Lett. 12, 334–350 (2009).
Article  Google Scholar 

19.
Buckley, L. B. et al. Can mechanism inform species’ distribution models? Ecol. Lett. 13, 1041–1054 (2010).
Article  Google Scholar 

20.
Brown, J. H., Gillooly, J. F., Allen, A. P., Savage, V. M. & West, G. B. Toward a metabolic theory of ecology. Ecology 85, 1771–1789 (2004).
Article  Google Scholar 

21.
Peck, L. S., Clark, M. S., Morley, S. A., Massey, A. & Rossetti, H. Animal temperature limits and ecological relevance: effects of size, activity and rates of change. Funct. Ecol. 23, 248–256 (2009).
Article  Google Scholar 

22.
Richard, J., Morley, S. A., Thorne, M. A. & Peck, L. S. Estimating long-term survival temperatures at the assemblage level in the marine environment: towards macrophysiology. PLoS ONE 7, e34655 (2012).
CAS  Article  Google Scholar 

23.
Dillon, M. E., Wang, G. & Huey, R. B. Global metabolic impacts of recent climate warming. Nature 467, 704–706 (2010).
CAS  Article  Google Scholar 

24.
Gillooly, J. F., Brown, J. H., West, G. B., Savage, V. M. & Charnov, E. L. Effects of size and temperature on metabolic rate. Science 293, 2248–2251 (2001).
CAS  Article  Google Scholar 

25.
Stuart-Smith, R. D., Edgar, G. J. & Bates, A. E. Thermal limits to the geographic distributions of shallow-water marine species. Nat. Ecol. Evol. 1, 1846–1852 (2017).
Article  Google Scholar 

26.
Janzen, D. H. Why mountain passes are higher in the tropics. Am. Nat. 101, 233–249 (1967).
Article  Google Scholar 

27.
Martin, T. L. & Huey, R. B. Why “Suboptimal” is optimal: Jensen’s inequality and ectotherm thermal preferences. Am. Nat. 171, E102–E118 (2008).
Article  Google Scholar 

28.
Pörtner, H.-O. et al. Cod and climate in a latitudinal cline: physiological analyses of climate effects in marine fishes. Clim. Res. 37, 253–270 (2008).
Article  Google Scholar 

29.
Sylvestre, E.-L., Lapointe, D., Dutil, J.-D. & Guderley, H. Thermal sensitivity of metabolic rates and swimming performance in two latitudinally separated populations of cod, Gadus morhua L. J. Comp. Physiol. B Biochem. Syst. Environ. Physiol. 177, 447–460 (2007).
Article  Google Scholar 

30.
Claireaux, G., Webber, D., Lagardère, J.-P. & Kerr, S. Influence of water temperature and oxygenation on the aerobic metabolic scope of Atlantic cod (Gadus morhua). J. Sea Res. 44, 257–265 (2000).
Article  Google Scholar 

31.
Björnsson, B. & Steinarsson, A. The food-unlimited growth rate of Atlantic cod (Gadus morhua). Can. J. Fish. Aquat. Sci. 59, 494–502 (2002).
Article  Google Scholar 

32.
Morley, S., Peck, L., Sunday, J., Heiser, S. & Bates, A. Physiological acclimation and persistence of ectothermic species under extreme heat events. Glob. Ecol. Biogeogr. 28, 1018–1037 (2019).
Article  Google Scholar 

33.
Edgar, G. J. & Stuart-Smith, R. D. Systematic global assessment of reef fish communities by the Reef Life Survey program. Sci. Data 1, 140007 (2014).
Article  Google Scholar 

34.
Grafen, A. The phylogenetic regression. Philos. Trans. R. Soc. Lond. Ser. B, Biol. Sci. 326, 119–157 (1989).
CAS  Google Scholar 

35.
Orme, D., Freckleton, R., Thomas, G., Petzoldt, T. & Fritz, S. The caper package: comparative analysis of phylogenetics and evolution in R. R. Package Version 5, 1–36 (2013).
Google Scholar 

36.
Halsey, L. G., Butler, P. J. & Blackburn, T. M. A phylogenetic analysis of the allometry of diving. Am. Nat. 167, 276–287 (2006).
Article  Google Scholar 

37.
Michonneau, F., Brown, J. W. & Winter, D. J. rotl: an R package to interact with the Open Tree of Life data. Methods Ecol. Evol. 7, 1476–1481 (2016).
Article  Google Scholar 

38.
Freckleton, R. The seven deadly sins of comparative analysis. J. Evol. Biol. 22, 1367–1375 (2009).
CAS  Article  Google Scholar 

39.
Freckleton, R. P., Harvey, P. H. & Pagel, M. Phylogenetic analysis and comparative data: a test and review of evidence. Am. Nat. 160, 712–726 (2002).
CAS  Article  Google Scholar  More

Radiation experiments
We experimentally quantified radiolytic hydrogen (H2) production in (i) pure water, (ii) seawater, and (iii) seawater-saturated sediment. We irradiated these materials with α- or γ-radiation for fixed time intervals and then determined the concentrations of H2 produced. Sediment samples were slurried with natural seawater to achieve a slurry porosity (φ) of ~0.83, which is the average porosity of abyssal clay in the South Pacific Gyre34. The seawater source is described below. To avoid microbiological uptake of radiolytic H2 during the course of the experiment, seawater and marine sediment slurries were pre-treated with HgCl2 (0.05% solution) or NaN3 (0.1% wt/vol). To ensure that addition of these chemicals did not impact radiolytic H2 yields, irradiation experiments with pure water plus HgCl2 or NaN3 were also conducted. HgCl2 or NaN3 addition had no statistically significant impact on H2 yields5,6,10.
Experimental samples were irradiated in 250 mL borosilicate vials. A solid-angle 137Cs source (beam energy of 0.67 MeV) was used for the γ-irradiation experiments at the Rhode Island Nuclear Science Center (RINSC). The calculated dose rate for sediment slurries was 2.19E−02 Gy h−1 accounting for the (i) source activity, (ii) distance between the source and the samples, (iii) sample vial geometry, and (iv) attenuation coefficient of γ-radiation through air, borosilicate, and sediment slurry. 210Po (5.3 MeV decay−1)-plated silver strips with total activities of 250 μCi were used for the α-irradiation experiments. For α-irradiation of each sediment slurry, a 210Po-plated strip was placed inside the borosilicate vial and immersed in the slurry. Calculated total absorbed doses were 4 Gy and 3 kGy for γ-irradiation and α-irradiation experiments, respectively.
The settling time of sediment grains in the slurries (1 week) was long compared to the time span of each experiment (tens of minutes to an hour for α-experiments, hours to days for γ-experiments). Therefore, we assumed that the suspension was homogenous during the course of each experiment.
H2 concentrations were measured by quantitative headspace analysis via gas chromatography. For headspace analysis, 30 mL of N2 was first injected into the sample vial. To avoid over-pressurization of the sample during injection, an equivalent amount of water was allowed to escape through a separate needle. The vials were then vigorously shaken for 5 min to concentrate the H2 into the headspace. Finally, a 500-μL-headspace subsample was injected into a reduced gas analyzer (Peak Performer 1, PP1). The reduced gas analyzer was calibrated using a 1077 ppmv H2 primary standard (Scott-Marrin, Inc.). A gas mixer was used to dilute the H2 standard with N2 gas to obtain various H2 concentrations and produce a five-point linear calibration curve (0.7, 2, 5, 20, and 45 ppm). H2 concentrations of procedural blanks consisting of sample vials filled with non-irradiated deionized 18-MΩ water were also determined. The concentration detection limit obtained using this protocol was 0.8–1 nM H2. Relative error was less than 5%. Radiation experiments were performed at a minimum in triplicate.
Sample selection and experimental radiolytic H2 yields, G(H2)
Millipore Milli-Q system water was used for our pure-water experiments. For seawater experiments, we used bottom water collected in the Hudson Canyon (water depth, 2136 m) by RV Endeavor expedition EN534. Salinity of North Atlantic bottom water in the vicinity of the Hudson Canyon (34.96 g kg−1) is similar to that of mean open-ocean bottom water (34.70 g kg−1)44,45.
The 20 sediment samples used for the experiments were collected by scientific coring expeditions in three ocean basins (expedition KN223 to the North Atlantic46, expedition KN195-3 to the Equatorial and North Pacific47, International Ocean Discovery Program (IODP) Expedition 329 to the South Pacific Gyre34, MONA expedition to the Guaymas Basin48, expedition EN32 to the Gulf of Mexico49, and expedition EN20 to the Venezuela Basin50). To capture the dominant sediment types present in the global ocean, we selected samples typical of five common sediment types [abyssal clay (11 samples), nannofossil-bearing clay or calcareous marl (2 samples), clay-bearing diatom ooze (3 samples), calcareous ooze (2 samples), and lithogenous sediment (2 samples)]. The locations, lithological descriptions, and mineral compositions of the samples are given in Supplementary Tables 1,  2,  3, and Supplementary Fig. 1. Additional chemical and physical descriptions of the sediment samples used in the radiation experiments can be found in the expedition reports for the expeditions on which the samples were collected34,46.
Energy-normalized radiolytic H2 yields are commonly expressed as G(H2)-values (molecules H2 per 100 eV absorbed)1. As shown in Supplementary Fig. 2, for all irradiated samples (pure water, seawater, and marine sediment slurries), H2 production increased linearly with absorbed α- and γ-ray-dose. We calculated G(H2)-values for each sample and radiation type (α or γ) as the slope of the least-square regression line of radiolytic H2 concentration versus absorbed dose (Supplementary Fig. 2). The error on the yields is less than 10% for each sample. G(H2)-values for each sample and radiation type (α or γ) are reported in Supplementary Table 3.
Although radiolytic OH• is known to react with dissolved organic matter51, total organic content does not appear to significantly impact radiolytic H2 production, since the most organic-rich sediment (e.g., Guaymas Basin and Gulf of Mexico sediment) did not yield particularly high H2 (Supplementary Table 3).
Calculated radiolytic production rates of H2 and oxidants in the cored sediment of individual sites
We calculated radiolytic H2 production rates (PH2, in molecules H2 cm−3 yr−1) for the cored sediment column at nine sites with oxic subseafloor sediment in the North Pacific, South Pacific, and North Atlantic; and seven sites with anoxic subseafloor sediment in the Bering Sea, South Pacific, Equatorial Pacific, and Peru Margin (see Supplementary Fig. 3 for site locations). For these calculations, we used the following equation from Blair et al.2:

$$P_{{mathrm{H}}_2} = {sum} A _{{mathrm{m}},i}rho left( {1 – {{upvarphi}} } right)E_i{mathrm{G}}({mathrm{H}}_2)_i$$
(1)

where i is alpha, beta, or gamma radiation; Am is radioactivity per mass solid; φ is porosity; ρ is density solid; (E_i) is decay energy; and ({mathrm{G}}({mathrm{H}}_2)_i) is radiolytic yield.
We calculated radiolytic oxidant production rates for these sediment columns from the H2 production rates. Because H2 production and oxidant production are stochiometrically balanced in water radiolysis [2H2O → H2 + H2O2], the calculated radiolytic H2 production rates (in electron equivalents) are equal to radiolytic oxidant production rates (in electron equivalents).
The in situ γ- and α-radiation dosages in marine sediment are, respectively, 13 and 15 orders of magnitude lower than the dosage used in our experiments. Because the measured G(H2) for pure water in our γ-irradiation experiment (dose rate = 2.19E-02 Gy h−1) is statistically indistinguishable from previously published G(H2) values at much higher dose rates (ca. 1.00E+3 Gy h−1)5, we infer that the γ-irradiation G(H2) value is constant with dose rate over five orders of magnitude. Similarly, our experimental pure water H2 yields following α-particle irradiation from a 210Po-source (dose rate of 2.55E+03 Gy h−1) are indistinguishable from the yield obtained by Crumière et al.6 [G(H2) = 1.30 ± 0.13] for air-saturated deionized water exposed to a cyclotron-generated He2+ particle beam at higher dose rate (dose rate 1.62E+05 Gy h−1). The close similarity in H2 yields obtained in both experiments implies that (i) radiolytic H2 yield from α-particle irradiation is identical to that from cyclotron-generated He2+ particle irradiation, and (ii) this yield is constant over a two-orders-of-magnitude range dose rate. Therefore, we use our experimentally determined α- and γ-irradiation G(H2) values for the low radiation dose rate found in the subseafloor. Because the G(H2) of β irradiation has not been experimentally determined for water-saturated materials, we assume that the G(H2) of β-radiation matches the G(H2) of γ-radiation for the same sediment types. In pure water, their G(H2) values differ by only 17%1. Because β radiation, on average, contributes only 11% of the total radiolytic H2 production from the U, Th series and K decay in marine sediment, these estimates of total H2 production differ by only 2–5% relative to estimates where the G(H2) of β radiation is assumed equal to that for pure water or for α radiation of the same sediment types.
To calculate H2 production rates for the entire sediment column at seven South Pacific sites and two North Atlantic sites, we measured downcore sediment profiles of U, Th, and K (i) 187 sediment samples from IODP Expedition 329 Sites U1365, U1366, U1367, U1368, U1369, U1370, and U137134,52, and (ii) 40 samples from KN223 expedition Sites 11 and 12 (ref. 46). Total U and Th (ppm) and K2O (wt%) for these sites are reported in the EarthChem SedDB data repository. We measured U, Th, and K abundances using standard atomic emission and mass spectrometry techniques (i.e. ICP-ES and ICP-MS) in the Analytical Geochemistry Facilities at Boston University. Sample preparation, analytical protocol, and data are reported in Dunlea et al.52. The precision for each element is ~2% of the measured value, based on three separate digestions of a homogenized in-house standard of deep-sea sediment.
To calculate H2 production rates for the sediment columns at North Pacific coring Sites EQP10 and EQP11 (ref. 47), we used radioactive element content data from Kyte et al.53, who measured chemical concentrations at high resolution in bulk sediment in core LL44-GPC3. Because Site EQP11 was cored at the same location as LL44-GPC3 (ref. 53) and the sediment retrieved at all three sites is homogeneous abyssal clay, we assume the radioactive element abundances measured in core LL44-GPC3 to be representative of Sites EQP10 and EQP11 (ref. 47). Calculated radiolytic H2 production rates for South Pacific sites are listed in Supplementary Table 4 and for North Atlantic and North Pacific sites in Supplementary Table 5.
For Bering Sea Sites U1343 and U1345 (ref. 54), sedimentary U, Th, and K content measurements are unavailable. Since sediment recovered at these two sites is primarily siliciclastic with a varying amount of diatom-rich clay, we use U, Th, and K concentration values reported for upper continental crust by Li and Schoonmaker for these Bering Sea sites55. Finally, we calculate downhole radiolytic H2 production rates for ODP Leg 201 Sites 1225, 1226, 1227, and 1230 (ref. 35). Sediment compositions for these sites include nannofossil-rich calcareous ooze (Site 1225), alternation of nannofossil (calcareous) ooze and diatom ooze (Site 1226), and siliciclastic with diatom-rich clay intervals (Sites 1227 and 1230). Because sedimentary U, Th, and K measurements are not available for Leg 201 sites, we used average U, Th, and K concentration values measured in North Atlantic46 and South Pacific Sites34,52 with corresponding lithologies.
We use isotopic abundance values reported in Erlank et al.56 to calculate the abundance of 238U, 235U, 232Th, and 40K from the measured ICP-MS values of total U, Th, and K concentration. We then converted radionuclide concentrations to activities using Avogadro’s number and each isotope’s decay constant2. We refer to Blair et al. for a detailed explanation of activities and radiolytic yield calculations2.
Calculation of global radiolytic H2 and oxidant production rates in marine sediment
We calculated global radiolytic H2 production in ocean sediment by applying Eq. (1) (ref. 2) globally. As with the rates at individual sites, we calculated global radiolytic oxidant production (in electron equivalents) from global H2 production and the stochiometry of water radiolysis [2H2O → H2 + H2O2].
Our global radiolytic H2 production calculation spatially integrates calculations of sedimentary porewater radiolysis rates that are based on (i) our experimentally constrained radiolytic H2 yields for the principal marine sediment types, (ii) measured radioactive element content of sediment cores in three ocean basins (North Atlantic46, North Pacific53, and South Pacific34,52), and (iii) global distributions of sediment lithology57, sediment porosity58, and sediment column length59,60.
To generate the global map of radiolytic H2 production, we created global maps of seafloor U, Th, and K concentrations, density, G(H2)-α values, and G(H2)-γ-and-β by assigning each grid cell in our compiled seafloor lithology map (Supplementary Fig. 4) its lithology-specific set of input variables (Supplementary Table 6). Because our model assumes that lithology is constant with depth, U, Th, and K content, grain density, and G(H2)-values are constant with depth.
The G(H2)-values (α, β, and γ radiation), radioactive element content (sedimentary U, Th, and K concentration), density, porosity, and sediment thickness are determined as follows.
Radiolytic yield [G(H2)] for α,β-&-γ radiation
Radiolytic yields for the main seafloor lithologies are obtained by averaging experimentally derived yields for the respective lithologies (Supplementary Table 6). We assume that G(H2)-β values equal G(H2)-γ values.
Sediment lithology
For these calculations of radiolytic chemical production, we generally used seafloor lithologies and assumed that sediment type is constant with sediment depth. For seafloor lithology, the geographic database of global bottom sediment types57 was compiled into five lithologic categories: abyssal clay, calcareous ooze, siliceous ooze, calcareous marl, and lithogenous (Supplementary Fig. 4). Some areas of the seafloor are not described in the database57. These include (i) high-latitude regions (as the seafloor lithology database extends from 70°N to 50°S)57 and (ii) some discrete areas located along continental margins (e.g., Mediterranean Sea, Timor Sea, South China Sea, Supplementary Fig. 4). We used other data sources to identify seafloor lithologies for these regions. We added an opal belt (siliceous ooze) in the Southern Ocean between 57°S and 66°S61,62. The geographic extent of this opal belt was based on DeMaster62 and Dutkiewicz et al.61. We defined the areas of the seafloor from 50°S to 57°S, from 66°S to 90°S, and in the Arctic Ocean as mostly composed of lithogenous material, based on (i) drillsite lithologies in the Southern Ocean [ODP: Site 695 (ref. 63), Site 694 (ref. 63), Site 1165 (ref. 64), Site 739 (ref. 65)], the Bering Sea and Arctic Ocean [International Ocean Discovery Program (IODP): Sites U1343 and U1345 (ref. 54), Site M0002 (ref. 66), ODP: Site 910 (ref. 67), Site 645 (ref. 68)] and between 50°S and 57°S [Deep Sea Drilling Project (DSDP): Site 326 (ref. 69), Ocean Drilling Program (ODP): Site 1138 (ref. 70), Site 1121 (ref. 71)], and Dutkiewicz et al.61.
In the North and South Atlantic, sediment type can be very different at depth than at the seafloor. For these regions, we departed from our assumption that sediment lithology is the same at depth as at the seafloor. Subseafloor lithologies at ODP Sites [1063 (ref. 72), 951 (ref. 73), 925 (ref. 74), and 662 (ref. 75)] and IODP Sites [U1403 (ref. 76) and U1312 (ref. 77)] indicate that sediment in the Atlantic Ocean basin is generally 30–90% biogenic carbonate content and detrital clay78, even where the seafloor lithology is abyssal clay57. Therefore, regions in the Atlantic Ocean described as abyssal clay in the seafloor lithology database57 were characterized as calcareous marl for our calculations (Supplementary Fig. 4). Because abyssal clay catalyzes radiolytic H2 production at a higher rate than calcareous marl, this characterization may underestimate production of radiolytic H2 and radiolytic oxidants in these Atlantic regions.
Radioactive element content
For four of the five lithologic types in our global maps (abyssal clay, siliceous ooze, calcareous ooze, and calcareous marl), we average U, Th, and K concentrations from sites in the North Atlantic46, North Pacific53, and South Pacific34,52. The average U, Th, and K concentration values are consistent with data reported in Li and Schoonmaker55 for the characteristic U, Th, and K content found in abyssal clay and calcareous ooze. For lithogenous sediment, we use U, Th, and K concentration values reported for upper continental crust by Li and Schoonmaker55. Lithology-specific radioactive element values are given in Supplementary Table 6 and used to calculate Am,i in Eq. (1).
Density
Characteristic density values for calcite, quartz, terrigenous clay, and opal-rich sediment were extracted from the Proceedings of the Integrated Ocean Drilling Program Volume 320/321 and are assigned to calcareous ooze, lithogenous sediment, abyssal clay, and siliceous ooze, respectively79.
Global porosity
For global porosity, we use a seafloor porosity data set by Martin et al.58 and accounted for sediment compaction with depth by using separate sediment compaction length scales for continental-shelf (0–200 m water depth; c0 = 0.5 × 10−3), continental-margin (200–2500 m; c0 = 1.7 × 10−3), and abyssal sediment ( >3500 m; c0 = 0.85 × 10−3)80,81. Once the porosity was 0.1%, the depth integration was halted.
Global sediment thickness
We calculated global depth-integrated radiolytic H2 production by summing the seafloor production rates over sediment depth in one-meter intervals (Fig. 3 in main text). Sediment thickness is from Whittaker et al., supplemented with Laske and Masters where needed82,83.
Ocean depth
For porosity calculations, water depths were determined using the General Bathymetric Chart of the Oceans84, resampled to a 5-arc minute grid, i.e. the resolution of the Naval Oceanographic Office’s Bottom Sediment Type (BTS) database “Enhanced dataset”57.
Dissolved H2 concentration profiles
H2 concentrations from South Pacific Sites U1365, U1369, U1370, and U1371, and the measurement protocol, are described in ref. 1. H2 concentrations from North Atlantic KN223 Sites 11, 12, and 15, and North Pacific Site EQP11 were determined using the same protocol and are posted on SedDB (see “Data availability”). The detection limit for H2 ranged between 1 and 5 nM H2, depending on site, and is displayed as gray vertical lines in Fig. 2 of the main text. H2 concentrations for Equatorial Pacific Site 1225 and Peru Trench Site 1230 were measured by the “headspace equilibration technique”, which measures steady-state H2 levels reached following laboratory incubation of the sediment samples85,86.
For comparison to these measured H2 concentrations, we use diffusion-reaction calculations to quantify what in situ H2 concentrations would be in the absence of H2-consuming reactions. The results of these calculations are represented as solid circles (•) in Fig. 2 of the main text. Temporal changes in H2 concentration due to diffusive processes and radiolytic H2 production in situ are expressed by Eq. (2):

$$frac{{partial {mathrm{H}}_2(x,t)}}{{partial t}} = frac{D}{{varphi F}}frac{{partial ^2{mathrm{H}}_2(x,t)}}{{partial x^2}} + P(x)$$
(2)

with
D: the diffusion coefficient of H2(aq) at in situ temperature
(varphi): porosity
F: formation factor
x: depth
Z: sediment column thickness
({mathrm{H}}_2): hydrogen concentration
P: radiolytic H2 production rate
t : time.
With constant radiolytic H2 production, P(x) = P with depth,
and at steady-state,

$$frac{{partial ^2{mathrm{H}}_2(x)}}{{partial x^2}} = – frac{{Pvarphi F}}{D}.$$
(3)

We integrate Eq. (3) over the length x twice,

$${mathrm{H}}_2(x) = – frac{1}{2}frac{{Pvarphi F}}{D}x^2 + Ax + B$$
(4)

where A and B in Eq. (4) are constants of integration. We use two boundary conditions to derive the value of these constants.
Boundary condition 1: concentration of H2 at the sediment-water interface, x = 0, is zero due to diffusive loss to the overlying water column.
Boundary condition 2: concentration of H2 at the basement–sediment-water interface, x = Z, is zero due to diffusive loss to the underlying basement.
With these boundary conditions, (A = frac{1}{2}frac{{Pvarphi F}}{D}Z) and B = 0
and

$${mathrm{H}}_2(x) = frac{1}{2}frac{{Pvarphi F}}{D}(xZ – x^2).$$
(5)

In cases where we expect radiolytic H2 production rates to significantly vary with depth due to changes in lithology, we adapted the boundary conditions and applied a two-layer diffusion model to account for this variation.
Calculation of Gibbs Energies for the Knallgas reaction
For H2 concentrations above the detection limits at South Pacific IODP Expedition 329 sites (Supplementary Fig. 5)34, we quantified in situ Gibbs energies (ΔGr) of the Knallgas reaction (H2 + ½O2 → H2O). In situ ΔGr values depend on pressure (P), temperature (T), ionic strength, and chemical concentrations, all of which are explicitly accounted for in our calculations:

$$Delta G_{mathrm{r}} = Delta G^circ _{mathrm{r}}left( {T,P} right) + 2.3,RT,{mathrm{log}}_{10}Q$$
(6)

where:
ΔGr: in situ Gibbs energy of reaction (kJ mol H2−1)
ΔG°r(T,P): Gibbs energy of reaction under in situ T and P conditions (kJ mol H2−1)
R: gas constant (8.314 kJ−1 mol K−1)
Q: activity quotient of compounds involved in the reaction.
We use the measured composition of the sedimentary pore fluid to determine values of Q.
For a more complete overview of in situ Gibbs energy-of-reaction calculations in subseafloor sediment, see Wang et al.87.
Calculation of organic oxidation rates (net rates of O2 reduction and DIC production)
We calculated the vertical distribution of net O2 reduction rates at nine sites where the sediment is oxic from seafloor to basement and the vertical distribution of DIC production rates at seven sites where the subseafloor sediment is anoxic (see Supplementary Fig. 3 for site locations). Dissolved O2 concentrations are from Røy et al.47 and D’Hondt et al.88. DIC concentrations are from ODP Leg 201 (ref. 35), and the Proceedings of the IODP Expedition 323 (Sites U12343, U1345)54 and IODP Expedition 329 (Site U1371 (ref. 34)).
The net rates are calculated using the MatLab program and numerical procedures of Wang et al.89, modified by using an Akima spline, rather than a 5-point running mean, to generate a best-fit line to the chemical concentration data. Details of the calculation protocol for O2 production rates and DIC production rates are respectively described in the supplementary information of D’Hondt et al.88 and in Walsh et al.90. The DIC reaction rates and their first standard deviations calculated for the seven sites are given in Supplementary Table 7.
To facilitate comparisons of radiolytic chemical rates to net DIC production rates, rates are converted to electron equivalents (2 electrons per H2, 4 electrons per O2, 4 electrons per organic C oxidized).
Estimation of sediment ages
We estimated sediment ages for Sites U1343 and U1343 using the sediment-age model of Takahashi et al.54, which is based on biostratigraphic and magnetostratigraphic data. Because detailed chronostratigraphic data are not available for the remaining sites (Equatorial Pacific sites (1225 and 1226), Peru Trench Site 1230 and Peru Basin Site 1231, South Pacific sites U1365, U1366, U1367, U1369, U1370, and U1371, North Pacific sites EQP9 and EQP10, and North Atlantic sites KN223-11 and KN223-12), we used the mean sediment accumulation rate for each of these sites (Supplementary Fig. 3) to convert its sediment depth (in meters below seafloor) to sediment age (in millions of years, Ma). Mean sediment accumulation rate was calculated by dividing sediment thickness by basement age91 (Supplementary Table 8). For Sites 1225, 1226, 1230, 1231, U1365, U1366, U1367, U1369, U1370, and U1371, sediment thickness was determined by drilling to basement34,35. For Sites EQP9, EQP10, KN223-11, and KN223-12, sediment thicknesses were determined from acoustic basement reflection data. More