Gene Networks in the Wild: Identifying Transcriptional Modules that Mediate Coral Resistance to Experimental Heat Stress.

Organisms respond to environmental variation partly through changes in gene expression, which underlie both homeostatic and acclimatory responses to environmental stress. In some cases, so many genes change in expression in response to different influences that understanding expression patterns for all these individual genes becomes difficult. To reduce this problem, we use a systems genetics approach to show that variation in the expression of thousands of genes of reef-building corals can be explained as variation in the expression of a small number of coexpressed "modules." Modules were often enriched for specific cellular functions and varied predictably among individuals, experimental treatments, and physiological state. We describe two transcriptional modules for which expression levels immediately after heat stress predict bleaching a day later. One of these early "bleaching modules" is enriched for sequence-specific DNA-binding proteins, particularly E26 transformation-specific (ETS)-family transcription factors. The other module is enriched for extracellular matrix proteins. These classes of bleaching response genes are clear in the modular gene expression analysis we conduct but are much more difficult to discern in single gene analyses. Furthermore, the ETS-family module shows repeated differences in expression among coral colonies grown in the same common garden environment, suggesting a heritable genetic or epigenetic basis for these expression polymorphisms. This finding suggests that these corals harbor high levels of gene-network variation, which could facilitate rapid evolution in the face of environmental change.

Introduction

The persistence of populations in the face of environmental variation in space and time is in large part determined by the ability of individual organisms to match parts of their phenotypes to their environments (Moran 1992; Kawecki and Ebert 2004; Marshall et al. 2010). When environments change, a variety of mechanisms play major roles in phenotype environment matches. Migration of species to more suitable habitat in response to environmental change has been observed in many changing ecosystems (Chen et al. 2011). Local adaptation and phenotypic plasticity have also been observed to play a major role in the regional persistence of populations faced with environmental change (Kozlowski and Pallardy 2002; Stillman 2003; Leimu and Fischer 2008; Hereford 2009; Sanford and Kelly 2011). Strong natural selection coupled with standing genetic variation may allow rapid adaptation in the face of environmental change (Barrett and Schluter 2008; Chevin et al. 2010; Pespeni et al. 2013; Orr and Unckless 2014). Plastic responses to environmental change may allow populations to persist long enough to adapt to changing conditions, provided that enough standing variation in adaptive traits exists at the limits of acclimation (Ghalambor et al. 2007). As widespread climate change increases in strength, maintenance of ecosystem function will in large part depend on the ability of ecologically important species to maintain a sufficiently close phenotypic match to their changing environment through some combination of these mechanisms (Barrett and Schluter 2008; Chevin et al. 2010; Somero 2010; Orr and Unckless 2014).

A substantial proportion of phenotypic variation both between and within species is generated by variation in gene expression (Wray et al. 2003). Populations in which variation in gene expression can be related to environmental history and stress tolerance can provide a mechanistic physiological framework for exploring the potential adaptive and acclimatory responses of populations to environmental change. For example, both constitutive and plastic variation in gene expression across environmental gradients has been repeatedly linked to population-level responses to environmental stress in widespread plant species (Swindell et al. 2007; Lasky et al. 2015). Yet, experimental studies of natural populations often identify transcriptional changes at thousands of genes as a function of season, diurnal timing, acclimation state, diet, acute environmental stress, or a host of other factors (Oleksiak et al. 2002; Jaenisch and Bird 2003). It is often difficult to determine which of these gene expression differences are most important in generating differences in organismal phenotypes, or when differences in the expression of individual genes are pivotal to organismal physiology. For example, we recently showed that thousands of genes show altered expression in reef corals in response to experimental heat stress (Seneca and Palumbi 2015) Determining the roles of these genes as coral reefs are faced with rapid environmental change is a huge challenge.

A great deal of recent progress in analyzing complex gene expression variation has come from studies using a systems genetics approach, which identifies groups of coexpressed genes with correlated responses across samples with genetic and environmental differences. These expression modules are thought to represent physiological and developmental units and have been shown to correlate with differences in physiological or morphological phenotypes. Coexpressed sets of genes can be characterized by examining them for enrichment of particular gene classes, and by testing them for associations with genetic or environmental influences. Overall, systems genetics approaches greatly simplify subsequent gene expression analysis and better reflect the underlying cellular physiology captured in gene expression studies. By reducing variation in transcripts from thousands of genes to a small number of functionally distinct transcriptional modules made up of sets of coexpressed genes, these approaches simplify hypothesis testing and reduce the number of independent statistical tests applied to highly intercorrelated individual gene expression profiles (Civelek and Lusis 2014).

