Global biogeography of microbes driving ocean ecological status under climate change

Domain: earth_science
Published: 2024-05-31
Editorial outcome: minor_revision

Abstract
Microbial communities play a crucial role in ocean ecology and global biogeochemical processes. However, understanding the intricate interactions among diversity, taxonomical composition, functional traits, and how these factors respond to climate change remains a significant challenge. Here, we propose seven distinct ecological statuses by systematically considering the diversity, structure, and biogeochemical potential of the ocean microbiome to delineate their biogeography. Anthropogenic climate change is expected to alter the ecological status of the surface ocean by influencing environmental conditions, particularly nutrient and oxygen contents. Our predictive model, which utilizes machine learning, indicates that the ecological status of approximately 32.44% of the surface ocean may undergo changes from the present to the end of this century, assuming no policy interventions. These changes mainly include poleward shifts in the main taxa, increases in photosynthetic carbon fixation and decreases in nutrient metabolism. However, this proportion can decrease significantly with effective control of greenhouse gas emissions. Our study underscores the urgent necessity for implementing policies to mitigate climate change, particularly from an ecological perspective.

Introduction
Microbial population of the oceans, the largest environment on Earth, reaches billions per liter of seawater, with an overwhelming diversity and complexity that is governed by both abiotic and biotic factors1. Previous research has provided insights into the diversity, structure and dynamics of microbial communities, which exhibit clear vertical stratification, regional variability, and temporal fluctuations in the oceans2–4. Marine microbes play a critical role in adding, removing, and transforming organic and inorganic materials from seawater, and their complex interactions drive the global biogeochemical fluxes of major elements, including carbon, nitrogen and sulfur5,6. Functional traits, including the ability in biogeochemical cycling, have been widely recognized as important perspectives for research on the ocean microbiome2. Global surveys using high-throughput nucleic acid sequencing have provided comprehensive information about the functional diversity, genomic potential and even transcriptomic activity of the microbial cycling of elements in the ocean7,8. In addition, the frequency of nutrient acquisition genes can be used as biomarkers of nutrient limitation in the ocean9,10.

Compared to other ecosystems on Earth, oceans potentially have the most complicated and changeable environmental conditions on broad spatial and temporal scales; as such, they are highly vulnerable to changes in climate from intensive anthropogenic disturbances, such as emissions of greenhouse gases11. These changing environmental conditions can eventually influence the structure of the microbial community and organismal interactions12,13. Previous research has demonstrated the impacts of anthropogenic climate change on the marine microbial community, including the reorganization of the main taxa12,14,15, loss of diversity16,17, and modification of primary productivity16–18. With recent advances in our understanding of ocean biogeography mediated by climate change, several studies have also identified the variable responses of marine microorganisms to climate change from low- to high-latitude regions and the poleward shifts of some specific taxa14.

Previous studies have raised questions regarding whether the intricate dynamics of microbial diversity, structure, and functional traits align spatially and how they collectively respond to climate changes over time. Although several earlier studies have explored these three ecological dimensions of microbial communities simultaneously7,19,20, they did not offer a systematic indicator for assessing the ecological status of marine microbial communities. In this study, we gathered a substantial dataset of 953 ocean samples from an extensive metagenome sampling project, namely, Bio-GO-SHIP21 (Fig. 1a and Supplementary Data 1), for taxonomic and functional annotations. For functional annotation, we compiled a database of biogeochemical marker genes associated with the core pathways of photosynthesis, carbon fixation, nitrogen metabolism, and sulfur metabolism (Supplementary Data 2). These marker genes are expressed by microbes and contribute to biogeochemical processes in the ocean7. By using the above dataset and database, we overviewed the variations in microbial communities in oceans and their links to environmental conditions. Subsequently, we constructed machine learning models for each microbial index and predicted their current distribution patterns in the global ocean. In this step, we also established ecological status, which serves as a composite representation of microbial communities, considering their functional traits, diversities, and structures. Finally, we harnessed machine learning techniques to elucidate how climate change might affect the future alteration of ecological status and to pinpoint the key drivers behind this transformation. This study is designed to bridge existing knowledge gaps in the biogeographic profiles of the global ocean under the influence of anthropogenic climate change. This topic holds significant importance in guiding management decisions and establishing effective policy objectives aimed at preserving the ecological integrity of our oceans. Furthermore, our comprehensive framework, encompassing ecological status, can be applied to other microbial ecosystems to assess and predict their resistance and resilience in the face of environmental or climate changes.

Results
Temporal and spatial variation in the ocean microbial communities
Bio-GO-SHIP21 is an international multidisciplinary project in which metagenomic samples are collected under standard pipelines without size fractionation (a detailed sampling protocol can be found in the ‘Methods’ section). Here, we collected 953 metagenomic samples from 8 cruises in Bio-GO-SHIP across 26 Longhurst Provinces22 that traveled from 2011 to 2020 (Fig. 1a and Supplementary Data 1). Notably, Longhurst Province was used to show the spatial variation in microbial profiles more clearly. Longhurst Province is a long-standing concept of biogeochemical partitioning in the global ocean23 and has also been used for spatial distribution patterns of chemical conditions24, animals25, protists26 and so on. Our results showed that the richness, Shannon index and relative abundance of the top five phyla and biogeochemical marker genes clearly varied among the different samples across temporal and spatial scales (Fig. 1b and Supplementary Figs. 1–6). Specifically, Bacteroidetes had the highest coefficient of variation, while richness had the lowest one. Besides, coefficient of variation of genes involved in photosystem I, photosystem II, the Calvin cycle, and sulfur oxidation were greater than those of genes involved in other functions (Fig. 1c).

To determine the links between environmental conditions and the microbial community, we obtained ten environmental factors for each sample (Supplementary Data 3) from the Geophysical Fluid Dynamics Laboratory Earth System Model version 4 (GFDL-ESM4)27. These environmental factors include temperature, salinity, partial pressure of carbon dioxide (pCO2), mixed layer depth, and concentrations of dissolved oxygen, nutrients (nitrate, phosphate, silicate, carbonate and iron); these factors have also been used in previous studies for evaluating the response of microbial communities to environmental change7,16,18. Consistent with microbial indices, these factors clearly varied among different samples (Supplementary Fig. 7) and contributed to temporal and spatial variations in both the structure and function of microbial communities (Fig. 1d, e).

Mapping the microbial profiles in the global ocean
The basic logic for our work is as follows: 1) Environmental variables change on broad spatial (e.g., different ocean regions) and temporal scales, and this progress can undoubtedly be influenced by the combination of natural processes and anthropogenic activities in the ocean28. Notably, natural processes (e.g., density and currents) can also be influenced by anthropogenic activities. 2) Changes in environmental conditions can subsequently alter the microbial profiles of the ocean1,5,8.

The tight links between environmental conditions and microbial communities have sparked questions regarding whether we can quantitatively predict microbial profiles in unknown oceans by using existing environmental data. To this end, we constructed a reliable regression model using machine learning for each microbial index (Supplementary Data 4). Our models simultaneously considered multiple environmental factors rather than a single stressor and were more reliable for predicting future conditions16,18. The results of hyperparameter tuning for the XGBoost random forest, linear regression, LASSO regression, and K nearest neighbor algorithms based on tenfold cross-validations showed that the random forest algorithm had the highest prediction performance (Supplementary Fig. 8 and Supplementary Data 5). Thus, we used the random forest algorithm with the best combination of hyperparameters to construct the final regression model for each microbial index (except richness because of the low R2).

Based on the regression models, we quantitatively predicted and mapped microbial profiles in the global ocean (Supplementary Figs. 9–11) under the current (2023) environmental conditions, which were obtained from GFDL-ESM4. We observed that maps of microbial profiles clearly divided the global ocean into several parts, indicating the feasibility of defining a comprehensive ecological index based on microbial profiles.

Definition of the ecological status of ocean microbial communities
A microbial community can be described by multiple ecological dimensions, such as diversity, structure, and functionality, as outlined above. These characteristics in turn determine the ecological status of communities in various habitats. In terrestrial ecosystems, Guerra et al. defined priority areas for soil nature conservation by synthetically considering soil biodiversity and ecosystem services29. Although some previous research has discussed the diversity, structure and functional traits of microbial communities simultaneously7,19,20, we still lack a systematic indicator to describe the ecological status of microbial communities in the global ocean. This situation hampers our ability to understand the comprehensive influence of anthropogenic activities on microbial communities in ocean systems. In this study, we defined the ecological status of an ocean microbial community not only as representing taxonomic and diversity changes, as was done in previous studies14,16, but also considering the biogeochemical potential (Supplementary Fig. 12).

