Phenotypic plasticity for improved light harvesting, in tandem with methylome repatterning in reef‐building corals

Acclimatization through phenotypic plasticity represents a more rapid response to environmental change than adaptation and is vital to optimize organisms' performance in different conditions. Generally, animals are less phenotypically plastic than plants, but reef‐building corals exhibit plant‐like properties. They are light dependent with a sessile and modular construction that facilitates rapid morphological changes within their lifetime. We induced phenotypic changes by altering light exposure in a reciprocal transplant experiment and found that coral plasticity is a colony trait emerging from comprehensive morphological and physiological changes within the colony. Plasticity in skeletal features optimized coral light harvesting and utilization and paralleled significant methylome and transcriptome modifications. Network‐associated responses resulted in the identification of hub genes and clusters associated to the change in phenotype: inter‐partner recognition and phagocytosis, soft tissue growth and biomineralization. Furthermore, we identified hub genes putatively involved in animal photoreception–phototransduction. These findings fundamentally advance our understanding of how reef‐building corals repattern the methylome and adjust a phenotype, revealing an important role of light sensing by the coral animal to optimize photosynthetic performance of the symbionts.


| INTRODUC TI ON
The modification of an organism's physical features (phenotype) through development and growth is affected by the interaction of gene expression (genotype) with environmental cues.This capacity for phenotypic plasticity allows organisms to optimize their physiological performance under different environmental conditions (Chevin et al., 2010;Nicotra et al., 2010;Torda et al., 2017).
While most organisms exhibit some degree of plasticity, the sessile condition of plants prevents movement to new environments if conditions become unfavourable.Consequently, they have evolved broad plasticity in their physical characteristics, such as leaf size or shape, root architecture or reproductive behaviours, to cope with changing environments and maintain optimal light harvesting (Borges, 2008).Animals generally exhibit far less plasticity than plants; except for reef-building corals.With similar life histories to plants, their colonies can display high levels of morphological plasticity.
Corals are modular, sessile organisms responsible for the net accumulation of calcium carbonate in coral reefs.The power to calcify is the result of animals acquiring photosynthetically fixed carbon through an obligate symbiosis with dinoflagellates (or microalgae) (family Symbiodiniaceae) (Falkowski et al., 1984;Gattuso et al., 1999;LaJeunesse et al., 2018;Trench, 1993).This means that corals, much like plants, make their living from light capture.
The metabolic integration is such that coral skeletons evolved to be efficient light collectors (Enríquez et al., 2017) with skeletal morphology adjusted in response to depth-dependent light availability (Malik et al., 2021).It is, therefore, not surprising that corals exhibit molecular signatures for perceiving and responding rapidly to changes in light availability.
While observed phenotypic plasticity is shaped by the interaction between genomes and environments, the role of epigenomes in this plasticity has captivated the interest of biologists.Phenotypic adjustments induced by environmental cues and gene expression may be influenced by chromatin factors like DNA cytosine methylation, a dynamic feature of many eukaryotic genomes, including plants, animals and fungi.DNA methylation is a process where methyl groups are added to cytosine bases of the DNA molecule and, in association with histone modifications, modify chromatin conformation (Buitrago et al., 2021;Hashimshony et al., 2003).
High-density methylation within promoter regions can silence genes, whereas lower-density intragenic methylation repatterning can influence alternative splicing of messenger RNA, leading to changes in an organism's phenotype (Bossdorf et al., 2010;Lev Maor et al., 2015).This is a reversible process influenced by environmental conditions, hence allowing phenotypic plasticity to occur (Verhoeven et al., 2010(Verhoeven et al., , 2016)).Moreover, methylation repatterning accompanies chromatin response to environmental changes without altering the DNA sequence and with the potential for heritable transmission.The repatterning is generally associated with other epigenetic effects such as histone modifications and changes in noncoding RNA.These methylome modifications can be assayed at single nucleotide resolution, providing the robust datasets required for identifying responsive underlying gene networks that could explain phenotypic adjustments (Hafner & Mackenzie, 2023;Kundariya et al., 2022;Sanchez et al., 2019;Sanchez & Mackenzie, 2016a, 2020;Yang & Mackenzie, 2020).
To examine the phenotype-to-methylome association, we conducted a reciprocal transplant experiment to induce light-mediated phenotypic responses in the reef-building Elkhorn coral Acropora palmata and investigated DNA methylation and transcriptional responses potentially responsible for plasticity.Extensive biometrics revealed changes in not only coral tissue pigmentation and metabolic rates but also skeletal morphology after 5 weeks.This skeletal remodelling was accompanied by intragenic methylome repatterning, discovered by signal detection with machine learning-based analysis.We further integrated differentially methylated (DMG) and expressed (DEG) gene datasets to elucidate how light responses integrate into gene regulatory networks controlling functional traits.By exploring the resulting hub genes and gene clusters, we were able to predict functional associations with observed phenotype changes and identify markers of plasticity in reef-building corals.Moreover, our results contribute to emerging evidence that epigenetics contribute to the machinery that can alter DNA structure during skeletal remodelling in metazoans.

| Light exposure conditions
We first characterized the natural light exposure within colonies of the coral Acropora palmata by measuring the incident light on both upperside and underside surfaces of in situ branches (base, mid-branch, tips).We estimated the vertical attenuation coefficient (K d ) of the water column and retrieved daily irradiance cycles from SAMMO light sensors (Meteorological and Oceanographic Monitoring Academic Service at UNAM).Light measurements were taken with a cosine-corrected quantum sensor (Diving-PAM, Walz, Germany) previously calibrated against a manufacturer-calibrated quantum sensor (LI-1400; LI-COR, USA).We identified two daily integrated PAR (Photosynthetically Active Radiation) conditions for in situ colonies that we also used for the experimental setup: HL surfaces represent fragments from upperside surfaces of branches exposed to ~20 mol quanta m 2 day −1 and LL surfaces represent fragments exposed to ~2 mol quanta m 2 day −1 (4-10% of upperside exposure).

| Coral sample collection and experimental setup
To induce phenotypic plasticity, we performed a reciprocal transplant experiment in the reef lagoon (30 m from the reef crest) in Puerto Morelos Reef National Park, Mexican Caribbean (Figure 2a).We sampled three colonies (representing three distinct multilocus genotypes, or genets, as detected with Standard Tools for Acroporid Genotyping STAGdb (Kitchen et al., 2020)) from a depth of 2-3 m (permits No. SGPA/DGVS/06960/17 and No. SGPA/DGVS/07846/17), each at least 300 m apart.We collected ~7 cm 2 fragments from HL (n = 42) and LL (n = 42) surfaces of branch surfaces keeping track of genet identification.We settled them in a reef-deployed PVC structure designed to simulate the light condition and colony position of source colonies.
We placed all coral fragments in their original light condition and orientation (i.e.HL facing up) and allowed them to heal and acclimate for eight weeks (when new growth was observed).Subsequently, we randomly divided HL fragments into control and treatment groups, and equal grouping was done for LL fragments.To induce phenotypic plasticity, coral fragments in treatment groups were flipped to the opposite light condition and position (i.e.HL fragments were flipped to a LL condition, while LL fragments were flipped to an HL position) (Figure 2a).This second acclimatory period was carried out for 5+ weeks.We estimated maximum excitation pressure determine successful acclimation to destination light condition (Iglesias-Prieto et al., 2004).We further compared Q m between experimental coral fragments (N = 84) and in situ colonies (30 colonies, N = 30 data points from upper and N = 30 from underside surface of branches) (Figure S3A).

