Environment-driven shifts in interindividual variation and phenotypic integration within subnetworks of the mussel transcriptome and proteome.

The environment can alter the magnitude of phenotypic variation among individuals, potentially influencing evolutionary trajectories. However, environmental influences on variation are complex and remain understudied. Populations in heterogeneous environments might exhibit more variation, the amount of variation could differ between benign and stressful conditions, and/or variation might manifest in different ways among stages of the gene-to-protein expression cascade or among physiological functions. Here, we explore these three issues by quantifying patterns of inter-individual variation in both transcript and protein expression levels among California mussels, Mytilus californianus Conrad. Mussels were exposed to five ecologically relevant treatments that varied in the mean and interindividual heterogeneity of body temperature. To target a diverse set of physiological functions, we assessed variation within 19 expression subnetworks, including canonical stress-response pathways and empirically derived coexpression clusters that represent a diffuse set of cellular processes. Variation in expression was particularly pronounced in the treatments with high mean and heterogeneous body temperatures. However, with few exceptions, environment-dependent shifts of variation in the transcriptome were not reflected in the proteome. A metric of phenotypic integration provided evidence for a greater degree of constraint on relative expression levels (i.e., stronger correlation) within expression subnetworks in benign, homogeneous environments. Our results suggest that environments that are more stressful on average - and which also tend to be more heterogeneous - can relax these expression constraints and reduce phenotypic integration within biochemical subnetworks. Context-dependent "unmasking" of functional variation may contribute to interindividual differences in physiological phenotype and performance in stressful environments.

INTRODUCTION

Interindividual phenotypic variation within a population is an important substrate for evolutionary processes, because it is typically underlaid by – or at least correlated with – genetic variation (Fusco & Minelli, 2010; Haldane, 1937; Hartl et al., 1997). However, this correlation between phenotypic and genetic variation can degrade under certain circumstances, such as during rapid environmental changes (Badyaev, 2005; Gerken et al., 2015; López‐Maury et al., 2008). Even within a single generation (i.e., before evolutionary processes can act), physiological plasticity expressed differently among individuals can change the apparent correlation between genetic variation and variation in phenotype (Crawford & Oleksiak, 2007; Gibert et al., 2019). Such plasticity‐driven shifts in phenotypic variation may alter the effects of selection, thereby influencing predicted responses to environmental change (Ghalambor et al., 2007; Guscelli et al., 2019; Oleksiak & Crawford, 2012). However, environment‐driven shifts in phenotypic variation remain understudied.

Here, we address three factors that might influence the magnitude and nature of within‐generation, interindividual variation (hereafter, simply “variation”) in molecular physiological phenotypes. First, variation may be elicited in habitats that are environmentally heterogeneous. This heterogeneity often manifests both in space (among individuals at any one time) and in time (individuals experience temporal fluctuations) (Dowd et al., 2015). For example, differences in shading, wave exposure, and community composition within the rocky intertidal zone drive heterogeneity in body temperatures of intertidal organisms over small to large spatial scales (Denny et al., 2011; Helmuth & Hofmann, 2001). Forest canopies are similarly heterogeneous, with very small‐scale variation in light availability, temperature, and nutrient availability (Lechowicz & Bell, 1991). Individuals within these environments acclimatize to microclimatic conditions that may differ substantially from those of their neighbors (Gleason et al., 2017). These processes can contribute to interindividual differences in phenotype (Mackay, 1980; Pincebourde et al., 2016). However, the effects of environmental heterogeneity on variation in underlying components of the physiological phenotype, such as patterns of transcript and protein expression, are largely unknown.

Second, the magnitude of physiological phenotypic variation might differ between environmental contexts (see figure 2c in Stearns & Kawecki, 1994), such as those that are benign or stressful on average (Tanner & Dowd, 2019). We define environmental stress as the infliction of macromolecular damage in the cells of an organism by environmental changes (Kültz, 2005). On one hand, environmental stress may be expected to mask phenotypic variation, resulting in a convergent “damage‐control” molecular phenotype (Kültz, 2005; Oleksiak & Crawford, 2012). This can be thought of as analogous to canalization of development in the face of genetic or environmental disturbances (e.g., Stearns & Kawecki, 1994). Alternatively, stress may be expected to unmask cryptic phenotypic variation, resulting in a diversity of responses to challenging environmental conditions that are then subject to selection (Ghalambor et al., 2015; Rohner et al., 2013). Physiologists have collected insufficient data to establish consensus for the relative importance of these competing masking and unmasking hypotheses (Gibert et al., 2019; Tanner & Dowd, 2019).

Third, the effect of stressful conditions on the magnitude of phenotypic variation also likely varies among physiological functions, levels of biological organization, and performance traits (Badyaev, 2005; Kingsolver & Buckley, 2017; Tanner & Dowd, 2019). For example, Tanner and Dowd (2019) found contrasting patterns of variation for different traits, including between global‐scale variation in the transcriptome and the proteome, within groups of intertidal mussels exposed to benign or thermally stressful treatments. More generally, traits fall along a spectrum of phenotypic variation. At one end, selection tends to canalize integrative traits with strong links to fitness, such as development time (Stearns & Kawecki, 1994) or genomic “modules” that generate morphological traits (Wagner & Altenberg, 1996). In contrast, other traits might be less constrained. For example, transcript levels and metabolic rate are reported to be highly variable among individuals experiencing similar or even identical environmental contexts (Lin et al., 2016; Steyermark et al., 2005). One plausible explanation for these discrepancies is that buffering mechanisms filter out the underlying “noise” to generate similar integrative phenotypes. In the context of the gene expression cascade – the series of events from initiation of transcription through protein translation – previous studies suggest that interindividual variation in transcript levels exceeds variation in levels of the corresponding proteins on a global scale (Hansen et al., 2018; Liu et al., 2016; Tanner & Dowd, 2019). There are several intermediate mechanisms (e.g., mRNA processing and turnover, preferential translation, or time lags between transcription and translation) that may lead to divergent variance patterns between transcriptome and proteome (e.g., Lackner et al., 2012; Liu et al., 2016; Maheshri & O’Shea, 2007; Vogel & Marcotte, 2012). Regardless of the mechanisms, discrepancies between these different representations of the molecular physiological phenotype must be considered when exploring connections between physiology and performance or fitness outcomes (Feder & Walser, 2005). We are unaware of studies that have simultaneously analysed the transcriptome and the proteome from the same individuals with the explicit goal of understanding how variation at each stage shifts in different environmental contexts. Even within a single step of the gene expression cascade, it is also likely that different functional units (we refer to these units within the transcriptome or proteome as subnetworks) are more or less prone to exhibiting phenotypic variation, a feature that would be concealed in global analyses such as those in Tanner and Dowd (2019).

