Why Environmental Biomarkers Work: Transcriptome-Proteome Correlations and Modeling of Multi-Stressor Experiments in the Marine Bacterium Trichodesmium

Culturing, Molecular Extractions, Sequencing and Proteomic Analyses

A depiction of the prior experimental design can be found in Fig. S1a and in Walworth et al. 2016 ²³ in addition to physiological and proteomic results. Detailed culturing²³, protein extraction²⁷, RNA extraction³¹, and analysis methods can be found in our previously published papers. Briefly, proteomic data was generated using data-dependent acquisition on a Thermo Fusion mass spectrometer and mapped to the IMS101 genome using SEQUEST with Proteome Discoverer. 1908 Trichodesmium proteins were detected across the sample set using a 0.3% false discovery rate and requiring 2 peptides per protein using Scaffold 3.0 (Proteome Software) and normalized spectral counts were reported (see Supplemental Dataset 1) and matched to transcript abundances. Spectral counts provided a robust relative abundance of the proteome, although future efforts could focus on peptide precursor or fragment ion peak areas which have broader dynamic range, although the former is subject to interference from sequence variants in metaproteomic samples. Transcript data was generated using Illumina Hi-Seq after rRNA removal yielding single end 50bp read libraries (see Supplemental Dataset 1).

Differential expression analysis

Differential expression data was from previously described datasets ²⁵.

Multivariate and Pairwise Analyses

Redundancy analysis was conducted on TMM-normalized transcript levels and nutrient concentrations using the vegan package ³². RDA ordinations were compared via the Procrustes test with the “protest” function in vegan using default settings.

Nonmetric multidimensional scaling (NMDS) was conducted using the “metaMDS” function from the vegan package with default settings (Bray-Curtis dissimilarities computed using Wisconsin double standardization) except for “autotransform=FALSE”. Hierarchical clustering with multiscale bootstrap resampling was conducted on Bray-Curtis dissimilarities calculated from TMM-normalized transcript levels using the pvclust package ³³ with the following settings: method.dist="euclidean", method.hclust="average", nboot=1000.

The “cor.test” function in “R” was used to calculate Pearson correlation coefficients on log2 fold changes (relative to the r380 condition; see below) of average TMM-normalized transcript levels vs log2 fold changes of average normalized spectral counts of the corresponding protein product detected in Walworth et al. 2016 ²³. Significance tests between mean correlation coefficients were carried out via the “perm” package³⁴ using a two-sample permutation test with Monte Carlo simulation and 2000 permutations.

Gene Ontology (GO) enrichment analysis

