s41467-022-32217-z.pdf

nature communications Article https://doi.org/10.1038/s41467-022-32217-z T h er o l eo fg e n ee x p r e s s i o na n ds y m b i o s i si n reef-building coral acquired heat tolerance Marie E. Strader 1,3 &K a t eM .Q u i g l e y2,3 Predicting how reef-building corals will respond to accelerating ocean warm- ing caused by climate change requires knowledge of how acclimation and symbiosis modulate heat tolerance in coral early life-history stages. We assayed transcriptional responses to heat in larvae and juveniles of 11 repro- ductive crosses ofAcropora tenuiscolonies along the Great Barrier Reef. Larvae produced from the warmest reef had the highest heat tolerance, although gene expression responses to heat were largely conserved by cross identity. Juvenile transcriptional responses were driven strongly by symbiosis– when in symbiosis with heat-evolved Symbiodiniaceae, hosts displayed intermediate expression between its progenitorCladocopium and the more stress tolerant Durusdinium, indicating the acquisition of tolerance is a conserved evolu- tionary process in symbionts. Heat-evolved Symbiodiniaceae facilitated juve- nile survival under heat stress, althoughhost transcriptional responses to heat were positively correlated among those hosting different genera of Symbio- diniaceae. These ﬁndings reveal the relative contribution of parental envir- onmental history as well as symbiosis establishment in coral molecular responses to heat in early life-history stages. Ocean temperatures are rising beyond the physiological limits of many marine taxa, driving rapid responses of organisms to counteract extinction either through range shifts, acclimation, or adaptation. Coral reefs globally are experiencing more frequent and extreme heat stress anomalies, which can induce short-term mortality and diminish recovery times, resulting in extreme reductions in population sizes. These changes drive community shifts to alternative stable states devoid of hard, reef-building corals and limit suitable substrate for coral recruitment, creating a negative feedback loop that further threatens coral reefs 1,2. Despite global population declines, thermally tolerant genotypes exist3,4, which may facilitate the rapid evolution of heat tolerance for applied conservation practices like assisted evolution 5. Genetic and environmental factors inﬂuence the adaptive and acclimation potential of organisms in response to global envir- onmental change. The relative roles of these factors depend on the environment experienced currently or in previous generations through either past selection pressures or acclimation 6.I nc o r a l sw i t h large population sizes spanning extensive ranges across temperature gradients, as with some species of Acropora in the Indo-Paci ﬁc, the potential for acclimation and adaptation to higher temperatures indicative of some future warming scenarios likely exists 7,8. However, the full picture of how each of these processes contributes to the future persistence of coral populations is yet to be fully understood, especially when the compounding effects of symbiosis on host phy- siology and responses to thermal stress in early life-history stages are taken into account. Genetic interventions to increase coral thermal tolerance are being considered worldwide to allow time for climate change action and mitigation. These methods include facilitating rapid acclimation, modifying symbiont communities (eukaryotic and prokaryotic), selectively breeding corals, and experimentally evolving their symbionts 5. The utilization of assisted gene ﬂow methods, a form of selective breeding, across temperature gradients on the Great Barrier Reef (GBR) can increase heat tolerance in corals 3,9 without severe trade-offs in growth, survival, and symbiont acquisition 10.F u r t h e r - more, when multiple approaches are combined, there is a Received: 17 June 2022 Accepted: 19 July 2022 Check for updates 1Department of Biological Sciences, Auburn University, Auburn, AL, USA.2Minderoo Foundation, Perth, WA, Australia.3These authors contributed equally: M a r i eE .S t r a d e r ,K a t eM .Q u i g l e y .e-mail: [email protected]; [email protected] Nature Communications| (2022) 13:4513 1 1234567890():,; 1234567890():,; synergistically positive increase in heat tolerance 9 resulting from increased adaptive genetic diversity in regions associated with heat tolerance 4. In this study, we use several assisted evolution methods synergistically to tease apart the role that gene expression plasticity (acclimation), and microbial community manipulation (symbiosis) have on the acquired heat tolerance of early life-history stages of a wide-spread, habitat-building, Indo-Paciﬁc coral species. Colonies of Acropora tenuis were selectively bred from ﬁve reef sites stretching over 6° of latitude on the GBR, resulting in 11 distinct crosses, including purebred and hybrid gametic crosses (Supplemental Fig. 1 but see 11 for full information on collection, coral spawning, reproduction, and cross design). Of these sites, three were sourced from the far North GBR (Curd (CU), Long Sandy (LS), and Sandbank 7 (SB)), with Curd Reef experiencing the highest mean sea surface temperature (SST) and two sites from the central GBR (Davies Reef (DR) and Backnumbers (BK)) (Supplemental Fig. 1). To quantify acclimation effects contributing to early life-history survival under heat stress, we characterized gene expression of both larvae and juveniles assayed under ambient con- ditions, heat stress, and under different symbiotic states. Transcrip- tional responses are a key determinant to survival and recovery during and after heat stress 12 and corals employ multiple transcriptional strategies in response to environmental change 13– 16.W ei d e n t i f yk e y genes associated with high larval survival under heat stress, which was speciﬁc to crosses produced from dams originating from Curd Reef, but report an overall conserved transcriptional response to heat across all crosses. Alternatively, host transcriptional responses to heat are functionally different in juveniles and are strongly associated with the symbiotic community. Juvenile symbiosis with a heat-evolved Clado- copium strain has host transcriptional pro ﬁles similar to juveniles hosting Durusdinium, and contributed to improved juvenile survival under heat stress. These results have important implications for understanding the molecular mechanisms behind heat tolerance acquisition and for the manipulation of symbiont communities and selective breeding in coral restoration. Results Larvae display conserved gene expression response to heat stress Larval gene expression (GE) was quanti ﬁed to assess if plastic responses in gene expression to heat stress varied depending on site of origin or parental identity. Larval survival under heat stress varied between crosses, with larvae produced from dams sourced from far Northern GBR sites exhibiting higher thermal tolerance (Fig. 1b). The cross with the lowest thermal tolerance (LSxSB) did not have any larvae survive the heat treatment (Fig. 1b, Supplementary Fig. 2). GE was examined in aposymbiotic larvae experiencing ambient conditions prior to the heat treatment (“pre”), larvae after exposure to simulated heat stress (35.5 °C for 56 hours, “post heat ”), and a simultaneous control temperature of 27 °C ( “post ambient”). Therefore, the “pre” larval treatment provided transcriptomic baselines of GE between genetic crosses while “post heat ” and “post ambient” comparisons show a baseline for cross-speci ﬁc heat responses without the con- founding effect of symbiosis found in the post-metamorphic phase. Using a principal coordinates analysis (PCoA), weﬁnd that GE patterns Fig. 1 | Larval gene expression responses to heat stress in genetically distinct crosses. a Discriminant analysis of principal components (DAPC) with heat treatment symbolized by color and cross ID symbolized by shape. The ﬁrst linear discriminant axis separates the heat treatments, 27 °C for “pre” (n =3 3 pools of 10 larvae each) and “post ambient ” (n = 33 pools of 10 larvae each), 35.5 °C for “post heat ” (n = 30 pools of 10 larvae each), with density plots showing the distribution of samples in each cross and treatment along LD1, ranked from low to high survival in heat stress. b Larval survival in heat and ambient treatments. c Eigengene expression for two WGCNA modules posi- tively, and signi ﬁcantly correlated with larval response to heat. For each cross, the origin of dams is denoted ﬁrst, and the origin of sires is second ( n = 3 pools of 10 larvae each per cross per treatment). Boxplots include the median values (center lines), upper and lower quartiles (box limits), 1.5× interquartile range (whiskers), and outliers (points). Article https://doi.org/10.1038/s41467-022-32217-z Nature Communications| (2022) 13:4513 2 in larvae were driven by treatment (“pre”, “post ambient”, “post heat”), explaining 29.2% of the variation in survival ( padonis < 0.001, Fig. 1a, Supplemental Fig. 2). Discriminant function of principal components (DAPC), a superior method compared to PCoA for partitioning the variance between populations 17, recapitulates PCoA patterns showing both treatment and cross-speciﬁc effects on larval GE. Using DAPC, we observe that“pre” treatment GE often lies midway between both of the post-treatment GE responses along DAPC-1 (Fig.1a). Weighted gene co- expression network analysis18 identiﬁed groups of genes (the grey60 module, 2069 genes, and the blue module, 1737 genes), which were positively correlated with larval heat treatment ( R 2 =0 . 8 4 ,p =3 e - 2 6 , and R2 =0 . 9 4 ,p = 2e-44, respectively.) (Fig.1c, Supplementary Fig. 3), but did not strongly correlate with larval survival under heat stress, which varied by cross (Fig. 1, Supplemental Fig. 3). This indicates a relatively conserved gene expression response to heat between crosses, despite differences in larval survival. The blue and grey60 modules also showed negative and positive correlations, respectively, with survival under ambient conditions (Supplementary Fig. 3). Gene Ontology (GO) categories enriched in genes within the blue module include “regulation of response to stress” and “nf-kappaB signaling”, while the grey60 module has GO enrichment of “protein folding” and “response to endoplasmic reticulum stress”, categories which are also enriched under a range of environmental stress responses 16 and implicated in the selective breeding for heat tolerance 4, suggesting their critical importance in a heritable response to heat tolerance acquisition (Supplementary Data 10). This cross-site breeding design allowed us to further clarify the impact of the local reef, and thus the environmental history, on the molecular responses to thermal stress in coral offspring. Cross identity explained 29.9% of the variation in larval GE (Supplemental Fig. 2, p adonis < 0.001), suggesting host genetics plays a vital role in larval heat tolerance. Larvae produced from maternal colonies from the warmest site, Curd Reef, which also exhibited high survival to heat stress (CUxBK, survival 81.25% ± 4.2 11,F i g .1b), demonstrated unique larval GE responses under elevated temperature. In the CUxBK cross, we observed an upregulation of thermal response genes (genes in the grey60 module) under“pre” and “post ambient” conditions, indicative of genetic assimilation (Fig. 1c) 19, a type of response that has been described in heat-tolerant corals sourced from shallow lagoons in other reef locations 20. A cross with the paternal contribution from Curd (LSxCU) exhibited intermediate GE responses to heat that were between crosses with Curd dams and the other crosses produced from cooler reefs (Fig. 1a), further indicating that individual colonies from Curd reef may be differentially adapted to heat stress given their northern, inshore location. Speci ﬁcally, the expression of 11 genes differentiated larvae with Curd dams in their responses to heat, including an nf-Kappa-B inhibitor (Supplementary Fig. 4). Additionally, modeling in DESeq2 of the “pre heat ” larval samples with survival under subsequent heat stress identiﬁed 438 genes predictive of larval survival in A. tenuis (Supplementary Data 8), which are further refer- enced throughout as “thermal tolerance” genes. GO categories “cell cycle processes” and “ribosome” were signiﬁcantly enriched in genes that were predictive of larval survival (Supplementary Data 5). The association of these genes in crosses produced from parents of Curd origin suggests that these genes may play a role in the local adaptation to the more extreme, inshore, warmer conditions experienced by this reef. Juvenile gene expression varies based on symbiont community While mechanisms elucidating acquired heat tolerance have been examined in aposymbiotic larvae 3,4 (Fig. 1,S u p p l e m e n t a r yF i g s .2– 4), no studies have yet resolved GE responses to heat or variable sym- biotic states in juvenile corals. Arguably, quantifying the role of population origin, plasticity, and symbiosis in this critical life-history stage is urgently needed as juvenile survival and sensitivity to envir- onmental stress are key for the conservation and restoration of coral populations globally. To address this, larvae from the aforementioned crosses were induced to metamorphose and then subsequently exposed to a range of cultured or wild Symbiodiniaceae treatments. Juveniles were infected with cultured Cladocopium goreaui (“C1”), Durusdinium trenchii (“D1a”), a heat-evolved species of Cladocopium goreaui (“SS1”) 21, or a diverse community of putatively heat-tolerant symbionts from the sediments collected from Curd (“SED”), with the dominant infection of each symbiont type conﬁrmed with transcripts mapping to each genus of Symbiodiniaceae (Supplementary Fig. 6). Juvenile Acropora uptake diverse communities of Symbiodiniaceae from their environment, a dynamic process comprised of a tax- onomically variable symbiont consortium that plays a substantial role in heat tolerance acquisition and survival in juveniles 10,22.A ta m b i e n t temperatures, juvenile host GE is mainly clustered by symbiont iden- tity although a signiﬁcant effect of the cross was also observed (Fig.2). A permutational multivariate ANOVA revealed that the cross explained 15.7% of the variance in juvenile GE under ambient conditions with symbiont treatment explaining 10.6%. In partitioning the variance using DAPC, we ﬁnd that two crosses (DRxDR and BKxBK) showed unique GE clustering (Fig.2) and despite close proximity of these two Fig. 2 | Juvenile gene expression responses to symbiosis.Discriminant analysis of principal components (DAPC) on gene expression of juveniles hosting different symbioses under ambient conditions. DAPC was performed on variance stabilized data (VSD) grouped by cross+symbiosis (a) and with the effect of cross removed by limma (b). For each cross, origin of dams is denoted ﬁrst, and the origin of sires is second. Article https://doi.org/10.1038/s41467-022-32217-z Nature Communications| (2022) 13:4513 3 sites on the offshore GBR, appear to have unique histories of selection and acclimation. Juveniles in the “SED” treatment hosted a highly diverse symbiont community, comprised of >300 taxonomic groups11, which may explain the pattern of high GE divergence in these juveniles compared to those hosting cultured symbionts. This emphasizes the inﬂuence of symbiont diversity in the modulation of host physiology, even during early ontogeny (Fig. 3, Supplementary Fig. 6). Moreover, juveniles in the ambient“SED” treatment exhibited the downregulation of genes associated with signal transduction, coenzyme transport, and secondary metabolite synthesis, and upregulation of energy produc- tion and conversion, and intracellular trafﬁcking, secretion, and vesi- cular transport (Fig.4a, Supplementary Data 9), compared to juveniles with cultured symbionts. These functional categories are a hallmark of host regulation of the symbiotic state 23, likely in the midst of the winnowing process24,25, and are most likely representative of juvenile host gene expression in situ on the reef. In juveniles, gene modulation in response to different symbiotic states included changes in genes associated with cell receptor signaling, cellular adhesion proteins, transcription, cellular stress responses, immunity, and the cell cycle (Fig. 4a, Supplementary Data 6). Juveniles hosting the SS1 strain, which was experimentally evolved to tolerate higher temperatures, upregu- lated genes associated with biological adhesion molecules, and reg- ulation of cell morphogenesis (Supplementary Data 6). These functional categories of genes have been previously implicated as signatures of selection and convergence in the evolution of vertical transmission in corals 26, suggesting that these pathways are speciﬁcally modulated by the host when the symbioses are perturbed. Larvae and juveniles have different responses to heat stress Under heat stress, juveniles hosting C1 experienced substantial mor- tality and are thus not represented in this study, while juveniles hosting SS1, experimentally evolved from C1, exhibited higher survival under heat stress. Survival of juveniles hosting SS1, D1a, and communities from the sediments all exhibited similarly high survival under heat stress 11 as well as signi ﬁcantly positively correlated functional responses to heat (Fig. 4d). While juveniles hosting SS1 at ambient conditions hosted a majority ofCladocopiumtranscripts, SS1 juveniles under heat stress all exhibited elevated proportions of Durusdinium- derived transcripts (Supplementary Fig 6). Therefore, the heat stress itself may have promoted a shift towards Durusdinium-derived tran- scripts in SS1 infected juveniles, potentially facilitating their thermal tolerance. Further, juveniles hosting SS1 showed positively signiﬁcant molecular responses to heat that are indicative of corals hosting D1a as seen in the KOG correlations between these two responses ( R 2 = 0.74 p = 8.5e-05, Fig. 4d), with stronger correlations to those hosting com- munities from the sediments ( R2 =0 . 8 p = 9.2e-06, Fig. 4d). Despite these overall functional similarities, key differences were observed that highlight the unique interactions between symbiont identity, symbiont community, and coral host origin that contribute to juvenile heat stress responses. For example, under heat stress, SS1 juveniles signiﬁcantly downregulated cell cycle control, cell division, and chro- mosome partitioning genes compared to the other two symbioses (Fig. 4d), indicative of a trade-off between growth and thermal toler- ance. Compared to the other symbioses, D1a juveniles downregulated signal transduction genes (Fig. 4d), genes likely to be involved in symbiosis establishment and maintenance. Heat-stress responses in juveniles hosting all symbiont treatments exhibited non-signiﬁcant correlations with the heat stress responses in larvae (Fig.4d). This underscores that heat tolerance is likely regulated by different molecular mechanisms in larvae and juvenile Acropora tenuis. Juveniles in all symbiont treatments showed strong upregula- tion of energy production and conversion genes compared to the larval heat stress response (Fig.4c). Larval heat tolerance also showed no functional correlations with either larval or juvenile responses to stress, suggesting different molecular mechanisms confer tolerance and response in early life-history stages. Discussion Heat tolerance may be regulated via several divergent evolutionary pathways, including genetic ﬁxed effects or via changes in gene expression 7,27,28. These mechanisms are likely linked, where the inter- play between geneticﬁxed effects and gene expression can be broadly connected by mechanisms like expression QTLs (reviewed in ref. 27). Here we use a unique experimental design to quantify the synergistic effects of acquired heat tolerance via gene expression plasticity and symbiosis in early life-history stages. We pro ﬁled gene expression in both aposymbiotic larval and symbiotic juvenile stages at pre-stress, ambient, and heat-stress treatments across 11 distinct genetic crosses Fig. 3 | Juvenile gene expression responses to symbiosis and heat stress. Dis- criminant analysis of principal components (DAPC) on gene expression of juveniles hosting different symbioses under ambient and heat conditions. DAPC was performed on variance stabilized data (VSD) grouped by cross+treatment+sym- biosis (a) and with the effect of the cross removed by limma (b) with density plots showing the distribution of samples in each symbiosis and treatment along LD1. Article https://doi.org/10.1038/s41467-022-32217-z Nature Communications| (2022) 13:4513 4 between and amongst northern and central GBR reef sites ofAcropora tenuis. We found that selective breeding between disparate sites revealed relatively conserved gene expression proﬁles underlying key differences in larval and juvenile heat tolerance between crosses, although mechanisms of heat tolerance differed substantially between these disparate life-history stages. Juveniles hosting SS1 under ambient conditions exhibited an intermediate heat response situated between D1a and their progenitor taxa C1 (Supplementary Fig. 5). These results suggest that the evolution of heat tolerance in the SS1 strain likely follows a conserved, or constrained, evolutionary pattern. In Symbio- diniaceae, given the recapitulation of molecular host responses to SS1 and D1a, but “residual” or ancestral similarity of SS1 to C1, it is likely that further experimental heat selection may be needed to completely convert the SS1 capacity for heat tolerance that is in line with the canonically stress-tolerant D1a. Together these data quantify, simul- taneously, the relative contribution of gene expression and symbiosis in critical early life-history stages and provide an applied foundation for optimizing genetics-based methods for the enhancement of coral survival. Methods Coral spawning, larval and juvenile experiments Reproductive Acropora tenuis coral colonies were collected across 6° of latitude and were spawned at the National Sea Simulator at the Australian Institute of Marine Science (AIMS) to produce purebred and hybrid crosses. The full details of reef locations (Supplemental Fig. 1), coral colony collections, spawning, and reproductive crossing can be found in 11,29. Larvae produced from these crosses were sampled for RNAseq (“pre” samples) prior to exposure at 0 hrs (directly before the heat ramp started) and cohorts of those larval cultures were subjected to control and heat stress for 36 hours (27 and 35.5 °C) in replicates of n =6( “post” exposure samples). Separate cohorts of each larval cross not used in heat trials were then induced to settle using sterilized crustose coralline algae as a settlement cue, exposed to the four symbiont treatments for symbiosis infection, and then exposed to control and heat stress conditions for 56 days (27 °C and 32 °C). The three cultured symbionts were grown under LED, providing illumina- tion set at 60 m −2 s−1 photosynthetically active radiation measured with a Li-Cor light sensor for a 14:10 light:dark photoperiod. Before inocu- lation, cells were measured using a Neubauer Haemocytometer (n =3 replicates/symbiont) to make up a ﬁnal cell density of 1 × 10 5 ml−1. RNA extraction and tag-seq All larvae and juveniles were preserved in RNAlater and stored at −20 °C until library preparation. Larvae and juveniles were both pre- served at the ambient and heat treatments once the heat treatment reached ~50% mortality. Ten pooled larvae and two juveniles were extracted per cross/treatment due to low RNA concentrations per Fig. 4 | Functional similarities in juvenile responses to hosting different sym- bionts and exposed to heat stress. Hierarchical clustering analysis of KOG enrichments in juveniles hosting different symbioses at ambient conditions (a)a n d host responses to heat in juveniles hosting different symbioses (c). Correlations of KOG delta rank values for each comparison under ambient conditions (b)a n d symbiosis responses to heat (d). KOG categories in bolded squares denote the signiﬁcance of Mann– Whitney U test at padj < 0.05. Exact Benjamini-Hochberg adjusted p values for each comparison are in Supplementary Data 9. Color scale in a and c are KOG delta rank values. Article https://doi.org/10.1038/s41467-022-32217-z Nature Communications| (2022) 13:4513 5 larva and juvenile. Extractions were performed using an RNAqueous kit (Ambion) and homogenized with an additional bead-beating step in 100 µl Lysis solution with ~6 glass beads (MPBio lysing matrix C) on the FastPrep (5.5 m/s for 30 sec × 2 times). Total RNA quantity was nor- malized to 375 ng for larval samples and 30 ng for juvenile samples using an RNA Pico Chip kit on the Agilent 2100 Bioanalyzer. Due to the high mortality of juveniles infected with C1 in the heat treatment, there was not adequate RNA to sequence these individuals. Samples were then selected and normalized based on RNA concentration and RIN (15 ng/µl for larvae, mean RIN = 6.6 and 2.5 ng/ µl for juveniles, mean RIN = 6.9) and sent to Genomic Sequencing and Analysis Facility at the University of Texas at Austin for custom library preparation and SR50 tag-based sequencing on the Illumina Hi-Seq 2500. Bioinformatics and statistical analysis Read processing of RNAseq reads employed established wrapper scripts for tag-seq ( https://github.com/z0on/tag-based_RNAseq) 30.I n summary, reads were trimmed of adaptor sequences, PCR duplicates were removed andﬁltered for quality using cutadapt (v.2.10) 31.F o rt h e larval dataset, sequencing yielded an average of 4,138,973 raw reads per sample (min: 351,509, max: 5,369,230, median: 4,133,164), which was reduced to a mean of 2,643,781 reads per sample after trimming (min: 247,979, max: 3,332,923, medi an: 2,686,320). For the juvenile dataset, sequencing yielded an average of 3,648,218 raw reads per sample (min: 1,688,621, max: 5,973,561, median: 3,619,747.5), which was reduced to an average of 1,695,570 reads per sample after trim- ming (min: 705,818, max: 3,285,162, median: 1,673,007.5) (Supple- mentary Data 1 and 2). Although recent work in corals using similar approaches has yielded deeper read depths (>35 M; see reference32), the tag-seq approach used here only sequences the 3 ’ end of the transcripts, resulting in substantially fewer reads required to assess gene expression, as the full assembly of full-length transcripts is not needed 33. To verify the appropriateness of mapping references, sequences were mapped to both the conspeciﬁc( A. tenuis34)a n dh e t - erospeciﬁc( A. millepora35) references. Overall, the results were similar, but the A. tenuis genome was selected based on the higher mapping efﬁciency (Supplementary Data 1 and 2). Cleaned reads were mapped to transcripts mined from the Acropora tenuis genome using Bowtie2 (v.2.2.9)36. The mapping ef ﬁciency of larval and juvenile samples averaged 63.4% and 54.4%, respectively (Supplementary Data 1 and 2). Gene counts were compiled using a custom perl script (samcounts.pl). Gene expression was analyzed using DESeq2 (v.1.30.1), implemented in R( v . 4 . 0 . 4 ) 37. Only highly abundant genes with >10 mean counts across samples were retained for statistical analysis: 11,192 genes for the larval data set, 9311 genes for the juvenile dataset. For the larval dataset, differentially expressed genes (DEGs) were identi ﬁed using a model that deﬁned ﬁxed factors of treatment ( “pre”, “post ambient”, “post heat”) and cross. Juvenile gene expression was analyzed in two dif- ferent ways. First, only samples under ambient conditions were con- sidered and a model was run specifyingﬁxed effects of symbiosis and cross. Second, all samples including juveniles exposed to heat were modeled with ﬁxed effects of treatment, symbiosis, and cross. Sig- niﬁcant DEGs were deﬁned by a false discovery rate adjustedp value of <0.05 (Supplementary Data 3 and 4). Stat values from each compar- ison were used for GO enrichment analysis using the GO_MWU pack- age (v.1.2) 38. This test involves a two-sided Mann – Whitney U test to identify functional enrichment of GO terms without having to set arbitrary cut-offs. Variance stabilized counts data from each model were output and used for additional multivariate analyses. Variance partitioning was determined using a principal coordinates analysis (PCoA) of Man- hattan distances and the signiﬁcance of each ﬁxed effect was assessed using the adonis function in R. Additionally, the variance associated with each ﬁxed effect was analyzed using a DAPC using the package adegenet (v.2.1.4) 39. Weighted gene co-expression network analysis (WGCNA18, v.1.69), implemented in R (v.3.6.1), was used to further identify groups, or modules, of genes whose expression signi ﬁcantly correlates with known factors in the larval dataset. A signed network was constructed using a soft threshold power of 18, a minimum module size of 30, and a module merging threshold of 10% dissimilarity. Module trait relation- ships are calculated using the cor() function in R, utilizing a Pearson correlation. GO enrichment analysis on gene modules was performed using Fisher’s exact test for genes present within the module or not using the GOMWU framework 38. For the juvenile dataset, a KOG enrichment analysis was performed. This analysis bins genes into 23 broad functional groups, reducing the complexity of functional annotations. The KOG_MWU function calculates delta rank values for these 23 broad functional groups using gene-level log2-fold changes generated from pairwise comparisons between treatments in DESeq2 40. For juveniles under ambient conditions, comparisons were between juveniles with a speciﬁcs y m b i o n tt y p e( C 1 ,D 1 a ,S S 1 ,o rS E D ) against all other samples at ambient (Fig. 2b). For the comprehensive juvenile dataset, heat responses were examined using pairwise com- parisons between ambient and heat for each symbiont type (C1 ambient vs. C1 heat, for example). As all the SS1 heat samples contained noticeable proportions ofDurusdiniumtranscripts, a comparison was performed 1) retaining all six of these samples and 2) removing those samples with approximately half of the transcripts derived from Dur- usdinium. Both comparisons yielded similar results (Supplementary Fig. 7). Additional comparisons were included in this analysis, includ- ing larval thermal tolerance ( “pre” samples modeled against survival under heat stress) and larval thermal response ( “post ambient” vs. “post heat”). Reporting summary Further information on research design is available in the Nature Research Reporting Summary linked to this article. Data availability ITS2 sequences are available at NCBI SRA at PRJNA720058 and RNAseq ﬁles are accessible through the NCBI GEO repository at GSE176051. Publicly available data used for our analysis include theAcropora tenuis genome: http://aten.reefgenomics.org/;a n dt h e Acropora millepora genome: https://www.ncbi.nlm.nih.gov/bioproject/767661. Code availability Analysis scripts are located at: https://github.com/LaserKate/AGF18_ RNAseq under Zenodo https://doi.org/10.5281/zenodo.682000741 (https://zenodo.org/badge/latestdoi/355383542). References