In this study, we attempt the complex task of combining these three lines of reasoning and extending them to analyses of functional subnetworks of the transcriptome and proteome. We use the dynamic rocky intertidal zone as a model system to explore these issues, quantifying variation in expression patterns of transcripts and their corresponding proteins. We focus on defences against thermal stress by investigating patterns of expression variation within groups of the sessile California mussel (Mytilus californianus Conrad) exposed to treatments that vary in temperature mean and variance (Jimenez et al., 2015). Importantly, the mean and interindividual range of mussel body temperatures in this heterogeneous habitat are tightly and positively correlated (Table 1 and references therein). Our treatments recreate this fundamental environmental reality, at the expense of distinguishing whether the mean, the variance, or only some combination ultimately drives the resulting patterns. To obtain a more nuanced assessment of context‐dependent shifts in variation, we reanalyse the expression data sets introduced in Tanner and Dowd (2019) in a more targeted fashion. Specifically, we assess variation within 19 subnetworks that span a range of physiological functions, including canonical stress‐response pathways and empirical, de novo coexpression clusters. Using these paired transcriptome and proteome data sets, we asked (1) whether the degree of constraint on relative expression levels within the 19 subnetworks changes consistently across environmental scenarios (does stress and environmental heterogeneity decrease the strength of correlations within subnetworks?), (2) whether patterns of phenotypic variation in expression in these subnetworks also shift in a consistent manner across environmental scenarios (is variation masked or unmasked in more thermally stressful treatments?), and (3) whether interindividual variation in expression levels within these subnetworks is consistently greater for transcripts than for proteins (is there evidence for buffering between these steps of the gene expression cascade?).

TABLE 1

Transcript and protein expression profiles of mussels were analysed following exposure to one of five treatments, which represented environmental contexts that vary in both the mean of daily maximum body temperature (T max) and the heterogeneity of body temperature among individuals (here quantified as the range of T max). Temperature data for the EFA and PFA sampling locations (i.e., the wave‐exposed and wave‐protected origin sites) were obtained over 24 days in July 2010 using iButton dataloggers implanted in silicon‐filled mussel shells (Denny et al., 2011). Data for POL and POH locations were obtained in July–August 2015 and 2016 for 21 and 29 days, respectively, using thermocouples implanted in live mussels (Miller & Dowd, 2017, 2019)

Treatment	Acronym	Origin site	Mean T max (°C)	Mean range of T max (°C)	Maximum range of T max (°C)	Sample size
Exposed field‐acclimatized	EFA	Wave‐exposed	22.2	3.7	11.7	9
Protected field‐acclimatized	PFA	Wave‐protected	27.4	6.4	12.7	9
Protected common garden	PCG	Wave‐protected	13.5 ± 1.0	Negligible	Negligible	5
Protected outplant low	POL	Wave‐protected	19.8, 16.7	4.5, 2.9	12.8, 15.8	8
Protected outplant high	POH	Wave‐protected	25.8, 24.2	7.0, 7.8	14.2, 14.0	10

MATERIALS AND METHODS

Treatments and tissue collection

Adult mussels (60–70 mm shell height) were exposed to one of five treatments at Hopkins Marine Station in Pacific Grove, CA (36.61788 N, 121.91678 W); these are the same individuals and expression data described in Jimenez et al. (2015) and Tanner and Dowd (2019). The treatments represent an ecologically relevant subset of possible combinations of body temperature mean and spatial variance (Table 1), but they do not capture all possible combinations. For example, field contexts with no spatial heterogeneity would be rare.

First, mussels were sampled directly from both a wave‐exposed (cool mean temperature) and a wave‐protected (warm mean temperature) origin site; these groups constitute the exposed field‐acclimatized (EFA) and protected field‐acclimatized (PFA) treatments, respectively (Table 1, Figure S1). We then used common garden and outplant approaches to manipulate the conditions experienced by groups of mussels originating only from the protected site to examine how the amount of variation in gene and protein expression changed with recent environmental experience. Shifts in variation of antioxidant capacities among mussels from the protected site were more pronounced in Jimenez et al. (2015), and it was not feasible to generate expression data for all samples. A group of mussels collected on the same date as the PFA treatment was placed in putatively benign, common garden conditions in the laboratory for one month. They were continuously submerged at constant and spatially homogeneous temperature of 13.5 ± 1°C in flow‐through acrylic aquaria supplied with seawater from Monterey Bay. Natural food in the seawater was supplemented daily with a commercial bivalve diet (Shellfish Diet 1800; Reed Mariculture). After 28 days, one subset of individuals from this group was sampled (protected common garden treatment, PCG), and the remaining individuals were returned to two field sites. Of those remaining mussels, half were placed back in the field at a low intertidal site (protected outplant low treatment, POL) characterized by a cool mean body temperature and moderate spatial heterogeneity. The other half were placed at a high intertidal site (protected outplant high treatment, POH) characterized by a warm mean temperature and high heterogeneity (Table 1). Mussels in the POH and POL treatments were sampled after 28 days in the field. Thus, treatments were sampled on three different dates (Figure S1).

We did not monitor the body temperatures of mussels. Given the consistent differences between the outplant high and low sites in mean and spatial heterogeneity of body temperature across years (Table 1), we assumed that mussels in this study would have experienced similar patterns.

Gill tissue was excised following each treatment, at a comparable point in the tidal cycle for all treatments just as tides were receding from the field sites. This standardization of sampling times was implemented to avoid potential impacts of rhythmic patterns of expression (Gracey et al., 2008). The samples analysed here constitute the “baseline” timepoint from Jimenez et al. (2015). Tissue was frozen immediately in liquid nitrogen and stored in a freezer at –80°C until subsampling for transcript and protein expression analyses.

Transcriptomics and proteomics methods

Methods for transcriptomics and proteomics were carried out as described in Tanner and Dowd (2019); this study represents a targeted analysis of the same data set (Tanner et al., 2022). Briefly, RNA was extracted from 0.04 g of gill tissue by grinding it in a mortar/pestle in liquid nitrogen with TRI reagent (Sigma), quantified by UV‐Vis spectroscopy (Implen Nanophotometer), and cleaned with the Qiagen RNeasy Mini Plus Kit following the manufacturer's protocol (Qiagen). Small fragment cDNA library construction, barcoding, and 150 bp paired‐end (PE150; 12 samples used to generate the de novo transcriptome) or 50 bp paired‐end (PE50; 49 samples from five treatments used to quantify expression levels) Illumina HiSeq 2500 rapid‐run sequencing was performed by the Genome and Cytometry Core at the University of Southern California. The average insert sizes were 350 bp (PE150) and 380 bp (PE50). Using CLC Genomics Workbench 9.0, sequences were trimmed for Illumina TruSeq adaptors, and quality controlled for ambiguity (maximum of two ambiguous nucleotides for PE150 reads and one ambiguous nucleotide for PE50 reads) and quality (limit 0.025). Reads shorter than 60 or 30 bp were discarded for PE150 and PE50, respectively.