As previously described in Walworth et al. 2016 ³¹, Gene Ontology (GO) mappings for Trichodesmium were downloaded from the Genome2D web server (http://genome2d.molgenrug.nl). The “phyper” function in “R” (R Core Team 2014) was used to test for significant enrichment of GO pathways and p-values were corrected with the Benjamin and Hochberg method ³⁵ using the “p.adjust” function (p ≤ 0.05). Finally, genes in enriched GO categories were manually checked.

RNA-Protein Model

A model simulating the biosynthesis of RNA transcripts and proteins was created in MATLAB (MathWorks) using parameters listed in Table 1. One-minute time steps over a duration of 3 days (4320 min) within a population of 100 cells was modeled. Expression was modeled with a signal parameter that varied between 0 and 1 representing the gene being turned on or off with a simple linear ramping function. Basal signal level (s₀) was set to 0.01, representing 1% of the population or 1 cell within 100 cells with gene expression to avoid zero denominator in fold-change estimates. In an alternate scenario representing metabolic functions in continuous use, a basal expression of 20% (0.2) was used. Rate constants for production of RNA and protein (k_RNA, k_protein), as well as their decay (k_-RNA, k_-protein) were obtained from experimental measurements in Escherichia coli and Saccharomyces cerevisiae where the degradation of each type of biomolecule was inhibited and measured ^36,37. Change in the inventories of RNA and protein biomolecules was estimated using simple production and loss terms (Equations 1 and 2):

Turnover time estimates of these biomolecules are, to our knowledge, unavailable for marine organisms. While they may be different, presumably slower based on observed growth rates, it seems a reasonable assumption that the large difference between RNA and protein decay would be maintained. Hence the trends observed here are likely to be consistent or further enhanced, particularly if the RNA degradation activity in marine microbes is slower than that of E. coli as might be expected.

This model has similarities to that of Moran et al. 2013 ³⁸, but is applied here to focus on any coordination between RNA and protein and how results affect fold change estimates. Moran et al. instead explored random triggers for transcript production using 10 random events per day of 1 transcript per event, and then simulated environmental stimuli by increasing to 3 transcripts per event. In contrast, the model in the present study estimates production of RNA and protein molecules based on their upstream signal (signal, transcripts) as described in equations 1-2, rather than assigning a constant number of seven protein copies per transcript. Moran et al. observed the total number of mRNA molecules in a typical marine bacterium was only ~200 per cell³⁸. These observations were incorporated into the model by allowing k_RNA to be 100 to enable signal input to be transcribed into 100 transcripts across a population of 100 cells, and then tuning the protein production k_Protein parameter to produce 6-7 proteins per transcript also based on Moran et al., 2013. Similarly Moran et al. used a RNA half-life of 1.5 minutes, whereas this model uses a more conservative RNA half-life of 6.8 minutes based on a tiling microarray study in E. coli³⁶.

Global comparison of CO₂-impacted, nutrient-limited molecular dynamics

To investigate shifts in molecular metabolism underlying Fe-limited, P-limited, and Fe/P co-limited metabolisms as they interact with predicted future CO₂ atmospheric conditions, we examined the global transcriptional output in each scenario from Walworth et al. 2018 ⁸ and their corresponding translational correlations with proteins identified in Walworth et al. 2016 ²³. Across all treatments, we detected transcription of ≥96% and translation of 37% (see below) of the IMG annotated genes. Instances where comparisons between transcripts and proteins were not possible due to one of those biomolecules not being detected were excluded. While not the goal of our study, future studies could apply deeper proteomic methods to explore correlative relationships in rarer biomolecules. Initial examination of the dataset in log2 fold change space revealed a cluster of points around the origin with a subset of data extending into the northeast and southwest quadrants implying a co-varying relationship between some transcript-protein pairs (Fig. 1). It is important to note that non-translated mRNAs can also potentially represent ‘irrelevant expression’, if these transcripts are not regulatory. Thus, there is much to be learned from future studies reconciling the constantly developing techniques of RNA and protein based ‘omics, with relatively few studies to date comparing multi-condition proteomes with concomitantly taken transcriptomes.

To this end, we first conducted exploratory analysis on the global transcriptional output without constraining for environmental variables (see below) in order to observe how replicates clustered without imparting treatment information (i.e. no assumptions). We used both nonmetric multidimensional scaling (NMDS) and hierarchical clustering with multiscale bootstrap resampling (replicates = 1,000) ⁴⁰. NMDS revealed consistent nutrient-limited patterns across replicates (Fig. S1b), which was further supported by hierarchical clustering (Fig. S1c) that generated two high-confidence clusters with AU-bootstrap p-values > 0.95 (i.e. rejection of the null hypothesis that the clusters do not exist at the 0.05 significance level). We labeled samples based on their environmental condition at the time of sampling with CO₂ concentration followed by the limiting nutrient. For example, 380-Fe signifies low-CO₂ and low iron. For replete conditions, we chose to put an “r” in front of the CO₂ concentration (i.e. r380 and r750). In these unconstrained analyses, the 380-Fe condition demonstrated greater inter-replicate transcriptional variation, although both replicates were indeed Fe-limited as evidenced by both reduced growth ²³ and several upregulated Fe-stress genes (Fig. S1b; see below and Walworth et al. 2018 ⁸). Interestingly, Fe-limited treatments (380- and 750-Fe) clustered with the nutrient replete ones (r380 and r750) supporting the notion that P-limitation evokes a more varied molecular response than Fe-limitation relative to replete nutrient regimes (Fig. S1c). This trend was further corroborated through hierarchical clustering of transcript levels of all differentially expressed (DE) genes (i.e., those DE in any treatment relative to r380; n = 1943; Fig. 2a). To better resolve sources of variation producing these groupings, we further clustered upregulated (n = 951; Fig. 2b) and downregulated (n = 1074; Fig. 1c) transcripts, which yielded strikingly different trends among Fe-limited and replete conditions. Upregulated genes exhibited significantly different (bootstrap confidence level > 0.95) clusters among Fe-limited and CO₂ regimes relative to those produced by clustering all DE genes while those of downregulated genes maintained similar groupings to those of DE genes. These patterns suggest that Fe-limited and replete clusters produced from all DE genes (Fig. 2a) may have been driven more so by downregulated transcriptional variation of the DE gene pool. Hence, mechanistic adaptations to high CO₂ and responses to Fe-limitation both seem to be more similar to each other than to molecular responses under P-limitation. Interestingly, both up- and downregulated P-limited clusters were statistically analogous to those produced by clustering all DE genes (Fig. 2a) and total genes (Fig. S1c), which further corroborates the notion that P-limitation evokes a highly conserved yet more varied transcriptional rearrangement than either high CO₂ or Fe-limitation.

Next, we employed constrained multivariate methods that attempts to explain variation based on the given environmental variables (e.g. RDA analysis) for both transcripts and proteins. When relating nutrient-limited proteome variation identified via redundancy analysis (RDA) in Walworth et al. 2016 ²³ (n = 1248; 37% of genome coding potential) to that of their corresponding gene transcripts (corresponding-RNA-RDA; Fig. 2d), similar nutrient-limited patterns across replicates were observed with most of the proportional variance explained by the environmental variables for both transcripts (adjusted R² = 0.56) and proteins (adjusted R² = 0.59), respectively. Moreover, a statistically significant correlation was observed (r = 0.75; P < 0.001) when comparing the corresponding-RNA-RDA and prot-RDA ordinations via permutational, least-squares orthogonal analysis (i.e., Procrustes analysis; Methods) indicating that, overall relative transcriptional variation among treatments was captured in the corresponding proteome variation (Fig. 2e).

Performing the same RDA analysis on the global transcriptional output (global-RNA-RDA; Fig. 2f; adjusted R² = 0.50; n = 4862) displayed a strikingly consistent ordination relative to that of both the prot-RDA (represented in Figure 1c in Walworth et al. 2016²³ and corresponding-RNA-RDA (Fig. 2d). In fact, when analyzing congruence of the global-RNA-RDA to both the prot-RDA ²³ and the corresponding-RNA-RDA, respectively, a robust, statistically significant correlation was observed relative to both the prot-RDA (r = 0.75, P < 0.001) and the correlated-RNA-RDA (r = 0.98; P < 0.001). Taken together, these analyses indicate that both transcriptional and translational variation inherent in the overlapping 37% of the detected proteome and 96% of the detected transcriptome are not only representative of one another across nutrient-limited conditions, but also of global transcriptional variation. In other words, this congruence implicates that the detected proteome variation is indeed representative of overall treatment-specific, transcriptional variation in IMS101 for Fe, P, Fe/P, and CO₂.

Another striking observation is that 380-Fe/P replicates were more closely associated to P-limited treatments (380- and 750-P) while 750-Fe/P replicates were more correlated with Fe-limited ones (380- and 750-Fe) for both ordination-based (RDA, Fig 2d,f; NMDS, Fig. S1b) and clustering analyses (Fig. 2a,b,c and Fig. S1c). Hence, these trends imply that regardless of the amount of variation explained by the given environmental variables, the global transcriptional profile of Fe/P co-limited metabolism is more similar to that of P-limitation when adapted to current CO₂ levels (380-Fe/P) but is more similar to Fe-limitation when adapted to high CO₂ (750-Fe/P).

From global to pairwise mRNA-protein correlations

Bacterial messenger RNA (mRNA) typically have much shorter half-lives than those of proteins and can be rapidly transcribed in response to quickly changing environmental conditions ⁴¹. For example, in an E. coli cell, many mRNAs are typically degraded within minutes, whereas their corresponding proteins can have longer lifetimes than the cell cycle ⁴². These data indicate that intracellular mRNA copy number at any instant may typically reflect recent transcriptional activity, while protein levels in the same instant can represent a relatively longer history of accumulated expression. Indeed, simultaneous measurements of mRNA and their corresponding protein concentrations have been shown to be uncorrelated within a single cell at a single time point, while mRNA levels integrated over many cells in a population generally correlate with protein concentrations ⁴². Different environmental variables invoking specific metabolic pathways can also impact overall mRNA-protein correlations ^1,2.

Unlike prior other microbial studies investigating mRNA-protein correlations across one or two treatments ², we examined correlations across eight different growth conditions. As in most organisms to date, IMS101 transcript levels modestly correlated to their corresponding protein abundances across all treatments (n = 1,246; squared Pearson correlation coefficient, r² = 0.36), suggesting that other forms of regulation might need to be invoked to explain the majority of variation ¹. Upon calculating Pearson correlation coefficients across all treatments for genes that were differentially expressed (DE) in at least one experimental condition relative to r380 versus those that were not called as differentially expressed (NDE) in any condition, we observed a significantly greater mean correlation coefficient for DE (r² = 0.35) vs NDE (r = 0.03) genes (p < 0.0005; two-sample permutation test with Monte Carlo simulation). This same robust trend was also observed when comparing average correlation coefficients for DE vs. NDE genes for each pairwise treatment comparison (Fig. S2). Hence, this greater correlation suggests that transcription and translation of environmentally responsive genes and their protein products may be more closely coupled on average than those that do not respond to changing environmental conditions.

Pathway-specific molecular relationships

To search for metabolic pathways harboring highly correlated transcript and protein abundances, we calculated Pearson’s correlation coefficients (PC) and kept genes retaining PC values ≥ 0.7 and p-value < 0.05. We then mapped genes onto their Gene Ontology (GO) pathways and searched for significant GO enrichment using the hypergeometric test with Benjamini-Hochberg correction (FDR ≤ 0.05; Dataset S1). Enriched GO pathways included translation, nitrogen fixation, oxidation-reduction, calcium ion binding, proteolysis, iron-sulfur cluster binding, and photosynthetic electron transport/light reactions. Hence, these data suggest that biosynthesis of these protein products increases with increasing transcript accumulation. Therefore, their abundances may be more representative of instantaneous cellular activity rather than those with low PC values, in which the latter may be a product of constitutive activity and/or reduced protein degradation. For example, well-characterized Fe- and P-stress genes exhibited high PC values in conjunction with genes involved in photosynthetic light reactions (Fig. 3). In all cases, the interaction of high CO₂ with limiting Fe irrespective of fluctuating P (750-Fe and 750-Fe/P) yields increases in both transcript and protein levels for all Fe-stress genes relative to their low CO₂/low Fe conditions (i.e. 380-Fe and 380-Fe/P). Hence, the expression of these Fe-stress biomarkers increases under limiting Fe in future ocean conditions suggesting their profiles to be robust for in situ profiling. The genes involved in photosynthesis light reactions (Fig. 3c) exhibit parallel trends, and could be related to mid-day sampling and increased transcription and translation of photosynthetic genes during the light period (Held et al., in press).. Interestingly, P-stress transcript and protein levels robustly respond to limiting P irrespective of Fe and CO₂, but increased CO₂ has little to no effect on P-stress molecular machinery (Fig. 3b). Additionally, these data suggest that different subunits of the same phosphate transport system (i.e. Pst) could respond differently to fluctuating P as transcript and protein levels of, for example, the PstB subunit increased more considerably than those of PstS under limiting P (Fig. 3b).

Conversely, GO-pathways specific to the NDE pool with low PC values involve general biosynthetic processes, magnesium ion binding, proteolysis, respiration (e.g., TCA cycle), nucleotide metabolism, inosine monophosphate (IMP) biosynthesis, and amino acid metabolism (Fig. 4). Of note, decoupling of TCA gene expression and enzyme abundance in yeast has also been previously observed ². Interestingly, the GO pathway assigned to photosystem I reaction center (GO: 0003989) experience significant enrichment of genes with low PC values (Dataset S1) while conversely, the overall photosynthesis light reaction GO pathway (GO:0019684) was enriched with genes with high PC values (Fig. 3c). Hence, it is important to distinguish between potentially constitutive vs. responsive components of large, multifaceted pathways with many components that can have widely differing roles (e.g., energy flow vs. structure). Overall, it seems the persistent activity of certain core pathways to maintain basic cellular homeostasis, irrespective of changing environmental conditions, may result in decoupling of transcription and translation on average (Fig. 4), relative to pathways that synthesize new transcripts to generate protein products specific to particular environmental triggers (e.g., Fig. 3).

RNA-Protein Modeling: Influence of biomolecule inventories on fold-change dynamics

In order to explore the potential causality of RNA-protein coupling and decoupling we generated a model for production and decay of both biomolecules (Fig. 5). Model parameters include previously published values for RNA and protein degradation based on studies of whole transcriptome and proteome degradation measurements from model organisms (E. coli ³⁶ and yeast ³⁷, respectively) as well as delays between RNA and protein expression observed in Prochlorococcus ⁴ (Table 1). Notably, RNA decay rates are approximately two orders of magnitude higher (~60-275 fold) than that of protein degradation. While differences in degradation rates occur between different specific transcripts and protein molecules, there is clearly a large difference between transcript and protein degradation rates, as expected based on the transient information transmission role of RNA versus the long-term functional roles of proteins. The choice of E. coli and yeast as organisms for obtaining decay parameters will likely vary from that of various marine microbes, with both degradation processes likely being slower within slower growing marine microbes such as Trichodesmium. Moreover, mis-folded and mis-metallated proteins are also known to be targeted for degradation by intracellular recycling processes. However, the large difference in cellular decay rates of RNA and protein rates is almost certainly a universal phenomenon, and hence we expect the general trends of temporal overlap between transcripts and protein inventories modeled here to be relevant even if future taxon-specific activities vary somewhat from those used here.

In particular, if the RNA degradation constant is lower in slower growing organisms due to lower abundance and/or efficiency of RNases, an increased temporal overlap between transcripts and proteins would result in greater opportunity for correlation between the two types of biomolecules. Additional parameters in the model include production of RNA and protein, as well as constants for temporal offsets between signal transduction and transcription (d_RNA), and transcription and translation processes (d_protein). For the latter, 3h was used, similar to the range observed in Prochlorococcus ⁴. Finally, the model includes initial condition parameters whose influence can be explored including basal expression levels (1% and 20%) and time settings for onset, duration, and decay of signals to initiate transcription. The notion of the transcriptional signal is generalized here, whether they occur from riboswitches, two component regulatory systems, sigma factors or other mechanisms. Moreover, the onset and decay of the regulatory signaling is also idealized here as a linear ramp (using a “long slow ramp” that repeats daily), although it could be further qualified to have positive feedbacks to increase sensory capabilities such as observed in increased expression of two component regulatory systems, e.g., phosphate regulatory systems ⁴³. The influences of post-translational modifications (PTMs) were purposely not included in this model in order to allow an examination of potential range of differences in RNA and protein fold-change expression levels prior to often invoked additional PTM regulatory controls. While PTMs likely contribute to protein regulation and enzyme activity in Trichodesmium and other marine microbes, those effects would likely not supercede the protein inventory dynamics measured and modeled here.

The RNA-Protein expression model results imply that changes in basal expression level can have a major influence on the extent of RNA-proteome coupling and decoupling observed experimentally. For example, when there is no significant basal expression (1% or s₀=0.01; where a non-zero value was chosen to allow fold-change ratio estimate) as would be expected prior to an environmental stimulus such as phosphate or iron stress, both RNA and protein signals respond with a large increase in inventory and fold-change (Fig. 6 and 7 left panels). In contrast, when there was a 20% basal level of expression (s₀=0.2), as expected for routinely used metabolic machinery, a significant inventory of protein product accumulates (Fig. 6 right panels), causing the fold change signal to become muted (Fig. 6). These results are also displayed as a movie format (Supplemental Movie M1) to display the temporal unfolding of these transcript-protein inventories and their influence on fold changes results. The fold change decrease is accentuated if calculated each day as the protein inventory accumulates (e.g., the choice of what to normalize expression to in the fold change estimate). Notably, proteins that respond to environmental stress are consistent with the lower basal expression model that yields correlations between RNA and protein.

In Trichodesmium, photosynthetic transcripts and proteins were also observed to be correlated (Fig 4c). While generally considered metabolic proteins with high inventories, this observed correlation can be particularly driven by iron limitation scenarios, which can remodel overall machinery by reducing the iron-rich photosystem I (PSI) in favor of PSII ²³. Iron stress transcripts and proteins like isiA can also concurrently increase in abundance solely in times of enhanced iron stress to form PSI-IsiA-chlorophyll–protein–antenna super-complexes ⁴⁴ while in other conditions like phosphorus stress, they show little correlation due to the lack of an environmental trigger. Additionally, diazotrophs have been shown to have strong diel cycles that affect both transcripts and proteins of photosynthetic and nitrogen fixation enzymes in order to promote intracellular oxygen protection and iron conservation efforts ^45,46. Hence, the choice of when to sample transcriptomic and proteomic experiments can also be consequential. The short-reported half-lives of RNA molecules cause a rapid decrease in transcripts if they are not continually being produced ³⁸. The expected large two order of magnitude difference in half lives of RNA and protein molecules results in a counterclockwise progression on RNA and protein space plots commonly used for assessing coherence between these biomolecules (Fig. 8). Choice of sampling time is typically estimated to allow sufficient time for RNA and protein biosynthesis, where if too early there is insufficient time for protein synthesis, and if too late RNA transcripts may have begun degradation ⁷. Notably, the magnitude of these counterclockwise traces is affected by the basal expression level, where low basal expression creates an easily observable pattern due to the high fold change values for transcripts and proteins, while high basal expression results in the circular trace to be compressed near the origin. This is indeed what many experimental datasets look like with metabolic proteins found near the origin in log2 transcript-protein space (Fig. 1), while genes that respond strongly to environmental cues have stronger signals that project into the northeast quadrant of these plots. The model results here may be overly conservative, as it seems likely that marine microbes with much slower reproduction rates than E. coli or yeast would also have lower concentrations of RNases and concurrent longer residence times of RNA, as mentioned above. If lower RNase activity were common among marine microbes like Trichodesmium, the overlapping time window for co-existence of RNA and proteins would be enhanced.

Based on this simple modeling approach, the (often negative) conclusions about the coherence of RNA and protein expression in microbes appear to have been clouded by use of relative units of fold change that discriminate against metabolic proteins with larger inventories. In these cases, there may be highly abundant proteins that change in inventory by 20% for example, which would be a large change in copy number, but a small change in fold change ratios. Use of absolute units for both transcripts and proteins and efforts to calibrate with standards to achieve copy number per cell estimates would resolve this problem. Importantly however for environmental applications, the successful use of environmentally responsive genes in both RNA and protein space can be explained by their large signals relative to basal expression that allows easy detection of large fold changes. In short, the lack of observed correlations of metabolic systems represented in transcripts and proteins can be explained, at least partially, by an over-reliance on ratios of data, and those observations should not detract from the statistically significant correlations of environmentally response genes in RNA-protein space and their application as research and monitoring tools in environmental microbiology and biogeochemistry.

While each experimentalist can determine which biomolecule(s) (proteins and/or transcripts) would complement their research program, these matters often default to practical decisions based on prior training or the availability of sequencing and mass spectrometry facilities. The results and model presented here could help to guide experimental design, and perhaps inspire future multi-omic studies. In circumstances of meta’omics where environmental samples containing highly diverse biological communities are measured, the measurement of transcripts is particularly useful when environmental Eukaryotic organisms with large genomes, as the transcripts allow for the generation of a sequence database without introns for peptide-to-spectrum mapping in data dependent acquisition approaches.

Together these empirical observations and model results validate the utility and provide a theoretical basis behind proteins that respond to environmental stimuli in marine microbes. Because organisms that comprise natural ecosystems are frequently limited in their productivity by nutritional elements, they often have evolved specific adaptive responses to nutritional scarcity. As such, these proteins and their mRNA transcripts represent ideal environmental biomarkers to detect controls on ecosystem productivity. Biomarkers for nutritional stress in marine phytoplankton and bacteria that have been identified are shown in Table 2, including those for nitrogen ^43,48-50, phosphorus ^5-6,51-53, iron ^{18-19, 54-61}, vitamin B₁₂^62,63, and zinc⁶⁴. Notably, proteins involved in adaptive responses in marine microbes are often distinct from those model organisms because of their very different environmental conditions. As a result, several of these marine biomarkers have only recently been discovered and it is almost certain that more remain to be discovered. In some cases, biomarkers for required elements such as cobalt do not appear to have adaptive responses (or their chemical species they respond to has not yet been identified)⁶⁵ and corresponding biomarkers perhaps due to the nutritional burden of such systems in highly streamlined genomes. With the natural environment undergoing rapid and unprecedented changes, the ability to observe and understand how the oceans are changing on a global scale is possible through the application of environmental biomarkers. Demonstration that RNA and protein do correlate for environmentally responsive genes and the theory behind when and why these correlations occur is valuable in increasing confidence for study of metaproteomic and other meta-‘omic studies to the global oceans ^66-70.