Schmitt, R. J., Holbrook, S. J., Davis, S. L., Brooks, A. J. & Adam, T. C. Experimental support for alternative attractors on coral reefs.Proc. Natl Acad. Sci. 116,4 3 7 2–4381 (2019).
Hughes, T. P. et al. Global warming transforms coral reef assem- blages. Nature 556,4 9 2–496 (2018).
Dixon, G. B. et al. Genomic determinants of coral heat tolerance across latitudes.Science 348,1 4 6 0–1462 (2015).
Quigley, K. M., Bay, L. K. & van Oppen, M. J. J. H. Genome-wide SNP analysis reveals an increase in adaptive genetic variation through selective breeding of coral.Mol. Ecol. 29,2 1 7 6–2188 (2020).
van Oppen, M. J. H., Oliver, J. K., Putnam, H. M. & Gates, R. D. Building coral reef resilience through assisted evolution.Proc. Natl Acad. Sci. 112,2 3 0 7–2313 (2015). 6 . G i e n a p p ,P . ,T e p l i t s k y ,C . ,A l h o ,J .S . ,M i l l s ,J .A .&M e r i l ä ,J .C l i m a t e change and evolution: disentangling environmental and genetic responses.Mol. Ecol. 17,1 6 7–178 (2008).
Torda, G. et al. Rapid adaptive responses to climate change in corals. Nat. Clim. Change 7,6 2 7–636 (2017). Article https://doi.org/10.1038/s41467-022-32217-z Nature Communications| (2022) 13:4513 6 8 . M a t z ,M .V . ,T r e m l ,E .A . ,A g l y a m o v a ,G .V . ,v a nO p p e n ,M .J .H .& Bay, L. K. Potential for rapid genetic adaptation to warming in a Great Barrier Reef coral. Plos Genet. 14, e1007220 (2018). 9 . Q u i g l e y ,K .M . ,R a n d a l l ,C .J . ,v a nO p p e n ,M .J .H .&B a y ,L .K . Assessing the role of historical temperature regime and algal symbionts on the heat tolerance of coral juveniles.Biol. Open 9, bio047316 (2020). 1 0 . Q u i g l e y ,K .M . ,W i l l i s ,B .L .&B a y ,L .K .M a t e r n a le f f e c t sa n dSym- biodinium community composition drivedifferential patterns in juvenile survival in the coral Acropora tenuis. R. Soc. open sci. 3, 160471 (2016).
Quigley, K. M. & van Oppen, M. J. H. Predictive models for the selection of thermally tolerant corals based on offspring survival. Nat. Commun. 13, 1543 (2022). 1 2 . J o v i c ,K .e ta l .T r a n s c r i p t o m er e silience predicts thermotolerance in Caenorhabditis elegans. BMC Biol. 17, 102 (2019).
Studivan, M. S. & Voss, J. D. Transcriptomic plasticity of mesophotic corals among natural populations and transplants ofMontastraea cavernosa i nt h eG u l fo fM e x i c oa n dB e l i z e .Mol. Ecol. 29, 2399–2415 (2020).
Cunning, R. & Baker, A. C. Thermotolerant coral symbionts mod- ulate heat stress‐responsive genes in their hosts.Mol. Ecol. 29, 2940–2950 (2020).
Kenkel, C. & Matz, M. V. Enhanced gene expression plasticity as a mechanism of adaptation to a variable environment in a reef- building coral. Nat. Ecol. Evol. 1,1 4( 2 0 1 6 ) . 1 6 . D i x o n ,G . ,A b b o t t ,E .&M a t z ,M .Meta-analysis of the coral envir- onmental stress response: Acropora corals show opposing responses depending on stress intensity.Mol. Ecol. 29, 2855–2870 (2020).
Jombart, T., Devillard, S. & Balloux, F. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations.BMC Genet. 11, 94 (2010).
Langfelder, P. & Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9,5 5 9 (2008).
Kelly, M. Adaptation to climate change through genetic accom- modation and assimilation of plastic phenotypes.P r o c .R .S o c .B 375,( 2 0 1 9 ) . 2 0 . B a r s h i s ,D .J .e ta l .G e n o m i cb a s i sf o rc o r a lr e s i l i e n c et oc l i m a t e change. Proc. Natl Acad. Sci. USA 110,1 3 8 7–1392 (2012).
Quigley, K. M., Alvarez Roa, C., Beltran, V. H., Leggat, B. & Willis, B. L. Experimental evolution of the coral algal endosymbiont, Cladocopium goreaui : lessons learnt across a decade of stress experiments to enhance coral heat tolerance. Restor. Ecol. 29, e13342 (2021).
Quigley, K. M., Willis, B. L. & Bay, L. K. Heritability of the Symbiodi- nium community in vertically- and broadcast spawning corals.Sci. Rep. 7, 8219 (2017). 2 3 . D a v y ,S .K . ,A l l e m a n d ,D .&W e i s ,V .M .C e l lb i o l o g yo fc n i d a r i a n - dinoﬂagellate symbiosis.Microbiol. Mol. Biology. Rev.76, 229–261 (2012).
Abrego, D., Willis, B. L. & van Oppen, M. J. H. Impact of light and temperature on the uptake of algal symbionts by coral juveniles. PLoS One 7,e 5 0 3 1 1( 2 0 1 2 ) . 2 5 . Q u i g l e y ,K .M . ,A l v a r e zR o a ,C . ,T o r d a ,G . ,B o u r n e ,D .G .&W i l l i s ,B .L . Co-dynamics of Symbiodiniaceaeand bacterial populations during the ﬁrst year of symbiosis withAcropora tenuisjuveniles.Microbiol. Open 9, e959 (2020). 2 6 . D i x o n ,G .B .&K e n k e l ,C .D .M o l e c u l a rc o n v e r g e n c ea n dp o s i t i v e selection associated with the evolution of symbiont transmis- sion mode in stony corals. P r o c .R .S o c .B286, 20190111 (2019).
Thomas, L. et al. Mechanisms of thermal tolerance in reef-building corals across a ﬁne-grained environmental mosaic: lessons from Ofu, American Samoa. Front Mar. Sci 4,( 2 0 1 8 ) .
Drury, C. Resilience in reef-building corals: the ecological and evolutionary importance of the host response to thermal stress. Mol. Ecol. 29, 448–465 (2020). 2 9 . Q u i g l e y ,K .M .e ta l .V a r i a b i l i t yi nﬁtness trade-offs amongst coral juveniles with mixed genetic backgrounds held in the wild. Front. Mar. Sci. 8, (2021).
Meyer, E., Aglyamova, G. V. & Matz, M. V. Proﬁling gene expression responses of coral larvae (Acropora millepora) to elevated tem- perature and settlement inducers using a novel RNA-Seq proce- dure. Mol. Ecol. 20,3 5 9 9–616 (2011).
Martin, M. Cutadapt removes adapter sequences from high- throughput sequencing reads.EMBnet J. 17, 10 (2011).
Voolstra, C. R. et al. Contrasting heat stress response patterns of coral holobionts across the Red Sea suggest distinct mechanisms of thermal tolerance. Mol. Ecol. 30, 4466–4480 (2021).
Lohman, B. K., Weber, J. N. & Bolnick, D. I. Evaluation of TagSeq, a reliable low-cost alternative for RNAseq.Mol. Ecol. Res 16, 1315–1321 (2016).
Cooke, I. et al. Genomic signatures in the coral holobiont reveal host adaptations driven by Holocene climate change and reef speciﬁc symbionts. Sci. Adv. 6, eabc6318 (2020).
Fuller, Z. L. et al. Population genetics of the coral Acropora mill- epora: Towards a genomic predictor of bleaching.Science 4674, 867754 (2020).
Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9,3 5 7–359 (2012).
Love, M. I., Anders, S. & Huber, W. Differential analysis of count data