The longer PE150 reads from 12 haphazardly selected mussels (n ≥ 2 per treatment) were used only to construct a de novo reference transcriptome. The Transcriptome Assembly Pipeline from the National Center for Genome Assembly and Support used SAMtools v1.9 and Bowtie v2.3.2 (Langmead & Salzberg, 2012; Li et al., 2009), as described in Tanner and Dowd (2019). The resulting reference transcriptome included 130,017 transcripts with a BUSCO score of 97.7% completeness, N50 of 2514 bp, and GC content of 35.27%.

Using the trinotate v3.1.1 pipeline (Bryant et al., 2017), the shorter PE50 reads (total of 49 individuals from five treatments) were then mapped to this reference transcriptome, transcript expression levels were quantified, and a normalized TPM expression matrix was generated from total reads per million using the cross‐sample normalization method, TMM. Also in the trinotate v3.1.1 pipeline, the transcriptome was annotated with gene ontology (GO) terms and mapped to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, using the NCBI nr database reduced to only Mollusca entries.

Label‐free proteomics data were collected in the laboratory of Dietmar Kültz at the University of California, Davis. Proteins were isolated from mussel gill tissue in liquid nitrogen and then trypsin‐digested into peptides as previously described (Kültz et al., 2013). Peptides were separated in a gradient of acetonitrile using an ultra‐performance liquid chromatography system (Waters) prior to injection into a Bruker Daltonics Impact II UHR‐QTOF mass spectrometer. Tandem mass spectrometry (MS1‐MS2) captured both parent peptides and associated fragmentation products in data‐dependent acquisition (DDA) mode. Expression levels were determined with the MaxQuant algorithm (Cox & Mann, 2008). Mussel proteins were matched to MS/MS fragment spectra for in silico predicted, tryptic peptides using four different search engines (Mascot 2.2.7, X!Tandem Alanine; PEAKSX, and Byonic 2.12). The gill‐specific reference transcriptome was translated to amino acid sequences for these searches. Protein IDs based on at least two unique peptides and meeting a protein level false‐discovery rate (FDR) of 1% and a peptide level FDR of 0.1% were considered valid. Due to tissue limitations, not all mussels were analysed for proteomics.

Ultimately, 41 mussels were represented in both the transcriptome and proteome data sets and were included in the analyses described below (Table 1). We focused on the 1519 genes shared between the transcriptomics and proteomics data sets. Of these genes, 1477 (97.2%) were annotated with GO and KEGG terms. Missing expression values were imputed using the local least squares method in the PCAmethods package from the Bioconductor suite in r (Stacklies et al., 2007); 34.56% of all protein expression values and 11.89% of all transcript expression values were imputed in this manner. Transcript and protein expression levels were centered and scaled using the r scale function to allow comparisons between the two data sets; this function was applied separately for each gene identifier across all individuals.

Canonical stress‐response pathways

For the analyses of phenotypic variation and constraint within subnetworks, we first focused on six KEGG pathways (Kanehisa, 2002) representing a variety of functions central to the cellular stress response (Kültz, 2005): glycolysis (20 genes), tricarboxylic acid (TCA) cycle (22 genes), apoptosis (20 genes), proteasome (19 genes), ubiquitin‐mediated proteolysis (4 genes), and mitogen‐activated protein (MAP)‐kinase signalling (18 genes). Numbers in parentheses indicate the genes within the data set that mapped to each pathway (see Table S1), using the KEGGREST package from the Bioconductor suite in r (Tenenbaum, 2019). Glycolysis and the TCA cycle were chosen for their roles in energy production (Calow & Forbes, 1998; Sokolova et al., 2012); apoptosis, the proteasome, and ubiquitin‐mediated proteolysis function in cellular death and protein damage responses (Kültz, 2003); and MAP‐kinase signalling cascades transduce environmental signals to the nucleus (Cowan & Storey, 2003; Whitmarsh et al., 1995).

Generation of empirical subnetworks from coexpression clusters

Empirical clusters, representing subnetworks of high covariance with no a priori assumptions regarding their function, were derived separately from coexpression matrices of transcript or protein levels. This approach can identify novel interactions within the transcriptome or proteome (Kenkel & Matz, 2017; Qiao et al., 2013; Yin et al., 2018), providing complementary information to canonical pathways. Mussel empirical de novo subnetworks were constructed from the set of 1519 genes and 41 individuals across all treatments in two steps. First, we used the graph_from_adjacency_matrix function in the igraph package for r, with a correlation threshold of Pearson's r > .6, to generate separate networks for the transcriptome and the proteome. These were exported to Cytoscape 3.7.1 using the RCy3 package from the Bioconductor suite in r (Gustavsen et al., 2019). Subnetwork “clusters” were identified using the algorithm MCODE (K‐core = 3, node score cutoff = 0.5) in Cytoscape and exported as lists of gene products for downstream analyses in r. These empirical subnetworks were then analysed for GO enrichment relative to the full data set of 1519 genes (unrestricted by category; includes biological process, molecular function, and cellular component), using a Fisher's exact test with a Benjamini‐Hochberg correction for multiple comparisons (Benjamini & Hochberg, 1995).

For the analyses of subnetwork‐wide variation described below, transcript and protein expression data for all individuals in a treatment were each mapped (in separate analyses) onto all empirical clusters generated from the transcriptome or the proteome. This process was then repeated separately for each treatment.

Subnetwork‐wide analyses of expression constraint using phenotypic integration metrics

Phenotypic integration refers to the covariation of related components of a complex phenotype (Matesanz et al., 2021). We co‐opted the analytical approach of Pavlicev et al. (2009), which was developed for morphological phenotypes, to quantify treatment‐dependent changes in subnetwork‐wide phenotypic integration of expression levels. We treated expression values within a subnetwork in a similar fashion to the treatment of individual components of phenotypically integrated morphological traits, such as a limb. In this framework, highly constrained subnetworks have expression levels that tend to rise and fall in a coordinated fashion (i.e., high correlation), much as bones of a limb grow proportionally. In contrast, a less‐constrained subnetwork would have expression patterns that vary independently among the genes involved (i.e., low correlation). Because only a subset of core stress response gene products within each subnetwork would be selectively transcribed and/or translated under environmental stress, we predicted a reduction in phenotypic integration within subnetworks under more stressful conditions.

To perform these analyses, treatment‐specific covariance matrices were generated separately for each transcript and protein subnetwork, followed by calculation of the corresponding eigenvalues. Genes with no variance within a treatment were discarded to eliminate singularity issues. The standard deviation (SD) and coefficient of variation (ICV) of eigenvalues within each subnetwork‐by‐treatment combination were calculated using the evolqg package in r (Melo et al., 2015). The SD is reported as the deviation from the expected integration value of a same‐sized matrix; therefore, SD is insensitive to the sample size. Empirically, tightly constrained expression within a subnetwork manifests as a large first eigenvalue (i.e., phenotypic variation is concentrated along a single dominant axis), smaller subsequent eigenvalues, and, therefore, a high SD across eigenvalues. For every treatment‐by‐data type combination, all transcript and protein covariance matrices were compared within a single analysis per subnetwork. To visualize the overall patterns of constraint, we determined the rank of SD across the 10 combinations of treatment and data type for each of the 19 subnetworks.

