Paternal heat exposure causes DNA methylation and gene expression changes of Stat3 in Wild guinea pig sons

Abstract Epigenetic mechanisms convey environmental information through generations and can regulate gene expression. Epigenetic studies in wild mammals are rare, but enable understanding adaptation processes as they may occur in nature. In most wild mammal species, males are the dispersing sex and thus often have to cope with differing habitats and thermal changes more rapidly than the often philopatric females. As temperature is a major environmental selection factor, we investigated whether genetically heterogeneous Wild guinea pig (Cavia aperea) males adapt epigenetically to an increase in temperature, whether that response will be transmitted to the next generation(s), and whether it regulates mRNA expression. Five (F0) adult male guinea pigs were exposed to an increased ambient temperature for 2 months, corresponding to the duration of the species' spermatogenesis. To study the effect of heat, we focused on the main thermoregulatory organ, the liver. We analyzed CpG‐methylation changes of male offspring (F1) sired before and after the fathers' heat treatment (as has recently been described in Weyrich et al. [Mol. Ecol., 2015]). Transcription analysis was performed for the three genes with the highest number of differentially methylated changes detected: the thermoregulation gene Signal Transducer and Activator of Transcription 3 (Stat3), the proteolytic peptidase gene Cathepsin Z (Ctsz), and Sirtuin 6 (Sirt6) with function in epigenetic regulation. Stat3 gene expression was significantly reduced (P < 0.05), which indicated a close link between CpG‐methylation and expression levels for this gene. The two other genes did not show gene expression changes. Our results indicate the presence of a paternal transgenerational epigenetic effect. Quick adaptation to climatic changes may become increasingly relevant for the survival of wildlife species as global temperatures are rising.