Here, we characterize variation in gene expression patterns in the context of heat tolerance by reef-building corals, which suffer a breakdown of their coral–dinoflagellate symbioses in response to small increases in temperature. Corals show high levels of gene expression variation between environments, and some populations of the same species are more stress tolerant than others (Oliver and Palumbi 2011; Barshis et al. 2013). Previous work has described conspecific populations with different levels of heat tolerance; transplants between them show that about half of the between-population phenotypic difference is due to acclimation between microclimates and half is due to constitutive differences between populations (Palumbi et al. 2014). Other work has shown that latitudinal differences in bleaching resilience are heritable in crossbreeding experiments and are related to differences in the expression of a large set of “tolerance-associated genes” (Dixon et al. 2015).

In this study, we take advantage of transplanted clones to characterize gene expression variation in natural populations of corals in two different common gardens experiencing different thermal regimes. We previously described an experiment in which 3 replicates of each of 20 colonies of the tabletop coral Acropora hyacinthus were reciprocally transplanted between adjacent back-reef lagoons with different temperature profiles (Seneca and Palumbi 2015). After a year, we measured gene expression profiles for each colony acclimated to each transplant site under normal environmental conditions, and after a standard heat stress that mimics strong coral bleaching effects. Previous analyses found that thousands of genes showed altered expression 5 and 20 h after heat stress (Seneca and Palumbi 2015).

Here, we show that these thousands of genes cluster into just a few coexpressed transcriptional modules with distinct functional enrichments. These modules show coordinated responses to environmental perturbations but also show consistent differences in expression between individual coral colonies. One such module includes a set of tightly coexpressed sequence-specific DNA-binding proteins including several transcription factors that are correlated with both colony-level and acclimatory variation in bleaching susceptibility. In addition, the existence of high levels of standing variation in induction of “bleaching modules” from individual to individual suggests that bleaching-related gene networks may harbor enough standing genetic variation to be capable of rapid evolution in response to environmental change.

Materials and Methods

Transplants and Experimental Heat Stress

The experiment discussed here was carried out in August 2011, and is described in full in a separate publication (Seneca and Palumbi 2015). Briefly, we reciprocally transplanted nubbins from 20 different colonies of the tabletop coral A. hyacinthus between 2 back-reef pools that experience different thermal regimes in Ofu, American Samoa. Corals in the more stressful highly variable (HV) pool regularly experience temperatures over 32 °C and as high as 35 °C and have high thermal tolerances compared with conspecifics found outside of the pool, while corals in the less stressful moderately variable (MV) pool rarely experience temperatures over 32 °C (Craig et al. 2001; Oliver and Palumbi 2011; Barshis et al. 2013; Palumbi et al. 2014). In December 2012, we subjected growing coral transplants to a standardized heat stress: We sampled transplants in each pool at 9:00 AM on the morning of each trial and placed them into custom-built stress tanks that allow us to tightly control water temperature over time. Each 6-l tank was exposed to standardized light conditions of 700 µE during daylight hours and was subject to a flow through of 5 l/h. During each trial, starting at 10:00 AM, the temperature of the stress tanks underwent a 4 h ramp from 29 °C to 35 °C, continued exposure to 35 °C for 1 h and then a return to 29 °C. We also subjected paired samples to consistent temperatures of 29 °C.

We sampled coral nubbins at 5 and 20 h into the standardized stress, representing time points immediately after heat stress and the following morning, when bleaching was apparent. We placed samples immediately into RNAlater at 4 °C, and transferred them to −80 °C after 24 h. Ultimately, we obtained samples from both pools and both experimental treatments for 16 coral colonies of origin. Additionally, we had samples from only the HV pool for one colony, samples from only the MV pool for two colonies, and an additional set of samples in the MV pool for one colony, as well as an additional set of samples in both pools for one colony (supplementary table S6, Supplementary Material online). At hour 20, we scored heat-stressed samples on a visual bleaching scale of 1 (no visible bleaching) to 5 (total bleaching).

RNA-Seq