Treatment‐specific analyses of subnetwork‐wide, multivariate variation

To complement the constraint metric, we evaluated variation among individuals on a subnetwork‐wide basis. This analysis first quantified how much individuals in a given treatment varied from one another in their overall expression patterns within a given subnetwork, and then asked if this level of phenotypic variation differed among groups in different treatments. We used a multivariate analog of a Levene's test, as applied to the full data set in Tanner and Dowd (2019). This modified approach compares deviations from the median rather than from the mean. It uses the underlying logic of the univariate median absolute deviation (MAD), a robust metric of variance that is calculated using the median of the absolute deviations of each individual's expression level () from the group median (). Multivariate dispersion for each treatment‐by‐subnetwork combination was calculated using the vegan package in r and the betadisper function (Oksanen et al., 2018). Distance from the treatment centroid (the median in multivariate space) was then calculated for each individual in a treatment; this distance is analogous to the deviation estimator in the univariate MAD equation. All pairwise comparisons of these distances among the five treatments were performed for each subnetwork and data type, using a modified t‐test; this resulted in 380 tests (19 subnetworks × 10 pairwise treatment comparisons × 2 data types). Significance was evaluated using the permutest function with 9999 permutations. Probabilities were corrected for multiple comparisons using the relatively conservative Benjamini‐Hochberg correction in the p.adjust function (Benjamini & Hochberg, 1995).

Analyses of gene‐by‐gene variation within subnetworks

The unmasking hypothesis predicts an overall greater amount of phenotypic variation in expression levels among individuals in more stressful (and, here, more heterogeneous) environments. However, it is also possible that the patterns of variation of individual gene products respond differently to the same environmental context because they perform different functions. For example, thermal stress might be expected to induce expression of heat shock proteins in all individuals albeit with discrepancies in the ultimate level of expression achieved, whereas other genes that play little or no role in the stress response might be more consistently expressed among individuals. To address these possibilities, we evaluated patterns of variation on a gene‐by‐gene basis for all genes within the 19 subnetworks. These comparisons followed the same logic as the multivariate variation calculations above. We calculated individuals’ absolute deviations from the treatment median expression level ; these deviations were used to evaluate pair‐wise differences in expression variation using a modified t‐test.

Three types of gene‐by‐gene comparisons of phenotypic variation were conducted, each with a different objective. First, we compared the magnitude of variation for each gene product among all five treatments (Table 1) within a given data type (transcript or protein); these comparisons identified treatments that masked or unmasked variation for that gene. Second, we compared variation between pairs of genes within each subnetwork‐by‐treatment combination (separately for each data type) to explore whether treatments had consistent effects on variation at specific nodes within each subnetwork. However, this second approach generated no significant results and will not be discussed further. Third, we compared variation between transcript and protein for each gene within each treatment to address disparities between these steps of the gene expression cascade. For each type of comparison, we corrected for multiple comparisons using the Benjamini‐Hochberg correction in the p.adjust function (Benjamini & Hochberg, 1995). Effect sizes were calculated using the effsize package to estimate Hedge's g, corrected for small sample sizes. Statistical analyses were performed in r v.3.6.2.

RESULTS

Mussels’ empirical subnetworks encompass diffuse cellular functions

Overall, coexpression structure in the proteome proved uninformative of patterns in the transcriptome, and vice versa. Accordingly, our analyses generated six de novo transcript‐based coexpression clusters and a nonoverlapping set of seven de novo protein‐based coexpression clusters (hereafter “transcriptome clusters” and “proteome clusters”, respectively) (Table S2). The average sizes of transcriptome and proteome clusters were 159.8 (range 5–410) and 20.9 (range 4–55) genes, respectively. Many of these empirical clusters contained more genes than mapped to the KEGG pathways (mean 17.2, range 4–22). There were 71 unique KEGG pathway genes that were shared with one or more empirical subnetworks (6.96% of all unique genes in the empirical subnetworks). One gene, pyruvate dehydrogenase E1, was shared among two KEGG and three empirical subnetworks.

The empirical subnetworks tended not to recapitulate known biochemical pathways or to represent consistent functions across the member genes. Only three of six transcriptome clusters and one of seven proteome clusters were significantly enriched in any GO terms (range 3–28 terms) relative to the full data set of 1519 genes (Table S3). The only clear case of functional enrichment was transcriptome cluster 2 (321 genes), which contained 49 ribosomal genes and was significantly enriched for ribosome and translation GO terms (Table S3). Otherwise, GO term enrichment occurred in small clusters (4–6 genes) in which one or two genes were associated with numerous GO terms, thus requiring caution in interpretation.

Benign conditions trend toward constrained coexpression patterns within subnetworks: Stressful and more heterogeneous environments tend to release this constraint

The pattern of constraint on expression of subnetwork components shifted considerably across treatments. Among the four treatment groups of mussels originating from the protected site, there was a clear trend toward higher constraint on relative expression levels for both transcript and protein in the two more benign treatments, protected common garden (PCG) and protected outplant low (POL) (Figure 1). Overall, the benign protected common garden treatment (PCG) exhibited the greatest constraint (the highest SD of eigenvalues) (Figure 1); this pattern also held for both the transcriptome and proteome. Across all 19 subnetworks, PCG transcript expression was the most constrained (mean rank of SD = 2.15 out of 10 combinations of treatment‐by‐data type), followed by PCG protein expression (mean rank 3.26). In contrast, protected outplant high (POH) protein expression was the least constrained (mean rank 8.68), and POH transcript expression was second least constrained (mean rank 7.44) (Figure 1). However, transcriptome and proteome patterns of constraint do not match for field‐acclimatized mussels originating from the relatively benign exposed site (EFA). In this treatment constraint on protein levels was more similar to the stressful PFA and POH treatments, whereas constraint on transcript levels occupied a median position relative to the other treatments. Moreover, there were some exceptions to these overall trends depending on the subnetwork type (Figure S2).

FIGURE 1

Common garden (PCG) conditions (blue) consistently induced the greatest constraint on relative expression levels of gene products across the 19 subnetworks analysed, whereas the outplant high (POH) treatment (red) had the lowest overall constraint. Box plots show the rank of the standard deviation (SD) of eigenvalues across the 10 treatment‐by‐data type combinations. Ranks were calculated from Tables S4–S6; each box represents 19 ranks, one for each of the subnetworks included in the analyses. Note that the y‐axis values are reversed, with high ranks at the top. High rank of SD corresponds with greater constraint on the relative expression levels of gene products within a given subnetwork, illustrated by the arrow to the right. Treatments colour‐coded as in Figure S1. Treatment codes: EFA, exposed field‐acclimatized; PCG, protected common garden; PFA, protected field‐acclimatized; POH, protected outplant high; POL, protected outplant low; POL, protected outplant low. Vertical line separates mussels originating from the wave‐exposed field site (left) from those originating from the wave‐protected field site (right)

