Challenges and Perspectives in the Epigenetics of Climate Change-Induced Forests Decline.

Forest tree species are highly vulnerable to the effects of climate change. As sessile organisms with long generation times, their adaptation to a local changing environment may rely on epigenetic modifications when allele frequencies are not able to shift fast enough. However, the current lack of knowledge on this field is remarkable, due to many challenges that researchers face when studying this issue. Huge genome sizes, absence of reference genomes and annotation, and having to analyze huge amounts of data are among these difficulties, which limit the current ability to understand how climate change drives tree species epigenetic modifications. In spite of this challenging framework, some insights on the relationships among climate change-induced stress and epigenomics are coming. Advances in DNA sequencing technologies and an increasing number of studies dealing with this topic must boost our knowledge on tree adaptive capacity to changing environmental conditions. Here, we discuss challenges and perspectives in the epigenetics of climate change-induced forests decline, aiming to provide a general overview of the state of the art.

Climate Change-Induced Tree Dieback and Forests Decline

Rising temperatures are exacerbating drought stress for trees at global scales, causing events of decline in several species (usually known as forest dieback; Allen et al., 2010). Climatic projections indicate that mean global temperatures are expected to rise at least 1.5°C until 2100 in response to anthropogenic emissions (IPCC, 2018). Such projections of warming indicate the duration and severity of droughts are expected to amplify, probably leading to widespread forest dieback and tree mortality, affecting tree species currently not assumed to be prone to drought stress, for instance in tropical or boreal forests, or already affected by drought stress, such as those of the Mediterranean forests (Allen et al., 2010). Even though several mechanisms have been described as physiological drivers of dieback, including hydraulic failure and carbon starvation (McDowell et al., 2016; Choat et al., 2018), we lack forecasting tools with sufficient predictive power, since the likelihood of drought-induced tree dieback and mortality is not straightforward (Trugman et al., 2021). This mismatch enhances the need to understand the processes that modulate contrasting adaptive capacity among trees within a given declining forest in order to be able to provide, for instance, adaptive management strategies. Epigenetic processes are known to be able to modulate the patterns of gene expression in plants (Rapp and Wendel, 2005), some being involved in physiological and morphological responses to drought (e. g., Zheng et al., 2013; Liang et al., 2014; Chwialkowska et al., 2016). As drought-induced tree mortality has far-reaching consequences that affect a wide-range of fields, from environmental conservation to climate change mitigation efforts, and it is expected to be more frequent under a climate-change scenario, we need to improve our knowledge about the degree to which epigenetic mechanisms are linked to mortality predictions (e. g., Sow et al., 2018; Amaral et al., 2020; González-Benito et al., 2020).

Trees are capable of coping with short drought periods by avoiding water loss. Thereby, several morphological and physiological adaptations have been developed by forest tree species, for instance variations in tracheids size (Carvalho et al., 2015), reduced leaf area, embedded stomatas, and quick stomata close (Sánchez-Salguero et al., 2015; Kijowska-Oberc et al., 2020). However, long-lasting and severe drought periods are predicted to occur more often, so these strategies might not be enough for their survival (IPCC, 2018).

Tree growth can be used as a marker to study forest tree response to environment alterations, as it is affected by drought stress. Dendrochronology studies have evinced growth decline during drought stress periods, which eventually leads to forest dieback (Camarero et al., 2015; Cailleret et al., 2017). Thus far, drought-induced forests dieback has been observed in several conifer species found in the Mediterranean region, like Abies pinsapo Boiss. (Linares et al., 2010), Abies alba Mill. (Macías et al., 2006; Camarero et al., 2011, 2015), and Cedrus atlantica (Endl.) Manetti ex Carrière (Linares et al., 2011). The Mediterranean basin is considered a climate change hotspot, being one of the most vulnerable areas in the world. As this region was an important refuge in the last glaciation, a large number of endemic and relict species can be found there. Hence, variations in the environment could result in death of these drought-sensitive species which will trigger a huge loss of biodiversity.

In addition, large-distributed forest tree species, like Pinus sylvestris L., often live under different environment conditions (Bose et al., 2020), so individuals must show variations that allow them to cope with those circumstances (Camarero et al., 2015; DeSoto et al., 2020). Studying not only endemic species but also largely distributed ones and identifying these differences could be key to understand why some individuals survive while others die (Cailleret et al., 2017).

Surprisingly, there is a lack of data about epigenomic changes in natural populations in spite of their potential role in adaptation to a rapidly changing environment. The most common studies compare individuals under drought conditions with regularly watered individuals, grown in greenhouses (e.g., Gourcilleau et al., 2010; Xu et al., 2018; Li et al., 2020). While this approach allows to accurately control the factors affecting the plants and makes the experiments replicable, it may lack realism, especially when working with long-lived perennial plants due to restriction of their rooting volume, which hinders the extrapolation of the results to natural populations (Gibson et al., 2004). Therefore, studying epigenetic variation in natural populations would be of great interest to understand how trees respond to climate change in the nature.

On summary, as said before, growth of forest tree species is usually slow and the upcoming climate variations will likely be faster than the speed of allele frequencies shift of forest trees, which emphasizes a likely role of epigenetic modifications as a more rapid way of adaptation.

Epigenetic Mechanisms Under a Climate Changing Scenario

Forest trees are sessile organisms which must cope with environment variations since they are not able to migrate to ensure their survival. Hence, they possess other abilities to compensate for their lack of mobility, such as phenotypic plasticity (Petit and Hampe, 2006). Phenotypic plasticity is the capacity to develop phenotypic modifications in response to multiple environmental factors, increasing the probability of surviving (West-Eberhard, 2008). This mechanism is the result of molecular modifications that trigger changes in gene expression. One of the most important processes involved in phenotypic plasticity is epigenetic regulation, which produces changes in genes expression without modifying the DNA sequence (Iwasaki and Paszkowski, 2014). In fact, phenotypic plasticity is one of the main biological processes that allows trees to respond to rapid variations, like climate change (Alberto et al., 2013). Climate change exacerbates drought stress for trees at global scales, causing stand-level growth decline and episodes of mortality (Allen et al., 2010; Anderegg et al., 2013). As drought-induced forests decline and mortality are reported at an unprecedented rate (Jump and Peñuelas, 2005), a particularly urgent scientific challenge is to understand epigenetic mechanics linking an environmental cue, here increasing drought stress, received by adult trees and its resulting effects on the phenotype of coming generations (Alberto et al., 2013).

In particular, conifers, a group of forest tree species which include ecologically and economically important genera like Pinus, Picea, and Abies, seem to be particularly drought-sensitive (McDowell et al., 2016; DeSoto et al., 2020) and prone to a high rate of epigenetic regulation due to their large genome sizes with high amounts of non-coding DNA (Yakovlev et al., 2012). The study of growth patterns in declining conifer forests suggests a lower capacity of resilience in stressed/dying trees. In fact, trees that died after drought were less resilient after previous droughts, relative to surviving conspecific individuals. Specifically, drought-related mortality risk in conifers has been associated with a low capacity to recover after drought (DeSoto et al., 2020), whereas the epigenetic mechanisms involved on this recovery mechanisms remain largely unknown. Here, we hypothesize that epigenetic mechanisms are likely related to the fate of trees that survive after drought, while some conspecifics died.

Three relevant epigenetic mechanisms have been described: DNA methylation, histone modifications and non-coding RNAs (ncRNAs). Of these, DNA methylation is the most studied one in plants, both model and non-model (Zhang et al., 2018). Two main reasons are responsible of this fact. In one hand, DNA methylation has been reported to be implicated in significant biological processes, like silencing of Transposable Elements (TEs), gene transcription regulation, plant development, and response to biotic and abiotic stress (Zhang et al., 2018). On the other hand, several genomic tools have been developed to easily study at a global level the methylations present in the DNA. Therefore, in this review we will pay special attention to this epigenetic mark.

Regarding histone modifications, these epigenetic marks are involved in chromatin compaction, playing a role in determining its state, which can be decondensed (euchromatin) or condensed (heterochromatin). Euchromatin, also known as open chromatin, is characterized by the presence of a high number of genes, which are usually susceptible to transcription as they are accessible to transcription machinery, while heterochromatin, also known as closed chromatin, is characterized by repetitive and silenced sequences (Amaral et al., 2020). Acetylations, methylations, phosphorilations, and ubiquitylations are some of the panoply of modifications that histones can be subjected to Bannister and Kouzarides (2011), creating patterns that are associated to different transcriptional states. Unfortunately, these modifications are numerous and challenging to analyze in trees as standard chromatin immunoprecipitation (ChIP) procedures are hard to apply to woody cells and tissues (Li et al., 2014). Thus, very few studies regarding histone modifications in trees have been published to date.

Non-coding RNAs are also rather relevant as an epigenomic mechanism, as they are mobile within and between cells and are not lost during cell divisions, unlike the majority DNA methylations and histone modifications (Yakovlev et al., 2012). They are classified by their length: if they have more than 200 nucleotides, they are called as long ncRNAs (lncRNAs) whereas if they have less than 200 nucleotides, they are called short ncRNAs. These short ncRNAs include micro RNAs (miRNAs), short interfering RNAs (siRNAs), and Piwi-interacting RNAs (piRNAs), and they are involved in genomic imprinting, X-chromosome inactivation, repression of gene transcription, and transposon repression (Wei et al., 2016). Therefore, given the elevated content of TEs in conifer genomes, ncRNAs stand as a major epigenetic mark that should be studied in depth.

Lastly, DNA methylation positions itself as a highly relevant epigenetic mark in conifers due to their large content in TEs and its role in silencing them. In addition, there are numerous studies that demonstrate its participation in the control of gene expression in a large number of traits, in plants in general and in conifers in particular, which will be addressed later. DNA methylation is a covalent change which consists in the addition of a methyl group to a cytosine residue, in most cases to its fifth carbon position, and can be, in some cases, inherited by the offspring. In plants, DNA methylation happens not only in CpG sites but also in CHG and CHH contexts, where H can be A, T, or C (Law and Jacobsen, 2010). CpG sites are the most frequent mark in/near protein-coding gene sequences (Ashapkin et al., 2020) especially of those with constitutive expression. Although this methylation is known to affect gene expression, its role is still unclear, as it can lead to both silencing and increased expression (Zemach and Ziberman, 2010). Regarding CHG and CHH sites, different studies have highlighted that they are mainly present in repetitive sequences and TEs, both abundant in tree genomes, and their methylation usually reduces or even blocks DNA transcription, allowing the maintenance of the genome integrity (Zemach et al., 2010a; Bewick and Schmitz, 2017).

It has been proved that DNA methylation is responsible of phenotypic variation in plants (Martin et al., 2009; Mougninot et al., 2021; Zhao et al., 2021). Considering that different types of stress, such as abiotic stress, can trigger methylation, eventually leading to variation of gene expression, this field of research has become a useful tool for understanding how tree species deal with the environment variations. Several studies have highlighted that epigenetics seems to act like a memory process to respond to different environmental stresses, which is called “stress memory” (Kinoshita and Seki, 2014). Thus, a new question arises: could it be possible to inherit this stress-induced DNA methylations? Most of the time, epigenetic marks are erased to start over in the progeny, however it is possible for this reset to be incomplete, allowing some marks to pass on to the offspring. Moreover, it is known that some germline cells come from somatic cells in plants, so epigenetic marks that appear in those cells could be passed on to the offspring. Hence, epigenetic marks can be transmitted to the offspring by meiosis in some cases. Studies carried out with Arabidopsis thaliana L. (Cervera et al., 2002; Boyko et al., 2010; Lang-Mladek et al., 2010) or Linaria vulgaris Mill. (Cubas et al., 1999) evidenced epigenetic modifications related to environment variations can be inherited. However, the mechanisms involved in this process are not well understood to date.

All in all, epigenetic modifications could be a rapid way to achieve short-term adaptations, since epigenetic marks bring new phenotypes which might affect fitness and, thus, the underlying evolutionary processes. Particularly, it could be key for forest trees survival, as they are sessile long-lived organism with long generation times, so the possible beneficial effects of natural selection on them are limited and slow and will soon be surpassed by the effects of climate change.

Techniques for Epigenomic Analysis

Over the years, different techniques have been developed to study DNA methylation. Their first requirement is high-quality DNA, which highlights the importance of selecting an appropriated method for sample conservations during and after each collection, as well as an extraction method as efficient as possible for each type of samples. Once the DNA is obtained, there is a wide variety of methods than can be used to analyze its methylation patters (Table 1). The most remarkable ones are shown below.

TABLE 1

Summary of the advantages and weaknesses of the most relevant DNA methylation techniques and their suitability for species without reference genomes.

Technique	Advantages	Weaknesses	Suitable for species without reference genome
MeDIP-seq	Low initial amount of DNA. Genome-wide scale.	No single base resolution. It does not discriminate methylation contexts. Relative levels of methylation.	No
MRE-seq	Reduction of the complexity of genomes. Single-base resolution.	Low coverage of the genome. Relative levels of methylation.	No
MSAP-seq	Reduction of the complexity of genomes. Useful with large-size genomes.	Low coverage of the genome. Loss of information of CHH and CHG contexts.	Yes
WGBS	It uses entire genomes. Single-base resolution. All contexts are represented. Absolute levels of DNA methylation.	Not suitable for large-size genomes.	No
RRBS	Reduction of the complexity of genomes. Single-base resolution. All contexts are represented.	Low coverage of the genome.	Yes

At the beginning, methods based on antibodies affinity to methylcytosine residues were used. One of them is the well-known methyl-DNA immunoprecipitation (MeDIP) technique, which captures the methylated DNA fragments with an antibody against methylated cytosines (Thu et al., 2009). Then a new approach of this technique, called MeDIP-seq, was developed (Zhao et al., 2014). The novelty of this method consisted in a final sequencing step based in next-generation sequencing (NGS). Both techniques require a previous DNA sonication step to obtain random DNA fragments. MeDIP-seq is frequently used in combination with MRE-seq (methylation sensitive restriction enzyme sequencing) (Maunakea et al., 2010). MRE-seq consists of a DNA fragmentation using several restriction enzymes which are methylation sensitive followed by a final sequencing step. Although MeDIP-seq and MRE-seq can be used independently, their combination can increase the accuracy of DNA-methylation studies (Li et al., 2015).

Subsequently, an improvement based on NGS allowed to modify the methyl-sensitive amplified polymorphism (MSAP) technique, giving rise to MSAP-seq. Originally, the MSAP technique appeared as a modification of the amplified fragment length polymorphism method (AFLP). This technique was described by Reyna-López et al. (1997) in fungi and then Xiong et al. (1999) adapted it to work with plants. MSAP starts with a genome fragmentation using restriction enzymes with differential sensitivity to cytosine methylation. MspI and HpaII are one of the enzyme isoschizomer pairs used in MSAP. Both recognize the same target 5’-CCpGG-3’, but HpaII will not cleave if the inner cytosine residue is methylated. Hence, this pair is used to identify differences in methylation marks between DNA samples in CpG contexts mainly. Following the cleavage step, a PCR amplification is carried out with the aim to obtain an overrepresentation of previously obtained fragments. Then, amplicons are separated and visualized via gel electrophoresis. MSAP-seq replaces this last sept by a sequencing one. This method was first described by Chwialkowska et al. (2016) in barley and was later improved to work with large genome species (Chwialkowska et al., 2017). It allows studies in non-model species with lack of a reference genome. However, restriction enzymes seem to cut less efficiently in CHH and CHG contexts than in CpG islands, which may lead to mistakes in methylation calling, which in turn might imply a loss in methylated states. This could be a remarkable weakness of this technique considering that the larger genome size is, the greater amount of cytosine methylations there is in the other contexts.

Finally, new approaches based on sodium bisulfite treatments have been developed to study DNA methylation, since bisulfite sequencing was developed in 1992 (Frommer et al., 1992). Bisulfite methods are useful to study genome-wide methylation and provide high single-base resolution. These techniques convert unmethylated cytosines to uracils, as methylated cytosines are protected against the deamination effect produced by the bisulfite treatment. After that, uracils are turned into thymines in a PCR amplification. Then, these changes can be easily identified by sequencing. To date, there are two bisulfite techniques: whole genome bisulfite sequencing (WGBS) and reduced representation bisulfite sequencing (RRBS), which is subdivided in three protocols based on RAD-seq (Baird et al., 2008) and genotyping by sequencing (GBS; Elshire et al., 2011) techniques: epiRADseq (Schield et al., 2016), bsRADseq (Trucchi et al., 2016), and epiGBS (van Gurp et al., 2016).

WGBS works with the entire genome, while the RRBS method targets specific regions. In RRBS methylation sensitive restriction enzymes are used to fragment the whole genome, and then DNA fragments are filtered by size. The characteristics of WGBS make it suitable for studies with model organisms with short length genome size and for species with a reference genome (Paun et al., 2019). On the other hand, the genome fragmentation carried out with RRBS methods constitutes a useful tool to reduce the complexity of genomes (Davey et al., 2011), which allows to work with large-size genomes, with or without a reference genome.

Although this review is focused on DNA methylation, here is a brief description of other techniques used to study the remaining epigenetic marks: histone modifications and ncRNAs.

One of the most common techniques to study in vivo histone modifications is chromatin immunoprecipitation (ChIP). First, samples are treated with formaldehyde to cross-link the chromatin interactions (protein–protein and protein–DNA). Then, small fragments are obtained by sonication or nuclease digestion. Next, specific antibodies are used to retain the DNA-histone modification complex of interest. Finally, different methods can be carried out to determine these regions such as PCR, real-time PCR, Southern blot, and Western blot. The development of NGS techniques has led to the addition of a final sequencing step to the ChIP protocol, creating the technique known as ChIP-seq. However, the presence of polysaccharide and phenolic compounds in plants may difficult these experiments, that could turn them into a challenge in some cases (Muhammad et al., 2020).

The most common technique used to study ncRNA sequences is RNA-seq (Wang et al., 2009). Then, sequencing data can be combined with RT-qPCRs to determine expression levels, or with fluorescence in situ hybridization (FISH) experiments to localize the ncRNAs in a cellular context, depending on the aim of the study. An in silico approach consists in looking for previously described or computationally predicted ncRNAs in the genomic region under study. Several databases of ncRNAs sequences are available online, such as CANTATAdb 2.0[1] (Szcześniak et al., 2015) which has over 200,000 lncRNAs from 39 plant species from which 2 are trees. As usual, forest trees are poorly represented in this kind of databases, mainly due to the lack of reference genomes.

In Search of a Reference

Forest trees span a huge range of genome sizes, including species with less than 1 Gb, such as the Populus genus, to the conifer giga-genomes, whose size varies from ∼6.5 to ∼37 Gb, being ∼15 Gb the most common 1C-value (Ahuja and Neale, 2005). Compared to the ∼0.6 Gb average of most angiosperms (Dodsworth et al., 2015), conifer genomes are remarkably large and, therefore, constitute a challenging work scenario for genomic and epigenomic studies. In addition, a great number of experimental techniques and bioinformatic pipelines are designed work with the 3.3 Gb-long human genome, making it difficult to adapt them to genomes which are on average five times bigger. For instance, techniques aimed at sequencing complete genomes, such as WGBS, might not be suitable for giga-genomes as optimal depth and coverage values will be hard to achieve.

To complicate matters further, approximately 75% of conifer genomes consists of highly repetitive DNA, mainly represented by non-coding DNA, such as mobile elements, along with rDNA repeat units (Ahuja and Neale, 2005). In addition, gene duplications have created large gene families, which are not necessarily found in certain regions but dispersed across the genomes (Kinlaw and Neale, 1997).

On the other hand, luckily, polyploidy is not a common phenomenon in conifers, and it has only been reported in the Cupressaceae family (Ahuja, 2005). All in all, however, the global genetic and, thus, epigenetic landscapes of these tree species do not provide a research-friendly environment.

While many plant genomes are being sequenced and assembled every year due to the improvement of next generation sequencing techniques and their decrease in cost (Mackay et al., 2012), few conifer reference genomes are available. The first gymnosperm giga-genome ever published was that of Picea abies (L.) H. Karst. (Nystedt et al., 2013). Since then, six more conifer genomes have been sequenced: Pinus taeda L. (Neale et al., 2014), Picea glauca (Moench) Voss (Warren et al., 2015), Pinus lambertiana Douglas (Stevens et al., 2016), Pseudotsuga menziesii (Mirb.) Franco (Neale et al., 2017), Larix sibirica Ledeb. (Kuzmin et al., 2019), and Abies alba Mill. (Mosca et al., 2019). The N50 statistic is used as an indicator for genome assembly quality. It is defined as the length of the contig that reaches 50% of the total genome length, being the contigs sorted by length in descending order. Therefore, the higher this value is, the better the assembly is, as the median contig size will be longer. The total length of the previously mentioned genomes ranges from ∼12 to ∼32 Gb, with N50 scaffold sizes between ∼5 to ∼2,500 Kb, which shows that further refinement is still needed. Nevertheless, over 45,000 genes have been successfully annotated in each of the cited examples.

The recent upsurge in popularity and number of users of third generation sequencing technologies, also known as long-read sequencing, could provide a substantial boost to the improvement of the quality of the existing reference genomes. Long reads, alone or in combination with short reads, have a high potential for assembling and re-assembling genomes (Amarasinghe et al., 2020). The reads obtained from third generation sequencing methods, such as PacBio single-molecule real-time (SMRT) and Oxford Nanopore, constitute a less fragmented representation of large genomes, facilitating their assembly in longer scaffolds. For instance, the quality of the reference genome of barley (Hordeum vulgare L.), which is not a forest tree but whose genome shares some characteristics with the conifer ones, such as being highly repetitive and ∼5 Gb-long, was superior for long-read assemblies than for the one based on short reads, obtaining higher contig N50 values (Mascher et al., 2021).

On the other hand, sequencing of transcriptomes is more extensive in conifers, due to their notably reduced size compared to whole genomes, which reduces costs and increases the probability of capturing the majority of the sequences. The One Thousand Plant Transcriptomes Initiative database[2] collects over 1,300 transcriptomes of green plants, of which 76 belong to 72 different conifer species, coming from diverse tissue samples, including leaves, branches, and young shoots. Hardwood Genomics Project[3] and TreeGenes[4] are two other examples of databases specialized in trees that include transcriptomes as well as genomes, annotations and genetic markers.

However, the information coming merely from transcriptome sequencing is clearly insufficient to address the study of epigenetic mechanisms in conifers. Some of the most common epigenomic methodologies, such as identification of differentially methylated regions (DMRs) and differentially methylated promoters (DMPs), need the information of the sequences surrounding the coding regions. Then, the transcripts identified in transcriptomic experiments are useless without a reference to provide their genomic context.

Regarding complete methylomes, it has been described that CpG and CHG methylation is positively correlated with genome size (Ausin et al., 2016; Takuno et al., 2016), and high methylation levels are also known to be characteristic of genomes with a high content of transposable elements, as they constitute the main target of DNA methylation (Baduel and Colot, 2021). The lack of methylation maps of conifers is therefore not surprising, given the vast amount of data that needs to be analyzed to create them. Withal, the Norway spruce (P. abies) DNA methylation map was successfully obtained (Ausin et al., 2016). As previously mentioned, the genome of this species is fully sequenced, which highlights the importance of having reference genomes in order to be able to conduct this type of high throughput epigenetic studies.

Bioinformatic Tools

Bioinformatic pipelines for the analysis of DNA methylation data obtained after bisulfite treatments consist of four main steps: quality check and read trimming, read mapping, methylation calling, and identification of differentially methylated regions (Wreczycka et al., 2017).

Quality check and trimming of raw reads can be performed with well-known software, including FastQC (Andrews, 2010), and Trimmomatic (Bolger et al., 2014) and Cutadapt (Martin, 2011), respectively. Generally, these steps include removing adapter sequences, filtering out low-quality reads and trimming reads of low-quality bases. Specific bisulfite-sequencing tools are also available, such as BSeQC (Lin et al., 2013), which evaluates quality and automatically trims technical biases that may interfere in accurate methylation calling. For example, it can detect 5’ bisulfite conversion failure and remove the reads affected by this bias, which causes artificially high methylation rates at that end of the reads.

Bisulfite conversion rates should also be checked, especially if non-CG methylation levels are the focus of the study. Non-CG context, as previously mentioned, play an important role in silencing repetitive sequences and TEs, which are highly abundant in conifer genomes, so their study could be of great interest. In particular, CHH contexts can be quite problematic due to their low methylation levels (Takuno et al., 2016), which require a high bisulfite conversion rate to ensure that they are not underestimated. Thus, a spike-in control, such as unmethylated lambda phage DNA, could be useful to check that the bisulfite conversion is working at optimal levels.

Read mapping is a particularly challenging step. After bisulfite treatment, unmethylated cytosines are converted to uracils and then to thymines in the following PCR amplification, while methylated cytosines are not affected. This creates a need for special reference genomes, modified to take into consideration these nucleotide changes. In addition to the increased mapping searching sequence, asymmetric cytosine to thymine alignments and multiple CpG heterogeneous methylation must be taken into consideration by the mapping software (Xi and Li, 2009). Nonetheless, several bisulfite mapping tools have been developed, such as BSMAP (Xi and Li, 2009) and Bismark (Krueger and Andrews, 2011). Recently, Grehl et al. (2020) compared the efficiency and impact of the mapping results on the identification of DMRs of eight mappers and concluded that these two are the best options, as they showed the highest precision with the shortest run time and lowest memory requirements, respectively.

Both tools also provide methylation level files showing methylation calling status for each cytosine position. The coverage files produced by the Bismark methylation extractor function can be directly used as inputs for DMR-finding programs, such as the popular edgeR package (Robinson et al., 2010), or can easily be edited to fit other tools’ input requirements, such as DSS (Park and Wu, 2016) and swDMR (Wang et al., 2015).

In addition to these tools that allow users to perform each step of the analysis separately, others combine all the procedures in one pipeline. Bicycle (Graña et al., 2017) and SMAP (Gao et al., 2015) are some of the recommended all-in-one software for WGBS and RRBS data analysis, respectively (Rauluseviciute et al., 2019).

Despite these apparently wide variety of tool choices, there are several issues to take into consideration when it comes to analyzing high-throughput bisulfite sequencing data. First, the need of a good quality reference genome to perform the read mapping step, which is not a frequent situation, as previously explained. Besides, some tools, like swDMR, require an annotation file as input, whose quality relies on the genome quality itself. For instance, if a DMR falls in a region that is not described, likely due to genome fragmentation caused by large number of scaffolds or small scaffold sizes, this DMR will provide little information regarding its biological meaning. In addition, the amount of Cs that would need to be analyzed will in all probability be quite high in comparison to other non-conifer species. Therefore, the computational requirements in epigenomic experiments on forest tree species are likely to be very demanding, particularly in terms of RAM memory usage, storage space and run time.

Researchers that choose non-bisulfite-based techniques are likely to face the above-mentioned problems derived from genome size and lack of references, as they also include a next-generation sequencing step. Nonetheless, as in the case of bilsulfite-based methods, there are many bioinformatic tools specifically designed to work with this kind of data. The MeDUSA pipeline (Wilson et al., 2012) and Batman software (Down et al., 2008) are suitable for MeDIP-seq data analysis, while the MSEQER pipeline (Chwialkowska et al., 2017) is designed to work with MSAP-seq data.

Previous Studies

As it is said above, studies of epigenetic marks in plants are mainly focused on DNA methylation (Kurdyukov and Bullock, 2016). Nevertheless, the role of histone acetylations in root growth and development (Ma et al., 2016) and in photosynthesis (Li et al., 2017) has been described in different species of poplar, as well as that of miRNAs in growth and phosphorus nutrition (Schönberger et al., 2016). It is worth noting that this kind of studies have been carried out on poplars because they are considered a forest tree model organism. As mentioned above, poplar has a small genome (∼0.5 Gb) with a low proportion of repeated sequences. This genome has been fully sequenced and annotated, so its genomic resources are of high quality. Except in this particular case, histone modifications and ncRNAs remain mostly undescribed in conifer species and represent a major challenge to advance in the knowledge of epigenomics in trees.

The first plant methylome published was that of A. thaliana L. (Cokus et al., 2008), and it was developed using a large-scale bisulfite sequencing experiment, while the first conifer one, of P. abies, did not come until 2016 (Ausin et al., 2016).

Over the years, some studies have focused on changes in DNA methylation during normal plant development (e.g., Bitonti et al., 2002; Zemach et al., 2010b; Perrin et al., 2020) and heteroblastic development (Hasbún et al., 2016). Others are about phenotypic plasticity and stress memory (e.g., Herrera and Bazaga, 2013; Le Gac et al., 2018). What it is common to the vast majority of them is the use of model species, like A. thaliana (e.g., Korotko et al., 2021; Markus et al., 2021), and important commercial species such as crops (e.g., Feng et al., 2010; Wang et al., 2010; Karan et al., 2012; Bhanu et al., 2020; Merce et al., 2020) or fruit trees (e.g., Avramidou et al., 2015; Badad et al., 2021; Rothkegel et al., 2021).

The knowledge about forest tree epigenetics is scarce. Nevertheless, some studies related with epigenetic memory during embryogenesis have been developed on forest trees of species such as P. abies (Kvaalen and Johnsen, 2008; Yakovlev et al., 2014, 2016), P. glauca (Webber et al., 2005), or P. sylvestris (Dormling and Johnsen, 1992), and more studies are coming (e.g., Lafon-Placette et al., 2018; Vanden Broeck et al., 2018).

Furthermore, studies have generally focused on the involvement of DNA methylation in abiotic stresses response and adaptation in forest trees and other plants (Table 2). Modifications in response to water have been found in the methylome of Populus species (Gourcilleau et al., 2010; Raj et al., 2011), as well as variations in the methylation pattern of Pinus radiata D. Don. in response to heat stress (Lamelas et al., 2020). Again, it is remarkable the huge amount of poplar (Populus spp.) studies which have been carried out over the years, thanks to the availability of a quality reference genome and its wide phenotypic variation.

TABLE 2

Brief description of some DNA methylation studies carried out with the main objective of studying plant species response to several abiotic stresses.

Plant species	Technique	Objective	Reference genome	References
Bruguiera gymnorhiza (L.) Lamk.	WGBS	Salt stress response	No	Miryeganeh et al., 2021
Ilex aquifolium L.	MSAP	Epigenetic variability	No	Herrera and Bazaga, 2013
Medicago sativa L. and Medicago spp	MSAP	Salt stress response	Yes	Al-Lawati et al., 2016
Morus alba L.	WGBS	Drought stress response	Yes	Li et al., 2020
Quercus ilex L.	MSAP	Drought stress response	No	Rico et al., 2014
Quercus lobata Née	RRBS	Environment response	Yes	Gugger et al., 2016
Eucalyptus globulus Labill.	MeDIP-seq, MRE-seq	Heteroblastic development	Yes	Hasbún et al., 2016
Populus simonii Carr.	MSAP	Various abiotic stresses	Yes	Song et al., 2016
Populus trichocarpa Hook.	WGBS	Drought stress response	Yes	Liang et al., 2014
Populus euphratica Oliv.	WGBS	Salt stress response	No	Su et al., 2018
Populus tremula L. x Populus alba L.	RNAi, HPLC, WGBS	Drought stress response	No	Sow et al., 2021
Pinus pinea L.	MSAP	Epigenetic variability	No	Sáez-Laguna et al., 2014

Studies in natural populations are of particular importance to understand how epigenetic modifications can contribute to adaptation to the new climatic conditions caused by climate change. One of these studies was carried out with P. sylvestris (Alakärppä et al., 2018), which is a widely distributed conifer, likely to show variations in its epigenetic marks as it can be found under quite different environmental scenarios, and a link between local adaptation and environment condition, and modifications in DNA methylation has been suggested. Another work with samples of natural populations of Populus simonii Carr. (Ci et al., 2016) was conducted to analyze population epigenetic structure, which showed a substructure population associated to different methylation patterns.

In the literature, several reviews about epigenetics in plants and forest trees can be found (Paun et al., 2019; Amaral et al., 2020; Ashapkin et al., 2020; Zhao et al., 2021). Most of them are focused on the epigenetic response to specific abiotic stresses, such as drought, salt, or heavy metals, and they are useful as an introduction to this challenging field. On top of that, this review should prove valuable to researchers who wish to design an experiment to study methylation, as it covers the whole process: from choosing an appropriate technique to performing the bioinformatic analysis of the data.

To sum up, epigenomics open a research field to identify some of the processes related with phenotypic plasticity and its study will bring us a better understanding of adaptive forest tree response to a changing environment. Fortunately, substantial efforts are being made to improve DNA methylation techniques to make them suitable and cost-effective for non-model species, specifically with trees or species with large genomes (e.g., Chwialkowska et al., 2017; Pereira et al., 2020; Werner et al., 2020).

Conclusion

It is expected that the growth rates of forest tree species become on average slow, while the upcoming climate change trends become faster than the rate of allele frequencies shift of forest trees. This emphasizes the putative role of epigenetic modifications as a fast adaptive mechanism to cope with current warmer and drier climate scenario.

Phenotypic plasticity is the principal strategy to compensate for trees’ lack of ability to move when it comes to dealing with rapid environmental variations, such as ongoing climate change. It is based on modifications of genetic regulation which depend on epigenetic mechanisms such as histone modifications, DNA methylation and non-coding RNAs. Despite of the importance of these processes, epigenomic studies in forest tree species are challenging and scarce, mainly due to their complex genomes and lack of references.

Throughout this review, we have gathered enough information to give a general overview about the key points for methylation experiments design when dealing with forest trees. We have also highlighted opportunities for future advances in this framework and the relevance of dedicating efforts to improve the techniques for epigenomic analysis and the bioinformatic tools to make them suitable for large-size genomes.

Author Contributions

IG-G and BM-C wrote the manuscript. All authors discussed the results and contributed to the final manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.