Hyperparameter tuned hierarchical clustering successfully divided the ocean area (0.5° × 0.5° spatial resolution, n = 44564) into seven types of ecological status (defined as ES1 to ES7) (Supplementary Data 6), which had significantly different microbial profiles from each other (Fig. 2a and Supplementary Data 7) and exhibited clear spatial variation (Fig. 2b). For example, most samples from ES4 existed in the tropical Atlantic Ocean and Indian Ocean. While ES7 dominated in the Antarctic area, ES2 was the main type of ecological status in the Arctic area. We also compared our results for 2023 with those for Longhurst Province22. The results showed that some Longhurst Provinces, such as the MONS and WTRA, had a single ecological status (Supplementary Fig. 13). However, the majority of Longhurst Provinces were comprised of multiple types of ecological statuses, indicating that the ecological status as defined in this study is a critical supplement for Longhurst Provinces.

The ecological status in this study could be used as a composite indicator to distinguish among ocean samples from a systematic perspective of the structure, diversity, and biogeochemical potential of the microbial community (Fig. 2c–e). Compared with the other ecological status groups, the ES4 group had the greatest abundance of Cyanobacteria, the lowest diversity, and the greatest functional potential for major biogeochemical processes. ES1 showed a median level of most microbial indices, including the median abundance for each dominant taxon, median diversity, and biogeochemical potential. Other types of ecological status all had low abundances of Cyanobacteria and were dominated by Proteobacteria; however, they can be distinguished from the perspective of diversity or biogeochemical potential. For example, ES3 and ES5 had no obvious difference in diversity (adjusted p value > 0.05 in pairwise comparisons performed by Dunn’s test, with |Cohen’s d| < 0.1), but they were significantly different in terms of biogeochemical progress. Compared with ES3, ES5 had greater functional potential in the phosphate acetyltransferase-acetate kinase pathway (Pta-Ack pathway) and Calvin cycle, but its functional potential in other pathways was much lower than that of ES3.

Overall, the ecological status clustered and defined here can be recognized as a composite description of ocean microbial communities, which we can easily use to construct a classification model and predict whole microbial alteration under changing environmental conditions. We applied machine learning classifiers using the same environmental factors in the regression model to predict the ecological status of the surface ocean (Supplementary Data 8). To reduce the bias resulting from different sample sizes for different types of ecological status, we performed random sampling for each type of ecological status (n = 3000). Compared to the original dataset, the randomly sampled dataset had the same composition of environmental factors (Supplementary Fig. 14) and correlations between them (Supplementary Fig. 15), and was then used as a training dataset in machine learning. And the remains were used as test dataset. The results of hyperparameter tuning for random forest, XGBoost, support vector machine, and logistic regression based on tenfold cross-validation indicated that XGBoost, with n_estimators=200, max_depth=9, and learning_rate=0.3, had the highest prediction accuracy (Supplementary Data 9). Thus, we used the XGBoost algorithm with the above combination of hyperparameters to construct the final classified model (accuracy = 99.09% on the training dataset and 98.31% on the test dataset). The high prediction accuracy of the machine learning models we constructed also confirmed the reliability of the ten environmental factors in classifying ecological status. The confusion matrix also confirmed the good performance of the final model in the classification of ecological status (Supplementary Fig. 16). Besides, phosphate was the most important environmental factor in the prediction of the ecological status of the ocean (Supplementary Fig. 17).

Changes in ecological status in the future ocean
Shared socioeconomic pathways (SSPs) induce a series of scenarios in which greenhouse gas (GHG) emissions (SSP1-1.9 (sustainability), SSP2-4.5 (middle of the road), SSP3-7.0 (regional rivalry), and SSP5-8.5 (fossil-fueled development without policy intervention))30 contribute to the prediction of future changes in microbial communities across diverse habitats14,16,29. A range of distinct outcomes of climate change were observed for ten environmental factors at the end of this century compared to the present under different climate scenarios (Supplementary Figs. 18, 19). In the SSP5-8.5 scenario, which has the highest level of anthropogenic climate change with immoderate fossil fuel development and GHG emissions30, the water temperature and pCO2 increased, but the concentrations of phosphate, carbonate and oxygen decreased in the global ocean. In some specific areas, the concentrations of iron (tropical Pacific Ocean) and nitrate (tropical Atlantic Ocean and Indian Ocean) also increased under the SSP5-8.5 scenario. These changes in environmental conditions incontrovertibly altered the ecological status at the end of this century (Fig. 3a, b and Supplementary Fig. 20).

Our prediction indicated that the ecological status would change in 32.44% of the surface ocean area by the end of this century compared to the present level, because of changes in environmental conditions under the SSP5-8.5 scenario (Fig. 3b). However, if we positively implemented policies to mitigate climate change and control the full range of GHG emissions, such as those outlined in the Paris Climate Agreement31, this proportion of changed areas in the surface ocean could effectively decrease (e.g., to 13.04% under the SSP1–1.9 scenario) (Fig. 3a). In addition, the proportions of different types of ecological status significantly changed in 2100 under the SSP5-8.5 scenario compared to the present scenario (Fig. 3c). Specifically, changes from ES1 to ES4 and from ES6 were common under the SSP5-8.5 scenario (Fig. 3d). However, environmental conditions under the SSP1-1.9 scenario only slightly altered the spatial distribution of each type of ecological status without obvious changes in their proportion in the global ocean (Fig. 3c and Supplementary Fig. 20), further confirming the necessity of implementing climate policies.

Using our regression models, we also focused on the detailed changes in diversity, key taxa, and biogeochemical potential between 2023 and 2100 under different climate change scenarios. We only considered ocean areas with changing ecological status (hereafter called “changed areas”) in this analysis. Overall, the detailed microbial profiles exhibited results similar to those for the ecological status: the SSP5-8.5 scenario had greater impacts on all microbial indices (Supplementary Figs. 21–24). Changes in environmental conditions under the SSP5-8.5 scenario profoundly increased the abundance of Cyanobacteria in 61.70% of the changed area, mainly in the low-latitude regions (Fig. 3e and Supplementary Fig. 22). We also observed an increase in the functional potential for photosynthesis and carbon fixation (e.g., photosystem I, photosystem II, and the Calvin cycle) (Fig. 3f). These results were consistent with previous studies showing that elevated CO2 and temperature under high GHG emissions (Supplementary Fig. 19) promoted the growth of marine phytoplankton32, primary production33 and carbon fixation34. The results also revealed that the abundances of Proteobacteria, Ascomycota, Bacteroidetes and Firmicutes exhibited poleward shifts (increased in high-latitude regions but decreased in low-latitude regions) under future climate scenarios, especially under the SSP5-8.5 scenario (Supplementary Fig. 22), which has been widely mentioned in previous studies12,14 as a thermal adaptation strategy for many microbes17. This poleward shift in the dominant taxa further resulted in considerable decreases in diversity (Supplementary Fig. 21), nitrogen metabolism (Supplementary Fig. 24c) and sulfur metabolism (Supplementary Fig. 24d) in low-latitude regions. However, poleward shifts in thermal adaptation cannot offset the negative impacts of future climate change on the biogeochemical potential of ocean microbes: more ocean areas will suffer from decreasing nitrogen and sulfur metabolism than from increasing nitrogen and sulfur metabolism (Fig. 3f).

Environmental drivers in the alteration of ecological status
In the future, anthropogenic climate change may influence ecological status and microbial profiles by modifying environmental conditions. Consequently, the effects of anthropogenic climate change on changes in ecological status across each grid point of the surface ocean result from a combination of factors: i) the impacts of anthropogenic climate change on environmental factors and ii) the repercussions of these environmental factors on the alteration of ecological status. Although we have already identified phosphate as the most critical factor for predicting ecological status, these factors alone are insufficient for accurately evaluating the individual contributions of each factor to the alterations in ecological status induced by climate change across surface ocean grid points.

Previous studies have reported shifts in the dominant drivers of plankton biogeography reorganization under climate change across different oceanic regions12,14. In line with this, we also quantified the relative contributions of the ten environmental factors to the changes in ecological status induced by climate change, employing the methodology from these prior studies12,14. Anthropogenic climate change primarily reshaped the ecological status by altering nutrient concentrations (including nitrate, phosphate, silicate, iron and carbonate; 45.79% of the changed area) under the SSP1-1.9 scenario (Fig. 4a), while the dissolved oxygen and carbonate contents acted as two major drivers of the changing ecological status under the SSP5-8.5 scenario (Fig. 4b).

Methods
Overview of ocean metagenome samples
We collected 953 metagenome samples of surface oceans from the large-scale metagenome sampling project Bio-GO-SHIP21, which was performed utilizing standard pipelines (Fig. 1a and Supplementary Data 1). The raw sequencing data were then downloaded from the NCBI-SRA database. Detailed information on the metagenomes is available in Supplementary Data 1, including the study ID/title, sample ID/title, sequencing strategy, laboratory, organization, isolation source, collection date, latitude, longitude, geological location and Longhurst Province code. Analysis with a bulk of metagenome samples from public datasets is becoming increasingly available for making new biological discoveries; however, a series of biases caused by distinct protocols of sampling (e.g., size fractionation) and DNA extraction should be noted and reduced before analysis. In this study, we only used the dataset from Bio-GO-SHIP21, which is a well-organized project and can minimize sample bias as follows (more details can be found in Larkin et al.21):

Sampling protocol. Taxonomic and functional traits are completely different for different fractionated microbes (for example, 0.22–3 μm for bacteria and 0–0.2 μm for viruses14), which may cause bias for samples with different size-fractionation schemes. In addition, samples from different depths cannot be compared if we only focus on the horizontal mapping of biogeographic patterns of ocean microbial communities (e.g., microbial profiles in different Longhurst provinces). For all the samples we used, researchers collected whole surface water (with a depth of ~3–7 m) into triple-rinsed containers and gently filtered it through a 0.22 μm pore size Sterivex filter (Millipore, Darmstadt, Germany). After filtration, the DNA was preserved with lysis buffer (4 mM NaCl, 750 μM sucrose, 50 mM Tris-HCl, 20 mM EDTA) and stored at −20 °C before extraction.

DNA extraction protocol. A previous study demonstrated that the method of DNA extraction had a considerable effect on the outcome of metagenomic analysis47. Fortunately, DNA extraction from each sample via Bio-GO-SHIP21 was performed via the same method. Sterivex filters were incubated with lysozyme and proteinase K and SDS buffer. DNA was then extracted from the Sterivex filters using TE buffer (10 mM Tris-HCl, 1 mM EDTA), precipitated in an ice-cold solution of isopropanol and sodium acetate, centrifuged and resuspended in TE buffer for 30 min. Afterwards, the DNA was purified using a genomic DNA Clean and Concentrator kit (Zymo Research Corp., Irvine, CA). Finally, DNA concentrations were quantified using a Qubit dsDNA HS Assay kit and a Qubit fluorometer (Thermo Fisher, Waltham, MA).

Annotation and calculation of the abundance of marker genes and taxa
For functional annotation, we compiled a database of biogeochemical marker genes associated with the core pathways of photosynthesis, carbon fixation, nitrogen metabolism, and sulfur metabolism (Supplementary Data 2). These marker genes can also be expressed by microbes and contribute to biogeochemical processes in the ocean7. We generated the nucleotide sequences for each marker gene from the Kyoto Encyclopedia of Genes and Genomes database (KEGG; release 109.0) and clustered them with 100% sequence similarity using CD-HIT48 (v4.7) to construct a database of biogeochemical marker genes.

We used FastQC (v0.11.5; https://github.com/s-andrews/FastQC) to qualify the raw data of the metagenomic samples and Trimmomatic49 (v0.36) to trim and filter the low-quality reads and obtain clean data for further analysis. The reads were annotated according to the biogeochemical marker genes in the database that we compiled using BWA50 (v0.7.13), a fast and accurate short-read alignment tool. We removed unmapped reads using Samtools (v1.3.1)51, a flexible tool for handling sequence alignment/map format, and then counted the number of mapped reads of cycling genes in each ocean sample using a script available on GitHub (https://github.com/ZhenyanZhang/Ecological_status)52. Additionally, we estimated the average genome size (AGS) by MicrobeCensus (v1.1.1)53 based on the abundance of 30 essential single-copy genes and calculated the genome equivalent for each metagenome sample as follows:

Genome equivalent = Library size (bp) / AGS (bp) (1)

where library size is the total number of sequences. Finally, the abundance of each cycling gene was normalized as reads per kilobase per genome equivalent (RPKG):

RPKG = Mapped reads / (Gene length (kbp) × Genome equivalent) (2)

The use of RPKG can improve the detection of differentially abundant genes between metagenome samples53. Taxonomic annotation was performed by Kraken254 (v2.1.2) with the default parameters, and the abundance of each taxon was also normalized by RPKG.

Links between microbial profiles and environmental factors
The coefficient of variation of each microbial index was defined as the standard deviation divided by the mean. The Longhurst Province, as a long-standing concept of biogeochemical partitioning in the global ocean22, was used to show the spatial distribution patterns of microbial communities more clearly in this study. We obtained the Longhurst province information for each sample by using ArcGIS (v10.8). The sampling time and site from GFDL-ESM4 were searched to obtain the environmental factors for each sample. GFDL-ESM4 from the National Oceanic and Atmospheric Administration (NOAA) provided the present and future values of all ten environmental factors under different scenarios and was also used in a previous study for predicting the reorganization of plankton biogeography under climate change14. These environmental factors include temperature, salinity, partial pressure of carbon dioxide (pCO2), mixed layer depth, and concentrations of dissolved oxygen, nutrients (nitrate, phosphate, silicate, iron) and carbonate, which have also been used in previous studies for evaluating the response of microbial communities to environmental change7,16,18. In total, we collected all ten environmental factors from 890 samples, and the other 63 samples collected before 2015 were not considered in further analysis. The Bray‒Curtis distances of the taxonomic and functional traits of the microbial community in the ocean samples were calculated based on the abundances of the genera and genes, respectively, using the ‘vegan’ R package55 (v2.5–7). To link the temporal and spatial variation of microbial communities with environmental factors, principal coordinate analysis with Bray‒Curtis dissimilarity of the taxonomic and functional traits of the microbial community and its links to environmental factors were calculated using the ‘vegan’ R package (v2.5–7) and visualized by the ‘ggplot2’ R package56 (v3.4.2).

Prediction of each microbial index in the global ocean
A series of regression models were used to predict the microbial indices in the ocean combined with tenfold cross-validation for better performance and less overfitting37. The ten environmental factors mentioned above were used as independent variables for machine learning in this study, while 26 microbial indices (5 dominant taxa, richness, Shannon indices and 19 biogeochemical pathway indices) were included as dependent variables for machine learning. In other words, we presumed that changing environmental conditions directly determine the microbial communities rather than the time and site of the samples. This premise ensures that we can predict microbial profiles and ecological status by using changing environmental data.

To identify the most suitable algorithm and related hyperparameters for machine learning, we first constructed regression models using different algorithms (XGBoost random forest, linear regression, LASSO regression, and K nearest neighbors) with different combinations of hyperparameters with tenfold cross-validation. In hyperparameter tuning for each algorithm, we used only 80% of the samples as the training dataset, and the remaining data were used as the test dataset. After hyperparameter tuning, the best hyperparameter combination for each algorithm and each microbial index was further evaluated by the prediction performance (represented as R2) on the overall dataset and test dataset. In total, random forest had the highest prediction performance for most of the microbial indices and was used to construct the final models (R2 on the test dataset ranged from 0.198 to 0.844), which were saved for further prediction. However, we failed to construct a promising regression model for richness with all combinations of the hyperparameters and algorithms (R2 on the test dataset ranged from −0.026 to 0.025). Thus, we did not consider richness in the following analysis and only used the Shannon index (R2 on the test dataset = 0.809) to represent diversity. After construction of the final models, we also determined the relative importance of each environmental factor in the prediction of each microbial index. All machine learning steps, including hyperparameter tuning, algorithm selection, model construction and feature importance determination, were performed with the ‘sklearn’ Python package (v1.1.3)57.

To delineate the microbial profiles in the current ocean, we obtained ten environmental factors in 2023 (average value in monthly data) from GFDL-ESM4 and predicted each microbial index with a 0.5° × 0.5° spatial resolution and one-year temporal resolution using the final models in the ‘sklearn’ Python package (v1.1.3). Maps of each microbial index in the global ocean were generated by the ‘Basemap’ Python package (v1.3.6).

Determining the ecological status of the ocean samples
To further determine how the diversity, structure, and biogeochemical potential of the microbial communities respond to climate change, we defined the ecological status as a composite description of the ocean microbial communities. We constructed a matrix that included 1) the abundance of the five dominant taxa in the ocean, which constitute >80% of the microbial community; 2) the Shannon index of the ocean microbial community, which is a critical index of diversity; and 3) the abundance of each pathway, which represents the biogeochemical potential of the microbial community. Before performing hierarchical clustering, a principal component analysis was conducted with the ‘sklearn’ Python package to reduce the dimensionality of the data and distance metrics) in the ‘sklearn’ Python package (v1.1.3). The silhouette score, an important internal evaluation index in hierarchical clustering, was used to evaluate clustering performance. Finally, the ‘ward’ clustering algorithm with Euclidean distance was chosen to generate seven clusters, which were defined as ecological status in this study. We then analyzed the detailed global distribution and microbial profiles for each ecological status.

Construction of a machine learning classifier for ecological status
To predict the ecological status in the ocean under future climate change scenarios, we also constructed a machine learning classifier for ecological status with ten environmental factors. The ten environmental factors above were used as independent variables for machine learning in this study, and the seven ecological statuses were included as dependent variables for the machine learning classifier. To reduce the bias resulting from different sample sizes for different types of ecological status, we performed random sampling (n = 3000) for each type of ecological status. Principal component analysis and correlation metrics between features were used to ensure the validity of the randomly sampled dataset (Supplementary Figs. 14, 15), which was then used as a training dataset in machine learning. The remaining dataset after random sampling was used as an independent test dataset.

To identify the most suitable algorithm and related hyperparameters for machine learning, we first constructed a classification model using different algorithms (random forest, logistic regression, support vector machine, and XGBoost) with different combinations of hyperparameters. Tenfold cross-validation was used for machine learning to ensure the performance of the random forest algorithm and avoid overfitting37. XGBoost with n_estimators=200, max_depth=9, and a learning rate=0.3 was then chosen because it was more accurate (average accuracy = 99.09% in tenfold cross-validation) than the other algorithms and hyperparameter combinations (average accuracy ranging from 60.64% to 98.28%). The final classification model was then built using the ‘sklearn’ Python package (v1.1.3). The confusion matrix was also used to confirm the good performance of the final model in the classification of ecological status.

Machine learning prediction
For predicting changes in ecological status and microbial profiles, we collected ten environmental factors used for machine learning with 0.5° × 0.5° spatial resolution and one-year temporal resolution for the surface ocean from GFDL-ESM427 under the SSP1-1.9 and SSP5-8.5 scenarios of climate change with different levels of trajectory of greenhouse gas (GHG) concentrations. We then predicted the changes in ecological status and microbial profiles in the surface ocean with the final models using the ‘sklearn’ Python package (v1.1.3).

Determining the contributions of the environmental factors
We used the methodology from Barton et al.12 to determine the contributions of environmental factors to the alteration of ecological status in the surface ocean induced by climate change. We calculated the probability (P) for each ecological status (n) under two conditions: i) all ten environmental factors (E) using the data from 2023 (present) and ii) all factors using the data from 2023 (present), except for factor e, for which we used the data for 2100 (future). We then calculated the relative contribution of factor e (RCe) to the alteration in ecological status induced by climate change at each grid point with a 0.5° × 0.5° spatial resolution using Eq. (3):

RCe = |P(n, E_future_e) - P(n, E_present)| / Σ |P(n, E_future_e) - P(n, E_present)| (3)

This equation first calculates the total changes in the probability of ecological status when only one environmental factor (for example, salinity) changes in the future. Then, such changes under different environmental factors were added together to calculate the relative contribution of one environmental factor (such as salinity). Finally, the dominant factor at each grid point was defined as the factor with the highest relative contribution at this point.