Across all 10 treatment‐by‐data type combinations, the KEGG subnetworks displayed the highest phenotypic integration (i.e., the most constraint on relative expression levels) in the protected common garden treatment (Table S4, Figure S2A). Transcript levels in PCG had the largest amount of variance among eigenvectors (overall mean SD = 0.66, mean ICV = 2.75; Table S4) in five of six pathways (Figure 2); the remaining pathway had transcript PCG as the second most constrained group. PCG also was the most constrained treatment for protein levels for five of the six KEGG pathways (Table S4). In contrast, protein levels for the thermally stressful POH and protected field‐acclimatized (PFA), as well as the relatively benign exposed field‐acclimatized (EFA) treatment, were among the least constrained groups for all pathways (POH: mean SD = 0.30, mean ICV = 1.47; PFA: mean SD = 0.31, mean ICV = 1.52; EFA: mean SD = 0.32, mean ICV = 1.54).

FIGURE 2

Subnetwork‐wide phenotypic integration metrics for transcript and protein expression levels in six KEGG pathways. Solid lines indicate transcript expression, dashed lines indicate protein expression. For each eigenvector the corresponding eigenvalue reflects the proportion of total variation in expression in that dimension. Subnetworks whose component gene products are tightly constrained in their relative expression levels across individuals have large first eigenvalues (for example, proteasome solid blue PCG treatment) and/or very small eigenvalues beyond the first few eigenvectors. For example, the PCG treatment's eigenvalues for transcript levels (solid blue lines) drop to zero rapidly with increasing eigenvector number, despite not always exhibiting the absolute largest first eigenvalue in all panels. Consequently, PCG transcript levels exhibit the largest SD among eigenvectors for all six pathways, indicating the most tightly constrained expression. By contrast, POH transcript levels (solid red lines) are less constrained across individuals in all panels, with larger eigenvalues at higher eigenvector numbers. Treatment codes: EFA, exposed field‐acclimatized; PCG, protected common garden; PFA, protected field‐acclimatized; POH, protected outplant high; POL, protected outplant low. The proteolysis pathway had few eigenvalues, because only four genes were represented in the data set

The degree of phenotypic integration shifted similarly among treatments within the empirical transcriptome clusters (Figure S2B). The PCG treatment exhibited the most constraint (i.e., the greatest SD among eigenvectors) for five of six transcriptome clusters (Figures S3A, B, D‐F; mean SD = 0.57, ICV = 7.47; Table S5). Three transcriptome clusters had the most constraint on relative expression levels for proteins in PCG (transcriptome clusters 1, 4, 5; Figures S3A,D,E) and two clusters had the most constraint for transcript levels in PCG (transcriptome clusters 2,6; Figures S3B,F). Similar to the KEGG subnetworks, transcript and protein expression from the thermally stressful POH treatment exhibited the two lowest levels of phenotypic integration for five of the six transcriptome clusters (1, 2, 4, 5, and 6) (Table S5).

Measures of constraint in the seven empirical proteome clusters revealed more complex patterns (Figure S2C; Figures S3G–M). Like the other subnetwork types, four proteome clusters had the most constraint for transcript expression levels in PCG (proteome clusters 1, 2, 3, 5; mean SD = 0.70, mean ICV = 3.80; Table S6). Similarly, in two proteome clusters PCG was the most constrained treatment for protein expression levels (proteome clusters 1 & 3; mean SD = 0.54, mean ICV = 1.68). The thermally stressful POH treatment tended to exhibit less constraint, although not as consistently as for the KEGG and transcriptome subnetworks. The POH‐protein group (proteome clusters 5 and 7), POH‐transcript group (Proteome Cluster 1), or both POH groups (proteome clusters 2, 4, and 6) fell among the three least constrained of the ten treatment‐by‐data type combinations for six of the seven proteome clusters (Table S6). However, Proteome Cluster 7 exhibited the highest constraint for transcript levels in the POH treatment (Figure S3M).

Treatment‐dependent shifts of multivariate phenotypic variation

Multivariate, subnetwork‐wide metrics of phenotypic variation increased only in the transcriptome and only in a stressful and heterogeneous treatment. All nine of the significant pairwise comparisons of multivariate MAD among treatments (out of 380 total comparisons) involved a larger magnitude of variation in the POH treatment for transcript levels (1.74 times greater MAD on average; Table 2; see Table S7 for uncorrected results). Multivariate MAD for transcript levels in POH was significantly elevated compared to the EFA treatment for five subnetworks, including glycolysis and proteasome pathways. Multivariate MAD in POH was elevated relative to the PCG treatment for transcriptome cluster 1 (Table 2). The POH treatment exhibited a greater magnitude of variation relative to PFA treatment for three subnetworks (Table 2).

TABLE 2

Significant pairwise comparisons of subnetwork‐wide, multivariate median absolute deviation (MAD) among the five treatments. Comparisons were limited to the same data type (transcript or protein levels). The treatment listed first (treatment 1; T1) exhibited greater multivariate MAD than treatment 2 (T2), interpreted as greater interindividual variation

Subnetwork	Data type	Treatment 1	Multivariate MAD (T1)	Treatment 2	Multivariate MAD (T2)	Adjusted p
Glycolysis	Transcript	POH	4.20	EFA	2.71	.03
Glycolysis	Transcript	POH	4.20	PFA	3.04	.03
Proteasome	Transcript	POH	4.55	EFA	2.35	<.001
Proteome cluster 1	Transcript	POH	6.21	EFA	4.43	.014
Proteome cluster 1	Transcript	POH	6.21	PFA	4.41	.014
Proteome cluster 4	Transcript	POH	2.23	EFA	1.12	.038
Transcriptome cluster 1	Transcript	POH	17.90	PCG	8.20	.003
Transcriptome cluster 1	Transcript	POH	17.90	EFA	9.97	.035
Transcriptome cluster 1	Transcript	POH	17.90	PFA	9.04	.005

Treatment codes – EFA, exposed field‐acclimatized; PCG, protected common garden; PFA, protected field‐acclimatized; POH, protected outplant high; POL, protected outplant low.

Expression levels of hundreds of genes exhibit greater variation in stressful treatments

Gene‐by‐gene comparisons of variation in expression levels among the five treatments revealed 222 unique genes (21.76% of the 1020 genes in the 19 subnetworks) for which transcript and/or protein variation significantly differed between any two treatments (Table S8). Overall, there were 57 significant pairwise comparisons within KEGG pathways (55 for transcript, 2 for protein) involving 15 unique genes (Table 3); 664 significant comparisons within empirical transcriptome clusters (646 for transcript, 18 for protein) involving 220 unique genes; and none within empirical proteome clusters (Table S8). Some differentially variable genes mapped to more than one subnetwork.

TABLE 3

