Role of gene body methylation in coral acclimatization and adaptation

Gene body methylation (GBM) has been hypothesized to modulate responses to environmental change, including transgenerational plasticity, but the evidence thus far has been lacking. Here we show that coral fragments reciprocally transplanted between two distant reefs respond with genome-wide increase or decrease in GBM disparity among genes. Surprisingly, this simple genome-wide adjustment predicted broad-scale gene expression changes and fragments’ fitness in the new environment. This supports GBM’s role in acclimatization, which may consist in modulating the expression balance between environmentally-responsive and housekeeping genes. At the same time, constitutive differences in GBM between populations did not align with plastic GBM changes upon transplantation and were mostly observed among FST outliers, indicating that they arose through genetic divergence rather than through transgenerational inheritance of acquired GBM states. One-sentence summary Genome-wide shifts in gene body methylation predict gene expression and fitness during acclimatization but do not contribute to epigenetic divergence between populations.


Abstract 21
Gene body methylation (GBM) has been hypothesized to modulate responses to environmental 22 change, including transgenerational plasticity, but the evidence thus far has been lacking. Here 23 we show that coral fragments reciprocally transplanted between two distant reefs respond with 24 genome-wide increase or decrease in GBM disparity among genes. Surprisingly, this simple 25 genome-wide adjustment predicted broad-scale gene expression changes and fragments' fitness 26 in the new environment. This supports GBM's role in acclimatization, which may consist in 27 modulating the expression balance between environmentally-responsive and housekeeping 28 genes. At the same time, constitutive differences in GBM between populations did not align 29 with plastic GBM changes upon transplantation and were mostly observed among F ST outliers, 30 indicating that they arose through genetic divergence rather than through transgenerational 31 inheritance of acquired GBM states. 32 33 GBM is a taxonomically widespread epigenetic modification the function of which 34 remains unclear (1, 2). GBM is bimodally distributed among genes: it is high in ubiquitously 35 expressed housekeeping genes and low in inducible genes (2, 3). Only the detrimental effect of 36 GBM is well understood: GBM causes hypermutability in protein-coding regions (4). Indeed, in 37 humans, GBM is the primary driver of deleterious parent-age-related mutations (5).  GBM is more consistent across fragments of the same original colony than GE, resulting in the 109 estimate of broad-sense heritability of 0.79 compared to 0.64 for GE (Fig. 1G, Fig. S3 and S4). In 110 addition, significant effects of origin and transplantation site were observed for GE and GBM, 111 both more pronounced for GE (Fig. 1 E and G). 112 Surprisingly, the change in GBM in response to transplantation from Orpheus to Keppel 113 consisted mainly in genome-wide reduction of disparity between high-and low-methylated 114 genes: highly-methylated genes became less methylated and low methylated became more 115 methylated (Fig. 2 A). This change was mirrored by less pronounced but clearly reciprocal GBM 116 adjustment in Keppel corals transplanted to Orpheus (Fig. 2 B, D) and recapitulated the 117 difference between fragments planted in their native environment (Fig. 2 C). Both exons and 118 introns underwent this change (Fig. S5). Absolute changes in methylation included both positive 119 and negative shifts ( Fig. S6 A, B), indicating that the change in relative GBM levels among genes 120 is not due to biased genome-wide methylation increase or decrease. Methylation levels of 121 intergenic regions and repeated elements remained relatively constant compared to genic 122 regions ( Fig. S6 C, D). 123 GE response was also reciprocal between transplantation directions but demonstrated 124 much broader range of fold-changes compared to GBM response (Fig. 2 D). Remarkably, gene 125 expression changes correlated with low-and high-methylated gene classes and were in the 126 opposite direction relative to GBM changes (Fig. 2 F, G). Although at the gene level the 127 correlation between GE and GBM was weak (but still highly significant, Fig. S7), it was very 128 strong when comparing broader functional groups of genes (Gene Ontology categories, Fig. 2 H  129 and Fig. S8). 130 To characterize GBM disparity among genes in each sample, we computed GBM class 131 difference ("GBMcd") as difference in mean log 2 (GBM) between high-and low-methylated gene 132 classes (Fig. 2 I). This measure aligned nearly perfectly with the first principal component of 133 GBM variation for the whole experiment, explaining 33% of total variation ( Fig. 1 F Fig. 3 A), and were positively rather than negatively 162 associated with differences in GE (Fig. 3 B). In addition, genes showing significant constitutive 163 differences in GBM also showed elevated F ST (Fig. 3 C) unlike genes undergoing significant plastic 164 changes (Fig. S9). This indicates that the mechanism giving rise to constitutive GBM changes 165 between populations is unrelated to the mechanism of plasticity and is most likely genetic 166 divergence. 167

Sample sizes 276
The experiment started with 15 colonies from each of the two sites, which were halved 277 and transplanted, resulting in a total of 60 fragments distributed across four study groups 15 278 samples each: two "natives" (KK and OO), and two "transplants" (OK and KO, Fig. 1C). Although 279 not all original samples were successfully analyzed using all approaches (Table S1), in the end 280 each of these four groups included 11 fragments representing unique coral genotypes and 281 paired by genotype between "home" and "away" groups that were analyzed for GBM, gene 282 expression, genotype, and fitness proxies. Adapter trimming and quality filtering reduced these the total read count 940 million, mean = 386 15 ± 0.64 SEM million reads per sample. The reference genome and annotations (version 1.1) 387 for Acropora digitifera (31) were downloaded from NCBI. Trimmed and filtered reads were 388 mapped to this concatenated reference using Bowtie2 (32). To ensure that using a reference 389 from an alternate species did not severely impair mapping, we compared the mean mapping 390 efficiency against the A. digitifera reference with that of a draft genome sequence for A. 391 millepora produced by David Miller and coworkers (James Cook University). Mean mapping 392 efficiency against the A. digitifera (78.4 ± 0.7% SEM) reference was only 5.2% lower than that of 393 A. millepora (83.6 ± 0.8% SEM), hence sequence differences between the two species did not 394 appear to substantially impair read mapping. Following mapping, PCR duplicates were removed 395 using Picard (https://broadinstitute.github.io/picard/). Mean duplication frequency was 14.3 ± 396 0.5% SEM. The BAM files were filtered to retain only highly uniquely aligned reads (mapping 397 quality 30, or 0.1% chance of non-unique alignment) using samtools (33), after which the 398 number of reads overlapping various genomic regions -exon, intron, genic (exon + intron), 399 repeated elements longer than 500b, and repeat-free intergenic regions at least 1 kb long and at 400 least 2 kb away from any gene -were counted using BEDtools (34). where, ! . is its distance along the discriminant axis between sample X and the mean DAPC 459 value for natives of its transplantation site, ! is the mean absolute distance for all transplants, 460 and s is the standard deviation of absolute distance for all transplants. Correlations of these 461 values with fitness characteristics of transplanted fragments (Fig. 4 E) were computed using 462 function envfit (R package vegan (42)) 463 464

Validation of MBD-seq results by bisulfite amplicon sequencing 465
To validate our MBD-seq results, we used targeted bisulfite sequencing. Genomic DNA 466 for this procedure was extracted and purified as described above. It should be noted that as 467 with the MBD-seq and Tag-seq, targeted bisulfite sequencing was performed from separate 468 isolations from the same or closely adjacent branch. We performed bisulfite conversion using an 469 EZ DNA Methylation-Gold kit (Zymo Research; cat. no. D5005). Thermocycler conditions for 470 conversion were 98°C for 10 minutes followed by 53°C for 4 hours. Primers for post-conversion 471 amplification were designed using Bisulfite Primer Seeker (Zymo Research; 472 http://www.zymoresearch.com/tools/bisulfite-primer-seeker). Primer sequences are given in 473 Table S1. We designed 13 primer sets to target coding sequences (excised from the A. digitifera

Analysis of genetic differentiation 496
The MBD-seq reads were also used to analyze genetic distances among the 22 497 sequenced colonies. For this procedure, we concatenated reads from clone pairs into single 498 files, and mapped the reads to the A. digitifera reference genome (31) using Bowtie2 (32). For 499 the subset of 12 samples where we sequenced both captured and flow-through fractions, both 500 sets of reads were concatenated for genotyping. Between-individual genetic distances for DAPC 501 analysis were determined as (1 -Identity-By-State) using ANGSD (44) with the following filter 502 settings: -uniqueOnly 1 -remove_bads 1 -minMapQ 20 -minQ 30 -baq 1 -minInd 10 -snp_pval 1e-503 2 -minMaf 0.1 (to summarize, only high-quality and uniquely-mapped reads were taken into 504 account; only sites covered in a minimum of 10 individuals were considered; the cutoff for the 505 p-value for the SNP being true was set at 0.01; and only alleles with frequency equal or 506 exceeding 0.1 were used). This resulted in genotyping information for roughly 507,000 507 polymorphic sites. The between-population F ST for genes passing 10% FDR cutoff for either for 508 origin or transplant effect and a random subset of 500 non-significant genes (rather than whole 509 genome, to speed up computation) were determined using ANGSD based on site frequency 510 spectra as priors, which were generated for every gene subset without using filters that would 511 affect representation of alleles of different frequencies (-uniqueOnly 1 -remove_bads 1 -512 minMapQ 20 -minQ 30 -baq 1 -minInd 4). 513 514

Determination of relative proportions of Symbiodinium clades 515
The MBD-seq reads were mapped to the concatenated reference including A. digitifera 516 genome and four Symbiodinium sp. trasncriptomes from four different genotypic groups 517 ("clades"). Transcriptomes for Symbiodinium clades A and B were from (45) and transcriptomes 518 for clades C and D were from (46). We then counted the relative proportions of reads producing 519 highly unique matches (mapping quality 40 or higher) to each Symbiodinium transcriptome, 520 using a custom perl script zooxType.pl. All corals were found to be dominated by  sharing the same first letter and the same number are clonal fragments from the same colony. 609 All of the 22 clone pairs except one (sample pair KK4 and KO4) clustered together. 610 number are clonal fragments from the same colony. All 24 available clone pairs clustered 618 together (six samples lacked data for clone pairs). 619  (Table S1).