We used Trizol/chloroform to extract RNA from coral samples, created libraries using Illumina’s TruSeq kit from poly-A+ RNA, and carried out multiplexed Illumina sequencing (12 samples per lane). We mapped reads to the A. hyacinthus (var surculosa) transcriptome (Barshis et al. 2013) using the Burrows–Wheeler Aligner (Li and Durbin 2009). Ultimately, 162,170,810 reads from 152 libraries were mapped to 33,496 coral contigs. We normalized expression counts with DESeq2 (Anders and Huber 2010). We carried out all further analyses discussed here on 15,737 contigs that averaged >5 normalized mapped reads per sample. We measured symbiont genotype proportions by mapping raw sequence reads to the Symbiodinium clade C and D 23S ribosomal genes and counting the proportion of reads uniquely mapping to each clade for each sample as described by Ladner et al. (2012). Gene expression counts and the reference transcriptome can be accessed at Dryad (doi:10.5061/dryad.hd922). Sequence data can be accessed at NCBI-Bioproject, Accession: PRJNA274410.

Gene Expression Analysis

All statistical analyses were carried out in R (R Development Core Team 2008). Among the 152 gene expression libraries that we sequenced, we analyzed the expression of 15,737 contigs with an average of 5 or greater normalized counts per sample (Seneca and Palumbi 2015). To examine these data for predictors of coral physiological responses to environmental stress, we used weighted gene coexpression network analysis (WGCNA) to identify transcriptional modules of coexpressed genes. This analysis identifies sets of genes that show coordinated variation across samples (Langfelder and Horvath 2008). The expression of these modules can be summarized as the expression of a single “eigengene,” calculated as the first principal component of the expression of all genes in a module across samples (Zhang and Horvath 2005). Genes can also be characterized with respect to their “module membership,” defined as the correlation between a single gene’s expression profile and a specific module’s eigengene (Zhang and Horvath 2005).

We used signed Pearson correlations with a weighting power of 6 and carried out blockwise identification of modules using a block size of 5,000.

We used ANOVA, as implemented in the built-in R package aov, to test the significance of colony, transplant environment, experimental treatment, time point, colony by environment, colony by treatment, treatment by environment, and colony by treatment by environment effects on eigengene expression. We used Pearson product-moment correlation tests, as implemented in the built-in package cor.test, to test for eigengene association with bleaching outcomes for each transplant at the 5 h and 20 h time points. We used DAVID for functional enrichment analyses (Huang et al. 2008). We calculated false discovery rate corrections using the Benjamini–Hochberg algorithm as implemented in the built-in p.adjust function (Benjamini and Hochberg 1995).

We quantified the amount of expression variation explained by module eigengene expression by calculating the sum of variance in normalized read counts across all contigs assigned to a module. We then calculated the sum of the residual variance of normalized read counts for all contigs assigned to a module after accounting for the expression of the module to which a given contig was assigned by WGCNA. We used the built-in R package lm to calculate residual variances. To facilitate the comparison of acclimatory differences in single genes and eigengene expression, we calculated acclimation intensities (the number of phenotypic standard deviations separating mean expression in the two transplant sites) for all contigs (Palumbi et al. 2014). We tested for enrichment of acclimatory genes in acclimatory modules using Fisher’s exact test, as implemented in the built-in R package fisher.test. We calculated the proportion of variance in plasticity explained by eigengene expression by calculating the sum of squares across all contigs explained by transplant location for each contig assigned to a module using the R function aov, and then calculating the residual sum of squares across these contigs explained by transplant location after accounting for eigengene expression.

Results

Variation in Expression across Thousands of Genes Can Be Explained by Variation in a Small Number of Coexpressed Transcriptional Modules

In our experiment, 86% of the analyzed contigs (13,596 of 15,737) were assigned to a transcriptional module that showed coordinated variation in expression across all samples, and 46% of the contigs had a module membership greater than 0.7 in at least one module (fig. 1A and B). To test whether this level of correlation would be expected in our data due to chance alone, we randomized gene expression data among samples in 100 separate trials, and re-ran the analysis on each randomized data set. With randomized data, across 100 trials, on average 39% (range 37–41%) of the genes were assigned to transcriptional modules, and on average only 2% (range 2.3–2.6%) had a module membership greater than 0.7.

F

