Functional Redundancy in Ocean Microbiomes Controls Trait Stability

Theory predicts that species extinctions within communities are less likely to cause the loss of specific functions when functional redundancy is high. Few experiments have empirically tested this prediction, however, especially in the context of microbial communities and at the landscape scale. In part, the lack of metrics to quantify functional redundancy in microbial ecosystems has prevented addressing this question. In a companion manuscript we proposed a quantitative metric for functional redundancy called Contribution Evenness (CE) that is optimized to reflect the sensitivity of microbially-catalyzed ecosystem functions to local species loss. Here, we use CE to predict the stability of marine microbial functions to species and transcription loss. Using transcriptomes deposited in the Ocean Microbial Reference Gene Catalog (OM-RGC.v2), a catalog of genes and transcripts sequenced by the TARA Ocean expedition, we quantified the functional redundancy for 4,314 KEGG Orthologs (KOs) across 124 marine sites. Functional redundancy was highly correlated with a latent variable consisting of four main ocean physiochemical parameters: oxygen and chlorophyll a concentrations, depth, and salinity. Functional redundancy was systematically higher at the poles than in non-polar regions. Simultaneously, regional β-diversity for individual functions was higher for functions with higher functional redundancy These observations provide evidence that higher functional redundancy indicates increased stability of microbial ecosystem functions on spatiotemporal scales consistent with surface ocean mixing. We suggest that future changes in ocean physiochemistry could likely influence this stability for functions with lower functional redundancy.

Introduction: 27 The ability to resolve complete genomes in microbial communities via contig binning and 28 single-cell genome sequencing has revealed that metabolic functions, or traits, are often present 29 in multiple taxa within a microbiome (1-4). The practical consequences of functional 30 redundancy, in terms of community composition and function, are not well understood for 31 microbial communities. Theoretical predictions, developed in the context of macroecosystems, 32 suggest that functional redundancy can buffer a community against function loss due to species notably climate change, may influence microbially-mediated biogeochemistry stability. 45 In a companion paper, we developed a metric for functional redundancy called 46 contribution evenness (CE) (14). CE uses genes or transcripts as traits (15,16) and measures 47 how evenly taxa contribute to the community-wide level of these traits (total gene abundance in 48 the case of genes, and total transcript abundance in the case of transcripts). CE is optimized to 49 reflect how stable microbially-catalyzed functions are to random species extinction or 50 transcription loss in communities. Here, we use CE to address three fundamental questions about 51 functional redundancy in marine microbial metacommunities: 52 1.) Does functional redundancy vary in ocean microbiomes, 53 2.) If so, does this variation relate to the physicochemical environment, and 54 3.) What ecological consequences are apparent from observed patterns in functional redundancy?

55
To answer these questions, we used KEGG orthology (KO) annotations of genes reported 56 from 180 sites sampled during the TARA Oceans expeditions, plus transcripts reported from a 57 subset of 124 of those sites (17). KEGG is an attractive ontology because the functions of KOs 58 are biochemically well-characterized and provide a reasonably large coverage of predicted genes 59 in marine microbiomes. The TARA Oceans dataset is ideal for this work because it encompasses 60 a diverse set of sites and marine ecosystems, which were sampled with generally standardized 61 methods (18). By applying CE to the TARA Oceans data set, we identified geographic variations 62 in functional redundancy, identified physicochemical factors relating to this variation, and 63 established that microbial functions were less stable in regions with lower functional 64 redundancy. These observations combined allowed us to identify some possible effects that 65 climate change may have on the resilience of microbial community function to local disturbance.  We sought to characterize spatial trends in the functional structure of microbial 100 communities. To do this, we performed redundancy analysis using the rda function from the R 101 package, vegan (21). In this analysis, the functional structure of microbial communities was 102 treated as a response variable and physio-chemical variables were treated as predictor variables.

103
For sites where a given KO was absent, the KO was assigned a value of zero. This situation 104 occurred ~13.6% and 18.5% of the time for gene and transcript profiles, respectively. Low 105 functional redundancy does not necessarily imply that a function is ecologically unimportant.

106
This is particularly true for traits with low functional redundancy that strongly correlate to depth (m), and silicate (mmol m -3 ). We chose these variables as they were in situ measurements 118 (18) and could scale with models that predict global ocean chemistry. We imputed missing data, 119 as removing incomplete cases can bias datasets (23). Random forest imputation was performed 120 using the missForest function, from the R package, missForest (24). Next, each environmental 121 variable was converted into a normal distribution using a boxcox transformation via the boxcox 122 function in the R package, MASS (34). Variables were then centered to have a mean of zero.

