DNA methylation differences between stick insect ecotypes

Epigenetic mechanisms, such as DNA methylation, can influence gene regulation and affect phenotypic variation, raising the possibility that they contribute to ecological adaptation. To begin to address this issue requires high-resolution sequencing studies of natural populations to pinpoint epigenetic regions of potential ecological and evolutionary significance. However, such studies are still relatively uncommon, especially in insects, and are mainly restricted to a few model organisms. Here, we characterize patterns of DNA methylation for natural populations of Timema cristinae adapted to two host plant species (i.e., ecotypes). By integrating results from sequencing of whole transcriptomes, genomes, and methylomes, we investigate whether environmental, host, and genetic differences of these stick insects are associated with methylation levels of cytosine nucleotides in CpG context. We report an overall genome-wide methylation level for T. cristinae of ~14%, being enriched in gene bodies and impoverished in repetitive elements. Genome-wide DNA methylation variation was strongly positively correlated with genetic distance (relatedness) but also exhibited significant host-plant effects. Using methylome-environment association analysis, we pinpointed specific genomic regions that are differentially methylated between ecotypes, with these regions being enriched for genes with functions in membrane processes. The observed association between methylation variation with genetic relatedness and the ecologically-important variable of host plant suggest a potential role for epigenetic modification in T. cristinae adaptation. To substantiate such adaptive significance, future studies could test if methylation has a heritable component and the extent to which it responds to experimental manipulation in field and laboratory studies.

Methods
Timema cristinae individuals from the selected populations were collected on the same date (25th April 2017) in the Californian spring using sweep nets, and kept in plastic containers at room temperature. Individuals were digitally photographed the following day under standard conditions, flash frozen using liquid nitrogen, and preserved at -80OC temperature. All procedures were performed to assure the methylation status was not considerably affected by variation in sampling conditions.

Half of each specimen’s body (cut longitudinally) was used to isolate its genomic DNA using DNeasy Blood and Tissue Kits (Qiagen). We included non-methylated cl857 Sam7 Lambda phage DNA (Promega Corporation) a spike-in in each sample (1% of the final volume). We submitted genomic DNA of one T. cristinae sample (individual 17_0015, ‘NO_BS.WGBS’ sample) for BS-seq, and as a control for the BS-treatment (i.e., sequenced without sodium-bisulfite treatment). The sodium-bisulfite treatment and high-throughput sequencing were performed by Biomedicum Functional Genomics Unit (FuGU, Helsinki). The libraries were sequenced using the Illumina NextSeq 500 system, with High Output 2 x 150 bp runs. In total, three flow cells with four lanes were run.

In addition, we performed RNA-seq of 18 individuals, which were the only ones that yielded enough material to perform further sequencing. The RNA extractions, library preparations and sequencing were performed by Genome Quebec. Total RNA for each individual was extracted from the remaining half of the specimens’ bodies using the Quiacube animal tissue kit and protocol. Libraries were multiplexed and sequenced on one lane of HiSeq4000 to obtain 150 base pair paired-end reads. More detailed information can be found in the folders for each particular analysis in the Online Supplementary Materials.