Gene expression variation across thousands of genes is related to variation in the expression of just a few transcriptional modules. (A) In this pairwise correlation heatmap, each point represents a pairwise signed Pearson correlation value (weighted at power 6; see Results) between two contigs across all samples, where the color denotes the strength and sign of the correlation. We used WGCNA to identify “transcriptional modules” or sets of genes which all show coordinated variation in expression across samples. Contigs are grouped into the modules identified by WGCNA, visible as large magenta blocks of highly coexpressed genes along the 1:1 line. The red and blue squares highlight two such modules, shown in greater detail in panel B. Here we have only shown contigs with a membership >0.7 in at least one module (see Results). (B) Genes in a given module show coordinated variation across samples. In this plot, each x-coordinate represents a sample, and each line represents the expression profile of a contig across all samples, colored by the module to which it was assigned. Samples are ordered by treatment condition, and by bleaching outcomes of the corresponding heat-stressed sample at 20 h within treatment condition. The overall expression of a module of coexpressed genes can be summarized by calculating an expression profile for the module eigengene, the first principal component of the expression profiles of all contigs in a module (see Results), shown as the white line for Module 10 and the black line for Module 12.

We found 23 modules ranging in size from 26 to 2,459 contigs (table 1). Eighteen of the 23 modules are enriched for at least one cellular function, and 7 were correlated with bleaching outcomes in at least one time point (table 1). Eigengene expression across these 23 modules explained 63% of variance in their associated contigs overall.

Table 1

Functional Enrichments and Bleaching Associations for Modules Identified by WGCNA

Module	Size	Bleaching Association		Molecular Function	Cellular Compartment	Biological Process
		5 h	20 h
1	2,459	0.3	0.62*	Transcription factor activity	Endosome	Programed cell death
2	1,430	−0.18	−0.36		Organelle lumen
3	1,368	−0.25	−0.31
4	1,355	−0.33	−0.67*		Intracellular organelle lumen	Ribonucleoprotein complex biogenesis
5	1,141	0.11	0.07
6	1,026	−0.17	0	RNA-directed DNA polymerase activity		RNA-dependent DNA replication
7	1,022	0.03	−0.15		Spliceosome	RNA processing
8	874	0.29	0.46			Positive regulation of molecular function
9	798	0.05	0		Cytosol
10	442	−0.66*	−0.62*	Extracellular matrix structural constituent	Extracellular region	Response to abiotic stimulus
11	384	−0.08	0.13	Calcium ion binding	Plasma membrane	Cellular component morphogenesis
12	277	0.55*	0.67*	Sequence-specific DNA binding		G-protein coupled receptor protein signaling pathway
13	162	−0.44	−0.48		Extracellular region
14	162	0.23	0.61*	Ion binding	Extracellular region
15	153	0.23	0.52*	Motor activity	Microtubule cytoskeleton	Microtubule-based movement
16	138	−0.27	−0.4		Intrinsic to plasma membrane
17	105	0.47	0.64*			Oxidation reduction
18	81	−0.11	−0.13	Structural molecule activity	Proteinaceous extracellular matrix
19	80	0.45	0.1	Structural constituent of ribosome	Ribosome	Translation
20	44	−0.33	−0.35	Sugar binding	Extracellular matrix
21	42	0.33	0.29
22	27	−0.24	0.05
23	26	0.02	−0.48

Note.—Module size denotes the number of contigs assigned to a module by WGCNA. Bleaching association values shown are Pearson correlations, and asterisks denote correlations with Bonferroni corrected P < 0.05. Functional enrichments are GO categories identified by DAVID reported at Benjamini–Hochberg adjusted P < 0.1 (see Results).

Transcriptional Modules Show Constitutive and Acclimatory Differences in Expression between Coral Colonies

Gene expression differences between coral colonies could be due to environmental differences or due to constitutive differences in gene expression between coral colonies with different genotypes (Palumbi et al. 2014). If environmental differences are more important, then gene expression variation among our transplants should be most strongly associated with the transplantation site of a given sample. Alternatively, if constitutive differences between colonies are more important, then transplants derived from different coral colonies should show consistent differences in expression across different transplantation sites. We found that about half (12 of 23) of modules, including the largest module (Module 1), varied significantly from colony to colony independent of transplant environment (ANOVA P < 0.05; supplementary table S1, Supplementary Material online). Fifteen of the 23 modules showed differences in eigengene expression between the HV and MV pools (ANOVA P < 0.05; fig. 2 and supplementary table S1, Supplementary Material online). Seven modules showed significant colony by environment interactions in eigengene expression (ANOVA P < 0.05; supplementary table S1, Supplementary Material online). Twelve modules showed significant colony by treatment by environment interactions, although all these higher order interactions accounted for less than 10% of variation in eigengene expression (ANOVA P < 0.05; supplementary table S1, Supplementary Material online).