Development of the ES_predictor
For better use of ecological status, we packed our prediction scripts and related files (such as the final classification model and Python environmental files for XGBoost) by using ‘pyinstaller’ Python package, and developed a software called ES_predictor, which will be licensed under a Creative Commons Attribution 4.0 International license (https://creativecommons.org/licenses/by/4.0/legalcode.en). ES_predictor is an easy-to-use tool developed for Windows users and can be used for determining the ecological status of surface ocean ecosystems with environmental variables. The software package, demonstrations, and instructional manuals can be obtained via https://doi.org/10.6084/m9.figshare.25627293.

Statistical analysis and visualization
Significant differences were identified using multiple methods. The adonis test (two-sided) with 999 permutations was performed with the ‘vegan’ R package (v2.5–7) to determine the temporal and spatial variation in both the structure and function of the microbial communities. The Kruskal–Wallis H test with pairwise comparisons (Dunn’s test) and Cohen’s d calculation were performed with the ‘scipy’ Python package (v1.13.0) to further evaluate the hierarchical clustering performance. The Friedman test with Nemenyi pairwise comparisons was performed with the ‘scipy’ Python package (v1.13.0) to determine the different performances of the five machine learning regression algorithms under the best hyperparameter combinations. All these statistical tests were two-sided, and the p values were adjusted by Bonferroni correction. All the maps were visualized on a world map using the ‘Basemap’ Python package with the same settings in latitude and longitude. Thus, for concision of the figures, we only defined the latitude and longitude in Fig. 1a, and the others were the same. Plots of the principal component analysis of ecological status were visualized using the ‘matplotlib’ Python package. Point plots, density plots, line plots, and box plots were visualized using the ‘ggplot2’ R package. Heatmaps were constructed using TBtools58 (v2.0.42). All schematic diagrams and elements in Fig. 3f and Supplementary Fig. 12 were created with BioRender.com released under a Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International license (https://creativecommons.org/licenses/by-nc-nd/4.0/legalcode.en). All other plots (e.g., pie plots and histograms) were generated using GraphPad Prism (v8.0.2).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This study focusses on a significant region of Earth addressing how global biogeographic pafterns of\nmicrobes are related to ecological funcfioning under scenarios of climate change. Previous work has\naddresses quesfions regarding the diversity, biogeography, temporality, funcfion and relafion to\nenvironment of global microbial populafions. Here, authors propose an indicator to assess the ecological\nstatus of marine microbial communifies. Authors worked with a vast dataset of 1433 ocean samples that\ncame from two metagenomic projects (Bio-GO-SHIP and BioGEOTRACES). Authors analyzed this data for\nfuncfional determinafions of biogeochemically relevant pathways, correlafing to species composifion.\nThis synthesis was compiled as the ecological status per species, accordingly to their funcfional diversity.\nUsing machine learning techniques, authors included global change scenarios to compute the changes in\ndistribufion, abundance, and presence of the microbial communifies analyzed. The approach is very\ninteresfing, and using machine learning to ring together datasets from metagenomics, climate models,\nbiogeography is clearly a fantasfic tool with great potenfial. This study corroborated the basin-wide\ndifferences in biogeography and funcfion of microbial communifies that have been previously reported\n(e.g. Tara Ocean) incorporafing finer points such as the relevance of potenfial elemental-cycling with\nlafitudinal change, which suggests the stability of funcfional traits and redundancy, which is relevant to\nunderstand.\nComments:\nFig S8. Check capfion-\nFig S11. Those variables define density and ocean regions- currents- hence, diversity pafterns will also be\nrelated.\nWould it be possible to dig deeper regarding the funcfional potenfial following the categories proposed\nby authors for C, N, P, S cycles? E.g. in Fig. 1, the relafive abundance is shown but not the genes or\nmetabolic features associated to each major element.\nAre the genefic markers for each element analyzed in flux-models to understand if metabolisms are\ncomplete? Or is the dataset organized in presence/absence of genes without correlafion to metabolic\nnetworks?\nThe delimitafion of oceanic status-provinces accounts for “The structure, diversity, and funcfionality of a\nmicrobial community can vary across the wide range of global environmental condifions which can be\ndetermined by the combinafion of natural processes and anthropogenic acfivifies in the oceans” lines\n182-183\nThe distribufion, structure and funcfion also vary temporarily and spafially. The datasets included in this\nstudy do not represent this. What implicafions does this have for the proposed model?\nAre the sampling campaigns comparable?\nWhere samplings done in similar moments?\nLines 245-246. The 245 increase of element-cycling potenfial under high GHG emission level scenario may\n246 result from the increase of element cycling within specific pathways.\nAlthough, Cyanobacteria are limited under current CO2ppm, hence under GHG emission scenarios, the\nincrease in photosynthefic potenfial would clearly affect the community structure and funcfion.\nLines 251-259. E.g. The Atlanfic basin is changing due to Sargassum blooms- how is this affecfing the\nproposed model? Is this considered?\nThe methods secfion needs to be refined since metabolic networks and hence funcfional classificafions\ncannot be stablished from presence/absence of gene surveys- there are network flow analysis needed to\nunderstand if the pathways are present or not, and their relafive abundance. This is fundamental since\necological status considered cycling potenfial.\nAlso, details are needed regarding the machine learning assumpfions used to infer and model ecological\nstatus. Authors menfion six variables, but these are not defined, explored or analyzed in this study.\nThe sampling protocol for the metagenomic datasets also needs further detail since authors state that\nbiases were accounted for but need explain the considerafions.", "reviewer_2": "This paper unifies metagenomes from two global survey studies, and measures several metrics along\ndiversity, composifion, and elemental cycling capability as inferred from a custom database of ‘element\ncycling genes’ for the C/N/P/S cycles. They then design a set of manual rules to divide these samples into\n‘ecological status’ clusters, and train a random forest classifier to predict the cluster label from\nenvironmental data.\nThe authors convey a passion for the urgency of understanding the effects of climate change on the\nmicrobial communifies that drive biogeochemical cycling in the ocean, and I commend them for this\ngoal. Unfortunately, they fall short of this objecfive in their metric design, and in many ways claim results\nthat are not supported by their data.\nI draw aftenfion to the parficular areas that can be improved below:\nAnalysis of differences across ocean basins:\n-The authors divided their samples into Atlanfic, pacific, and Indian Ocean regions (line 108). I don’t see\nwhy. What’s the mofivafion? Especially given the lackluster stafisfics that follow:\n--The differences they highlight in Fig 1 do not seem supported by their data. The quarfiles in their\nbarcharts all overlap, and the abundance differences are minor, especially considering the small catalog\nof genes they’re considering.\n--Fig 1C-E should report t-test p-values, if these differences are in fact significant, to support their\nconclusions here.\nFig S2 – the authors claim there is a correlafion with specific species, but their maximum R2 value is\n0.155 and most are ~0.01. That is not a notable correlafion, even if the p value is significant. The data is\ntelling us the opposite of what they claim in the paper (line 83).\nSame story with the correlafion with geographic distance. R2 of 0.02 is not a posifive correlafion. (line\n106)\nThe authors then go on to say that there are lafitudinal pafterns (paragraph starfing line 122), but since\nthe Indian Ocean is all low lafitudes this could all be an arfifact of the different ocean basins. What do\nthe authors think is driving this, lafitude or geographic region? How could they differenfiate between the\ntwo given the confoundedness of their data?\nThey also sfipulate lafitudinal differences are due to temperature but don’t show this (line 125). Support\nor move to discussion.\nThere are several fimes where the authors make statements that are befter suited to interpretafion in\nthe discussion than the results secfion, e.g. lines 128-129.\n•\nMy main objecfion to this paper is the formulafion of the ecological status metric. This metric is\nconstructed completely manually, and is very suscepfible to the confounding factors discussed above.\n--There are so many clustering approaches you could use to be quanfitafive here, then you could see\nhow to interpret clusters based off the individual features.\nThe authors make claims about lafitudinal vs temp/salinity differences in their interpretafion for climate\nchange, after poinfing out lafitudinal differences and region differences, and make no effort to\ndisentangle these potenfially confounding factors. This is really necessary if you want to make claims\nabout climate change.\nI like the way the authors interpret Figure 3, where they break out individual contributors in their\necological status clusters to see how they are projected to change. However, this makes me quesfion the\npremise of the paper. You could just build separate models to predict each individual feature, and\ncompare predicfions of each model under different climate scenarios. This would also enable you to use\nquanfitafive values rather than arbitrary categories that have been constructed from a confinuous\nvariable, as has been done here.\nThe authors offer no suggesfion on how to interpret the meaning of their ecological status metric. What\ndoes it mean to have high diversity and low cycling potenfial and be dominated by cyanobacteria? Does\nthis have any impact on the funcfion of these ecosystems? I think this is impossible to deduce as the\nmetric is constructed, and makes me quesfion the ufility of such a metric.\n--Related to this are quesfions about the meaning of the cycling potenfial genes and how they should be\ninterpreted. A lot of N acquisifion genes and a lot of N remineralizafion genes would have very different\neffects on N cycling and funcfion of ecosystems, but would be weighted the same by this metric, I\nbelieve.\nI again commend the authors for an ambifious effort, despite my serious concerns about the validity of\ntheir metric and interpretafion of their approach. I do not think this paper is ready for publicafion, but I\nhope that by revisifing their approach to clustering, metric design, and data interpretafion they can\nachieve their research objecfives in a future paper."}, "author_response": "Response to Reviewers\nReviewer #1 (Remarks to the Author):\nThis study focusses on a significant region of Earth, addressing how global\nbiogeographic patterns of microbes are related to ecological functioning under\nscenarios of climate change. Previous work has addresses questions regarding the\ndiversity, biogeography, temporality, function and relation to environment of global\nmicrobial populations. Here, authors propose an indicator to assess the ecological status\nof marine microbial communities. Authors worked with a vast dataset of 1433 ocean\nsamples that came from two metagenomic projects (Bio-GO-SHIP and\nBioGEOTRACES). Authors analyzed this data for functional determinations of\nbiogeochemically relevant pathways, correlating to species composition. This synthesis\nwas compiled as the ecological status per species, accordingly to their functional\ndiversity. Using machine learning techniques, authors included global change scenarios\nto compute the changes in distribution, abundance, and presence of the microbial\ncommunities analyzed. The approach is very interesting, and using machine learning to\nring together datasets from metagenomics, climate models, biogeography is clearly a\nfantastic tool with great potential. This study corroborated the basin-wide differences\nin biogeography and function of microbial communities that have been previously\nreported (e.g. Tara Ocean) incorporating finer points, such as the relevance of potential\nelementalcycling with latitudinal change, which suggests the stability of functional\ntraits and redundancy, which is relevant to understand.\nResponse: Thank you for your positive evaluation of the interest, relevance and utility\nof our study and your insightful comments, which have greatly improved our work.\nAccording to your detailed comments, we have made substantial revisions to our work\nas follows:\n1) We investigated the functional potential of metabolic networks, including the\nenzymes and reactions related to each core pathway involved in photosynthesis,\ncarbon fixation, nitrogen metabolism, and sulfur metabolism. To this end, we\nrevised our database to include genes that can unambiguously be used as\nmarker genes for specific pathways (Supplementary Data 2). A previous study\ndemonstrated that these marker genes can be expressed by microbes and\ncontribute to biogeochemical processes in the ocean1.\n2) We analysed all the environmental variables used in machine learning. We\nobtained data on ten environmental factors, including temperature, salinity,\npartial pressure of carbon dioxide (pCO ), mixed layer depth, and\n2\nconcentrations of dissolved oxygen and nutrients (nitrate, phosphate, silicate,\ncarbonate, and iron). These factors have also been used in previous studies for\nevaluating the response of microbial communities to environmental change1-3.\nIn the revised manuscript, we analysed their variations in our dataset\n(Supplementary Fig. 7), links to microbial profiles (Figs. 1d, e) and changes\nunder future climate change scenarios (Supplementary Figs. 18 and 19).\n3) We revised our dataset of metagenomic samples and considered only samples\nfrom Bio-GO-SHIP4, which is an international multidisciplinary project. We\nalso provided the sampling protocol and DNA extraction protocol for the\nmetagenomic dataset and explained how they can minimize sample biases and\nmake the sampling campaigns comparable (‘Methods’ section; lines 400-420).\n4) We provided more details in the ‘Results and Discussion’, such as the impacts\nof CO on Cyanobacteria and photosynthetic carbon fixation (lines 259-262).\nversion: the ecological status of approximately 32.44% of the surface ocean may\nundergo changes from the present to the end of this century, assuming no policy\ninterventions, and this proportion could decrease significantly with effective control of\ngreenhouse gas emissions (Fig. 3a, b).\nAbove all, substantial revisions according to your comment greatly refined our\nwork and made our story much more solid and clearer. Thank you again for your\nvaluable comments.\nComments:\nFig S8. Check caption-\nResponse: In the revised manuscript, we checked the information of all the figures (e.g.,\nlegends and captions) and ensured that they clearly show the results of our work. Thank\nyou for making us notice this.\nFig S11. Those variables define density and ocean regions- currents- hence, diversity\npatterns will also be related.\nResponse: We fully agree with your comment. The environmental variables define the\ndensity, ocean regions, and currents, as you mentioned. However, the basic logic for\nour work is as follows (lines 130-136): 1) Environmental variables change on broad\nspatial (e.g., different ocean regions) and temporal scales, and this progress can\nundoubtedly be influenced by the combination of natural processes and anthropogenic\nactivities in the ocean5. Notably, natural processes (e.g., density and currents) can also\nbe influenced by anthropogenic activities. 2) Changes in environmental conditions can\nsubsequently alter the ecological status (defined by multiple microbial profiles) of the\nocean6-8. In this work, we focused on the latter and directly obtained environmental data\nfrom GFDL-ESM4, which has already considered the influence of natural processes\nand anthropogenic activities on environmental variables under different scenarios.\nUsing these ten environmental variables, we analysed the distribution patterns of the\ndiversity, structure and biogeochemical potential of ocean microbial communities.\nWould it be possible to dig deeper regarding the functional potential following the\ncategories proposed by authors for C, N, P, S cycles? E.g. in Fig. 1, the relative\nabundance is shown, but not the genes or metabolic features associated to each major\nelement.\nResponse: Thank you for pointing this out. We dug deeper into the functional potential\nof detailed pathways rather than categories in all analyses of the revised manuscript,\nincluding their spatial and temporal variations, definitions of ecological status and\nfuture predictions.\nAre the genetic markers for each element analyzed in flux-models to understand if\nmetabolisms are complete? Or is the dataset organized in presence/absence of genes\nwithout correlation to metabolic networks?\nResponse: Thank you for raising this question, which is fundamental for our analysis\nof microbial profiles and the ecological status of the ocean. In the revised manuscript,\nwe curated a set of biogeochemical marker genes in the core pathways involved in\nphotosynthesis, carbon fixation, nitrogen metabolism, and sulfur metabolism. We also\nprovided detailed information on the enzymes and reactions encoded by these genes in\nmetabolic networks (Supplementary Data 2; Fig. 3f).\nThe delimitation of oceanic status-provinces accounts for “The structure, diversity, and\nfunctionality of a microbial community can vary across the wide range of global\nenvironmental conditions which can be determined by the combination of natural\nprocesses and anthropogenic activities in the oceans” lines 182-183\nThe distribution, structure and function also vary temporarily and spatially. The datasets\nincluded in this study do not represent this. What implications does this have for the\nproposed model?\nResponse: We agree with your comments that the distributions of structure and\nfunction vary temporarily and spatially. We added the related analysis to the revised\nmanuscript. The results showed that the richness, Shannon index and abundance of the\ndominant phyla and biogeochemical marker genes clearly varied among the different\nsamples across temporal and spatial scales (Fig. 1b and Supplementary Figs. 1-6),\nwhich resulted from the variation in environmental conditions (Fig. 1d, e and\nSupplementary Fig. 7), confirming the tight links between environmental conditions\nand microbial communities. Based on these results, we constructed machine learning\nmodels that considered environmental factors as independent variables and microbial\nprofiles and ecological status as dependent variables. In other words, we presumed that\nchanging environmental conditions directly determine the microbial communities\nrather than the time and site of the samples. This premise ensures that we can predict\nmicrobial profiles and ecological status by using changing environmental data.\nAre the sampling campaigns comparable?\nResponse: Yes, they are. We revised our dataset of metagenomic samples and\nconsidered only samples from Bio-GO-SHIP4, which is an international\nmultidisciplinary project. We also provided the detailed sampling protocol and DNA\nextraction protocol for the metagenomic dataset and explained how they can minimize\nsample biases and make the sampling campaigns comparable (‘Methods’ section; lines\n400-410).\nWhere samplings done in similar moments?\nResponse: In the revised manuscript, we collected metagenomic samples from 8\ncruises in Bio-GO-SHIP that travelled from 2011 to 2020 (Fig. 1a and Supplementary\nData 1). Therefore, the samplings were not performed at similar times. Based on the\ntight links between environmental conditions and microbial communities, we\nconstructed machine learning models that considered environmental factors as\nindependent variables and microbial profiles and ecological status as dependent\nvariables. In other words, we presumed that changing environmental conditions directly\ndetermine the microbial communities rather than the time and site of the samples. We\nobtained the environmental factors for each sample according to their detailed sampling\ntime and site to reduce the uncertainty caused by data collection.\nLines 245-246. The 245 increase of element-cycling potential under high GHG\nemission level scenario may 246 result from the increase of element cycling within\nspecific pathways.\nAlthough, Cyanobacteria are limited under current CO2ppm, hence under GHG\nemission scenarios, the increase in photosynthetic potential would clearly affect the\ncommunity structure and function.\nResponse: Thank you for pointing this out. We added pCO as one of the key\n2\nenvironmental variables for altering ecological status and detailed microbial indices.\nChanges in environmental conditions (including dramatic increases in pCO ) under the\nthe changed area, mainly in the low-latitude regions (Fig. 3e and Supplementary Fig.\n22). We also observed increasing functional potential for photosynthesis and carbon\nfixation (e.g., photosystem I, photosystem II and the Calvin cycle) (Fig. 3f). These\nresults were consistent with previous studies showing that elevating CO and\n2\ntemperature under high GHG emissions promoted the growth of marine phytoplankton9,\nprimary production10 and carbon fixation11.\nLines 251-259. E.g. The Atlantic basin is changing due to Sargassum blooms- how is\nthis affecting the proposed model? Is this considered?\nResponse: Thank you for this suggestion. The occurrence of Sargassum blooms in the\nAtlantic basin is increasing, similar to the occurrence of extreme events in global\nclimatic systems. Both of these factors induce complexity and uncertainty in future\nenvironmental conditions, which play a decisive role in our prediction of ecological\nstatus. Unfortunately, we cannot obtain available data for Sargassum blooms in the\npresent or future ocean. We mentioned this limitation at the current stage in the\ndiscussion (lines 358-361). However, given the urgency of maintaining global\necological integrity and addressing climate change crises12, we must operate under a\ndegree of uncertainty13. This study provides up-to-date information to help scientists\nand policymakers collaboratively address these global crises.\nThe methods section needs to be refined since metabolic networks and hence functional\nclassifications cannot be stablished from presence/absence of gene surveys- there are\nnetwork flow analysis needed to understand if the pathways are present or not, and their\nrelative abundance. This is fundamental since ecological status considered cycling\npotential.\nResponse: Thank you for raising this question, which is fundamental for our analysis\nof microbial profiles and the ecological status of the ocean. In the revised manuscript,\nwe curated a set of biogeochemical marker genes in the core pathways involved in\nphotosynthesis, carbon fixation, nitrogen metabolism, and sulfur metabolism. We also\nprovided detailed information on the enzymes and reactions encoded by these genes in\nmetabolic networks (Supplementary Data 2; Fig. 3f).\nAlso, details are needed regarding the machine learning assumptions used to infer and\nmodel ecological status. Authors mention six variables, but these are not defined,\nexplored or analyzed in this study.\nResponse: We apologize for the lack of these important details. We analysed all the\nenvironmental variables used in the revised manuscript. We obtained data on ten\nenvironmental factors, namely, temperature, salinity, partial pressure of carbon dioxide\n(pCO ), mixed layer depth, and concentrations of dissolved oxygen and nutrients\n2\n(nitrate, phosphate, silicate, carbonate, and iron). These factors have also been used in\nprevious studies for evaluating the response of microbial communities to environmental\nchange1-3. In the revised manuscript, we analysed their variations in our dataset\n(Supplementary Fig. 7), links to microbial profiles (Figs. 1d, e) and changes under\nfuture climate change scenarios (Supplementary Fig. 18 and 19).\nThe sampling protocol for the metagenomic datasets also needs further detail since\nauthors state that biases were accounted for but need explain the considerations.\nResponse: We revised our dataset of metagenomic samples and considered only\nsamples from Bio-GO-SHIP4, which is an international multidisciplinary project. We\nalso provided the detailed sampling protocol and DNA extraction protocol for the\nmetagenomic dataset and explained how they can minimize sample biases and make\nthe sampling campaigns comparable (‘Methods’ section; lines 400-410).\nReviewer #2 (Remarks to the Author):\nThis paper unifies metagenomes from two global survey studies, and measures several\nmetrics along diversity, composition and elemental cycling capability as inferred from\na custom database of ‘element cycling genes’ for the C/N/P/S cycles. They then design\na set of manual rules to divide these samples into ‘ecological status’ clusters, and train\na random forest classifier to predict the cluster label from environmental data.\nThe authors convey a passion for the urgency of understanding the effects of climate\nchange on the microbial communities that drive biogeochemical cycling in the ocean,\nand I commend them for this goal. Unfortunately, they fall short of this objective in\ntheir metric design, and in many ways, claim results that are not supported by their data.\nI draw attention to the particular areas that can be improved below:\nResponse: Thank you for your appreciation of our research goal and your valuable\ncomments, which greatly improved our work and made our manuscript more sound and\nclearer. We have reanalysed all the data and rewritten the manuscript according to your\ncomments as follows:\n1) We used hyperparameter-tuned hierarchical clustering rather than manual\nclustering to cluster and define the ecological status in the revised manuscript.\nThe silhouette score, an important internal evaluation index in hierarchical\nclustering, was used to evaluate clustering performance. Finally, the ‘ward’\nclustering algorithm with Euclidean distance was chosen to generate seven\nclusters (Supplementary Data 6), which were defined as ecological status in\nthis study. The Kruskal–Wallis H test with pairwise comparisons (Dunn's test)\nand Cohen's d were performed to further evaluate the clustering performance\n(Supplementary Data 7). After that, we analysed the detailed global\ndistribution and microbial profiles for each ecological status.\n2) We constructed a series of reliable regression models using machine learning\nfor each microbial index to quantitatively predict and map detailed microbial\nprofiles in the global ocean (Supplementary Data 5). Additionally, to predict\nthe overall alterations in ecological status in the future ocean, we constructed\na classification model for ecological status (Supplementary Data 9). Based on\nthe tight links between environmental conditions and microbial communities,\nwe constructed these machine learning models by considering environmental\nfactors as independent variables and microbial indices and ecological status as\ndependent variables. Our models simultaneously considered multiple\nenvironmental factors rather than a single stressor and were more reliable for\npredicting future conditions.\n3) We carefully checked all the statistical analysis in the revised manuscript and\nprovided the p values or original data. Specifically, we calculated Cohen's d to\nevaluate the extent of difference in microbial indices between different types\nof ecological status (Supplementary Data 7).\n4) We revised the discussion to avoid overstating our results. In addition, we\ndivided the ‘Results and Discussions’ section into two sections in the revised\nmanuscript. This makes our descriptions more precise and clearer.\nDespite the above revisions, we still obtained similar results to those in the\nprevious version: the ecological status of approximately 32.44% of the surface ocean\nmay undergo changes from the present to the end of this century, assuming no policy\ninterventions, and this proportion could decrease significantly with effective control of\ngreenhouse gas emissions (Fig. 3a, b).\nAbove all, substantial revisions according to your comment greatly refined our\nwork and made our story much more solid and clearer. Thank you again for your\nvaluable comments.\nAnalysis of differences across ocean basins:\n-The authors divided their samples into Atlantic, pacific, and Indian Ocean regions (line\n108). I don’t see why. What’s the motivation? Especially given the lackluster statistics\nthat follow:\n--The differences they highlight in Fig 1 do not seem supported by their data. The\nquartiles in their barcharts all overlap, and the abundance differences are minor,\nespecially considering the small catalog of genes they’re considering.\n--Fig 1C-E should report t-test p-values, if these differences are in fact significant, to\nsupport their conclusions here.\nFig S2 – the authors claim there is a correlation with specific species, but their\nmaximum R2 value is 0.155 and most are ~0.01. That is not a notable correlation, even\nif the p value is significant. The data is telling us the opposite of what they claim in the\npaper (line 83).\nSame story with the correlation with geographic distance. R2 of 0.02 is not a positive\ncorrelation. (line 106)\nResponse: Thank you for pointing out these mistakes in the analysis of differences\nacross ocean basins. In the revised manuscript, Longhurst Province was used to show\nthe spatial variation in microbial profiles more clearly, as it is a long-standing concept\nof biogeochemical partitioning in the global ocean14 and has also been used for spatial\ndistribution patterns of chemical conditions15, animals16, protists17 and so on. In\naddition, we carefully checked all the statistical analyses in the revised manuscript and\nprovided the p values or original data. Specifically, we calculated Cohen's d to evaluate\nthe extent of difference in microbial indices between different types of ecological status\n(Supplementary Data 7).\nThe authors then go on to say that there are latitudinal patterns (paragraph starting line\n122), but since the Indian Ocean is all low latitudes this could all be an artifact of the\ndifferent ocean basins. What do the authors think is driving this, latitude or geographic\nregion? How could they differentiate between the two given the confoundedness of\ntheir data?\nThey also stipulate latitudinal differences are due to temperature but don’t show this\n(line 125). Support or move to discussion.\nResponse: We completely agree. Latitudes and ocean regions are two spatial factors\nfor the variation in microbial indices and environmental conditions, but their\ncombination increases the confoundedness of the analysis. Different latitudes have\nhighly distinct patterns of terrestrial systems (e.g., the Southern Hemisphere has more\nocean regions than does the Northern Hemisphere; the Indian Ocean all has low\nlatitudes), which influences the analysis of latitudinal patterns of microbial profiles.\nThus, we did not consider latitude in our revised manuscript.\nThere are several times where the authors make statements that are better suited to\ninterpretation in the discussion than the results section, e.g. lines 128-129.\nResponse: Thank you for your valuable suggestions. In the revised manuscript, we\ndivided the ‘Results and Discussions’ section into two sections. However, potential\nreasons for some descriptions of our analysis and results need to be explained\nimmediately to validate our analysis and results. Thus, we also included the part of\nsupplementary descriptions of our results in the ‘Results’ section. In the ‘Discussion’\nsection, we focused on the main findings in our work and their ecological significance.\nIn addition, we discuss the limitations of our work in the current stage and propose\nsome feasible suggestions for future studies. Thank you again for making our\ndescriptions more precise and clearer.\nMy main objection to this paper is the formulation of the ecological status metric. This\nmetric is constructed completely manually, and is very susceptible to the confounding\nfactors discussed above.\n--There are so many clustering approaches you could use to be quantitative here, then\nyou could see how to interpret clusters based off the individual features.\nResponse: Thank you for your valuable suggestion. We used hyperparameter-tuned\nhierarchical clustering rather than manual clustering to cluster and define the ecological\nstatus in the revised manuscript. This reversion makes the ecological status more\nmeaningful and greatly reduces subjectivity. The silhouette score, an important internal\nevaluation index in hierarchical clustering, was used to evaluate clustering performance.\nFinally, the ‘ward’ clustering algorithm with Euclidean distance was chosen to generate\nseven clusters, which were defined as ecological status in this study. The Kruskal–\nWallis H test with pairwise comparisons (Dunn's test) and Cohen's d were performed\nto further evaluate the clustering performance (Supplementary Data 7). After that, we\nanalysed the detailed global distribution and microbial profiles for each ecological\nstatus.\nThe authors make claims about latitudinal vs temp/salinity differences in their\ninterpretation for climate change, after pointing out latitudinal differences and region\ndifferences, and make no effort to disentangle these potentially confounding factors.\nThis is really necessary if you want to make claims about climate change.\nResponse: We apologize for the unclear description. The basic logic for our work is as\nfollows (lines 130-136): 1) Environmental variables change on broad spatial (e.g.,\ndifferent ocean regions) and temporal scales, and this progress can undoubtedly be\ninfluenced by the combination of natural processes and anthropogenic activities in the\nocean5. Notably, natural processes (e.g., density and currents) can also be influenced\nby anthropogenic activities. 2) Changes in environmental conditions can subsequently\nalter the ecological status (defined by multiple microbial profiles) of the ocean6-8. This\nlogic suggests that environmental conditions (such as temperature/salinity differences)\ncan naturally differ in different ocean regions; however, climate changes caused by\nintense anthropogenic activities reshape them.\nIn the revised manuscript, we analysed variations in environmental conditions in\nour dataset (Supplementary Fig. 7) and their changes under future climate change\nscenarios (Supplementary Figs. 18 and 19). Our results also confirmed their\ncontributions to the temporal and spatial variations in the richness, Shannon index and\nabundance of the dominant phyla and biogeochemical marker genes of ocean microbial\ncommunities (Fig. 1d, e and Supplementary Fig. 7).\nI like the way the authors interpret Figure 3, where they break out individual\ncontributors in their ecological status clusters to see how they are projected to change.\nHowever, this makes me question the premise of the paper. You could just build separate\nmodels to predict each individual feature and compare predictions of each model under\ndifferent climate scenarios. This would also enable you to use quantitative values rather\nthan arbitrary categories that have been constructed from a continuous variable, as has\nbeen done here.\nResponse: Thank you for your insightful suggestion, which made our analysis much\nmore sound. Following your suggestion, we constructed a series of reliable regression\nmodels using machine learning for each microbial index to quantitatively predict and\nmap detailed microbial profiles in the global ocean.\nA microbial community can be described by multiple ecological dimensions, such\nas diversity, structure, and functionality, as outlined in our regression models. These\ncharacteristics in turn determine the ecological status of communities in various habitats.\nIn terrestrial ecosystems, Guerra et al. defined priority areas for soil nature conservation\nby synthetically considering soil biodiversity and ecosystem services29. Although some\nprevious research has discussed the diversity, structure and functional traits of microbial\ncommunities simultaneously7,19,20, we still lack a systematic indicator to describe the\necological status of microbial communities in the global ocean. This situation hampers\nour ability to understand the comprehensive influence of anthropogenic activities on\nmicrobial communities in ocean systems. Thus, we defined the ecological status of an\nocean microbial community not only as representing taxonomic and diversity changes,\nas was done in previous studies14,16, but also considering the biogeochemical potential\nof the ocean. In addition, our predictions under different climate change scenarios also\ndemonstrated that alterations in ecological status could be effectively summarized and\nrepresented by changes in detailed microbial profiles, including poleward shifts in the\nmain taxa, increases in photosynthetic carbon fixation and decreases in nutrient\nmetabolism (Figs. 3e, f and Supplementary Figs. 21-24). Thus, ecological status can be\nrecognized as a more convenient and comprehensive index for evaluating the influence\nof anthropogenic activities on microbial communities in ocean systems.\nThe authors offer no suggestion on how to interpret the meaning of their ecological\nstatus metric. What does it mean to have high diversity and low cycling potential and\nbe dominated by cyanobacteria? Does this have any impact on the function of these\necosystems? I think this is impossible to deduce as the metric is constructed, and makes\nme question the utility of such a metric.\nResponse: Thank you for pointing this out. We used hyperparameter-tuned hierarchical\nclustering rather than complete manual clustering to cluster and define the ecological\nstatus in the revised manuscript. This revision made the ecological status more\nmeaningful and greatly reduced the subjectivity.\nA microbial community can be described by multiple ecological dimensions, such\nas diversity, structure, and functionality, as outlined in our regression models. These\ncharacteristics in turn determine the ecological status of communities in various habitats.\nIn terrestrial ecosystems, Guerra et al. defined priority areas for soil nature conservation\nby considering soil biodiversity and ecosystem services29. Although some previous\nresearch has discussed the diversity, structure and functional traits of microbial\ncommunities simultaneously7,19,20, we still lack a systematic indicator to describe the\necological status of microbial communities in the global ocean. This situation hampers\nour ability to understand the comprehensive influence of anthropogenic activities on\nmicrobial communities in ocean systems. Thus, we defined the ecological status of an\nocean microbial community not only as representing taxonomic and diversity changes,\nas was done in previous studies14,16, but also considering the biogeochemical potential\nof the ocean.\nOur predictions under different climate change scenarios also demonstrated that\nalterations in ecological status could be effectively summarized and represented by\nchanges in detailed microbial profiles, including poleward shifts in the main taxa,\nincreases in photosynthetic carbon fixation and decreases in nutrient metabolism (Figs.\n3e, f and Supplementary Figs. 21-24). Thus, ecological status can be recognized as a\nmore convenient and comprehensive index for evaluating the influence of\nanthropogenic activities on microbial communities in ocean systems.\nFor better use of ecological status, we provide an easy-to-use tool called\nES_predictor (available at https://doi.org/10.6084/m9.figshare.25627293), which can\nbe used for determining the ecological status of surface ocean ecosystems with\nchanging environmental variables under diverse anthropogenic pressures depending on\nthe research goals. For example, if researchers have already measured the ten\nenvironmental variables in areas under different fishing pressures, they can easily\nevaluate the comprehensive impacts of fishing on microbial ecology by our tools\nwithout metagenomic sequencing, which is costly and requires bioinformatics.\nConsidering the amount of data currently available and in the future, this tool will also\nbe continually refined. (lines 343-352)\n--Related to this are questions about the meaning of the cycling potential genes and how\nthey should be interpreted. A lot of N acquisition genes and a lot of N remineralization\ngenes would have very different effects on N cycling and function of ecosystems, but\nwould be weighted the same by this metric, I believe.\nResponse: We completely agree. We investigated the functional potential of the\nmetabolic networks in the revised manuscript, including the enzymes and reactions\nrelated to each core pathway involved in photosynthesis, carbon fixation, nitrogen\nmetabolism, and sulfur metabolism (Supplementary Data 2). This approach was useful\nand invaluable for overviewing and predicting the genetic potential of microbes in\nbiogeochemical processes of the global ocean and their response to climate change.\nHowever, different genes may affect biogeochemical progress to different extents,\nas you mentioned, and genetic potential is not completely linked to actual\ntranscriptomic activity, metabolomic composition or element cycling. This limitation\nof our work can be solved by using multidisciplinary methods, such as transcriptomics,\nmetabolomics and chemometrics. But at present, these approaches still lack\ninternational projects with standard sampling protocols. We mentioned this limitation\nat the current stage in the discussion (lines 373-383).\nI again commend the authors for an ambitious effort, despite my serious concerns about\nthe validity of their metric and interpretation of their approach. I do not think this paper\nis ready for publication, but I hope that by revisiting their approach to clustering, metric\ndesign, and data interpretation they can achieve their research objectives in a future\npaper.\nResponse: Thank you again for your appreciation of our research goal and your\nvaluable comments, which greatly improved our work and made our story increasingly\nsound. We have reanalysed all the data and rewritten the manuscript according to your\ncomments. We hope our revisions can convince you of the value of our study and can\nlead you to recommend its acceptance. Many thanks."}, {"round": 2, "reviewer_comments": {"reviewer_1": "In this study, Qian et al.,focus on a significant region of Earth addressing how global biogeographic pafterns of microbes are related to ecological funcfioning under scenarios of climate change. I want to congratule authors for taking in to account all concerns raised during the first round of review. Their arficle is now more solid and clear. This is fundamental since authors make predicfions of biogeochemical funcfion shift under different scenarios of Global Climate Change, hence, the verificafion of metabolic pathways and methods was fundamental. The approach is very interesfing, and using machine learning to bring together datasets from metagenomics, climate models, biogeography is clearly a fantasfic tool with great potenfial.", "reviewer_2": "The authors have extensively addressed the comments in my previous review, and I think the revised manuscript is much stronger. I sfill think the ecological status clusters are hard to interpret in terms of funcfional effects of changes in distribufions due to climate change, but this is a limitafion for the field as a whole that they sufficiently address in the discussion.\nThe only thing I want to point out is there is a sentence fragment in describing the datasets used for training the random forest classifier at lines 538-540, \"Principal component analysis and correlafion metrics between features were used to ensure the validity of the randomly sampled dataset (Supplementary Figs. 14 and 15), which was then used as a training dataset in machine learning and used as a test dataset.\" Given the very high accuracy they should doublecheck they used an independent test set, it isn't clear from this descripfion.\nAside from that minor comment, I think this manuscript is much improved and I have no other concerns to raise before publicafion."}, "author_response": "Reviewer #1 (Remarks to the Author):\nIn this study, Qian et al., focus on a significant region of Earth addressing how global\nbiogeographic patterns of microbes are related to ecological functioning under\nscenarios of climate change. I want to congratulate authors for taking in to account all\nconcerns raised during the first round of review. Their article is now more solid and\nclear. This is fundamental since authors make predictions of biogeochemical function\nshift under different scenarios of Global Climate Change, hence, the verification of\nmetabolic pathways and methods was fundamental. The approach is very interesting,\nand using machine learning to bring together datasets from metagenomics, climate\nmodels, biogeography is clearly a fantastic tool with great potential.\nResponse: Thank you for your time and effort in handling our manuscript. Your\nprevious comments have greatly improved our work.\nReviewer #2 (Remarks to the Author):\nThe authors have extensively addressed the comments in my previous review, and I\nthink the revised manuscript is much stronger. I still think the ecological status clusters\nare hard to interpret in terms of functional effects of changes in distributions due to\nclimate change, but this is a limitation for the field as a whole that they sufficiently\naddress in the discussion.\nResponse: Thank you for your time and effort in handling our manuscript. We\nappreciate your previous comments, which greatly refined our work and made our story\nmuch more solid and clearer. In this study, we propose a new concept, ecological status,\nfor ocean biogeography by integrating diversity, biogeochemical potential, and\ndominant microbial taxa. Ecological status can be recognized as a more convenient and\ncomprehensive index for evaluating the influence of anthropogenic activities on\nmicrobial communities in ocean systems.\nThe only thing I want to point out is there is a sentence fragment in describing the\ndatasets used for training the random forest classifier at lines 538-540, \"Principal\ncomponent analysis and correlation metrics between features were used to ensure the\nvalidity of the randomly sampled dataset (Supplementary Figs. 14 and 15), which was\nthen used as a training dataset in machine learning and used as a test dataset.\" Given\nthe very high accuracy they should doublecheck they used an independent test set, it\nisn't clear from this description.\nResponse: Thank you for pointing this out. We revised this sentence to “Principal\ncomponent analysis and correlation metrics between features were used to ensure the\nvalidity of the randomly sampled dataset (Supplementary Figs. 14 and 15), which was\nthen used as a training dataset in machine learning. The remaining dataset after random\nsampling was used as an independent test dataset.” (lines 538-540)\nAside from that minor comment, I think this manuscript is much improved and I have\nno other concerns to raise before publication.\nResponse:Thank you again for your valuable comments."}]

Decision letter