the DESeq2 package. https://doi.org/10.1101/002832(2014). 3 8 . W r i g h t ,R .M . ,A g l y a m o v a ,G .v ,M e y e r ,E .&M a t z ,M .V .G e n e expression associated with white syndromes in a reef building coral, Acropora hyacinthus. BMC Genomics 16,3 7 1( 2 0 1 5 ) .

Jombart, T. Adegenet: A R package for the multivariate analysis of genetic markers.Bioinformatics24,1 4 0 3–1405 (2008).
Strader, M. E., Aglyamova, G. V. & Matz, M. V. Red ﬂuorescence in coral larvae is associated with a diapause‐like state. Mol. Ecol. 25, 559–569 (2016).
Quigley, K. M., Strader M. E., Github Code from “The role of gene expression and symbiosis in reef-building coral acquired heat tol- erance.” https://doi.org/10.5281/zenodo.6820007(2022). Acknowledgements We would like to thank the Traditional Owner communities who allowed us to collect northern colonies from their Sea-Country, in particular the Lama Lama Traditional Owners in the far northern Great Barrier Reef. We thank Sara Bell for performing larval and juvenile RNA extraction and normalization, and the University of Texas Gen- ome Sequencing and Analysis Facility for performing the library preparation and RNA tag-based sequencing. The Alabama Super- computer provided computational resources for data analysis. We would also like to thank Carlos Alvarez-Roa, Sam Noonan, Veronique Mocellin, and the GBRLegacy team for their help in collecting the Acropora tenuis colonies from the far north, and the coral collection teams from SeaSim at AIMS for collecting coral colonies from the central reefs. Corals were collected under permits G12/35236.1 and G18/41667.1. We would also like to thank Line Bay, Magena Marzonie, Josephine Neilsen, Gina Matthews, and Charlie Nicholas for their help with coral spawning and their help in sampling and measurements of juvenile photographs. M. Strader was supported by NSF-OCE #1935308 to MES. K. Quigley was supported by the Australian Insti- tute of Marine Sciences. Article https://doi.org/10.1038/s41467-022-32217-z Nature Communications| (2022) 13:4513 7 Author contributions K.M.Q. designed, performed, and managed the experiments. M.E.S. and K.M.Q. analyzed the data and wrote the paper. Competing interests The authors declare no competing interests. Additional information Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s41467-022-32217-z. Correspondence and requests for materials should be addressed to Marie E. Strader or Kate M. Quigley. Editorial Note This manuscript has been previously reviewed at another journal that is not operating a transparent peer review scheme. The manuscript was considered suitable for publication without further review at Nature Communications. Reprints and permission information is available at http://www.nature.com/reprints Publisher’s note Springer Nature remains neutral with regard to jur- isdictional claims in published maps and institutional afﬁliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. © The Author(s) 2022 Article https://doi.org/10.1038/s41467-022-32217-z Nature Communications| (2022) 13:4513 8