F

Acclimation to different environments is accomplished in part by changes in modular expression patterns. We calculated acclimation intensities (the number of phenotypic standard deviations separating mean expression in the two transplant locations) for both single genes and for module eigengenes. Strongly positive acclimation intensities correspond to increased expression in the more stressful HV pool relative to the MV pool. Strongly negative acclimation intensities correspond to decreased expression in the more stressful HV pool relative to the MV pool. Large points represent eigengenes, and small points represent single genes. Module eigengenes with significant acclimatory effects (ANOVA P < 0.05) and the genes within these modules are shown in red for eigengenes that showed increased expression in the HV pool relative to the MV pool and in blue for eigengenes that showed decreased expression in the HV pool relative to the MV pool.

Of 787 contigs in the 5% tail of greatest acclimatory differentiation (i.e., highest acclimatory intensity, see Methods) between environments, a very high fraction was assigned to the most strong acclimatory modules; 479 (61%) of these contigs were assigned to modules 5 or 7, which made up only 14% of the contigs analyzed but showed the largest effects of transplant environment on eigengene expression (supplementary table S1, Supplementary Material online). This enrichment is much greater than that expected by random chance (Fisher’s exact test, P < 10 × 10−16). Eigengene expression accounted for 66% of acclimatory variance across all contigs assigned to modules. Previous analyses have found that many coral genes show acclimatory differences in expression (Palumbi et al. 2014; Bay and Palumbi 2015); here we find that individual genes with the strongest responses to environmental differences are highly concentrated in just a few coexpressed transcriptional modules.

Transcriptional Modules Are Enriched for Different Molecular Functions

We tested the gene complement of each transcriptional module for enrichment for genes involved in different molecular and cellular functions using DAVID. Eighteen of the 23 modules show significant enrichment for at least one cellular function (table 1). For instance, Module 1 was highly enriched for endosomal (Benjamini–Hochberg adjusted P = 3 × 10−7) and apoptotic (Benjamini–Hochberg adjusted P = 0.006) proteins, and was strongly induced under heat stress at the 5 h time point (supplementary fig. S2, Supplementary Material online). Coral colonies showed differences in expression of this module in both control and heat stress treatments (table 1). Module 6 was highly enriched for RNA-directed DNA polymerase activity (Benjamini–Hochberg adjusted P = 6.6 × 10−7) and showed HV expression across samples, but minimal overall differences between treatments (supplementary fig. S2, Supplementary Material online) or colonies (table 1). Closer inspection of the contigs in this module showed that it contained several contigs annotated as transposon (12 contigs) and reverse transcriptase (13 contigs) proteins, suggesting that it represents transposon activity or viral interactions. Module assignments, expression patterns, functional enrichments, and the effects of clonal variation for all 23 modules are summarized in table 1, supplementary figure S2, Supplementary Material online, and supplementary tables S1, S3, and S4, Supplementary Material online.

Two Transcriptional Modules Predict Bleaching Outcomes

Our 5 h time point occurs before A. hyacinthus corals show substantial bleaching. Instead, bleaching occurs only overnight and does not increase later if corals are held at normal temperatures. Nevertheless, the expression of two modules (Modules 10 and 12; table 1) at the 5 h time point was strongly correlated with bleaching outcomes 15 h later (fig. 3). Module 10 was enriched for extracellular matrix proteins (Benjamini–Hochberg adjusted P = 0.0002) and especially collagens (Benjamini–Hochberg adjusted P = 0.0002) and was negatively correlated with bleaching. Module 12 was positively correlated with bleaching and was enriched for sequence-specific DNA-binding proteins such as transcription factors and zinc-finger proteins (Benjamini–Hochberg adjusted P = 0.04; table 1 and supplementary tables S3 and 4, Supplementary Material online). These DNA binding proteins were particularly enriched for E26 transformation-specific (ETS) binding domains (Benjamini–Hochberg adjusted P = 0.01), and included high module membership (>0.7) putative orthologs of the ETS-family transcription factors ERG, SPDEF, ETV6, ETS97D, and ELK1.

F

Two transcriptional modules have an early relationship to bleaching outcomes. Eigengene expression profiles for Modules 10 and 12 in the 5 h heat stress treatment were correlated with variation in bleaching outcomes observed in the corresponding 20 h heat stress samples. Boxplots show median eigengene expression and interquartile range for each level of bleaching observed at 20 h. Points represent eigengene expression levels for individual transplants at 5 h corresponding to a given level of bleaching at 20 h.

Corals transplanted into the more stressful HV pool induced Module 12 less in response to a common heat stress, but colonies also showed constitutive differences in response to heat stress across transplant sites (fig. 4 and supplementary table S1, Supplementary Material online). In contrast, Module 10 expression at the 5 h time point was strongly correlated with bleaching outcomes, but it did not show strong effects of transplant site and did not vary significantly among coral colonies (supplementary table S1, Supplementary Material online). Neither module showed significant effects of colony by transplant environment effects, but both modules showed weak but significant effects of colony by treatment by environment effects (supplementary table S1, Supplementary Material online). High variation in heat stress expression of Module 10 may therefore be explained by higher order interactions between genotype and environment, but further work is needed to discern the nature of control of this key bleaching module.

F

Influences of colony differences and environment on bleaching module expression under heat stress. In this reaction norm plot, each line segment represents 5 h heat stress expression of the Module 12 eigengene for a single coral colony transplanted into our two different transplant sites—the more stressful HV pool and the less stressful MV pool (see Methods). Coral colonies showed constitutive differences in Module 12 eigengene expression regardless of transplant location, but corals transplanted to the MV pool showed consistently higher levels of Module 12 eigengene expression.

Modules 10 and 12 predict bleaching before it is visually apparent, but only after heat stress. As a result, we also sought evidence for gene clusters that could predict bleaching before imposition of heat stress. Such predictive clusters could help identify resilient colonies before bleaching events. However, none of our clusters in control samples was correlated with bleaching after acute heat stress (supplementary table S1, Supplementary Material online).

In addition, we sought evidence for a role of symbiont type in expression and bleaching. Many previous studies have reported an important role for coral symbionts in determining colony-level bleaching resilience (Baker 2001; Baker et al. 2004). Rapid shifts in symbiotic associations from less tolerant to more tolerant symbionts could allow coral to withstand sudden environmental changes (Baker 2001). Our transplantation experiment allowed us to test whether symbiont associations played a major role in determining colony-level differences in bleaching resilience, and whether shifts in symbiotic associations allowed corals to cope with transplantation. Colony-level associations with symbiont clades C and D were stable between transplants from a common colony of origin, and were not significantly associated with variation in bleaching in this experiment (supplementary fig. S5, Supplementary Material online).

Discussion

Coral gene expression variation across thousands of genes can be grouped into differences in the expression of a relatively small number of modular coexpression clusters. The largest module includes genes seen regularly in coral bleaching studies, including proteins involved in apoptosis. Different stress-responsive modules include clusters of genes with more specific cellular functions, including extracellular matrix proteins and transcription factors. These two gene classes were especially common in Modules 10 and 12, respectively. Expression of genes in these modules is highly correlated with colony bleaching that occurred the day after we measured their expression. Thus, change in expression in these modules anticipates coral bleaching.

How changes in these gene expression patterns serve to modulate bleaching is unknown. Transcription factors that dominate Module 12 and increase with bleaching might be responsible for myriad downstream gene expression shifts. Extracellular matrix genes that dominate Module 10 and decrease with bleaching might reduce the connection among coral gastrodermal cells in order to facilitate exocytosis of the symbiont. They might also alter cytoskeletal anchoring of the symbiont within the coral host cell to allow expulsion. Discovery of these modules provides better understanding of the cellular processes of bleaching, and suggests potential bioassays of bleaching likelihood that enhances previous studies of single gene expression patterns. For example, assaying the expression level of the Module 12 ETS-family transcription factors in wild corals experiencing heat anomalies could potentially give insight into the stress level of those corals, and allow wildlife managers to assess how close they are to inducing a bleaching response.

Acclimation and Adaptation via Changes in Modular Gene Expression Patterns

Our data suggest that both acclimatory and genetic variation in bleaching resilience in this species is mediated in part by differences in the induction of bleaching-related gene networks. For most modules, expression varies widely among colonies, suggesting a genetic or persistent epigenetic effect on overall expression. However, because we followed individual coral colonies after they were fragmented and transplanted to more and less thermally stressful microclimates, we could also compare gene expression in the same genotype across different environments. Several modules showed roughly equal contributions of colony differences and acclimation on their expression (supplementary table S1, Supplementary Material online). This is similar to the expression patterns reported by Palumbi et al. (2014) based on single gene analyses. In the present case, the analysis suggests that acclimation is in part accomplished not by independent changes in gene expression among many genes one at a time, but through regulation of coexpressed gene modules.

Like many other modules, Module 12 showed both colony-level and acclimatory differences in eigengene expression. However, the other early bleaching module, Module 10, was an exception to this pattern. Module 10 did not show strong colony or acclimatory effects, but it did show a weak colony by treatment by environment effect. Although these estimated colony and environmental interactions account for very little of the variation in expression in this module (about 1%; supplementary table S1, Supplementary Material online), nevertheless expression of Module 10 genes was negatively associated with bleaching. One possible explanation for this is that Module 10 responds to some environmental factors that vary strongly within pools and not much between them. Depth and current flow are possible candidates, but further work is needed to ascertain the nature of control of the extracellular matrix proteins and other components of Module 10.

The Transcriptional Response of Corals to Heat Stress

The genes that we identified as responding to heat stress in this experiment are broadly similar to the genes reported by other studies of coral responses to heat stress. Apoptosis genes have frequently been associated with the coral heat stress response (Desalvo et al. 2008; Voolstra et al. 2009; Bellantuono et al. 2012). In addition, growth, cell division, and metabolism are also often implicated in transcriptomic studies of coral environmental stress responses (Desalvo et al. 2008; Portune et al. 2010; Kenkel et al. 2013; Bay and Palumbi 2015). A growing number of studies have found that variation in the expression of oxidative stress and extracellular matrix genes may be particularly associated with bleaching resilience in coral hosts (Barshis et al. 2013; Dixon et al. 2015; Seneca and Palumbi 2015). In this experiment, we previously identified genes involved in innate immunity, apoptosis, extracellular matrix formation, and cytoskeletal processes to be highly enriched among stress responsive genes (Seneca and Palumbi 2015).

Where the present systems genetic analysis differs from other single gene approaches is in the way it takes advantage of variation between many different coral colonies in their responses to a common stress, using this variation to group many genes into distinct transcriptional modules. We use this modular structure to discover distinct functional enrichments, environmental responses, and relationships to physiology (in this case, bleaching outcomes) for these distinct sets of genes. Using this approach, we find that the most abundant classes of stress responsive genes (e.g., apoptosis genes) are not the most strongly related to variation in bleaching outcomes. Instead, the expression of other, less abundant classes of genes (e.g., ETS-family transcription factors) is strongly related to differences in coral bleaching within a single population. Other experiments have identified these genes as responsive to heat stress in coral (Desalvo et al. 2008; Polato et al. 2013), but where single gene analyses often overlook the importance of less abundant classes of genes, our systems genetic analysis shows these genes to occupy a place of potentially pivotal importance in coral bleaching gene networks.

System Genetics Applied to Gene Expression in the Wild

Systems genetics approaches have made great progress in identifying gene regulatory networks involved in generating complex phenotypes. They have been particularly useful in better understanding heritable disease, and have revealed underlying causes of physiological and developmental variation between strains of Drosophila and mammals (Ayroles et al. 2009; Brawand et al. 2011). Here we demonstrate the application of these methods to questions about physiological and developmental variation in wild populations of ecologically and economically important nonmodel organisms (Filteau et al. 2013). By identifying physiological and developmental units in the form of transcriptional modules first, and characterizing their functional enrichments and relationships to genetic, environmental, and physiological factors later, we can gain insights into the potential roles of genes with variable expression and can better understand the regulatory architecture of complex traits like coral bleaching resilience.

The presence of a module made up of highly expressed transposon-related transcripts with highly variable expression between samples was an unexpected source of variation in this population. The activation of silenced transposons by environmental stress has been observed across diverse array of organisms, including corals (Grandbastien 1998; Capy et al. 2000; de la Vega et al. 2007; Beauregard et al. 2008; Desalvo et al. 2008). This activation can disrupt cellular function, but the mutagenic effects of transposition can also generate novel genetic variation, including novel stress-induced gene expression patterns, potentially facilitating adaptation to stressful environments (Beauregard et al. 2008). A deeper understanding of the causes and consequences of transposon activation in coral populations might clarify the relative importance of different effects of transposon activation in shaping population responses to environmental change.

Supplementary Material

Supplementary figures S2 and S5 and tables S1, S3, S4, and S6 are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/).