KEGG pathway genes identified as significantly different in univariate MAD for transcript or protein levels between two or more treatments. In all cases, the first listed treatment exhibited higher inter‐individual variation than the second. Bold font indicates genes for which transcript level variation in the POL treatment exceeded that in POH (see text). Two genes (phosphoenolpyruvate carboxykinase, cytosolic [GTP]‐like and dihydrolipoyl dehydrogenase, mitochondrial) were represented by two isoforms in the data set; both isoforms had the same significant treatment comparisons. All but two significant treatment comparisons were for transcript levels; the two significant comparisons for protein levels are denoted with asterisks

Gene description	Significant treatment comparisons
Triosephosphate isomerase	EFA‐POL*, POH‐EFA, EFA‐PFA*, POL‐POH, POH‐EFA
Isocitrate dehydrogenase	POH‐EFA, POH‐PFA
Heat shock cognate 71 kDa protein	POH‐EFA, POH‐PCG, POH‐PFA, POL‐EFA, POL‐PFA
cdc42 homologue	POH‐EFA, POH‐PCG, POH‐PFA
Succinate dehydrogenase (ubiquinone) iron‐sulphur subunit, mitochondrial‐like	POH‐EFA
Phosphoenolpyruvate carboxykinase, cytosolic (GTP)‐like (two isoforms)	POH‐PCG, POH‐EFA, POH‐PFA, POL‐POH
Heat shock protein 70	POH‐EFA, POH‐PCG, POH‐PFA, POL‐POH
Fructose‐1,6‐bisphosphatase 1‐like isoform X2	POH‐EFA, POH‐PFA, POL‐POH
Pyruvate carboxylase, mitochondrial‐like isoform X1	POH‐EFA, POH‐PFA, POL‐POH
Serine/threonine‐protein kinase PAK 3‐like isoform X3	POH‐EFA, POH‐PCG, POH‐PFA, POL‐POH
Ubiquitin‐like modifier‐activating enzyme 1	POH‐PFA, POL‐POH
Septin‐2‐like isoform X2	POH‐POL
Filamin‐A‐like isoform X1	POH‐EFA, POH‐PCG, POH‐PFA, POL‐POH
Dihydrolipoyl dehydrogenase, mitochondrial (two isoforms)	POH‐EFA, POH‐PCG, POH‐PFA, POL‐POH
Phosphoglycerate kinase	POH‐EFA, POH‐PCG, POH‐PFA, POL‐POH

The unmasking hypothesis predicts a predominant pattern of enhanced variation in expression levels in more stressful conditions. Across all subnetwork types, 650 of the 721 significant pairwise treatment comparisons (90.2%) involved the stressful and heterogeneous POH treatment, 494 of which (76.0%) indicated higher variation in expression in that treatment (178 relative to EFA, 123 relative to PCG, 178 relative to PFA, 15 relative to protected outplant low [POL]). The mean effect sizes for those genes with higher variation in the POH treatment were 1.21 ± 0.30 [SD] and 1.36 ± 0.41 in KEGG pathways and transcriptome clusters, respectively.

Closer examination of the gene‐by‐gene patterns of phenotypic variation revealed three additional observations. First, variation in the POL treatment exceeded that in POH in a considerable number of cases. Specifically, POL variation exceeded POH variation for transcript levels of nine unique KEGG pathway genes (bold in Table 3). These nine genes primarily play roles in energy metabolism or protein homeostasis. Within transcriptome clusters, 118 of 174 unique genes (67.8%) that were identified as differentially variable in at least one treatment comparison involving POH exhibited elevated variation in POH relative to the same three treatments (PCG, EFA, and PFA). Yet, for each of these 118 genes the magnitude of variation for expression in POL exceeded POH (see shaded rows in Table S8). Nearly all of these 118 genes differed in variation for transcript levels; only one instance involved shifts in protein‐level variation for a putative NAD(P)H: quinone oxidoreductase (Table S8). These 118 genes were associated most commonly with processes related to protein homeostasis (GO:0006457 protein folding, 14 genes; GO:0042026 protein refolding, 5 genes; GO:0006886 intracellular protein transport, 4 genes), redox balance (GO:0045454 cell redox homeostasis, 5 genes; GO:0055114 oxidation‐reduction process, 5 genes), and cytosolic carbohydrate metabolism (GO:0006094 gluconeogenesis, 3 genes; GO:0006096 glycolytic process, 3 genes), all of which are components of the cellular stress response.

Second, a small subset of 15 genes in transcriptome clusters 1–3 exhibited higher MAD for one of the more benign treatments relative to the POH treatment (Table S9). Seven instances were between PCG and POH (mean effect size = 1.40 ± 0.39), and eight were between EFA and POH (mean effect size = 1.19 ± 0.22). In addition to the antioxidant glutathione S‐transferase, this group included several genes involved in chromatin (re)organization. The only protein‐level example in this group was ribose‐5‐phosphate isomerase, which contributes to nucleic acid/nucleotide biosynthesis via the nonoxidative branch of the pentose phosphate pathway.

Third, surprisingly there were no genes for which MAD values for transcript and protein expression levels shifted significantly in the same direction between the same two treatments. One gene had opposing patterns for MAD of transcript and protein levels between the same two treatments; within transcriptome cluster 1, ribose‐5‐phosphate isomerase had elevated variation in EFA relative to POH for protein expression, whereas transcript variation shifted in the opposite direction (Table S8).

Little evidence for buffering of variation in expression levels between transcriptome and proteome

We also compared variation between expression levels of transcripts and their corresponding proteins to search for a systematic bias toward elevated variation at the transcript level, as would be expected if there were buffering between these steps of the gene expression cascade. In our subnetwork‐wide analyses, we found considerable evidence refuting the notion that molecular phenotypes are more constrained at the protein than at the transcript level. Specifically, within KEGG pathways, transcript relative expression levels were more constrained than the corresponding protein levels in 29 of 30 (96.7%) pathway‐by‐treatment combinations (Figure S2A, Table S4). Similarly, transcript expression within empirical transcriptome and proteome clusters was more constrained for the large majority of cluster‐by‐treatment combinations (20 out of 30 [66.7%], and 22 out of 35 [62.9%], respectively) (Figures S2B,C; Tables S5 and S6).

Similar patterns emerged in gene‐by‐gene comparisons; although relatively few, significant comparisons overwhelmingly demonstrated greater phenotypic variation at the protein than at the transcript level. When we pooled across all 41 mussels, none of the 1519 genes in the data set exhibited a global difference in variation between transcript and protein levels. Inspection of treatment‐specific patterns revealed 12 instances of variation for protein expression levels exceeding that for transcript levels; these were largely concentrated in the POH and POL treatments (Table 4). In KEGG pathways, four unique genes (3.9% of KEGG genes analysed) exhibited significantly greater variation for protein levels (mean effect size = 2.25 ± 0.23, Table 4). Within empirical transcriptome clusters, protein variation was higher than transcript variation for seven genes (0.46% of all empirical subnetwork genes analysed; mean effect size = 2.94 ± 0.86), including two in POH and five in POL (Table 4). There were no instances in the empirical proteome clusters.

TABLE 4

Genes with significant differences in the magnitude of interindividual variation between protein and transcript expression levels within a treatment group. Protein variation is higher than transcript variation in each case

Subnetwork	Treatment	Gene description	Effect size	Adjusted p
Transcriptome cluster 1	POH	Adenylate kinase 8‐like	2.42	.020
Transcriptome cluster 1	POL	Fructose‐1,6‐bisphosphatase 1‐like isoform X2	3.53	.005
Transcriptome cluster 1	POL	Uncharacterized protein LOC110455829	3.20	.005
Transcriptome cluster 1	POL	Major egg antigen‐like	2.60	.012
Transcriptome cluster 1	POL	Endoplasmic reticulum chaperone BiP	3.61	.012
Transcriptome cluster 2	POL	Heat shock protein 70	3.80	.028
Transcriptome cluster 5	POH	Calmodulin 3b (phosphorylase kinase, delta)	1.39	.032
Glycolysis	POH	Pyruvate dehydrogenase E1 component subunit beta, mitochondrial‐like	2.58	.001
MAP‐kinase signalling	EFA	Filamin‐A‐like isoform X1	2.17	.017
MAP‐kinase signalling	PFA	Heat shock protein 70	2.06	.031
TCA cycle	PFA	Dihydrolipoyl dehydrogenase, mitochondrial	2.04	.035
TCA cycle	POL	Dihydrolipoyl dehydrogenase, mitochondrial	2.38	.014

Treatment codes – EFA, exposed field‐acclimatized; PFA, protected field‐acclimatized; POH, protected outplant high; POL, protected outplant low.

Protein levels also tended to exhibit higher variation among individuals than transcript levels when variation was averaged across all of the genes within a subnetwork. Specifically, in 72 of the 95 subnetwork‐by‐treatment scenarios analysed, protein levels exhibited higher mean variation (quantified as mean of the scaled univariate MAD values for the subnetwork) relative to transcript levels (Figure S4, black bars in A[all], B[i–v], and C[i–iii, v, vii]). The majority of the 23 instances of higher mean transcript variation occurred in two treatments. Specifically, six subnetworks had higher average transcript than protein variation in the POH treatment (Figures S4A[i, iii, iv], C[i, iv, vii]), and 10 subnetworks had higher average transcript than protein variation in the PCG treatment (Figures S4A[iii, v], B[iii, iv, v, vi], C[ii, iii, iv, vi]).

DISCUSSION

Three complementary methods of analysing interindividual phenotypic variation within subnetworks of the gill transcriptome and proteome of adult M. californianus tell a broadly consistent story. Environmental stress in combination with environmental heterogeneity is associated with greater magnitudes of variation in expression levels, as well as with a breakdown in phenotypic integration (i.e., reduced correlation of expression levels) within subnetworks of the transcriptome and proteome. These results generally support the unmasking hypothesis introduced above, and they provide a more nuanced assessment than global trends for these same data sets at the transcriptome‐wide and proteome‐wide scale (Tanner & Dowd, 2019). Other physiological measurements conducted on the same individuals analysed here indicate that variation in antioxidant capacities and organismal metabolic rate is similarly unmasked in more stressful and heterogeneous treatments (Jimenez et al., 2015). Context‐dependent patterns of phenotypic variation have only infrequently been explored for physiological traits, but behavioural studies demonstrate robustly that stress acts to expose variation (Roche et al., 2016; Westneat et al., 2015). Such findings further illustrate that the relationships between phenotypic variation and the surrounding environment are dynamic and complex (Canino‐Koning et al., 2019; Chen et al., 2019; Oleksiak & Crawford, 2012; Takahashi, 2018).

Notably, these trends appear only in some subnetworks of the mussel transcriptome and proteome, and for only some of their component genes. Inconsistent patterns of variation among genes and/or functions necessitate caution when attempting to generalize the effects of environmental conditions on phenotypic variation. Patterns of variation also manifest uniquely in these data at the transcript and protein stages of the gene expression cascade, with little overlap. Finally, the apparent release of expression constraints within subnetworks under stressful conditions implicates specific, labile subnetwork “nodes” (see Table 3) as possible contributors to variation in physiological phenotypes and performance.

Relative expression levels within mussel subnetworks are more constrained in benign than in stressful conditions

Common garden conditions consistently exhibited high levels of phenotypic integration, with high correlation of expression levels within the subnetworks analysed here (Figure 1). Mussels in this putatively most benign treatment were continuously submerged, under constant cool temperature. The concentration of most of the variance in each subnetwork into one eigenvalue (Figure 2 & Figure S3) implies relatively tightly coordinated expression mechanisms across genes (Hanisch et al., 2002; van Waveren & Moraes, 2008). Individuals did vary in their absolute expression levels in common garden, but expression levels of the genes within their relatively constrained subnetworks tended to rise and fall together. In contrast, we consistently found the least constraint on subnetwork‐wide expression levels in more stressful and heterogeneous treatments, particularly in the protected outplant high treatment.

This pattern of reduced subnetwork‐wide phenotypic integration under stressful conditions appears consistent with broadly conserved patterns in the stress response. Specifically, because environmental stress is known to selectively induce expression of a core suite of defence genes (Kültz, 2005), relative expression levels of most other genes must decrease by comparison. The logical prediction, supported by many of our observations, is decreased integration among the components of subnetworks that contain those core stress response genes when compared to more benign conditions.

The more difficult, still unresolved question is whether relative capacities to upregulate only core stress response genes truly differentiate individuals in their abilities to cope with environmental stress. This conjecture appears commonly in the literature, for example in discussion of the role of molecular chaperones in heat resilience (e.g., Chen et al., 2018; Moseley, 1997). Alternatively, differences among individuals in the ability to sustain expression of a complementary suite of genes outside of the core stress response may ultimately predict interindividual differences in physiological outcomes (Dixon et al., 2020; Wang et al., 2013; Zhao et al., 2016); this alternative rests on the assumption that nearly all individuals sufficiently activate the core stress response. To date, we have insufficient longitudinal expression or performance data on our study organism to fully evaluate these alternatives, but our empirical results tend to support the latter explanation.

Specifically, the empirical coexpression clusters generated from our data sets encompass a diffuse set of biological functions. For example, transcriptome cluster 2 was significantly enriched in numerous GO terms related to ribosome structure and translation, a pattern driven by genes not frequently implicated in the core cellular stress response. Yet, this large coexpression cluster also contained genes representing a variety of classical core stress response functions: antioxidants (e.g., thioredoxin, glutathione S‐transferase), protein stabilization and turnover (e.g., several proteasome subunits, heat shock protein 70), and energy metabolism (e.g., cytochrome c oxidase subunits, succinate dehydrogenase, glyceraldehyde‐3‐phosphate dehydrogenase) (Abele & Puntarulo, 2004; Buttemer et al., 2010; Tomanek, 2015). Indeed, nearly all of the genes in the six KEGG pathways – which were chosen for their known or putative roles in the stress response – were also represented in the de novo empirical subnetworks. Yet, these KEGG genes were distributed across the empirical subnetworks and interspersed with a variety of other functions (Table S2). Further investigation of these apparent linkages among diverse functional groups is needed under different environmental conditions to expand our knowledge of the many facets of the stress response.

Furthermore, it is likely that the covariance structures and composition of coexpression subnetworks themselves shift among environmental contexts (Adamo, 2012; Fait et al., 2020; Guy et al., 2008). With larger sample sizes than employed here, an even more nuanced approach would be to regenerate the empirical subnetworks using only individuals from a single treatment condition and to then compare the composition of subnetworks across treatments.

Stressful and heterogeneous conditions unmask phenotypic variation in expression, both subnetwork‐wide and on a gene‐by‐gene basis

Overwhelmingly across the 19 subnetworks, expression levels in the protected outplant high (POH) treatment tended to be more variable among individuals, particularly when compared with common garden and field‐acclimatized conditions. We anticipated that analyses of multivariate or univariate MAD would largely mirror subnetwork‐wide trends from the phenotypic integration analysis, such that MAD would be higher in less constrained treatment‐by‐data type groups. In line with this prediction, a pattern of high MAD in POH treatment was evident for both entire subnetworks and their individual component genes. Over 91% of all significant comparisons among treatments for individual genes revealed higher variation in the POH or POL treatments. Many of these genes play roles in the cellular stress response (e.g., molecular chaperones, oxidative stress response proteins) and/or cytosolic (i.e., anaerobic) energy metabolism. In contrast with this overall pattern, a much smaller subset of genes had greater variation in benign conditions. This latter group of genes appear largely associated with chromatin organization and transcriptional regulatory functions.

The contrasting patterns of variation in expression for stress/metabolism and transcription are reminiscent of previous findings in M. californianus (Gracey et al., 2008). That study identified a “cell‐division” cluster of genes for which expression levels oscillated between high and low tide. Notably, this cell‐division cluster was anticorrelated with expression of a “metabolism” gene cluster (Gracey et al., 2008). These oscillations were larger in magnitude at a high intertidal site. Taken together, the current results and those of Gracey et al. (2008) highlight the complex, disparate regulation of different cellular functions in mussels.

The present study captures patterns of variation primarily using field treatments, but further work is needed to disentangle the environmental factors that are responsible for these patterns. For example, at the POH field site prolonged periods of air exposure allow individuals to reach high body temperatures, but an increase in the magnitude of heterogeneity in body temperature accompanies the increased mean (Denny et al., 2011; Miller & Dowd, 2019). Interestingly, variation of individual genes was greater in the POL treatment than in POH more frequently than the reverse pattern occurred. Mussels at the POL site can also reach high mean body temperatures, albeit usually for a shorter duration during low tide. Furthermore, heterogeneity of mussel body temperature is also substantial at the POL site (Miller & Dowd, 2019; Table 1). As noted above, our experimental design cannot distinguish between an elevated mean temperature, elevated thermal heterogeneity, or only the two in combination as the driver of variation in expression levels. However, given the positive correlation between mean and heterogeneity for body temperatures of our study organism (Miller & Dowd, 2019), it remains reasonable to conclude that more thermally stressful environmental conditions will tend to coincide with elevated variation in expression levels among neighbouring mussels.

Expression variability is higher for proteins than for transcripts in some cases

Overall, we found little evidence for buffering of variation between the transcript and protein stages of the gene expression cascade. Rather, our results implicate translational control as a likely mechanism generating variation among individuals at the protein level (Gawron et al., 2014; Lackner et al., 2012; Vasquez et al., 2014). Although transcriptome‐wide variation exceeded proteome‐wide variation in the whole data set (Tanner & Dowd, 2019), we observed only the opposite pattern in the current, targeted analyses. At both the subnetwork‐wide and individual‐gene scales, these instances of higher phenotypic variation for protein expression relative to transcript expression were more prevalent in more stressful and heterogeneous treatments. These seemingly contradictory results might be reconciled in part by the fact that the 19 subnetworks encompassed only a fraction of the total transcriptome/proteome.

Given the larger percentage of imputed values in the protein expression data set, these patterns should be interpreted with caution. However, we suspect that our imputation method artificially underestimates the magnitude of variation in protein expression levels. Therefore, additional data are unlikely to systematically reverse the direction of the observed differences between protein and transcript variation.

We cannot attribute the nearly complete disconnect between patterns of variation for transcript and protein levels to any specific mechanism, in part because the data capture only a single snapshot in time. Importantly, there are several molecular mechanisms that might introduce temporal lags between the transcriptome and proteome, including recent discoveries in mRNA processing and turnover and the role of micro‐RNAs in regulation of expression (Baek et al., 2008; Greenbaum et al., 2003; Liu et al., 2016; Mata et al., 2005). Additional work incorporating a time‐course of sampling is needed to examine the persistence of these observed patterns. In the mussel system, phenotypic variation could arise from cumulative environmental effects of benign/stressful treatment, single acute stress events immediately before sampling, or some combination. Some available data, including measurements of antioxidant capacities on the same individuals analysed here (Jimenez et al., 2015), suggest that the magnitude of physiological variation might not change over acute environmental stress episodes (Tanner & Dowd, 2019). We attempted to control for any such effect in our design by sampling individuals during periods when the tide was just receding and before stressful conditions would be encountered on that day.

Other aspects of our results further highlight the disparities between the mussel transcriptome and proteome, perhaps most notably the fact that for no single gene did we observe transcript and protein expression variation shifting in the same direction between two treatments. This surprising result demands further study to examine whether it is generalizable across organisms and types of stressful environmental scenarios.

CONCLUSION

The relationships among environmental heterogeneity, environmental stress, and phenotypic variation are complex, as illustrated here by context‐dependent and subnetwork‐dependent patterns of variation in transcript and protein expression levels. Stressful and environmentally heterogeneous conditions tended to unmask interindividual variation and reduce phenotypic integration within mussel subnetworks. Additional work is needed in the M. californianus system to further disentangle the roles of chronic and acute stress from that of environmental heterogeneity per se. Studies that expand on the instantaneous “snapshots” of variation captured here and that link variation in expression levels to performance metrics that are reasonable proxies for fitness should prove very informative. Nonetheless, the approaches implemented here can identify specific subnetworks or individual genes for which the environment strongly influences the magnitude of interindividual variation. Because such variation is crucial for evolutionary processes (Gould, 1985; Whitehead & Crawford, 2006), deeper exploration of interindividual variation in biochemical and physiological traits is essential. Circumstances arising from global change, such as an increasing frequency of stressful thermal events (Dowd & Denny, 2020; Stillman, 2019), make understanding of these relationships critical to forecasting future biological responses.

AUTHOR CONTRIBUTIONS

Lani U. Gleason and W. Wesley Dowd designed the research, Lani U. Gleason generated transcriptome data, Richelle L. Tanner analysed data, Richelle L. Tanner and W. Wesley Dowd wrote the manuscript.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.