Introduction
Past epigenetic studies have mainly focused on medical and developmental aspects in humans and model species, while studies on ecologically relevant traits in nonmodel species are still scarce (Richards 2006;Pertoldi and Bach 2007;Bossdorf et al. 2008). However, to cope with environmental factors, (wild) animals need to adapt to environmental changes, for example in habitat and temperature (Kilvitis et al. 2014). The regulation of gene expression is fundamental for immediate adaptation processes in the same generation. In addition, the inheritance of responses to experienced changes (adapted traits) is fundamental for long-term adaptational memory and thus evolutionary processes. The mechanism regulating gene expression and conferring such immediate and inherited adaptation is "epigenetic response" (Jablonka and Raz 2009). This response is achieved by epigenetic marks, including DNA methylation, histone modifications, polycomb proteins, and siRNAs, which are controlling the accessibility of the chromatin to the transcriptional machinery. The best-studied epigenetic modification is DNA methylation (e.g., Boyes and Bird 1991;Tate and Bird 1993;Jones and Takai 2001;Bird 2002). In mammals, DNA methylation is mostly achieved by the covalent binding of a methyl (-CH 3 ) group to the 5-carbon site of a cytosine, forming 5-methylcytosine (5mC). It mainly occurs at single cytosine-phosphateguanine dinucleotides (CpG) of which about 80% are methylated (mCpG). CpG methylation occurring in the promoter region of a gene is often strongly associated with gene repression (He et al. 2000;Deaton and Bird 2011). Thus, this epigenetic mechanism can regulate the gene activity, without changing the DNA sequence underneath. The most intriguing aspect of this mechanism is that the epigenetic pattern is dynamically changing in response to a changing environment in the same generation ("immediate response" or "epigenetic plasticity") and can also be transmitted to the next generation ("inherited response" or "transgenerational epigenetic plasticity"). Although the exact mechanisms of transgenerational epigenetic inheritance and the differences between maternal and paternal effects are still unclear, the phenomenon itself is progressively confirmed (Weaver et al. 2004;Champagne 2008;Carone et al. 2010;Ng et al. 2010).
To ensure epigenetic plasticity through generations, DNA methylation patterns must either survive epigenetic reprogramming during early embryogenesis (Haaf 2006) or must be re-established afterward. Such "survivors" are parentally imprinted genes that are important in early embryonic development (Reik et al. 2001). However, recent studies showed that nonimprinted genes may also "survive" the reprogramming process (Borgel et al. 2010;Wei et al. 2014). Due to the close intrauterine relationship, as well as maternal care after birth, maternal effects on the offspring are more evident than paternal effects. As a result, researchers mainly focused on maternal effects and transmission (e.g., Weaver et al. 2004;Wolff et al. 1998), whereas paternal effects have been widely neglected (e.g., Carone et al. 2010;Dias and Ressler 2014). An advantage of examining paternal effects is that effects during pregnancies and behaviorally induced effects can be excluded (in many mammal species, males are not actively involved in rearing the offspring), simplifying detection of molecular mechanisms (Curley et al. 2011). Thus, to achieve a comprehensive understanding of transgenerational epigenetic effects, paternal effects need to be studied too, as well as their role in adaptation.
Besides, studying paternal effects in wild species is especially crucial, because in most wild mammal species, males are the dispersing sex and thus have to cope with more rapid changes in habitats and temperature, while females are often philopatric. This also holds true for the Wild guinea pig, Cavia aperea, which is living in harem structures with one dominant male, who defends several females from other males (Asher et al. 2008). The nondominant, roaming males thereby need to adapt quickly to new habitats and temperatures before finding accepting female(s). Climate changes have always impacted on the evolution of species, and ambient temperature is one major external selection factor animals have to cope with. The species long-term survival depends on their shortterm responses. The ability to convey heritable phenotypic plasticity to descendant generations will likely contribute to fitness increase and thus to improved chances of survival.
We here investigated the potential paternal contribution to an epigenetically adaptational response to increased ambient temperature in a wild species. Therefore, as described in detail in a recent publication (Weyrich et al. 2015), we temporally exposed adult male Wild guinea pigs to an increased ambient temperature and let them mate to the same females before and after the heat exposure. We then analyzed whether their sonsone group sired before, the other after the fathers' heat treatmentshowed altered DNA methylation patterns compared with the situation prior to the heat exposure (experimental setup in Fig. S1). To study the epigenetic responses to heat, we here focused on the main metabolic, heat-producing and thus thermoregulatory organ, the liver. Mammals, such as C. aperea, are homeostatic; that is, they compensate external thermal fluctuations (thermoregulation) to maintain internal temperature (homeostasis). In endotherms, thermoregulation is a crucial and complex process to ensure optimal energy usage at a given temperature (Sonna et al. 2002;Boyles et al. 2011). It is regulated by many genes of which the Signal Transducer and Activator of Transcription 3 gene (Stat3) (also named: Acute Phase Response Factor) is crucial for maintaining temperature homeostasis (Sonna et al. 2002). Stat3 becomes activated by various cytokines and regulatesdepending on the tissue typethe expression of genes essential for systemic thermoregulation, embryogenesis, and immune response (Durant et al. 2010;Qi and Yang 2014). Due to its role in many different physiological processes, Stat3 is an important target for investigation.
In this study, we aimed to examine whether heat exposure of fathers would affect the methylation patterns in their F1 sons (measured in comparison with F1 sons sired either before (F1 C ) or after paternal heat exposure (F1 H )) and whether these altered patterns would influence gene expression. With our approach, we address both (1) the paternal contribution of epigenetic transmission and (2) understanding the epigenetic adaptability in the context of globally rising temperatures. Both may contribute to a fitness increase of the offspring and would thus play an important role in evolution.