| Structural traits
Phenotypic traits were measured based on parameters describing coral morphology and physiology after 13+ weeks of experiment (Table 1) (n = 12 per group condition for non-invasive techniques, and n = 9 per group condition for invasive techniques such as host total protein content, chlorophyll a content and symbiont cell counts) (Figure S1).Morphological features were described by polyp density (number of polyps per area) and by corallite height (mm), defined as the vertical distance between the corallite base and the top of the theca.Corallite height was delineated into "height classes" C 1 (0-1.5 mm), C 2 (>1.5-3 mm) and C 3 (>3 mm) to further estimate polyp

Traits
Parameters Units

O2
Maximum gross photosynthetic rate per area (P max ) Pressure over Photosystem II (Q m ) Dimensionless

TA B L E 1 Table of terms, definitions
and units for the structural, optical and photosynthetic parameters used to describe phenotypes.
density at each height class.The projected area of each fragment was estimated with photography and used to normalize physiological and morphological parameters.
To describe the structural and optical properties of the tissue, we measured reflectance (R) of the intact coral tissue as [De 675 = log (1/R 675 )] and estimated absorbance at 675 nm (De 675 ) (Enríquez et al., 2005;Vásquez-Elizondo et al., 2017).In addition, we estimated symbiont density, Chlorophyll (Chl) a and c from each fragment (Enríquez et al., 2005;Iglesias-Prieto & Trench, 1994;Vásquez-Elizondo et al., 2017).Briefly, tissue extractions were carried out using an air gun and filtered seawater (FSW).Slurries obtained from this method were subsequently homogenized at low temperature (Tissue-Tearor Homogenizer BioSpec Inc., USA) and centrifuged.The resulting pellet was re-suspended in filtered seawater (FSW) and preserved for symbiont cell counts (counted in a haemocytometer after the addition of 200 μL of iodine preservation solution), and for Chl concentration (extracted with acetone/dimethyl sulfoxide 95:5 vol/vol).Chla and c concentrations (ρ pigment content per projected surface area in mg Chl m −2 ) were estimated spectrophotometrically (three reads per sample) with a modular spectrometer (Flame-T-UV-VIS; Ocean Optics Inc., USA) using the equations described by Jeffrey and Humphrey (1975).
Downstream calculations of optical properties were performed by integrating the parameters detailed above (Enríquez et al., 2005;Scheufen et al., 2017).Symbiont density and Chla were normalized to calculate Chla per symbiont cell (Ci in pg Chla sym −1 ).The specific absorption coefficient of Chla (a* Chla ), a descriptor of the light absorption efficiency of the holobiont, was estimated as [a* Chla = (De 675 /ρ)•ln(10)] (Enríquez et al., 2005).Other calculations included the specific absorption coefficient of symbionts in hospite (a* sym ), a descriptor of the light absorption efficiency of symbionts, and the host mass-specific efficiency (a* M ), a descriptor of the light absorption efficiency per host mass, a descriptor that quantifies the potential benefits returned to the host, from the capacity of the symbiosis to collect solar energy (Falkowski et al., 1985;Scheufen et al., 2017).

| Coral tissue sampling and nucleic acid extraction
After 13+ weeks of experiment (Table 1), coral tissue from each fragment (n = 12 per group condition) was split into two samples and preserved, in both 95% ethanol and RNAlater (Ambion, Life Technology) and stored at −80°C until processing.Genomic DNA (gDNA) was extracted from 32 A. palmata fragments (n = 8 per group condition; 16 fragments per genet).We used the DNeasy Blood and Tissue Kit (Qiagen, Switzerland), as per the manufacturer's protocol.
Total RNA was extracted from the same 32 A. palmata fragments used in WGBS as previously noted.The tissue from each fragment was homogenized in TRIzol reagent (Ambion, Life Technology) before centrifugation with chloroform for 15 min at 12,000 g at 4°C.
The aqueous phase was then isolated and cleaned using Qiagen RNeasy Mini kit (Qiagen), as per the manufacturer's protocol with an additional on-column DNase treatment using RNase-Free DNase Set

| Methylome sequencing and data processing
Sequencing libraries were prepared with an average sequenc- The expected bisulphate conversion efficiency for this method is >99%.While Lambda DNA was not spiked-in to estimate nonconversion and mis-sequencing rate, we can assume that only CpG methylation context occurs in corals (Liew et al., 2018;Trigg et al., 2022); hence, non-CpG methylation should be close to zero.
We detected an average of 13.96% (SE = 0.07) methylation in CpG context and a neglectable 0.54% (SE = 0.01) in non-CpG contexts (CHG and CHH, where H is A, T or C), which resulted in an inferred bisulphite conversion efficiency of ~99.5%.

| Data processing
Methylation analysis was performed using Bismark 020.0 (Krueger & Andrews, 2011).Briefly (Figure S2), FASTQ files were qualityfiltered, and adapter sequences were trimmed using Trim Galore 0.5.0 (Krueger, 2015).A bisulphite-converted reference genome file was generated using Bismark-Bowtie2 algorithm, and the epig- Genes overcoming constraints 1-5 and displaying significant difference between control and treatment comparison according to likelihood ratio test (LRT) derived by the anova function from stats packages were identified as DMGs (Figure S2).A detailed description of how to define and compute DMPs and potential DMGs is included in the Methyl-IT vignettes and the package manual (Sanchez, 2020).
Acropora palmata genome annotation file (GSE245494) was used to annotate genome features.

| Methyl-IT downstream analysis (Methyl-IT. Utils)
A HC was performed to provide an initial estimation of the number of possible groups and information on their members.The effectivity of HC depends on the experimental dataset, the metric used, and the glomeration algorithm applied.Ward's agglomeration algorithm was used as it seems to perform much better on biological experimental datasets than the other available algorithms (e.g.UPGMA, UPGMC) (Ward, 1963).

| RNA sequencing and data processing
RNA libraries for 2 × 150 bp paired-end sequencing were prepared using the NEBNext Ultra II Nondirectional Library Prep Kit with polyA selection (New England Biolabs, Inc.).Samples were run on one plate of the Illumina NovaSeq platform.Illumina universal adapters and reads below PHRED of 22 were trimmed using Cutadapt 3.5 (Martin, 2011).Filtered reads were mapped to the Acropora palmata genome (JAOVVL01) using the RNA-seq aligner STAR (2.5.3a) with read count data generated by the -quantMode GeneCount parameter.Reads were verified using the generated BAM files for input into htseq count.

| Differential gene expression analysis
Gene count normalization and differential expression analysis were performed using DESeq2 3.12.0.Significant Differentially Expressed Genes (DEGs) were determined via pairwise comparison among control and treatment groups (HL vs. HL→LL, LL vs. LL→HL) and genet, with a false discovery rate (FDR)-adjusted p-value of <.05 (Figure S2).

| Network-associated responses
Based on available data, we chose gene-gene interaction networks (predicted protein-protein interaction networks), as they are undirected, and the graph is non-sequential (X affects Y, but we do not know how) (Le Novere, 2015).A key feature for the biological interpretation of graph properties are hubs.Small-degree nodes (non-hubs) are the most abundant, but high-degree nodes or hubs, although less frequent, have a much higher number of interactions (Dennison & Alberte, 1982), and its loss will cause a major breakdown of the network and be lethal to an organism (Albert, 2005;Jeong et al., 2001;Said et al., 2004).

| Weighted gene correlation network analysis (WGCNA): Agnostic network analysis of DMG and DEG datasets
To understand the interaction between the change in methylomes (n = 32) and transcriptomes (n = 32), and how genes contributed to the change in phenotype, we performed a WGCNA 1.71 (Langfelder & Horvath, 2008) (Figure S2).We generated a dataset (gene list) comprising the subsets of (1) all DMGs (from Methyl-IT 0.3.1.2) and ( 2) DEGs (from DESeq2 3.12.0)that presented at least one DMP (from Methyl-IT 0.3.1.2) after the change in light treatment.This dataset included outputs from all groups (HL, HL→LL, LL, LL→HL) and genets to discover general patterns of gene contribution.Additionally, a binary annotation was included to keep track of DMG, DEG and both DMG-DEG.As a result, each sample was represented as vector of 3272 genes/coordinates, where each coordinate was given by the sum of HD at each DMP on the given gene.
We first performed an HC applying Ward's agglomeration algorithm to provide an initial estimation of the number of possible groups and information on their members.Methyl-IT function pcaLDA was used to perform a PCA and a PCA + Linear Discriminant Analysis (LDA).
Based on the cumulative proportion of variance, the PC1 and PC2 car-

| Network analysis from curated databases: Predicted biological networks from DMGs
To identify the biological meaning of potential relationships among DMGs, we inferred gene interaction networks from stringApp in Cytoscape 3.8.2(Figure S2).The associations in the string database provide known and predicted protein-protein associations data for many organisms, including both physical interactions and functional associations, by integrating available experimental data and pathways from curated databases (Doncheva et al., 2019;Szklarczyk et al., 2017).We used only our detected DMGs (without network expansion) as input in string protein query (Swiss-Prot hit name) to retrieve an arbitrarily long list of nodes and interactions.This approach is generally used to retrieve string networks from proteomics and transcriptomics studies (Doncheva et al., 2019;Szklarczyk et al., 2017).
The best hit for baseline networks was reached with string query for Homo sapiens.We recognize that cross-species knowledge transfer is quite challenging because the phylum cnidaria diverged from Bilateria 550 million years ago and may have fundamentally different genetic architectures.Also, as species diverge, protein functions change and are re-purposed through divergent and convergent evolution, and genetic interactions are often rewired (Fan et al., 2019).Parallel to this, in network-based approaches, most predicted interactions for each species are not experimentally verified.
Despite these limitations, the best hit in string query was still Homo sapiens, perhaps, due to the presence of conserved stress response, conserved pathways (e.g.extrinsic and intrinsic apoptotic pathways, ion trafficking system) and (predicted) gene products between early and late branching metazoans at the molecular level (Bhattacharya et al., 2016;Courtial et al., 2017;Davy et al., 2012;Drake et al., 2013;Ottaviani et al., 2020).Furthermore, genes that emerged from our data and clustering analyses predicted gene interactions associated to coral biological processes that aligned well to the measured phenotypes, suggesting our analytical approach was plausible.
To identify hub genes in the retrieved baseline network, we performed a statistical analysis incorporating extended centrality measures from CentiScape 2.2 App (Eigenvector of centrality).We chose Eigenvector centrality because this attribute ranks nodes by taking into consideration not only the number of interactions of a node (degree) but also the centrality of the interactions that it is connected to.Then, we assigned methylation signal and Eigenvector of centrality as attributes to the nodes (Sanchez & Mackenzie, 2020).
In networks, a protein with a very high Eigenvector is a protein inter-

| Predicted biological networks from DEGs
To identify the biological meaning of potential relationships among DEGs, we inferred gene interaction networks from stringApp in Cytoscape 3.8.2 with the same workflow as for DMGs interaction networks (Figure S2).The output dataset generated in DESeq2 3.12.0 was imported to this baseline network to assign gene regulation (up-regulated or down-regulated) as a node attribute.

| Predicted biological networks from DMG and DEG datasets
To further explore the association between DMGs and DEGs, we integrated DEGs into DMGs network datasets (Figure S2).The integration was done at the cluster level (after clusterMakerApp independent analyses) with the criteria that a DEG cluster be selected if it contained at least one gene also identified as DMG.The output datasets generated from Methyl-IT 0.3.1.2and DESeq2 3.12.0were imported to this baseline network to assign node attributes.With this approach, we were able to enhance networks by adding new attributes to the nodes: DMG, DEG and both DMG-DEG.We maintained the attributes Eigenvector of centrality, methylation signal (from Methyl-IT 0.3.1.2) and gene regulation (up-regulated or down-regulated from DESeq2 3.12.0).A new clustering was performed based on Eigenvector of centrality ranks and 1st, 2nd and 3rd neighbours of high-ranked nodes.
To identify over-represented functions in the DMGs-DEGs integrated clusters, we performed an NEA as explained before.Because key coral biological processes emerged from the new clustering, we further explored potential key regulators and candidate genes involved in lightmediated phenotypic plasticity of structural traits in corals.

| Corals exploit intra-colonial environmental differences through plasticity
The branching coral Acropora palmata is the dominant reef builder on shallow, wave-exposed Caribbean reefs.Colonies exhibit treelike morphologies with strong intra-colonial light gradients ranging from 3% to 100% of irradiance just below the sea surface (E s ).In our experiment, we quantified intra-colonial phenotypic plasticity by measuring traits from High Light (HL) surfaces (fragments from upperside surfaces of branches) exposed to 70% of E s and Low Light (LL) surfaces (fragments from underside surfaces of branches) exposed to 3-7% of E s (Figure 1, Table 1).Through the examination of structural, optical and physiological traits of HL and LL fragments we identified two distinct phenotypes.HL phenotypes had significantly more polyps per area (Figure 1f), more C3 corallites (height class >3 mm) per area (Figure 1g), total host protein (Figure 1h) and symbiont cell density (Figure 1i).Surfaces with taller corallites can favour the formation of internal light gradients, increase levels of pigment self-shading and reduce the proportion of polyp-surface exposed to the external HL levels (Enríquez et al., 2017;Ow & Todd, 2010).In contrast, LL phenotypes showed less polyps per area, a flat surface with a majority of C1 corallites (height class 0-1.4 mm) that can facilitate the lateral spread of light.
We expected these differences in skeletal features to affect how corals collect and utilize light for colony growth.To disentangle this, we used algal symbiont density, chlorophyll a (chla), host soluble proteins and in vivo light absorption of the intact coral tissue (Table S1) to describe light absorption efficiency of HL and LL phenotypes.We estimated three optical traits: a* Chla (m 2 mg Chla −1 ), which describes the holobiont's efficiency to absorb light (Enríquez et al., 2005), a* sym (m 2 sym −1 ) that describes in hospite light absorption efficiency of the algal symbionts, and a* M (cm 2 mg protein −1 ) indicative of the potential return for the host (mass) of the energy absorbed (Falkowski et al., 1985;Scheufen et al., 2017) (Table S1).We detected less algal cell densities (Figure 1i) but more chla per cell in LL phenotypes (Figure 1j), and opposite traits in HL phenotypes.This resulted in equal chla concentration in both HL and LL phenotypes (Figure 1k).
These findings contradict the general assumption that more light availability induces less pigmentation in photosynthetic organisms (Dubinsky & Schofield, 2010) and confirm the ability of coral skeletal features to rewire the algal light environment.Consequently, a* Chla showed that both phenotypes of the coral colony are equally efficient in absorbing light (Table S1), a response reached by adjusting skeletal features.
Acropora palmata exploits a wide range of light environments, with its algae species (Symbiodinium "fitti") (Baums et al., 2014) adjusting to the strong intra-colonial light gradient.These colonies fine-tune structural traits: algal density, chla density and skeletal morphology (Figure 1).The mechanism minimizes "pigment packing" in underside surfaces in response to LL conditions, resulting in light collectors as efficient as those of HL surfaces.Although photosynthesis is more active on upper side surfaces (Figure 1n), the ratios of photosynthesis/respiration (P/R = 3), photosynthetic efficiency (α) (Figure 1l) and minimum quantum requirement (Table S1) were similar on both sides of branches, indicating an optimization of photosynthesis.These observations further highlight the ability of A. palmata to optimize light absorption and utilization through plasticity as a central strategy to compensate for the strong intracolonial light gradient and maximize colony productivity for growth.

| Light-induced phenotypic plasticity
Acropora palmata frequently reproduce asexually via branch fragmentation, a result of physical disturbance (i.e.waves and storms) (Baums et al., 2006).Branches are often turned upside down when they land on the benthos.Fragmentation thus induces strong and rapid changes in light regimes, where survival is dependent on their successful acclimatization to the new light conditions.Presumably, upper and underside branch surfaces interchange phenotypes through acclimation to new light regimes.We took advantage of this life history trait to induce plasticity by altering light exposure in a reciprocal transplant experiment (see 'Section 2', Figure 2a).Table S1).Interestingly, a* Chla , the holobiont's efficiency to absorb light, indicated that fragments may continue changing pigmentation and skeletal features to fully optimize performance (Figure 2d,e).Nonetheless, most significant changes were observed in skeletal features (taller corallites per area), polyp density and the balance between chla density and symbiont density (Figure 2, Table S1).Significant differences were found among the four group conditions (R = 0.336, p < .001, Figure S5), suggesting a response driven by light-mediated phenotypes and not the genet (R = 0.188, p > .05).

| Light-mediated methylome repatterning
Following induced phenotypic plasticity, we investigated DNA methylome response to coral group conditions (n = 8 per group condition).
We identified CpG context methylation (~14%) to be higher in the A.
To associate the observed coral phenotypic plasticity with highresolution DNA methylome variation, we used a signal detectionmachine learning approach (Methyl-IT 0.3.1.2) (Sanchez, 2020), designed to discriminate methylation signal induced by environmental variation at individual cytosine positions (Sanchez et al., 2019;Sanchez & Mackenzie, 2020;Yang et al., 2020).We assessed geneassociated, differentially methylated positions (DMPs) with no regard to methylation density, context or directionality (hypo/hypermethylation) (Figure 3c).Parallel analysis of DMP variation within control samples allowed discrimination of treatment-associated DMPs and their classification on the basis of hierarchical clustering (HC) and principal component analysis (PCA), which enabled an unbiased view of methylome repatterning (Figure 3d, Figure S7).We detected significant separation of control and treatment samples, indicating that light-mediated methylome modifications were driving the first two principal components (PC1 and PC2 carried 76% of the total sample variance).Genes with the strongest discriminatory power from PCscores in PC1 were associated with cell cycle, extracellular matrix (ECM), regulation of transcription and transduction, and biomineralization (Table S2).
DNA methylation as a chromatin feature that influences DNA stability can be dynamic and rapidly responsive to environmental change.Under experimental conditions, spontaneous natural methylation or "background methylation" may occur.This stochasticity has complicated the discrimination of treatment-associated signal and the understanding of phenotypic adjustments with methylome modifications (Hafner & Mackenzie, 2023;Sanchez & Mackenzie, 2023;Yang & Mackenzie, 2020).As the significance of individual cytosine variations in phenotypic responses becomes more evident, innovative methods, such as MethyIT (Sanchez et al., 2019), have incorporated these features into machine learning to distinguish treatment-associated differential methylation.The approach treats methylation data as probability distributions, permitting variation within multiple control samples to be subtracted from the treatment datasets to discriminate treatment-specific signal and the related methylation regulatory machinery.Machine learning then permits validation of treatment association with over 98% confidence.

Further validation of this approach is accomplished, in models like
Arabidopsis, with incorporation of mutations in the RNA-directed DNA methylation pathway (Kundariya et al., 2022).However, the approach is especially valuable in non-model systems where changes cannot be confirmed with targeted mutation(s).The association between phenotype change and treatment-associated methylome modification is informative in understanding the underlying molecular features of the phenotypic change by directly identifying responsive gene networks (Hafner & Mackenzie, 2023;Kundariya et al., 2022;Sanchez et al., 2019;Sanchez & Mackenzie, 2016a, 2016b).

| Underlying molecular features of the phenotypic change: Responsive gene networks from DMG and DEG datasets integration
Significant coral phenotypic changes and gene-associated methylome repatterning warranted deeper investigation to uncover functional relationships with gene expression.We identified 861-2255 DMGs (Methyl-IT 0.3.1.2) and 1334-6479 DEGs (DESeq2 3.12.0.) (Figure S6).We further used a network-based approach to integrate these datasets, which provided us with a collection of nodes and edges representing putative gene interactions (Albert, 2005).
Since network-based analyses can be influenced by the available annotation for a given species, we first performed a Weighted Gene Correlation Network Analysis (WGCNA 1.71) (Langfelder & Horvath, 2008) of coral gene expression and methylome modifications (Figure S2), with visualization and statistical analysis in Cytoscape 3.8.2(Sanchez & Mackenzie, 2020).
The HC (Figure 4a), PCA (PC1 and PC2 carried 92% of the total sample variance) and a linear discriminant analysis (LDA) (Figure S7) to reduce dimensionality, showed a classification of samples separating controls (purple) and treatments (red) groups regardless of genet or destination light treatment.Furthermore, methylometranscriptome-derived gene network information revealed an extraordinary number of hub genes (Albert, 2005;Sanchez & Mackenzie, 2020) likely to be integral to morphologic plasticity in corals.General biological processes included visual and sensory perception, growth and immunity, including carriers, transporters and receptors.Two main cluster categories were identified (Figure 4b).
A Type I sub-network (Figure 4c) showed genes with strong genegene interaction (edges with strongest correlation weights), denoting genes with similar contribution to the change in phenotype.The main sub-network under this category was enriched in ECM gene products, collagen-like domains, signalling activity, cell-cell adhesion and EGF domains, putatively associated with soft tissue growth and biomineralization (Figure 4c).139/199 genes in this sub-network were found to be differentially methylated, highlighting the influence of methylation on changes in coral skeletal growth.Interestingly, we We further aimed to investigate gene interactions in available experimental data by consulting pathways from curated databases (Doncheva et al., 2019;Szklarczyk et al., 2017) to identify predicted gene networks (see 'Section 2'; Figure S2).We performed an initial screening of network-based analysis from methylome repatterning.
With the set of only DMGs as input in stringApp (Cytoscape 3.  (Stadler & Richly, 2017), are essential in DNA methylation patterning (Huck-Hui & Bird, 1999).We performed a separate screening from differential gene expression.We used the set of only DEGs to run the same network analyses separately (Figure S9).However, we found significantly lower values of node Eigenvector of centrality (Figure S10).This limited the identification of hubs from gene expression data alone.
The identification of key regulators in DMG sub-networks led us to investigate integrated DMG and DEG networks in available experimental data and pathways from curated databases (Doncheva et al., 2019;Szklarczyk et al., 2017).We detected sub-networks of hub genes contributing to phenotypic changes and associated cellular processes, including inter-partner recognition and phagocytosis, regulation of host-symbiont biomass, and calcification (reviewed in (Davy et al., 2012)).Sub-networks aligned with these processes, in- and calcification sub-networks (Figure 5b).Previous work has elucidated the role of light in the coral animal symbiont-independent metabolism.It was found that enhanced coral calcification but not photosynthesis occurs under blue light exposure (Cohen et al., 2016), and evidence suggests a light-mediated electrical potential in coral epithelia (Taubner et al., 2019).Moreover, Cnidaria constitutes the earliest branching phylum containing a well-developed visual system.Some box jellies (Cubozoa) have advanced camera-type eyes with photo receptor cells similar in sensitivity to vertebrate eyes (Kozmik et al., 2008).These findings suggest that coral animal photoreception-transduction is implicated in changes in growth patterns and skeletal features, although the role of photoreceptionassociated proteins is yet to be explored.
Consistent with significant changes in coral skeletal features for the optimization of light harvesting, we identified a coral growth and biomineralization sub-network (Figure 5b, 86 nodes and 200 edges) that revealed NEA categories related to ECM proteins (KW-0272, 18 enriched genes, FDR < 0.001), cell-cell adhesion (KW-0130, 14 enriched genes, FDR < 0.001) and EGF domains (KW-0245, 13 enriched genes, FDR < 0.001) as most prominent.From the subnetwork of hub genes, integrin and spectrins interacted with glycoproteins and lipoproteins ( DMG PTK2, DMG-DEG THBS1 and DMG LRP5) involved in coral biomineralization (Drake et al., 2013;Gutner-Hoch et al., 2017;Hemond et al., 2014).Thrombospondin has previously been suggested for its role in biomineral remodelling (Mummadisetti et al., 2021).Consistent with this network topology, an important transcriptional interaction identified was SOX9, a transcription factor (TF) with a role in skeletal development.Associated nodes involved TFs DEG HIF1A, DEG HOPX and DEG FOXL1, with the receptor protein DEG NOTCH2 up-regulated in this transcriptional interaction (independent of the destination light condition).All hub genes interacted with a strongly connected cluster of collagen-like domains that was mainly transcriptional and up-regulated when transplanted to HL, while down-regulated when transplanted to LL. Collagen plays a structural role in the skeletal organic matrix (SOM), and presence of SOM in calcifying organisms appears to be a prerequisite for the formation and growth of biominerals (Allemand et al., 1998;Mummadisetti et al., 2021).They represent candidates for reverse genetics.
Our analysis provides evidence of the association of genic methylation repatterning with programmed changes in phenotype.One interpretation of these outcomes is that dramatic shifts in gene expression are accompanied by methylation changes that stabilize local chromatin to help re-establish homeostasis.However, environmental changes also induce non-random changes in gene body methylation that influence alternative splicing activity to modulate gene expression and phenotype (Ausin et al., 2012;Lev Maor et al., 2015;Yang et al., 2014;Zhang et al., 2020).While our data do not provide information on the regulation of gene expression, the integration of methylation and transcriptional information makes significant inroads in the identification of networks underpinning coral phenotypic plasticity and provides a roadmap for studies of other non-model organisms.As methylome repatterning is a signature of chromatin reorganization and 3D architecture, integrating datasets enhances the potential to identify essential interactions and likely signatures of gene network coordination that may not be seen in DEG datasets alone (Ouyang et al., 2020).Our integrated analyses offered the potential to predict phenotype at the gene network level and further postulate a light-mediated sequential response triggered by animal sensing of initial light stimulus to the change in skeletal morphology (Figure 6).
Our study measured the remarkable phenotypic plasticity exhibited by reef-building corals, emphasizing its pivotal role in optimizing light harvesting mechanisms.Likewise, it highlights their ability to coordinate an epigenomic response involving an extensive methylome modification.Through a sophisticated interaction of environmental, genomic and epigenomic processes, these corals emerge as true ecosystem engineers, strategically adjusting their morphological features to exploit light gradients.This plasticity is a key driver in maximizing both colony and reef growth.

E TH I C S S TATEM ENT
This article does not contain any studies with human participants or animals where formal consent is required, and ethics approval for this study was not required.

B EN EFIT-S H A R I N G S TAT E M E N T
Benefits generated include a research collaboration developed with scientists from the countries providing genetic samples, all collaborators are included as coauthors, the results have been shared with local and the broader scientific community, and the research disentangles previously unknown aspects of the biology of a critically endangered species.More broadly, our group is committed to international scientific partnerships, as well as institutional capacity building.Other benefits accrue from the sharing of our data and results on public databases as described above.

(
Qiagen).To maximize concentration of eluted RNA, the same 35 μL of RNAse-free water was twice passed through the spin column for the final isolation step.Concentration and purity (A 260/280 ) were analysed via NanoDrop ND-1000 spectrometer and quality assessed (RIN > 7) with Agilent 2100 Bioanalyzer (Agilent Technologies).Total RNA was sent to Admera Health (New Jersey, USA) where concentration and purity were re-analysed (Qubit® dsDNA BR Assay Kit on a Qubit® 2.0.Fluorometer).
library insert size of 450 bp and according to TruSeq® DNA Methylation Kit protocol, Sample Preparation Guide (Illumina Inc., USA).Briefly, standard reaction mix consisting of 130 μL of the CT Conversion Reagent and 20 μL of each DNA sample were used for bisulphite conversion (EZ DNA Methylation-Gold™ Kit) in a thermal cycler (Eppendorf® Mastercycler® Pro S).After the incubation period, bisulphite-converted DNA was purified following the protocol of EZ DNA Methylation-Gold™ Kit.The bisulphite-converted sequencing library was enriched following the TruSeq® DNA Methylation Kit protocol.Libraries were validated and quantified with qPCR.Sequencing was performed in Illumina Hi-Seq platform with pair-end reads of 2 × 150 and a mean depth coverage of 30× (Illumina Inc., USA).The estimated number of passed filter reads per sample was 110-120 million paired-end reads (55-60 M reads in each direction).
enome library sequenced data were aligned to the Acropora palmata genome (JAOVVL01).Methylation information was then extracted from the output SAM files, and resulting genome tracks were used for the visualization and reporting of downstream differential methylation calculations.Methylated and unmethylated read counts for all cytosines across the genome in the CG, CHG and CHH context were obtained from census files.The approach used for downstream analysis was based on the identification of Differentially Methylated Positions (DMPs) and Differentially Methylated Genes (DMGs) using R package (Methyl-IT 0.3.1.2) (Sanchez, 2020) (FigureS2).Briefly, methylation count (COV) files were read into R to calculate Hellinger Divergence (HD, a variable used to measure methylation level divergence) by using the pool of methylation counts for control samples as reference.Potential DMPs (pDMPs) were estimated based on critical values of HD α=0.05 for each sample from the best fitted probability distribution model; in this case, a three-parameter gamma distribution model.Final DMPs were estimated from the set of pDMPs by calculating the optimal cut-off threshold for HD based on Youden index.Generalized linear regression analysis (generalized linear model, GLM) was applied to test the difference between group DMP counts (among control and treatment groups; HL vs. HL→LL, LL vs. LL→HL) for selected genomic features.The fitting algorithm approaches provided by glm and glm.nb functions from the R packages stat and MASS were used for Poisson, Quasi-Poisson and negative binomial models with logarithmic link.The "countTest2" function in Methyl-IT was used to implement the selected model.The following parameters were applied to identify significant DMGs (Yang et al., 2020): (1) the minimum DMP count per bp on gene body: CountPerBp = 2.5; (2) a minimum count per sample (on average) in at least 5 DMPs in one group: minCountPerIndv = 5; (3) a maximum coefficient of variance for each group: maxGrpCV = 1; (4) minimum value of the logarithm of fold changes: Minlog2FC ≤1; (5) p-value cut-off: pvalCutOff = .01;(6) p-value adjustment was performed by Benjamini-Hochberg method: pAdjustMethod = "BH".Parameters 1-3 are addressed to prevent spurious DMGs, which cannot be rejected by the generalized linear regression algorithms.
ried 92% of the total sample variance and could split the samples into meaningful groups.PC-scores for each gene (gene-score) indicate the discriminatory power in the clustering (control vs. treatment) and its genomic/epigenomic contribution (in terms of proportion of the whole phenotypic variance) to the change in phenotype.Genes PC-scores (from 14 PCs) were then used to build the pairwise correlation matrix for the WGCNA.Kendal's tau correlation was selected since it is better at detecting nonlinear behaviours and is more conservative than Pearson's correlation.The resulting weighted correlation matrix was then constructed as a network in the R package WGCNA(Langfelder & Horvath, 2008) and exported as edge list (interactions with weights) and node list files with assigned modules into Cytoscape 3.8.2 for visualization.Each entity of the dataset is a (gene) node, and two nodes are connected if their correlation or distance reaches a threshold (here set to 0.4).Network topology included gene discriminatory power (gene-scores) and a measure of how similarly they contributed to this classification (weight from correlation).The correlation network was analysed and visualized in Cytoscape 3.8.2(Sanchez& Mackenzie, 2020).
acting with several important proteins (regulating them or being regulated by them), thus suggesting a central super-regulatory role or a critical target of a regulatory pathway.We used Eigenvector as parameter to perform k-means clustering algorithm (clusterMakerApp in Cytoscape 3.8.2.) (Sanchez & Mackenzie, 2020) to identify clusters of hub genes.To identify over-represented functions in the large set of DMGs, we performed NEA.Enriched terms were retrieved as UniProt KnowledgeBase (kw) categories in String Enrichment App in Cytoscape 3.8.2(Sanchez & Mackenzie, 2020).
Coral fragments from three colonies representing three distinct genets were manipulated, so that HL phenotypes (n = 21) and LL phenotypes (n = 21) experienced unchanged light fields, and treated fragments were switched to the opposite light condition HL→LL (High Light to Low Light, n = 21) and LL→HL (Low Light to High Light, n = 21) (Figure 2a).Within five weeks, treated coral fragments significantly adjusted their phenotype.The acclimation was gradual (tracked by visual inspection), and transplants F I G U R E 1 Phenotypic plasticity of Acropora palmata in response to light availability.(a) Acropora palmata colonies have a strong intracolonial light gradient.The branching morphology exhibits modules (polyps) exposed to direct sunlight (HL surfaces) and modules growing in the shade (LL surfaces).(b) Morphological skeletal features of the branch cross-section showing the transition from upperside to underside of the branch.(c) HL surfaces and (d) LL surfaces show distinct skeletal morphology, (e) with corallites significantly taller on the surface exposed to HL. (f-n) Phenotypic traits of HL (n = 21) and LL surfaces (n = 21) from three genets.Centre lines show the median and centre squares the mean; box limits indicate the 25th and 75th percentiles; whiskers extend 1 time the interquartile range.For all panels, ***p < .001;**p < .01;*p < .05;ns p >.05, two-tailed, unpaired Student's t-test.(f) Polyp density (# polyps cm −2 ), (g) Density of corallites larger than 3 mm in height (# polyps cm −2 ), (h) soluble host protein (mg protein cm −2 ), (i) symbiont density (# sym cm −2 ), (j) Chla per symbiont cell (Ci, pg Chla sym −1 ), (k) Chla density (mg Chla m −2 ), (l) photosynthetic efficiency (μmol O 2 μmol quanta), (m) respiration rate (μmol O 2 m −2 s −1 ), (n) maximum photosynthetic rate (μmol O 2 m −2 s −1 ).became increasingly similar in morphology to coral surfaces of the destination light condition (Figure 2c, Figure S4), comparable to what is observed when coral colonies are transplanted along depth gradients (Malik et al., 2021).Pressure over photosystem II (Q m ), metabolic rates and visual growth indicated that corals acclimated successfully to the destination light conditions (Figure 2b,c, F I G U R E 2 Induced phenotypic plasticity with reciprocal transplants.(a) Schematic representation of the experimental design.Control fragments from HL (n = 7 per genet) and LL (n = 7 per genet) remained unchanged, while treatment fragments (HL→LL, n = 7 per genet; LL→HL, n = 7 per genet) were manipulated in a reciprocal transplant that altered their light exposure by ~80%.After 13 weeks, light phenotypes were described, and tissue was collected for genomic and epigenomic analyses.(b) Fold change of main phenotypic traits showing the acclimatory mechanism to the destination light condition after 5+ weeks of transplant.(c) Visual inspection of one genet after 5+ weeks showing the change in corallite height and density.(d, e) Changes in optical traits based on specific absorption coefficients, a* Chla which describes the holobiont's efficiency to absorb light and a* sym , which describes in hospite light absorption efficiency of the algal symbionts.
Phenotypic plasticity was induced by modifying approximately 80% of available light, equivalent to about 18 mol quanta m −2 day −1 .These changes, although substantial, are representative of conditions frequently encountered by coral species throughout their life cycle.The optimization of whole colony metabolic performance is achieved through adaptations at the individual module (polyp) level, primarily for resource acquisition.Similarly, in plants, phenotypic plasticity emerges as a local response, often observed in shaded branches, which serves to enhance light capture and utilization for growth (De Kroon et al., 2005).Despite the distinction that corals are colonial animals, where each polyp (module) is an individual, coral photoacclimation should consider the integrated response to local condition (within the colony) experienced by local modules.Furthermore, modular plasticity may represent an evolutionary trait subject to selection pressures, akin to what has been proposed for plants (De Kroon et al., 2005).

F
I G U R E 3 DNA Methylation context and light-mediated methylome repatterning.(a) Proportion of methylated cytosines (n = 8 per group condition) was highest in the CpG context and neglectable in CHG and CHH contexts.Pie charts show the number of cytosines (×10 6 ) in each context.(b) Mean methylation levels of all cytosines were highest at genic regions; 2 kb upstream of Transcription Start Site (TSS), and 2 kb downstream of Transcription End Site (TES) are shown.Methylation levels were computed, divided into 60 bins, and plotted by genet and group condition.(c) Number of DMPs per group conditions identified by Methyl-IT, with centroid of control groups used as reference.DMPs were always higher in treatments than control samples.Two A. palmata genets are shown for comparison.(d) Hierarchical clustering of DMPs in genic regions classified by Hellinger Divergence.Classification of samples separated control (purple) and treatment (red) samples regardless of genet or destination light treatment.observed low gene-score values in this cluster, meaning that individually each gene carries small proportion of the whole phenotypic variance.Seemingly, these genes were co-methylated to give rise to a quantitative cumulative effect on phenotypic variation among treatments, leading to the change in growth pattern, but with low discriminatory power individually.F I G U R E 4 Responsive gene networks integrating methylome and transcriptome data: WGCNA.The data were prepared by combining DMG and DEG datasets into one large dataset.To estimate the initial number of possible groups we performed a (a) hierarchical clustering, which showed a classification of samples separating controls (purple) and treatments (red) groups regardless of genet or destination light treatment.(b) Whole network of gene-gene interactions.(c) Type I sub-network showing genes with strongest gene-gene interactions (edges with strongest weights from correlation but low gene-score), denoting genes that have similar contribution to the change in phenotype (n = 199 genes).(d, e) Type II sub-networks of hub genes showing strong interactions and loadings (highest gene-scores) denote hub genes with strongest contribution to the change in phenotype (discriminatory power of treatments from controls).(f-h) Top 10-20 genes based on gene-score in each sub-network.The coloured line between genes represents weight values from correlation matrix, low weight values (yellow) to high weight values (purple), node colour indicates if DMG (yellow), both DMG-DEG (blue), DEG (grey).Sub-networks Type II represented critical genes with the strongest discriminatory power, denoting hub genes with strongest contribution to the change in phenotype (Figure 4d,e).Accordingly, we considered these loci to be strong candidates for biomarker identification.Annotated hub genes included DEG SLC7A2, DEG SLC22A13, DEG SLC17A6B, DEG PANP, DEG IFI30, DEG TRIM71, DEG MELC2, DEG ADAMTS18, DEG CTRC, DMG HECTD4, DMG LIPK, DEG TNR, DMG C0H691, DMG TRPV6, DEG HSP-16.2,DEG SUSD2, DEG COL6A5, DEG BTBD2, DEG CASR DMG GPR133, DMG RABL6, DEG PRSS27, DEG SSTR5, DMG MPDZ, DMG SPTAN1, DEG ABHD4, DEG COL12A1, DEG SPDEF, DEG HES4-a, DEG PTK7, DEG ASIC3 DEG MFSD6, DEG COL11A1, DMG CHMP6, DEG PRDM6, DMG FGFR1, DMG TMPRSS15, DMG TMPRSS11D and DEG B3GAT3.
8.2), we built a baseline network and performed clustering analyses (with network centrality measures).Independent of the destination light condition and genet, enriched function categories were posttranscriptional protein variants (kw-0621 polymorphism, kw-0025 alternative splicing, kw-0597 phosphoprotein), cytoskeleton proteins (kw-0206), calcium-dependent proteins (kw-0106 calcium), growth-associated proteins (kw-0131 cell cycle, kw-0965 cell junction) and phagocytosis (kw-0966 cell projection) (Figure S8).A common feature of main candidate hub genes (Figure S8, ATR, PIF1, FANCD, SIRT1) was DNA damage repair response (DRR), a mechanism regulated by chromatin conformation.Chromatin remodelling complexes, at the core of DRR dependent of the destination light conditions and genet, and generally targeting the same enriched categories, pathways or protein families (Figure 5, Figures S11-S14).Network topology and node hierarchy based on Eigenvector of Centrality suggested a general DMG-DEG interplay targeting key regulators (hubs) and triggering a cascade of responses, evident when source-sink nodes (Figure 5a CEP290; Figure 5b THBS1, PTK2, SOX9; Figure 5c RAB1A; Figure 5d ANAPC1) interacted with strongly connected clusters.Cells have mechanisms to detect specific environmental signals to transduce and trigger the appropriate responses.We identified a symbiont-independent photoreception-phototransduction subnetwork (Figure 5a, 51 nodes and 139 edges) potentially involved in animal sensing of a light stimulus (p-Cluster I).Network Enrichment Analysis (NEA) yielded mainly cilium assembly (GO.0060272, 19 enriched genes, FDR <0.001).The putative hub gene DMG-DEG SPTAN1 interacted with DMG CEP290 and a transcriptional non-hub cluster (USH2A, CRB1, EYS, MYO3A) involved in phototransduction reception.Interestingly, this sub-network further interacted with growth Coral holobiont biometrics showed the strong regulation of algal symbiont density by light-mediated adjustments in skeletal features.These changes are achieved by phagocytosis-exocytosis-endocytosis with associated innate immune recognition(Davy et al., 2012).We identified associated sub-networks with enriched categories related to vesicle/vacuole-mediated transport (GO.0016192, 41 enriched genes, FDR <0.001) (Figure5c, 78 nodes and 174 edges).The core sub-network of hub genes (Figure5c) interacted with spectrins, ankyrins and key regulators of intracellular membrane traffic GTPases-RAB.Previous studies have emphasized pattern recognition receptors (PRR) as key players in symbiosis establishment(Davy et al., 2012).Main PRR identified are endocytosis mediator C-type lectin domain family member ( DMG-DEG MRC1) and Toll-like and Toll/interleukinlike receptors ( DEG TLR6 and DEG TLR1).Furthermore, we identified a sub-network associated to immune system responses (NEA, HSA-1280218, 18 enriched genes, FDR <0.001) (Figure 5d, 30 nodes and 195 edges), where hub genes ( DMG-DEG POLA1, DMG ATR, DMG PIF1, DMG FBXO18) interacted with the strongly connected component through an E3 ubiquitin ligase that targets cell cycle regulatory proteins for degradation ( DMG ANAPC1).The targeted cluster was mainly transcriptional (both up-and down-regulated), and the few DMGs were also E3 ubiquitin ligases ( DMG-DEG HECTD1, DMG HERC1, DMG HUWE1).4 | CON CLUDING REMARK S Whole colony plasticity in the branching coral Acropora palmata resulted from integrating of modular responses to local variation in light availability.The plasticity of local modules, HL and LL phenotypes, was quantified by transplanting fragments to contrasting light environments, where significant changes in skeletal features F I G U R E 5 Responsive gene networks integrating methylome and transcriptome data: curated databases.Main sub-networks of hub genes retrieved from integration of DMGs and DEGs.(a) Photoreception-phototransduction (network from StringApp without attributes: https:// versi on11.strin g-db.org/ cgi/ netwo rk.pl? netwo rkId= sqcU0 gyKux2Y).(b) ECM proteins, cell-cell adhesion and EGF domains associated with soft tissue growth and calcification (network from StringApp without attributes: https:// versi on11.strin g-db.org/ cgi/ netwo rk.pl? netwo rkId= 8RQ Kx Pbg9zzZ).(c) Vesicle/vacuole-mediated transport, Ca 2+ metabolism and cytoskeletal protein binding associated with symbiont trafficking (network from StringApp without attributes: https:// versi on11.strin g-db.org/ cgi/ netwo rk.pl? netwo rkId= Cefq2 PjoZN5R).(d) Innate immune response associated to inter-partner recognition (network from StringApp without attributes: https:// versi on11.strin g-db.org/ cgi/ netwo rk.pl? netwo rkId= pZerN p9HZxM0).Larger nodes indicate key regulators or a critical target of a regulatory pathway.The line between genes represents interactions.Node colour indicates if DMG (yellow), both DMG-DEG (blue), DEG (grey).Font size represents methylation (signal density variation from Methyl-IT) and font colour up (red) -down (blue) regulation.Genet 1 LL to HL are shown for interpretation.promptedthe optimization of light absorption and utilization to maximize metabolic outputs and growth.Recognizing this modular concept of plasticity is essential in the face of the heritability of colony-level traits versus local module-level traits.Animal colonies consisting of many modules may remain coherent entities where colony traits have the evolutionary potential to respond to natural selection(Simpson et al., 2020).Remarkably, variants may arise locally, and single modules may have some evolutionary potential(Vasquez Kuntz et al., 2022).Local modular plasticity was accompanied by significant methylome and transcriptome modifications.Our data showed a significant light-mediated change in coral morphology, a phenotypic adjustment reflected in molecular signatures of changes in growth, including soft tissue growth and biomineralization.These observations suggest that symbiotic corals have acquired the capability of effecting an epigenomic response that incorporates whole methylome repatterning and is associated with changes in coral morphology.This interpretation aligns with previous wholegenome views that focused on the function of DNA methylation in phenotypic plasticity(Liew et al., 2018).The approach allowed us to integrate genome-wide DNA methylation with gene expression datasets into meaningful biological networks.Our comprehensive network analysis based on interactions from correlation matrix networks and gene interaction from curated databases for predicted networks offered a powerful approach to identifying potential markers of plasticity in the interaction of DMGs and DEGs.Annotated genes COL6A5, SPTAN1, PTK7, FGFR1, FBXO30, P4HA2, TNR, NID1, ITGAX, ANK3, RCHY1 and GBF1, and the transcription factor SHOX, were key regulators identified in sub-networks of hub genes in both analyses.These genes are involved in regulating F I G U R E 6 Predicted model for light-mediated phenotypic plasticity of structural traits in the branching coral Acropora palmata based on key regulators from DMGs-DEGs integrated networks.A significant change in the light environment activates photoreception mechanisms to detect cues and transduce information within cells (symbionts, cytoskeleton, extracellular matrix -ECM and nucleus-Nu are labelled).This activates signalling pathways to control growth, both soft tissue and skeletal growth; and in parallel, to initiate cellular transport related to symbiont recognition and changes in symbiont population densities (network from StringApp without attributes: https:// versi on11.strin g-db.org/ cgi/ netwo rk.pl? netwo rkId= uh6Y1 lbNXqJR).the cell cycle, ECM, vesicular trafficking, immune regulation, regulation of transcription and transduction, and biomineralization.
Corals optimize growth through morphological plasticity by means of genomic and epigenomic processes.Comprehensive biometrics and machine learning applied to methylation and transcriptional information identified key methylome-transcriptome networks containing a large number of hub genes that appear integral to inducing changes in biomineralization, polyp formation and ultimately colony shape.Moreover, phenotypic change corresponds with irradiance and substantiates the importance of light sensing in the ecology and biology of corals.These findings make significant inroads in the identification of epigenetic processes underpinning phenotypic changes and provide a roadmap for studies of non-model organisms.AUTH O R CO NTR I B UTIO N S R.I.P., I.B.B., S.A.M. and K.G.C. conceived and designed the experiments.K.G.C. performed the experiments, integrated and analysed the data.I.B.B., S.A.M., R.S. and X.Y. contributed to genomic-epigenomic analysis and interpretations.S.A.M. and R.S. to Methyl-IT and biological network analyses.R.I.P. and S.E. to phenotypic descriptions interpretation.I.M.R. assisted in fieldwork and phenotypic descriptions.T.M. assisted in methylation computational work.C.C.O.assisted in laboratory and gene expression computational work.K.G.C and R.I.P. wrote the manuscript.K.G.C., S.E., I.B.B., S.A.M. and R.I.P. reviewed and edited the manuscript.Funding acquisition: R.I.P., I.B.B. and S.A.M.All authors approved the manuscript.ACK N O WLE D G E M ENTS We thank UNAM and ICMyL for staff support, facility assistance, SAMMO data.PASPA-DGAPA to S.E. for sabbatical at PSU-Biology.M.K. Durante for assistance in laboratory and computational work.Collection permits No. SGPA/DGVS/06960/17 and No. SGPA/ DGVS/07846/17 to SE, CITES permit MX90225 to SE. FU N D I N G I N FO R M ATI O N This project has received funding from The Pennsylvania State University SEED Grant program awarded to I.B.B., R.I.P. and S.A.M.; from NSF grant (NSF OCE-1537959) to I.B.B., NIH grant (NIH R01 GM134056-01) to S.A.M., and The Pennsylvania State University Start-up to R.I.P.All data are available in the main text or the Appendix S1.All data generated by this study were deposited to the National Center for Biotechnology Information (NCBI) (GSE245494).Custom codes used for methylation analysis with Methyl-IT are available at genomaths.github.io.
Osinga et al., 2012).As a downstream calculation, the minimum quantum requirement of photosynthesis ( −1 ; photons absorbed to evolve one O 2 molecule) was estimated based on the linear regression of photosynthetic rates in the light-limited region of PE curves. −1 is defined in terms of the light being absorbed and used to drive photosynthesis or photosynthetic utilizable radiation (PUR) [ −1 = 1/(a*A PUR ); where A PUR is the average absorptance for