123
Last, the significance of each variable as well as the first two canonical axes was verified using 124 the anova.cca function in the R Package, vegan (21).

125
The first canonical axis derived from the redundancy analysis was used to predict mean 126 transcript functional redundancy from the individual metatranscriptomes (sites) using OLS 127 regression. This model was compared to an OLS regression using salinity, nitrate, phosphate, 128 oxygen, chlorophyll A, depth, and silicate as predictors of mean transcript functional 129 redundancy. Similar to before, missing data were imputed using missForest (24), variables were 130 transformed with a boxcox transformation (25), and variables centered so the mean distribution 131 was zero prior to regression. The best predictor subset was determined using the regsubsets 132 function from the R package, leaps (26). The criteria defining the best model was minimizing Here, we used a diversity order q=0, or the richness to calculate diversity. As such, γ-diversity 154 was either 0 (regionally absent) or 1 (regionally present), and α-diversity was the proportion of 155 local communities with the trait. Since β-diversity is not defined for traits completely absent 156 from a region, the β-diversity calculation simplified to: which is simply the inverse proportion (f) of local communities in a region with a given function. 160 We used CE (14) to measure functional redundancy for KEGG ortholog (KO)-annotated 161 genes from 180 and 124 TARA Oceans metagenomes and metatranscriptomes, respectively (17).

162
CE uses genes (CE with respect to genes) or transcripts (CE respect to transcripts) as traits (15,163 16) and measures how evenly species in a community contribute to these traits as a whole. After 164 calculating CE for metagenomes (Fig. 1A) and metatranscriptomes (Fig. 1B), we found that KOs

173
Although the nature of these observations is similar, the ecological implication is different. Earth's oceans using the derived coefficient from our OLS regression of mean functional 205 redundancy versus the first canonical axis of the RDA. Our model predicts that functional 206 redundancy is highest in polar regions and near river outflows, and lowest in subtropical gyres 207 (Fig. 3A). Variance was highest in polar regions, coastal regions, and river outflows (Fig. 3B).

208
The transition from high to low functional redundancy between polar and non-polar latitudes was 209 consistent with previously reported ecological boundaries for ocean microbiomes, where the 210 transition from non-polar to polar latitudes corresponded to compositional changes in 211 metatranscriptomes and metagenomes (17).

213
In the companion paper, we presented a numerical simulation demonstrating a 214 relationship between CE and trait stability to random extinctions. Here, we sought to establish 215 whether there is empirical evidence that functional redundancy increases trait stability in ocean 216 microbiomes, as previously hypothesized (31). Specifically, we hypothesized that if higher levels 217 of functional redundancy do indeed increase trait stability, then the sets of functions present in 218 nearby communities should be more similar when functional redundancy is higher. To test this 219 hypothesis, we compared functional redundancy at each epipelagic site to the β-diversity (27) of 220 each function (coded as present or absent) in a local region consisting of the site in question plus 221 the nine nearest sites. In this context, lower β-diversity indicates that a larger portion of local 222 communities shared a given function (equation 3). Since microbial communities exhibit spatial 223 autocorrelation (32), these ad-hoc local regions serve as a reasonable alternative to the ideal data 224 set, which would consist of subsequent metagenomes taken from a parcel of water over time. 225 We observed a negative monotonic relationship between functional redundancy and beta 226 diversity for every KO for every unique region (Fig. 4). This qualitative evidence is consistent 227 with the idea that functional redundancy maintains trait stability on a regional level. Notably, 228 CE>0.01 corresponded to functions becoming almost entirely stable within their respective 229 regions (Fig. 4). This threshold suggests a large portion of functions should be stable across 230 Earth's ocean microbiome (Fig. 3). Nonetheless, select functions do appear vulnerable to 231 regional instability. Notably, nitrogen metabolism functions tend to have CE<0.01 (Fig. 1B).

232
Climate models and historical observations indicate that major ocean currents are slowing (33), 233 oxygen minimum zones are expanding (34), and there will be future changes in regional primary 234 productivity (35). With respect to oxygen and primary productivity, these variables are important 235 in our factor model and suggest that future changes in these parameters might disproportionately 236 influence the stability for lower functional redundancy functions such as nitrogen energy 237 metabolism.