Materials and Methods
Animal care and treatment All husbandry and experimental procedures were approved of by the German Committee of Animal Welfare in Research (permit no. V3-2347Research (permit no. V3- -35-2011. Wild guinea pigs originating from Argentina and Uruguay (Asher et al. 2008) were obtained from F. Trillmich (University of Bielefeld) and housed at the IZW field station in Niederfinow, Germany. Wild guinea pigs are polyoestric. Male cavies were born in mid-November 2010 and females in April-May 2011. Mating took place in January/February 2012 [control] and in September 2012 [heat]. Animals were housed in shorttunnel-connected indoor-outdoor-enclosures (1.3 m 2 ) under natural photoperiod and temperature conditions. All animals were fed guinea pig pellets (Altromin Spezialfutter GmbH & Co. KG, Lage, Germany). Water and hay were provided ad libitum, and supplementary apples, peppers, or carrots were given daily. Vitamin C was added to the drinking water once a week. During a two-month period, the indoor cages (80 9 80 cm 2 ) were placed on a heating plate (Candor GmbH Leipzig) which heated the floor to a temperature of 30°C (covering a full cycle of spermatogenesis) (Fig. S1). The slightly cooler edges of the heating plate were fenced off by a mesh, installed in 5 cm distance from the cage wall, reducing cage size to 70 9 70 cm 2 . The experiment to assess the influence of environmental stressors was carried out with five adult males kept exclusively indoors (60 days, 24 h; 30°C).
Each of the five males mated with two females before (control) and after heat exposure (heat) by introducing them to the female's cage. After an observed mating, males were transferred back to their own cage. Males were mated with the same two females before and after heat exposure. Mating resulted in 16 F1 sons from the control mating (F1 C N = 16) and 18 from the mating after heat exposure (F1 H N = 18). Seven days after birth, whole livers (the main thermoregulatory and metabolic organ) of F1 sons were harvested and individually sequenced (N = 34).

DNA methylation level analysis and group comparison
The bisulfite conversion rate was calculated as the number of mapped nonmethylated CpGs divided by the total number of mapped CpGs. After bisulfite treatment and alignment to the reference sequence, Cs in a read that mapped to a C in the reference were assumed to have been methylated cytosines (mCs). Thymines that mapped to a C position were regarded as formerly unmethylated cytosines (Cs), respectively (converted to U by bisulfite and substituted by subsequent PCR). Accordingly, the methylation level of each cytosine position was calculated as the number of reads mapping to this position and carrying a C divided by the number of reads carrying either C or T at this position.
Methylation level was calculated by the equation: We compared methylation states of sons before (control) and after the fathers' heat treatment. Sites with coverage of ≤5 per sample per group were excluded from the analysis. Significance of differences was determined using Fisher's exact test. Strong hyperand strong hypomethylation were defined as absolute methylation difference greater than 30% (Gu et al. 2010).

RNA extraction and cDNA synthesis
RNA was extracted from 20 mg of liver tissue of the sacrificed guinea pigs using the innuSPEED Tissue RNA kit (Analytik Jena, Germany) and peqGOLD DNase I Digest Kit (Peqlab, Erlangen, Germany). RNA extracts (elution with 50 lL H 2 O) were stored at À80°C. RIN was determined as an objective indication of RNA quality (Schroeder et al. 2006) using the Agilent 2100 bioanalyzer and a RNA 6000 Pico Kit following the manufacturer's guidelines. RIN values for all samples were between 5.9 and 8.3 (SD = AE0.7). RNA samples were diluted to a final concentration of 25 ng/lL in a volume of 22 lL. The Revert Aid First Strand cDNA Synthesis Kit (Thermo Scientific, Munich, Germany) and its associated protocol were applied, including 100 mM Oligo(dT) 18 . The obtained cDNA was frozen and stored at À20°C.

Primer design and verification
Primers (Table 1) were designed based on the C. aperea whole genome sequence (Weyrich et al. 2014) using the Geneiousâ software. For maximum quantification accuracy, primers producing amplicons of 120-200 bp were chosen with an optimal annealing temperature (T a ) of 55°C and a length of 20+/À 2 nucleotides. Theoretical T a values were calculated with the CFX Manager Software (Bio-Rad GmbH, Munich, Germany). The cDNA template for PCR testing was freshly synthesized from randomly selected samples. Negative controls were set up for each primer pair. A positive PCR control was prepared using the Stat3 primer pair. For visualization of PCR products, a 1.5% agarose gel was prepared (1.5 g agarose GOLD â (Peqlab) + 100 mL 1x TBE buffer, 3 lL GelRed TM (Biotium, Cologne, Germany)). 5 lL PCR product was added to 1 lL loading dye (6x orange DNA loading Dye; Fermentas, St. Leon-Rot, Germany). The 50bp ladder (O'GeneRuler TM 0.1 lg/lL, Thermo Scientific) was used for size determination. PCR reaction conditions for primer pairs were optimized. Amplicons were verified in size using agarose gel electrophoresis as well as by Sanger sequencing (Sanger et al. 1977) using the BigDye â Terminator v3.1 Cycle Sequencing Kit and BigDye â XTerminator TM Purification Kit (Life Technologies, Darmstadt, Germany). Capillary electrophoresis was performed using the ABI Prism 3130xl Genetic Analyzer and associated Collection software v. 3.0 (Life Technologies). Sequences were analyzed using the Sequence Analyser â v. 5.2 (Life Technologies) and BioEdit v. 7.2.5 software (Hall 1999). NCBI BLAST search finally confirmed sequence specificity.

Real-time quantitative PCR
Quantitative PCR experiments were performed on the CFX96 Real-Time System C1000 Thermal Cycler equipped with CFX Manager software (both Bio-Rad). As standard curve, a 1:6 mixture of all 34 samples with five serial dilutions (in steps of 1:10) allowed calculation of run efficiency (E). Samples, standards, and controls were run in duplicates and prepared using the SsoFast Eva-Green TM Supermix on white optical plates. All qPCR runs were prepared by the same analyst and performed on the same platform to minimize inter-run variation, thus minimizing distortion of results (Derveaux et al. 2010).

Real-time quantitative PCR data analysis
Run quality was assessed based on the efficiency (E) and r 2 values generated by the software. E-values > 90% and r 2 > 0.98 (no between-cycle variation in efficiency) were accepted. If the standard deviation between duplicates was > 0.4, samples were repeated. Efficiency-corrected expression values (Q) (Pfaffl et al. 2002) were calculated from the software-generated cycle threshold values (Ct) according to the equation: Q = E ÀCt Relative quantification relies on the normalization of the expression values for the target genes with reference genes, accounting for intraindividual and technical variation. Reliable normalization requires the use of several reference genes showing stable expression for the tested tissue and treatment. Normalization of the Stat3 expression values was performed using the software GeNorm to define a combination of three reference genes resulting in the most stable expression levels (Weyrich et al. 2010). Eventually, we used Hmbs (efficiency E = 110.4%), B2 m (E = 108.8%), and Gapdh (E = 90%) ( Table 1). For each sample, a normalization factor (NF) was calculated based on the geometric mean (GM) of the Q values for the three reference genes (Qr) (Vandesompele et al. 2002): Normalization factor was then applied to expression values for the target genes using the efficiency-corrected DCt method (Pfaffl et al. 2002), giving the normalized expression values (NE) used for statistical analysis: With Q1 = efficiency-corrected expression value for target gene 1, NF1 = normalization factor for sample 1.

Statistical analysis
We applied Fisher's exact test for methylation levels and a paired Wilcoxon signed-rank test to test for significant differences between F1 C and F1 H for normalized expression values using R v.3.1.1 (http://cran.r-project.org/bin/ windows/base/). A paired Wilcoxon signed-rank test was also applied to investigate a potential maternal effect; wherefore, we first averaged the expression values of the sons per father-mother combination (n = 10). The plots were generated in R with ggplot2 v1.0.0 and are displayed as whisker plots with 1.5*IQR (http://cran.r-project.org/ web/packages/ggplot2/citation.html). Normalized expression values were also tested for significant associations with specific methylation changes using Pearson's test in R.

Methylation changes
After comparing CpG-methylation levels of the control group with those of the heat group (F1 C vs. F1 H ), we detected significant differences in several genomic loci (P < 0.05; Fisher's exact test), indicating paternal epigenetic transgenerational inheritance. We then focussed on the three genes with the greatest abundance of significant mCpG-differences (occurring in promoters as well as in coding sequences) between control and heat group: the thermoregulatory gene, Signal Transducer and Activator of Transcription 3 (Stat3); the Cathepsin Z gene (Ctsz), a lysosomal cysteine proteinase with function in tumorigenesis; and the Sirtuin 6 gene (Sirt6), a gene with function in intercellular and epigenetic regulation (Fig. 1). We investigated methylation levels for each CpG in F1 C (N = 16) and F1 H sons (N = 18) (Table S1).
In the Stat3 gene, six CpGs were differentially methylated, one in the promoter and five in the coding region ( Figs. 1 and 2). While five of the six CpGs showed greater methylation levels in the F1 H group, one CpG had a higher methylation level in the F1 C group (Fig. 2).
Similar to the Stat3 methylation pattern, in Ctsz, hypermethylation of CpGs was also mostly detected in the F1 H group (in 10 of the 11 CpGs). Heat exposure of fathers led to the hypomethylation of CpG 11 in F1 H . Interestingly, this CpG was the only one located in the promoter region. In Sirt6, all differentially methylated CpGs were located in the gene's promoter region. Four of these CpGs were hypomethylated and two were hypermethylated in F1 H individuals relative to control individuals, reflecting varying levels of methylation.

Gene expression analysis
Measuring expression levels in both F1 C and F1 H using quantitative PCR (qPCR) revealed significant differences in gene expression of Stat3, but not of Sirt6 or Ctsz ( Fig. 3; Table S1). Stat3 mRNA levels were significantly lower in F1 H than in F1 C (mean F1 C = 0.635, SD = AE0.233; mean F1 H = 0.376, SD = AE0.289; paired Wilcoxon singed-rank test P = 0.01099). By applying a paired Wilcoxon signed-rank test to the averaged expression values of the sons per father-mother combination (n = 10), we verified treatment specificity and excluded maternal effects on gene expression (P = 0.03125). These results showed that DNA methylation in F1 sons was paternally influenced, indicating transgenerational epigenetic plasticity.

Linking Stat3 DNA methylation and gene expression
Focussing on Stat3, we applied a Pearson correlation to investigate the association between paired mRNA expression values and the methylation levels at each of the six single CpGs (Figs. 1 and 2, Table 2). There were clear, yet weak associations between methylation differences at each of the tested CpG positions (CpG 1 to CpG 6) and gene expression changes. Differential methylation at position 1 and 3 had the strongest (negative) association with gene expression (r 1 = À0.371 and r 3 = À0.418), followed by CpG 4 and CpG 6 (r 4 = À0.167, r 6 = À0.168).
[Correction added on 18 March 2016, after initial online publication on 28 February 2016. Citation and placement of tables 1 and 2 and reference to 'Weyrich' in the Abstract were incorrect and have been updated in this version.]

Main findings
Our results demonstrate a paternal effect on DNA methylation patterns of the next generation after exposure to increased temperature. The investigation of three genes revealed that one of them Stat3, a key thermoregulation gene involved in numerous biological processes was not only differentially methylated between F1 C and F1 H siblings, but also significantly differentially expressed (P < 0.05). In F1 H, Stat3 was hypermethylated in comparison with F1 C and gene activity reduced.

Regulation of expression
Our findings add to the evidence (Ng et al. 2010;Wei et al. 2014) for paternally inherited changes in DNA methylation patterns and their impact on gene activity. Although the exact mechanism remains still unknown, the genomic location of CpGs at which methylation changes take place seems crucial, because promoter methylation has often been strongly associated with gene repression (He et al. 2000;Deaton and Bird 2011), while intragene methylation can have either a gene silencing or an activation function (Hahn et al. 2011;Jjingo et al. 2012). Some studies describe one CpG/mCpG position as being sufficient for gene regulation (e.g., Ng et al. 2010), others suppose that stretches of differentially methylated CpGs, so-called differentially methylated regions (DMRs), are the main regulatory units (Radford et al. 2014). Besides direct effects, DNA methylation can inhibit transcription also indirectly via a methyl-binding (MBD) protein (Boyes and Bird 1991) and eventually interact with other epigenetic modifications, such as histone modifications and RNAi (Wei et al. 2014).

Stat3 -Correlation results
In our study, Stat3 differential CpG methylation for CpG 1-4 and CpG 6 was negatively associated with the gene's expression level (all hypermethylated in F1 H ) and positively associated (albeit weakly) for CpG 5 (hypermethylated in F1 C ), which was the only one located in a CG-poor region ( Fig. 2A). The strongest associations were detected for CpG 3 and CpG 1, both of which were located in the transcribed region of the gene. The site located in the promoter region, CpG 6, had only a weak negative correlation to gene expression. This was unexpected because methylation of mainly promoter regions had been shown to influence gene expression (Deaton and Bird 2011;He et al. 2011). However, in case of Stat3, the CpGs located in the transcribed region of the gene are likely candidates for having caused the expression differences observed. CpG methylation in transcribed regions is associated with both activation and deactivation of gene expression: in cells with low proliferation rates such as neural cells, it is linked to gene deactivation. In fastdividing cells such as liver cells, it is associated with gene activation and thus with higher expression levels (Hahn et al. 2011;Jjingo et al. 2012). The mechanism, by which intragenic methylation within the coding region increases a gene's activity, is currently poorly understood. It has often been suggested that it prevents uncontrolled transcription and promotes transcriptional elongation (Hellman and Chess 2007;Aran et al. 2011). It may also well be that the different CpGs act collectively, which would indicate a more complex interplay between the methylations/demethylations at different positions. Correlations between specific DNA methylation changes and gene expression changes were weak. A reason might be the limited number of animals studied. However, despite the Wild guinea pig being a heterogeneous mammal species with phenotypic, genetic, epigenetic, and transcriptional variability and an average genetic dissimilarity among animals of~0.17% (Weyrich et al. 2015), we found a significant response across individuals in the Stat3 gene. This response across individuals of Stat3 indicates an ecological relevance. It also identifies Stat3 to be a suitable candidate for studies on globally changing temperatures.
The differentially methylated CpGs in the Stat3 gene accompanied with gene expression changes (lower expression of Stat3 in F1 H liver) suggest that the paternal experience transmitted via epigenetic transmission might contribute to a fitness increase of the offspring.
The gene product of Stat3 itself is also interacting with epigenetic factors. In T cells, STAT3 activates transcription of genes in co-regulation with histone 3 lysine 4   trimethylation (H3K4me3), a marker with known transcription-activating function (Durant et al. 2010). We cannot fully exclude an age effect on DNA methylation as F0 males grew~7.5 months older until the second mating after heat exposure. Studies on monozygotic human twins have shown that young twins (3 years) are indistinguishable in their DNA methylation patterns, whereas older twins (50 years) show substantial methylation changes (Fraga et al. 2005). Methylation patterns remain relatively constant for long periods, and those age-related DNA methylation changes were observed after 15-20 years (Fraga et al. 2005, Fig. S5). Even when taking the shorter life span of the Wild guinea pigs into account, significant age-related DNA methylation changes are not expected within a period of 7.5 months during adulthood.

Ctsz and Sirt6
Ctsz and Sirt6 were not affected in their gene activity, even though methylation had changed between F1 C and F1 H . This was particularly surprising for Sirt6 where all six differentially methylated CpGs were located in the promoter region.
However, recent studies have demonstrated that DNA methylation does not always correlate with gene expression (Huang et al. 2014), which thus indicates a more complex regulatory mechanism. We hypothesize that the lack of expression changes might reflect a regulatory mechanism, in which an increase in exposure time is positively correlated to an increase in the number of CpGs to be altered in their methylation in regulatory regions. Reaching this threshold of differentially methylated CpGs may alter gene expression (dosage effect). Following our hypothesis, in case of Ctsz and Sirt6, methylation changes were thus not strong enough to be of biological meaning. Changes might also occur later in life or in case the sons would again be exposed to temperature increase. Other reasons might be the age of the animals at the time of sampling. As the exact molecular and mechanistic basis for these observations remains largely unclear, further studies will be needed to address these questions.

Comparison with other studies
As this is the first study examining the molecular effects of paternal heat stress on subsequent generations in a nonmodel mammalian species, comparisons with other studies in the field of paternal epigenetic transgenerational inheritance are limited. In rats, paternal high-fat diet increased the risk of obesity and diabetes in female offspring (Ng et al. 2010), and low-protein diet influenced hepatic gene expression and metabolic function in both female and male offspring (Carone et al. 2010). Behavioral conditioning in male mice changed gene expression and neurophysiology in male and female offspring (Moore et al. 2013).

Evolutionary relevance
The evolutionary significance of epigenetic transgenerational inheritance of ecologically relevant traits is a subject of debate (Lachmann and Jablonka 1996;Takai and Jones 2002;Root and Schneider 2006;Moore et al. 2013). While the underlying mechanisms have not been fully understood, our results suggest that acquired responses to heat exposure are transmitted to the next generation in a nonmodel species. Its transmission to the next generation and this seemingly targeted response to heat is the prerequisite for epigenetically mediated microevolution. For a comprehensive assessment of the evolutionary significance of epigenetic transgenerational inheritance, long-term studies including several subsequent generations will be needed.

Outlook
Thermoregulation is of great ecological and evolutionary importance. Transgenerational epigenetic inheritance may be of great advantage to adapt to rapidly as well as persistently rising temperatures (Lachmann and Jablonka 1996;Root and Schneider 2006). In times of globally rising temperature, this mechanism might become even more important for species survival.

Supporting Information
Additional Supporting Information may be found in the online version of this article: