Unique features of the m⁶A methylome and its response to drought stress in sea buckthorn (Hippophae rhamnoides Linn.)

ABSTRACT

In plants, recent studies have revealed that N6-methyladenosine (m⁶A) methylation of mRNA has potential regulatory functions of this mRNA modification in many biological processes. m⁶A methyltransferase, m⁶A demethylase and m⁶A-binding proteins can cause differential phenotypes, indicating that m⁶A may have critical roles in the plant. In this study, we depicted the m⁶A map of sea buckthorn (Hippophae rhamnoides Linn.) transcriptome. Similar to A. thaliana, m⁶A sites of sea buckthorn transcriptome is significantly enriched around the stop codon and within 3ʹ-untranslated regions (3ʹUTR). Gene ontology analysis shows that the m⁶A modification genes are associated with metabolic biosynthesis. In addition, we identified 13,287 different m⁶A peaks (DMPs) between leaf under drought (TR) and control (CK) treatment. It reveals that m⁶A has a high level of conservation and has a positive correlation with mRNA abundance in plants. GO and KEGG enrichment results showed that DMP modification DEGs in TR were particularly associated with ABA biosynthesis. Interestingly, our results showed three m⁶A demethylase (HrALKBH10B, HrALKBH10C and HrALKBH10D) genes were significantly increased following drought stress, which indicated that it may contributed the decreased m⁶A levels. This exhaustive m⁶A map provides a basis and resource for the further functional study of mRNA m⁶A modification in abiotic stress.

KEYWORDS:

• Sea buckthorn
• N6-methyladenosine
• drought stress
• m⁶A modification
• m⁶A demethylase

Introduction

As a rising novel epigenetic modification, N6-methyladenosine (m⁶A) of mRNA has been widely studied in mammals and plants [1–3]. As one component of a methyltransferase complex, N6-adenosine methyltransferase like 3 (METTL3) were considered as a methyltransferase of m⁶A, which identified in Saccharomyces cerevisiae, Drosophila melanogaster, Homo sapiens and Arabidopsis thaliana [4–7]. Expect, methyltransferase complex also contained methyltransferase-like 14 (METTL14) and Wilms tumour 1-associated protein (WTAP), which combined with METTL3 to catalyse m⁶A RNA methylation in mammals [8,9]. All these N6-adenosine methyltransferases can affect organism growth and development, including embryo, shoot, root, leaf and floral [2]. On the contrary, past researches revealed that FTO and ALKBH5 were confirmed as two m⁶A RNA demethylases in mammalian development [10]. Then, thirteen m⁶A demethylases (ALKBH family) are also studied in Arabidopsis [11,12]. For example, AtALKBH9B can remove m⁶A modification of single-stranded RNA in vitro. In addition, several RNA-binding proteins (RBPs), which play a role in recognizing m⁶A marks, have been identified in animals and Arabidopsis. All these findings reveal the reversible RNA modification can regulate functions of its target mRNA.

In recent researches, transcriptome-wide m⁶A methylation have been identified in Arabidopsis, rice and maize [3,13–15]. For instance, the m⁶A content of mRNA varies across different tissues in Arabidopsis. As an N6-adenosine methyltransferase, MTA is important for shoot, root, leaf and floral development of Arabidopsis [2]. An m⁶A eraser protein (ALKBH10B) can influence floral transition by regulating transcript levels of SPL3, SPL9 and FT [11]. Then, m⁶A reader proteins (ECT2, ECT3 and ECT4) can also regulate plant organogenesis [16]. With the development of m⁶A-seq, m⁶A modification was found in many species in response to environmental stress [17–21]. For example, the m⁶A can enhance the salt tolerance of Arabidopsis by regulating the stability of salt-tolerant transcripts [22]. Similarly, m⁶A-seq reveal that N6-methyladenosine can also regulate the salt tolerance of sweet sorghum [23]. In response to heat stress, the reader ECT2 can relocate to stress granules in the cell of Arabidopsis [24]. As a grave global ecological problem, drought stress has a great impact on plant development. In Arabidopsis, levels of m6A methyltransferase, m6A demethylase and m6A-binding proteins were not significantly modulated by drought stresses. However, drought stresses can up-regulate the expression level of m⁶A demethylase in maize [19]. It implied the different roles of mA in different species in the response to drought stress.

As a deciduous shrub, sea buckthorn (Hippophae rhamnoides Linn.) are widely grown and cultivated in Asia and Europe [25–27]. Now, sea buckthorn has been widely used to protect the environment because of its light-tolerant and drought-tolerant properties. The physiological responses of sea buckthorn seedlings under drought stresses have been studies on previous researches, which showed that epigenetic modification and post-translational modification may play important roles in sea buckthorn to drought stress [28]. In this study, to further explore the roles of mA in plants response to drought stress, we show transcriptome wide mA profiling in sea buckthorn seedlings and leaf under CK and TR condition. We found that mA is significantly enriched around the stop codon region and within 3ʹUTRs in sea buckthorn, which is similar to those in A. thaliana. m⁶A deposition and mRNA levels have a positive correlation, which indicates a regulatory role of m⁶A in plant gene expression. In general, as compared with the normal sample, we detected 13,287 DMPs within mRNAs. Intriguingly, the DEGs harbouring DMPs were involved in many important biological pathways associated with photosynthesis.

Results

m⁶A methylation is abundant in sea buckthorn mRNA

We used the m⁶A-targeted antibody and performed the MeRIP-seq to obtain the transcriptome-wide m6A map of the whole H. rhamnoides plants. We obtained almost 30 million reads per library (input and IP libraries). After quality control of raw data, 90% clean reads were uniquely mapped into the sea buckthorn reference genome. exomePeak2 software package was used to identify the m⁶A peaks between the IP and control samples. We totally detected 18,102 m⁶A enriched peaks in sea buckthorn leaf samples, which have an average length of 331 bp. The minimum length of peaks was 50 bp. Based on these results, we found that the H. rhamnoides transcriptome contains 0.7 m⁶A peaks per 1,000 nucleotides, which are similar to those obtained in mammals. And, most genes contain one or two m⁶A peaks (Figure 1(a)).

Figure 1. Overview of m⁶A methylome in sea buckthorn. (a) Number of genes with different m6A peak numbers. (b) Accumulation of m⁶A-IP and input reads of transcripts. Each transcript contained three parts: 5ʹ UTRs, CDS and 3ʹ UTRs. (c) The m⁶A peak distribution in different genic regions. (d) The m6A peak distribution along a metagene. Enrichment scores are calculated as (n/N)/p; n, number of peaks belonging to each category; N, number of total peaks; p, proportion of each category within the genome by length. (e) Distribution of m6A peaks across sea buckthorn chromosomes. (f) The RRACH and DRACH conserved sequence motifs of m6A peak regions

We investigated the m⁶A peaks distribution according to gene annotations. Among them, most reads of m⁶A peaks were localized around 3ʹUTRs (Figure 1(b)). This distribution tread was similar to those in Arabidopsis. The result of the locations of m⁶A on transcripts showed that m⁶A peaks are abundant around the stop codon (43.87%), followed by the coding regions (20.99%), start codon (19.09%) and then 3ʹ UTRs (13.72%), which consistent with the distribution of reads (Figure 1(c)). Then, we also observed that m⁶A peaks is significantly enriched around the stop codon (enrichment score > 25) by segment normalization of each gene portion (Figure 1(d)). Whole-genome density of m⁶A modification peaks revealed that m⁶A peaks have different distribution modes across chromosomes (Figure 1(e)). Similarly, m⁶A peaks have a higher density at the telomeric ends compared with other regions. Further analysis showed that the average density of chromosome 9 was the highest whereas chromosome 8 had the lowest average density (Figure 1€). These distribution modes among chromosomes indicated that m⁶A modification may have a relationship with chromatin state. Based on the analysis of m⁶A consensus sequence, we found that 15,387 (94.39%) m⁶A peaks have at least one RRACH motif in the CK group and 10,873 (94.15%) m⁶A peaks have at least one RRACH motif in the TR group (Fig. 1F). The same sequence motif, which is observed in eukaryotic mRNA, may be essential for the formation of m⁶A modification.

Characterization of mRNAs containing m⁶A modification

Past researches revealed that m⁶A is critical for plant development. To further analyse the functional roles of m⁶A modification in H. rhamnoides, we performed gene ontology (GO) enrichment analysis of genes containing m⁶A modification using the enrichment tool at HRIA. We found that most of these genes are involved in metabolic biosynthesis (Figure 2(a)). To determine whether the m⁶A distribution patterns of genes have a relationship with its GO categories, we divided genes into three subgroups (PeakStart: m⁶A peaks around start codon, PeakStop: m⁶A peaks around stop codon, PeakBoth: m⁶A peaks around both start and stop codons) based on the distribution of mA peaks. GO enrichment analysis of each subgroup showed that mA-containing genes in PeakStart and PeakBoth didn’t have enrichment terms. But, mA-containing genes of PeakStop exhibit high enrichment of metabolic process, signal transduction and RNA binding (Figure 2(b)). Certainly, we identified a series of metabolic-related genes carrying mA peaks. For instance, Chalcone Synthase (HrCHS), Chalcone isomerase (HrCHI) and flavanone 3-hydroxylase (HrF3H) contained mA peaks around stop codon, which plays important roles in flavonoid biosynthesis (Figure 2(c)). Phytoene synthase (HrPSY), which also contained mA peaks at stop codon, play an important role in phytoene biosynthesis. The results suggest that mA mRNA methylation may regulate metabolites biosynthesis in plants.

Figure 2. Enrichment analysis of genes with m⁶A peaks. (a) GO enrichment of genes with m6A peaks. MF: Molecular Function, CC: Cellular Component, BP: Biological Process. (b) KEGG enrichment of genes with m⁶A peaks. (c) Examples of genes with m6A peaks at stop codon

m⁶A methylation level is decreased in sea buckthorn under drought stress

To determine the change of m⁶A methylation in sea buckthorn response to drought stress, we profiled the m⁶A methylomes of the sea buckthorn leaf under CK and TR treatments. We performed the comparison of m⁶A profiles of CK and TR and identified 16,302 and 11,549 m⁶A differential methylated peaks (DMPs) in CK and TR, revealed a significant decrease in TR compared to CK. The Venn diagram identifies 7,496 common m6A peaks both in CK and TR (Figure 3(a)). Different analysis showed that 13,212 m⁶A methylation peaks were considered as DMP, including 6609 up-regulated DMS and 6603 down-regulated sites (Figure 3(b)). Then, we detected 9,596 transcripts that have differentially m⁶A methylation peaks (DMPs) in TR compared to CK. In detail, 6609 up-regulated DMPs contained 2174 up-regulated genes, 2932 down-regulated genes and 96 unaltered genes. 6603 up-regulated DMP contained 2106 up-regulated genes, 3077 down-regulated genes and 100 unaltered genes. More than 40% DMPs are located in stop codon regions (Figure 3(c)). Further analysis of the whole transcriptome of CK and TR, we found that m⁶A modification can certainly affect genes expression levels (Figure 3(d)). Therefore, we performed gene ontology (GO) enrichment analysis of genes containing DMPs to explore functional characteristics of DMPs in the context of genic location. The GO enrichment analysis of up-regulated DMP contained genes showed that the biological processes of blue light signalling pathway (GO:0009785), electron transport chain (GO:0022900) and photosynthesis, light reaction (GO:0019684) were significantly enriched (Figure 3(c)). The GO enrichment analysis of up-regulated DMS contained genes showed that the cellular component of the cell wall (GO:0005618) was significantly enriched. But, the GO enrichment analysis of down-regulated DMPs contained genes showed that the biological processes of protein ubiquitination (GO:0016567), nuclear-transcribed mRNA catabolic process (GO:0000956), cell wall modification (GO:0042545) were significantly enriched (Figure 3(e)). Otherwise, KEGG analysis reveal that up- and down-regulated DMPs contained genes were enriched metabolic pathway (Figure 3(f)).

Figure 3. Overview of m⁶A methylation of sea buckthorn leaf under CK (control) and TR (drought) treatment. (a) Numbers of group-common and group-specific m⁶A peaks between CK and TR. (b) volcano analysis of differentially m6A methylation peaks. (c) Distribution of differentially m⁶A methylation peaks across the mRNAs. (d) Cumulative distribution of mRNA expression changes between CK and TR for m⁶A-modified genes and non-target genes. (e) Gene ontology enrichment terms for the genes with DMPs between CK and TR. (f) KEGG enrichment pathways of the genes with DMPs between CK and TR

The m⁶A pathway modulates genes responses to drought stress

To explore the roles of m⁶A methylation in sea buckthorn response to drought stress, we profiled the modification between normal and drought stress treatment leaf samples. We performed co-differential analysis to further associate m⁶A modification and gene expression. Finally, 541 DMGs were identified as significant co-differential genes with DMPs, including 228 up-regulated genes and 313 down-regulated genes (Figure 4), Supplementary Tables 1 and 2). Previous studies also showed that m⁶A modification has a great change during plant development and under abiotic stress. Furthermore, GO and KEGG analysis research revealed that up-regulated genes were significantly enriched in abscisic acid metabolic process, response to ethylene, fatty acid metabolic process, negative regulation of abscisic acid-activated signalling pathway (Supplementary Table 3). For example, CRTISO (Sph_LG3G002144) plays an important role in response to drought stress by affecting ABA biosynthesis in Arabidopsis. In addition, CYP707A2 (Sph_LG2G000371) encode ABA 8′-hydroxylases have hyper- and hypo- m⁶A peaks in TR samples, which decrease the ABA levels in Arabidopsis (Figure 4). The methylation levels of m⁶A peaks of CRTISO mRNA have a two-fold increment in TR samples compared to CK samples. However, the transcripts level of CRTISO in TR samples was higher than in CK samples. RMA1H1 (Sph_LG9G000969) was hypermethylated in TR, which can be induced by dehydration and conferred strongly enhanced tolerance to drought stress in Arabidopsis plants. Together, all these results indicate that m⁶A methylation can regulate the expression level of ABA-related genes to candidate the response to drought stress in sea buckthorn.

Figure 4. Analysis of DMP-related gene and DEGs between sea buckthorn leaf under CK (control) and TR (drought) treatment. DEGs: differentially expression genes, DMP: differentially m⁶A methylation peaks. (a) Venn map of DMP related gene and DEGs. (b) Gene ontology enrichment terms of the DEGs with DMPs. (c) Examples of DEGs with DMPs at stop codon. (d) Expression level of DEGs with DMPs at stop codon

The m⁶A pathway modulates proteins responses to drought stress

In addition, we determined whether m⁶A affects the protein level in the drought response process by characterizing the m⁶A methylome and proteome in CK and TR groups. First, we selected 9,596 genes with DMPs, which didn’t have change at transcripts level during drought stress. Next, we focused on whether the m⁶A level of selected 9,596 genes is associated with its protein levels (Supplementary Table 4). Using the proteomic data that MALDI‐TOF‐TOF to deeply profile the leaf proteome under control and drought stress, we found that the m⁶A level and protein expression level of some genes have a correlation between TR and CK (Supplementary Table 5). For example, we observed a decrease in m⁶A methylation in the 3′ UTR and 5′ UTR of the Photosystem II Subunit (PsbS), Oxygen-evolving enhancer protein (PSBO), ribulose bisphosphate carboxylase/oxygenase activase (RCA1) are which correlated with a decreased protein level in TR, but no significant change at mRNA levels. Contrarily, aminomethyltransferase (AMT) has increased m⁶A methylation level in its mRNA and increased protein level in TR compared with CK. All these results suggested m⁶A modification may regulated protein levels of key genes to response to drought stress in plants.

ALKBH10B may contribute to the decreased m⁶A methylation level in TR samples

Genome-wide m⁶A methyltransferases, demethylases and binding proteins in sea buckthorn were identified and constructed a phylogenetic tree using ML methods (Figure 5(a)). And, chromosome distribution analysis showed these genes were widely distributed in 12 chromosomes (Figure 5(b)). Based on the RNA-seq data, we quantify the expression abundance of m⁶A related genes in the leaves of sea buckthorn under control, drought stress and rehydration treatment. Interestingly, three ALKBH10B homologs (HrALKBH10B, HrALKBH10C and HrALKBH10D) were significantly upregulated by drought stress and downregulated by rehydration (Figure 5(c)). Meanwhile, the m⁶A abundance in samples under drought stress was significantly lower than that in samples under control treatment, suggesting global m⁶A hypomethylation induced by drought stress. According to the expression level of m⁶A demethylase genes (Figure 5(d)), we considered that the m⁶A demethylase genes most likely induced decreased m⁶A methylation during drought stress in sea buckthorn.

Figure 5. Expression analyses of sea buckthorn m⁶A methyltransferases, demethylases and binding proteins. (a) Phylogenetic analysis of sea buckthorn m⁶A RNA methyltransferases (left), demethylases (middle) and binding proteins (right). (b) Chromosome distribution of m⁶A methyltransferases, demethylases and binding proteins genes. (c) Expression level of sea buckthorn m⁶A methyltransferases, demethylases and binding proteins genes in leaf under CK, TR and RE. Arabidopsis genes (AtMTA, AtFIP3, AtVIR, AtHAKAI, AtALKBH9A, AtALKBH9B, AtALKBH9C, AtALKBH10A, AtALKBH10B, AtALKBH6, AtECT2, AtECT3, AtECT4, AtECT5) were used as a comparison. (d) Relative mRNA levels of sea buckthorn m⁶A methyltransferases, demethylases and binding proteins genes in leaf under CK and TR

Discussion

Recently, several research studies have illustrated the changes and roles of m⁶A methylation in the eukaryotic, including mice, human, Arabidopsis, rice, maize and tomato [3,14,19,29–33]. Since the finding of methyltransferase and demethylase proteins, m⁶A was considered as a dynamic and reversible RNA modification [34,35]. For example, m⁶A methyltransferases play an important role in plant development and regulate shoot stem cell fate of Arabidopsis [16]. It is been shown that m⁶A modification can dynamically regulate cellular responses to heat shock, ultraviolet light and oxidative stress in mammal [36–38]. However, the study about m⁶A modification in plant under drought stress was lacking. In this research, we used the sea buckthorn to present the dynamic m⁶A modification in response to drought stress. We found that m⁶A is of high abundance in leaf tissue of sea buckthorn. This resource indicated that mA modification may play a functional role in sea buckthorn leaf.

In these sea buckthorn leaf mA-seq data, we discovered that the features of the mA distribution in sea buckthorn are similar to those in Arabidopsis 3. m⁶A in sea buckthorn transcriptome is significantly enriched at stop codon. These results indicate that distribution features of the m⁶A were conserved in plants. In mammals, the RRACH motif is the most common motif of m⁶A methylation peaks, which has also been found in Arabidopsis and maize [3,15,19]. A recent study about motif analysis of all m⁶A modification peaks of sea buckthorn showed that most peaks have high enrichment in canonical RRACH motif and DRACH motif. These are common between sea buckthorn and Arabidopsis likely indicated that the conservative of these two m⁶A methylation motifs between two species, which can be clearly clarified by further functional analysis of m⁶A-related proteins. m⁶A modification normally is enriched stop codons and 3′ UTR to impact mRNA expression level. According to the distribution of m⁶A modification peaks, we found some genes with m⁶A peaks around their stop codon were involved in the metabolic process, including flavonoid and porphyrin and chloroghyII.

A recent study revealed that m⁶A modification may candidate plant response to abiotic stress. Here, we found a great change in m⁶A modification of sea buckthorn leaf under drought stress, indicating that m⁶A may be considered as an important role in drought responses. Importantly, transcripts with m⁶A modification have a great change compared with transcripts with non-mA modification, which showed that m⁶A modification can regulate transcripts level during drought stress. This result has also been found in two maize strains (B73 and Han21) [19]. Furthermore, enrichment analysis revealed that DEGs with m⁶A methylation was associated with some important biological processes, implying the role of m⁶A modification in regulating expression of gene-related drought response. As an important gene involved in ABA biosynthesis pathway, a homolog of CRTISO (Sph_LG3G002144) were hypermethylated in leaf under drought stress, which plays an important role in Arabidopsis under drought stress [39]. Compared with leaf under control, CYP707A2 (Sph_LG2G000371) have hyper- and hypo-m⁶A peaks in TR samples, which can decrease the ABA levels in Arabidopsis [40]. All these results indicated that m⁶A may contribute to drought response by ABA biosynthesis in sea buckthorn leaf. In addition, RMA1H1 (Sph_LG9G000969) were increased and gained m⁶A peaks in TR. Overexpression RMA1H1 conferred strongly enhanced tolerance to drought stress in transgenic Arabidopsis plants [41]. So, m⁶A modification may regulate the expression of key genes to enhance the drought stress tolerance in sea buckthorn.

Drought response is governed by genetic and epigenetic regulation in the plant. In this study, we firstly contrasted a transcriptome-wide m⁶A modification map of sea buckthorn and showed the role of m⁶A modification in response to drought stress. These results provide resources for further uncovering m⁶A modification functions in abiotic stress response and plant metabolism.

Materials and methods

Plant material

Plant growth and analyses generally followed the procedures described in past research [42]. Two-year-old sea buckthorn plants were cultured in a greenhouse with controlled light, temperature and humidity conditions (16 h dark, 8 h light; 25°C under day, 20°C under dark, 60–70% relative humidity) for three months. Then, sea buckthorns plants were treated with water deficit conditions for 10d (25–30% field capacity). The control group was treated with 75–85% field capacity by watering (CK). The leaf, stem and root tissue were collected and mixed as seedlings sample. And, leaf under control (CK) and water deficit condition (TR) were collected and immediately frozen in liquid nitrogen.

RNA extraction and library construction

We extracted total RNA from sea buckthorn seedlings sample, leaf under control condition (CK), leaf under water deficit condition (TR) using Trizol (Invitrogen) following a published protocol. Then, we tested the quality of extracted total RNA by Nanodrop and gel electrophoresis. Polyadenylated RNA was extracted and fragmented to 200nt using RNA Fragmentation Reagents (Ambion). Fragmented RNA and m⁶A antibody were incubated in IP buffer for 3 h at 4°C. We eluted the mixture with elution buffer and treated the eluted RNA with RNasin following the manufacturer’s instructions. According to published protocol, we constructed immunoprecipitated RNA and input RNA library using mRNA Sample Prep Kit (Illumina) for each sample. Thus, we performed m⁶A-seq and input RNA-seq with NovaSeq 6000 by paired-end strategy.

Analysis of m⁶A-seq

According to the past research [20,33,43], we analysed sequence data using bioinformatic software. Firstly, we used FastQC to remove the adapter and low-quality data (Q < 25) of raw data from IP RNA-seq and input RNA-seq. Using HISAT2 software [44], we mapped the clean data from the input and IP-sequenced samples to the sea buckthorn reference genome (http://hipp.shengxin.ren/). Based on the background information (input RNA-seq data), we used R package (exomePeak2) to identify m⁶A peaks (P-value < 1e-5). Following past researches [15,20], we estimated the distribution of m⁶A in the different parts of the gene and transcripts. The conserved sequences of identified m⁶A peaks were found by MEME Suite (MEME-chip) and HOMER [45]. The de novo motifs were uncovered by MEME-chip software with setting the motif length (3–8bp and 8–15bp). Then we used Homer software to establish the enrichment significance of RRACH and DRACH in sequences of m⁶A peaks. We used methods in past research to study the distribution pattern of the m⁶A peaks [20,23]. PeakAnnotator (version 2.0) was performed to annotate m⁶A peaks to the sea buckthorn annotation file [46]. The m⁶A site differential algorithm was used for the identification of differentially methylated peaks (DMPs) between the CK and TR sample with a criterion of P-value < 0.05 and fold change ≥ 2.

Analysis of RNA-seq data

Previous studies reveal that the input sequencing reads in the m⁶A-seq can be used for RNA-seq analysis. Briefly, tophat2 and cufflinks were used for the reads mapping and assembly of each sample. We calculated the fragments per kilobase of exon per million mapped fragments (FPKM) as Gene expression level using Cuffdiff [47]. The FPKM were used to identify differentially expressed genes (DEGs) with a cut-off criterion (log2FC ≥ 1 and P-value < 0.05). The P values of each DEG were also adjusted to control the false discovery rate (FDR) by Benjamini and Hochberg’s approach. The Gene Ontology (GO) and KEGG pathway enrichment analysis of m⁶A modified genes and DEGs with DMPs was performed based on the tool in HRIA. The significance values of GO and KEGG pathways with P-value < 0.05 were considered significantly enrichment terms.

Identification, phylogenetic and expression analysis of the m⁶A methyltransferases, demethylases and binding proteins genes in sea buckthorn

Based on the sequence of the mA methyltransferases, demethylases and binding proteins, we identified their homologs in sea buckthorn by BLAST homology searches. The reference genome sequence and annotations of sea buckthorn were downloaded from the Hippophae rhamnoides Information Archive (http://hipp.shengxin.ren). MEGA 7.0 were used for constructing a maximum likelihood phylogenetic tree of m6A methyltransferases, demethylases and binding proteins with 1000 bootstrap replicates. The expression level of all these genes in CK, TR and RE groups was calculated.

m⁶A-IP-qPCR

According to past research [33], we performed m⁶A-IP-qPCR for selected genes. Briefly, RNA fragmentation buffer (10 mM Tris-HCl, pH 7.0 and 10 mM ZnCl2) was used to fragment intact mRNAs into ~ 300 nucleotide-long fragments by incubation at 95°C for 25 s. We purified fragmented mRNAs using standard ethanol precipitation and resuspended with 250 μL of DEPC-treated water. Then, 50 μL of fragmented mRNAs and 2.5 μg of anti-m6A polyclonal antibody (Synaptic Systems) were incubated at 4°C for 2 h in IP buffer. 5 μL of fragmented mRNAs was taken as the input control. The mixture and 10 μL of Dynabeads protein-A (Invitrogen) were subsequently incubated at 4°C for 2.5 h. After washing with high-salt buffer, the m6A-containing fragments and 6.5 mM N6-methyladenosine were incubated in 100 μL of IP buffer at 4°C for 2 h. Then, both the m6A-containing mRNAs and the input mRNAs were submitted to quantitative RT-PCR. The primers were listed in Supplementary Table 6. The value for the immunoprecipitated sample was normalized against that for the input. Relative m6A enrichment in a specific region of a transcript was calculated using the 2^−ΔΔCt method. The experiment was performed with three biological replicates.

Measurement of gene expression levels using real-time qPCR

We used the plant RNA extraction kit (Tiangen, DP441) to extract total RNAs from the sea buckthorn leaves. The extracted RNAs were used for reverse transcription reaction by the PrimeScript™ RT reagent Kit (Takara, RR036A). The synthesized cDNAs were then employed as templates for PCR amplification using the TB Green Fast qPCR Mix (RR430S) and StepOne Plus Real-Time PCR System (Applied Biosystems) in the following program: 94°C for 10 min, followed by 40 cycles of 94°C for 15 s and 60°C for 30 s. Transcription levels of selected genes were quantified by the cycle threshold 2^−ΔΔCt method. Sea buckthorn 18S rRNA gene was applied to normalize the expression values. The primers for PCR amplification are listed in Supplementary Table 6. The experiment was conducted with three biological replicates.

Acknowledgments

This work was supported by the Fundamental Research Funds of CAF (CAFYBB2020SY006) and the Joint Funds of the National Natural Science Foundation of China (U2003116).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by the Fundamental Research Funds of CAF [CAFYBB2020SY006]; the Joint Funds of the National Natural Science Foundation of China [U2003116].