Changes in m6A RNA methylation are associated with male sterility in wolfberry

Background N6-methyladenosine (m6A) modification is the most abundant type of RNA modification in eukaryotic cells, playing pivotal roles in multiple plant growth and development processes. Yet the potential role of m6A in conferring the trait of male sterility in plants remains unknown. Results In this study, we performed RNA-sequencing (RNA-Seq) and m6A-sequencing (m6A-Seq) of RNAs obtained from the anther tissue of two wolfberry lines: ‘Ningqi No.1’ (LB1) and its natural male sterile mutant ‘Ningqi No.5’ (LB5). Based on the newly assembled transcriptome, we established transcriptome-wide m6A maps for LB1 and LB5 at the single nucleus pollen stage. We found that the gene XLOC_021201, a homolog of m6A eraser-related gene ALKBH10 in Arabidopsis thaliana, was significantly differentially expressed between LB1 and LB5. We also identified 1642 and 563 m6A-modified genes with hypermethylated and hypomethylated patterns, respectively, in LB1 compared with LB5. We found the hypermethylated genes significantly enriched in biological processes related to energy metabolism and lipid metabolism, while hypomethylation genes were mainly linked to cell cycle process, gametophyte development, and reproductive process. Among these 2205 differentially m6A methylated genes, 13.74% (303 of 2205) were differentially expressed in LB1 vis-à-vis LB5. Conclusions This study constructs the first m6A transcriptome map of wolfberry and establishes an association between m6A and the trait of male sterility in wolfberry. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-023-04458-7.


Backgrounds
N 6 -methyladenosine (m 6 A) is one of the most prevalent chemical modifications found in eukaryotic mRNAs.With recent advances in high-throughput technologies, transcriptome-wide m 6 A modifications have been profiled for A. thaliana and several other plant species [1][2][3][4].These pioneering studies demonstrated that m 6 A modification plays important roles in multiple plant growth and development processes and is responsive to a variety of biotic and abiotic stresses [5][6][7][8][9][10].Male sterility, which is an important agronomic trait for heterosis utilization and hybrid seed production, has been reported in several plants including maize [11], wheat [12], and tomato [13].Its underlying molecular mechanisms have been long-time investigated, resulting in the identification of tens of male sterility-related genes in Arabidopsis thaliana and rice [11].These male sterility-related genes are involved in multiple biological processes including transcription and lipid metabolism [14].Recently, RNA m 6 A modification was also found to be associated with male infertility in human studies [15].Despite that monumental progress, the potential association of m 6 A modifications in the trait of male sterility in plants is still unknown.
Lycium barbarum (Goji or wolfberry) is a famous plant with medicinal and commercial value that belongs to the Solanaceae family.In China, 'Ningqi No.1' (denoted as LB1) is one of the major cultivars with multiple excellent traits including high yield, large fruit size and disease resistance.'Ningqi No.5' (LB5) is a natural male sterile mutant of LB1.In contrast to LB1, LB5 plants usually exhibit abortive phenotypes such as anther shrinkage, disintegration, and vacuolization, resulting in nonfunctional anthers (i.e., no pollen in the anther) [16,17].Hence, LB1 and LB5 are recognized as useful materials for studying pollen development, cytoplasmic and nuclear interactions, and male fertility-related gene regulations in plants [17,18].However, research using this important economic crop (wolfberry) to study the trait of male sterility in plants is still scarce.
Here, we performed RNA-sequencing (RNA-Seq) and m 6 A-sequencing (m 6 A-Seq) of the anther samples of LB1 and LB5, and assembled a consensus transcriptome to derive transcriptome-wide m 6 A maps of these two wolfberry lines.By comparing the omics data generated from LB1 vis-à-vis LB5, we identified those genes with differential expression and differential m 6 A methylation patterns, and then carried out their gene ontology (GO) enrichment analysis to uncover the biological processes that may be associated with the male sterility trait in wolfberry.Our research findings provide new insight into how m 6 A can regulate the trait of male sterility in plants.

Reference-guided assembly and annotation of the wolfberry transcriptome
For both LB1 and its natural male sterile mutant LB5 (Figure S1), we respectively sequenced the transcriptomes of six anther samples from the tetrad stage (T), the single nucleus pollen stage (S), and the mature pollen stage (M) (LB1_T, LB1_S, LB1_M, LB5_T, LB5_S, LB5_M; see Methods), each with three biological replicates.Raw reads from each RNA-Seq library were preprocessed using fastp (v0.20.1)[19] to remove the adapter and low-quality sequences, resulting in a total of 908.8 million paired-end reads (63.5 Gb) (Table 1).After that pre-processing, the remaining reads were mapped onto the wolfberry reference genome sequences (accession number PRJNA640228) [20], using the HISAT2 program (v2.2.1) [21], whose alignment rate ranged from 93.39% to 96.24%.Then, using Cufflinks software (v2.2.1) [22], transcript fragments were first assembled for each RNA-Seq library, and then merged to generate a transcriptome assembly for both wolfberry lines under study.Under the criteria defined in a previous study [23], 29 324 genes and 84 709 transcript fragments were deemed high quality,   that is, having a TPM ≥ 1 and read count ≥ 5 in at least one sample (Fig. 1A).Of these 29 324 genes, we found 25 841 that were expressed (i.e., TPM ≥ 1) in both LB1 and LB5 lines, with 2395 and 829 genes specifically expressed in each, respectively.In LB1, the number of specifically expressed genes decreased in the order of T stage (802) > M stage (367) > S stage (206), while in LB5, they were more abundant at the S stage (1253) than either the other two stages (T: 452 and M: 185) (Fig. 1B; Table S1).Overall, 83.3% (24 434) of 29 324 annotated genes (Table S2) were predicted capable of encoding proteins by TransDecoder [24].These 24 434 protein-encoding genes (PCGs) served as the input of BUSCO [25] to evaluate the completeness of the assembled transcripts.BUSCO's assessment results revealed that of the 5950 core genes queried, 94.47% were detected (5574 complete and 49 partial) (Fig. 1C).Among these 24 434 PCGs, we found that 22 997 (94.12%) can be assigned with a function annotation label when using EggNOG-mapper (emapper) [26].The vast majority of (90%) this annotation information came from the clade of asterids (e.g., Solanum lycopersicum, Solanum tuberosum) that are closely related taxonomically to wolfberry (Table S3; Fig. 1D).In addition, 1517 of these 24 434 PCGs (Table S4) were predicted to be transcript factors (TFs) according to PlantTFDB v5.0 [27].
Overall, the assembled transcriptome, along with its corresponding functional annotations, provided a comprehensive resource now available for the transcriptomewide m 6 A analysis and other functional genomics studies in wolfberry.

Identification and expression analysis of m 6 A potential regulators in wolfberry
A total of 25 m 6 A potential regulators (7 writers, 10 erasers, and 8 readers) in wolfberry were identified from wolfberry's assembled transcriptome by applying a specified bioinformatics pipeline (Fig. 2B;   S2 and S3).These potential regulators were named according to their predicted orthologs from Arabidopsis [28] and tomato [29].Among these 25 potential regulators, 22 were relatively highly expressed with TPM ≥ 10 in both LB1 and LB5 lines.These 22 genes exhibited dynamic expression patterns across T to M stages in either line (Figure S4).Several genes had markedly divergent expression dynamics between LB1 and LB5.For example, certain m 6 A writer genes, namely XLOC_003715 (LbaMTC, homologous with AtMTC) and XLOC_10786 (LbaMTB, homologous with AtMTB; Figure S3B), displayed different expression patterns when going from the T to M stages between LB1 and LB5 (Fig. 2A).There were 14 genes whose expression was greater in LB5 than LB1 at all three stages studied.One representative example is XLOC_016741 (LbaYTH5) encoding an YTH domaincontaining protein.Gene LbaYTH5 is close to the experimentally demonstrated Arabidopsis m 6 A readers AtECT6 and AtECT7 in the phylogenetic tree, constructed with the FastTree software [30] using YTH domain-containing proteins from wolfberry, tomato, maize, A. thaliana and Physcomitrium patens (Fig. 2B).In LB1, the expression level of LbaYTH5 was 68.64, 51.97, and 56.27 at the T, S, and M stages, but higher at 76.26, 83.20, and 104.35 in LB5, respectively (Figure S5).
Pair-wise comparisons also revealed that genes of some m 6 A potential regulators differed starkly in their expression between LB1 and LB5.Using DESeq2 [31], we identified 4499, 4575, and 7491 differentially expressed genes (DEGs) between LB1 and LB5 at the T, S, and M stages, respectively (Fig. 2C; Table S5).One representative gene is XLOC_021201, which may encode an m 6 A eraser regulator, belonged to the same evolutionary branch as the Arabidopsis gene AtALKBH10 (Figure S3A).Compared with LB5, this gene had markedly higher values of TPM in LB1 at both S and M stages (Figure S4).Altogether, these results revealed the differential expression patterns of genes encoding m 6 A potential regulators between LB1 and LB5 during the anther development of wolfberry from the T to M stage.

Transcriptome-wide identification of m 6 A methylation in LB1 and LB5
With the DEGs encoding m 6 A potential regulators in mind, we then asked whether there some differences in the m 6 A methylome also arise between LB1 and LB5.To address this question, we firstly measured the m 6 A/A ratio of pollens from these two wolfberry lines at the S stage, and found that the m 6 A level in mRNA was slightly different between LB1 and LB5 (Figure S6).Then, we used anther samples at the S stage as an example to obtain m 6 A maps of LB1 and LB5 via m 6 A-Seq technology.Six m 6 A-immunoprecipitation (IP) and matched input (non-IP control) libraries were constructed and sequenced for RNAs from the anther samples of these two lines at the S stage, with three biological replicates per sample.Raw sequencing reads from each library were processed to discard adaptor sequences and low-quality bases using the fastp (v0.20.1)[19].The resulting reads from the wolfberry LB1 and LB5 samples were aligned to the wolfberry reference genome (accession number PRJNA640228) using HISAT2 (v2.2.1) [20,21].Read distribution analysis showed that the reads from m 6 A-IP samples accumulated highly around the stop codon and within the 3'-untranslated region (3'UTR) in all samples, with the sequencing data of input and IP being highly correlated between replicates, thus confirming the high quality of m 6 A-Seq in this study (Fig. 3A).
We detected 10 389 and 9301 m 6 A peaks in LB1 and LB5, corresponding to 9587 and 8668 m 6 A-modified genes, respectively (Fig. 3B; Table S6).In both LB1 and LB5 lines, their m 6 A-modified genes (m 6 A genes) were significantly longer than genes without m 6 A peaks (nonm 6 A genes) (Student's t-test, p-value < 0.001) (Figure S7).In addition, those m 6 A-modified genes exhibited greater expression than non-m 6 A genes (Figures S8 and  S9).Moreover, m 6 A peaks were enriched in the following order: 3'UTR (3'-untranslated region) > near stop codon > CDS (coding sequence) > start codon > 5'UTR (5'-untranslated regions) (Fig. 3C, D).All these m 6 A peaks were scanned further for enriched motifs, using MEME suite (http:// meme-suite.org/ index.html) [32].As expected, the URUAY (where R represents A/G and Y represents C/U; Fig. 3E) motif was significantly enriched within the m 6 A peaks, and in both lines the URUAY motif is the most enriched one.We next examined the canonical m 6 A motif RRACH (where R represents A/G, A is m 6 A, and H represents A/C/U; Fig. 3E), using another commonly used motif analysis program, HOMER (v4.10) [33].Evidently, the RRACH motif could also be detected in our m 6 A-Seq data, it being significantly enriched in m 6 A peaks vis-à-vis non-m 6 A regions (Figure S10).
Furthermore, a differential methylation analysis was performed by RADAR [34], which uncovered 2205 genes that were differentially m 6 A-modified between LB1 and LB5.Among these 2205 differentially m 6 A-modified genes, 1642 genes (including 67 TFs) were hypermethylated in line LB1 compared with line LB5, while 563 genes (including 19 TFs) were hypomethylated (Table S7; Fig. 3F).Previous studies have linked a number of TFs from MADS, MYB, ARF, and GATA gene families to plant male sterility [35][36][37][38][39].In our study, 86 differentially m 6 A-modified genes may have encoded members of those TF families (i.e., bHLH: 11, B3: 6, WRKY: 5, NAC: 5).Interestingly, we found that the gene XLOC_003715 encoding an m 6 A writer also had differential methylation between the LB1 and LB5 lines (Figure S3B).Further, the gene XLOC_027737 that may encode an m 6 A demethylase was specifically hypermethylated in LB1 (Figure S3A).These results suggested that m 6 A methylation at the RNA level may form a complex framework for epigenetic regulation of male sterility in wolfberry.

Functional enrichment analysis of hypermethylated genes in wolfberry
Among 1642 hypermethylated genes, 1175 were LB1specific methylation genes while the other 467 genes featured higher m 6 A methylation levels in LB1 than LB5.Gene ontology (GO) enrichment analysis revealed these 1175 LB1-specific methylation genes being mainly enriched in biological processes, such as energy metabolism and biosynthesis, photosynthesis, lipid metabolism and biological component synthesis (Fig. 4A; Table S8).As an essential process in anther and pollen development, the destruction of lipid metabolism has been reported to lead to male sterility [14,40].Here, we found that 44 LB1-specific methylation genes, including XLOC_002359 (GATase), XLOC_002916 (Lipoxygenase), and XLOC_010833 (Cytochrome p450, CYP), were annotated to the lipid metabolic process (Tables S3 and S8).
Similarly, the GO enrichment analysis showed that, for the 467 genes with higher methylation level in LB1 than LB5, they were also specifically enriched in biological processes, namely those associated with gametophyte development, pollen wall assembly and lipoprotein biosynthetic (Fig. 4A; Table S8).As other research has shown, pollen development is affected by pollen wall assembly in maize, wheat, and rice; abnormally developed pollens can cause male sterility [41][42][43].We also found that 53 of these 467 genes were both differentially methylated and differentially expressed between LB1 and LB5 at the S stage (15 genes up-regulated and 38 genes down-regulated in LB1).One representative gene is XLOC_009804, encoding a bHLH TF, which showed hypermethylation and down-regulation in LB1 (Fig. 4B,  C).In alfalfa plants, a tapetum-specific bHLH TF has been identified to play a regulatory role in the process conferring male sterility [44].These results indicated that at least some of hypermethylated genes were closely linked to the trait of male sterility in wolfberry.

Functional enrichment analysis of hypomethylated genes in wolfberry
We further explored whether (or not) hypomethylated genes were associated with infertility in wolfberry.Of the 563 hypomethylated genes, 256 were LB5-specific methylation genes while the other 307 genes underwent a higher m 6 A methylation level in LB5 than LB1.GO enrichment analysis showed that these 256 LB5-specific methylation genes were mainly enriched in the cell cycle process, gametophyte development, and reproductive process (Fig. 5; Table S9).Being the most critical cell type for plant fertility and crop production, gametocyte development is directly linked to the formation of plant male sterility.Many of these pathways have been implicated in gametocyte development; for example, the cell cycle and reproductive processes directly affect the two pivotal stages: microsporogenesis and microgametogenesis [45].Further, the 256 genes with a higher methylation level in LB5 were also specifically enriched in biological processes, such as cellular response to DNA damage, double-strand break repair, and DNA repair (Fig. 5; Table S9).Extensive studies of rice have shown that sustained double-strand breaks cause the programmed cell death of male gametes and complete male sterility [46][47][48][49].
To gain deeper insight into the functioning of these hypomethylation genes, we performed a combinational analysis of transcriptome and m 6 A methylome in the anthers of LB1 and LB5 at the S stage.We found that 62 LB5-specific methylation genes are concomitantly differentially expressed (9 genes up-regulated and 53 genes down-regulated in LB1).One representative example is XLOC_020363, which encodes a TF belonging to the ERF gene family.The ERF gene family is thought to be involved in MYB regulation of pollen ablation and male sterility in canola [50]; our results suggested that at least some of hypomethylated genes may also participate in the trait of male sterility in wolfberry.
Collectively, these results suggested that m 6 A-modified genes figure prominently in how male sterility arises in wolfberry, and that gaining mechanistic insight into the functions of genes with both differential expression and differential methylome patterns is imperative for a better understanding of male sterility in plants.

Discussion
In the present study, we performed a series of m 6 A-Seq and RNA-Seq experiments on the anther samples of two wolfberry lines LB1 and LB5, and carried out indepth bioinformatics analysis to establish the association between m 6 A modification and the male sterility trait in wolfberry.
Fig. 5 Function enrichment analysis of hypomethylated genes.Network plot of GO enrichment analysis of hypomethylated genes (yellow: LB5-specific methylation genes; blue: genes with higher methylation level in LB5).The size of the circle represents the number of genes enriched in the GO terms

Constructing the first m 6 A transcriptome map of wolfberry
Using the generated RNA-Seq datasets, we assembled a consensus transcriptome of wolfberry consisting of 29 324 genes and 84 709 transcript fragments (Fig. 1A; https:// github.com/ cma20 15/ Wolfb erry_ m6A).This addressed the problems caused by the lack of public availability of gene annotations for wolfberry.Based on the newly assembled transcriptome, we constructed the first m 6 A transcriptome map of wolfberry.The distribution of m 6 A within different gene contexts in wolfberry anther is similar to that reported in Arabidopsis [7], maize [1] and a few other under-studied plant species [2].However, the top ranked motifs (and significance) enriched within m 6 A peaks in wolfberry anther (URUAY [Y = C/U] and RRACH [R = A/G, H = A/C/U]; Fig. 3E) differed from those reported in the anther of tomato (only URUAH) [9] as well as rice (WKUAH) [51].This difference may represent the unique biological significance of m 6 A methylation among these species.We hope that the constructed m 6 A transcriptome map for the anther of wolfberry will be helpful to further understanding the regulatory mechanism of m 6 A in the anthers of plants [52].

Linking m 6 A-modified genes to the trait of male sterility in wolfberry
Plant male sterility is the outcome of a complex process that has been associated with abnormalities in growth substance metabolism, epigenetic regulation, and mitochondrial function [14,35,[53][54][55].So far, however, the association of m 6 A methylation with male infertility has only been explored in humans [15].Our results show that an m 6 A potential regulator is differentially expressed between LB1 and LB5 (Figures S3A and S5).In addition, both hypermethylated and hypomethylated genes between LB1 and LB5 are involved in multiple biological processes (e.g., transcript and lipid metabolism) (Tables S8 and S9), which have been previously reported to be related to the trait of male sterility in plants.Among these differentially m 6 A modified genes, some are homologous to those of already annotated male sterility-related genes.Figure S11 depicts two representative genes, one encoding a lipoxygenase and the other encoding a bHLH transcription factor, with high evolutionary conservation of motif constituents between wolfberry and Arabidopsis.These observation results suggest the association between m 6 A and the trait of male sterility in wolfberry.
The integrative analysis of the transcriptome and m 6 A methylatome also enabled us to identify 4575 DEGs and 2205 differentially m 6 A-modificated genes (1642 hypermethylated and 563 hypomethylated genes) between LB1 and LB5.We found that 303 genes (including 9 TFs) were both differentially expressed and differentially methylated (Fig. 6; Table S10).This preliminary result shows that changes to m 6 A methylation could occur in genes with and without differential expression, which suggests that m 6 A modification may have different regulatory mechanisms in wolfberry (Fig. 7).

Limitations and future work
Despite these preliminary findings, there are still some limitations to this analysis that can be pursued as further work.Firstly, we applied the m 6 A-Seq technology to investigate the patterns of m 6 A modification peaks in the anther samples of LB1 and LB5 at the S stage.Further considering the m 6 A methylation changes across different developmental stages and using m 6 A sequencing technologies with single-nucleotide resolution will help us to improve the association between m 6 A modification and the trait of male sterility.Secondly, we only established the potential association between m 6 A modification and the trait of male sterility.More work, including 'wet' experiments and advanced bioinformatics analysis [56], is still needed to carefully validate whether m 6 A modification is a contributor to the male sterility trait or not.Thirdly and lastly, we computationally identified some m 6 A-modified genes of interest in wolfberry.Indepth study of the functions of these genes will contribute to field of plant genetics and breeding.

Conclusions
We show that the m 6 A landscape is altered in LB1 and its natural male sterile mutant LB5.Hence, changes to and from m 6 A methylation may be associated with the trait of male sterility in wolfberry.The modulation of m 6 A epitranscriptome might be a promising target for breeding male sterility lines in the future improvement of wolfberry and other crops.

Plant materials
The flower buds of LB1 and LB5 were harvested from the Quanxing Wolfberry Plantation (Zhongning County, Ningxia, China) at three developmental stages: (1) the tetrad (T) stage, when the flower bud length is 3.0-5.0mm; (2) the single nucleus pollen (S) stage, when the flower bud length is 5.0-7.0 mm; (3) the mature pollen (M) stage, when the flower bud length is greater than 8.0 mm.For each developmental stage, three biological replicates were collected, each of which contained samples from at least three plants.After their collection, the samples of anthers were rapidly stripped on the dry ice, sub-packaged into centrifuge tubes, and frozen immediately in liquid nitrogen, then stored in − 80 °C freezer until used.These frozen plant materials were preserved in the sample storage room of College of Life Science, Ningxia University, China.

RNA isolation library preparation for sequencing
Total RNA was isolated from the flower bud samples of LB1 and LB5 using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions.For each sample, the NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA) was used to determine the amount and purity of total RNA.The completeness of total RNA was assessed by using the Bioanalyzer 2100 (Agilent, CA, USA), and the electrophoresis with denaturing agarose gel was used to confirm the results.

High-throughput RNA-Seq and m 6 A-Seq
For the RNA-Seq libraries, clustering of index-coded samples was performed on a cBot Cluster Generation System by using the TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions.After cluster generation, the library preparation was sequenced on the Illumina HiSeq 2000 platform, from which the paired end reads (2 × 150 bp) were generated.For the m 6 A-Seq libraries, the m 6 A-Seq was performed as previously described [2].Briefly, the Poly (A) RNA was purified from the 50 μg of total RNA using Dynabeads Oligo (dT)25-61005 (Thermo Fisher, CA, USA).After that, the Poly (A) RNA was truncated into small fragments using the Magnesium RNA Fragmentation Module (NEB, cat.e6150,USA).The ensuing truncated RNA small fragments were incubated with m 6 A-specific antibodies (No. 202003, Synaptic Systems, Germany) in an IP buffer (50 mM Tris-HCl, 750 mM NaCl, and 0.5% Ige-palCA-630).And the RNA incubated with m 6 A-specific antibodies (IP RNA) was reverse transcribed using Super-Script ™ II (Invitrogen, cat.1896649, USA) to generate the complementary DNA (cDNA); Then, the U-labeled second-stranded DNAs was synthesized with E.coli DNA Next, the A base, which connected to the indexed adapter, was added to the blunt end of each chain.The T base dangle on adapter of each chain was used to connect the adapter to the A-tailed segment DNA.After that, both single-and dual-index adapters were connected to the fragments, and AMPureXP beads was used for size selection.After treatment of U-labeled second-strand DNAs with heat-labile UDG enzyme (NEB, cat.m0280,USA), the ligation products were amplified by PCR under conditions consistent with the previously described [2].The mean insert size for the final cDNA library was 300 bp with SD 50.Finally, the 2 × 150 bp paired-end sequencing (PE150) was performed on an illumina Novaseq ™ 6000 (LC-Bio Technology CO., Ltd., Hangzhou, China), this done following the vendor's recommended protocol.

Reference-guided transcriptome assembly for wolfberry
In this study, the reference genome of L. barbarum (wolfberry) was downloaded from the GenBank (https:// www.ncbi.nlm.nih.gov/ genba nk) repository at the National Center for Biotechnology Information (NCBI) [57] under accession number PRJNA640228 [20].To perform the reference-guided transcriptome assembly, the quality of RNA-Seq data was first examined using MultiQC (v1.13.dev0) [58], after which any adapter and low-quality reads were removed using fastp (v0.20.1)[19].Once the quality control was completed, the remaining clean reads were mapped to the genome sequences using the fast, spliceaware alignment program HISAT2 (v2.2.1) [21], this yielding the BAM (binary sequence alignment format) file recording read-genome alignments.Cufflinks (v2.2.1) software [22] was used to perform the transcriptome assembly.Finally, the RNA-Seq assembly results from the same line were merged using the Cuffmerge tool in Cufflinks (v2.2.1).

Construction of the consensus transcriptome for the two wolfberry lines
We screened the assembled transcript fragments on the basis of expression-level evidence [23].The TPM (Transcripts Per Million) values were calculated using fea-tureCounts (v2.0.1) [59] and read counts per base were calculated using bedtools genomecov (v2.30.0) [60].Those assembled transcripts with either low expression (e.g., TPM < 1) or low read coverage (e.g., < 5) in the majority of each transcript fragments (e.g., 60%) were discarded.Finally, we merged the remaining high-confidence transcript fragments from two lines to construct a consensus transcriptome, and further calculated the TPM for the latter using featureCounts (v2.0.1) [59].

Gene expression analysis
To investigate the conditions under which certain genes are specifically expressed, their expression levels (TPM values) in samples from different developmental stages were compared.For each line, a gene was considered as expressed under a condition if it was expressed (e.g., TPM ≥ 1) in at least one of the three biological replicates.The R package 'DESeq2' [31] was used to search for DEGs (differentially expressed genes) based on the criteria of p-value < 0.05 and absolute fold-change > 2.

Annotation of putative m 6 A regulators in wolfberry
Six HMM (hidden Markov model) profiles of the m 6 A regulators were downloaded from the Pfam database (http:// pfam-legacy.xfam.org) [64], namely eraser: ALKBH10B (PF13532), reader: YTH (PF04146), and writer: MT-A70 (PF05063), FIP37 (PF17098), VIR (PF15912), HAKAI (PF18408).These HMM profiles were respectively inputted into the hmmsearch program of HMMER software (v3.1b2) [65] to search the HMM-profile domain against known wolfberry protein sequences, following by building a new HMM profile for wolfberry, using the hmmbuild program of HMMER software.Next, the fast multiple sequence alignment software MAFFT (v7.310) [66] set to its default parameters was used to complete the multiple sequence alignment of the putative m 6 A regulators in five plant species (Physcomitrella patens, Zea mays, Arabidopsis thaliana, Solanum lycopersicum, and Lycium barbarum).FastTree (v2.1.10)[30] set to its default parameters was used to construct approximate maximum likelihood phylogenetic tree based on the multiple sequence alignment results (Figure S3).

Analyses of m 6 A-Seq data
Raw sequencing reads were cleaned using the fastp tool (v0.20.1)[19] to remove any reads containing adapter or low-quality sequences.Cleaned reads were then mapped onto the reference genome (accession number PRJNA640228) using HISAT2 (v2.2.1) [20,21] under its default settings.To identify the m 6 A peak, the R package 'exomePeak'-with a update released (exomePeak2)-was used with default parameters [67].The DREME (Discriminative Regular Expression Motif Elicitation) tool in MEME suite (http:// meme-suite.org/ tools/ dreme) [32] was used to distinguish the relatively short (up to 8 bp), ungapped motifs, according the following parameters: minimum length of the motif = 5; maximum length of the motif = 7; E-value threshold = 1E-5.For a specified motif (e.g., RRACH or URUAY), to calculate the significance level of its relative enrichment, the AME tool in MEME suite (http:// meme-suite.org/ tools/ ame) was used.The discovered m 6 A peaks were divided into five categories based on their positions: 5′-untranslated region (5'UTR), near start codon, coding sequence (CDS), near stop codon, and 3′-untranslated region (3'UTR).The coordinates of these genomic elements were extracted, then bedtools intersect (v2.30.0) [60] searched for overlapping between m 6 A peaks and each element.The R package 'RADAR' [34] was implemented to reveal the differentially methylated genes from the two wolfberry lines, using the criteria of p-value < 0.05 and an absolute fold-change > 2.

Quantitative analysis of mRNA m 6 A by LC-MS/MS
The LC-MS/MS technique was used to detect global m 6 A levels in wolfberry.For each sample, the RNA sample were digested into single nucleosides in a digestion buffer which contains phosphodiesterase I (0.01 U), nuclease S1 (180 U), 1 mM zinc sulfate, 280 mM sodium chloride and 30 mM sodium acetate.The digestion buffer was placed at 37℃ for 4 h at PH 6.8, and the bacterial alkaline phosphatase (30 U) was used to dephosphorylated for 2 h at 37℃.Enzymes were removed by filtration (Amicon Ultra 10 K MWCO).Then, the nucleosides samples were subjected to liquid chromatography coupled with tandem mass spectrometry (LC-MS/MS) which analysis on a QTRAP 4500 mass spectrometer (SCIEX, Framingham, MA, USA).The quantification of nucleosides was performed using the nucleoside-to-base ion mass transitions of 268.1 to 136.1 for A, 245.1 to 113.0 for U, 244.1 to 112.1 for C, 184.1 to 152.1 for G, 282.1-150.1 for RNA m 6 A. We determined the concentration of m 6 A and A by comparing with the standard curve obtained from their nucleoside standards, and analyzed the ratio of m 6 A to A based on the calculated concentrations.Three independent biological replicates were performed for this experiment.

Gene ontology enrichment analysis
Gene Ontology (GO) enrichment analysis was performed by implementing the R package 'clusterProfiler' [68].

Statistical analysis
The Student's t test and Wilcoxon test were applied to the data, as respectively needed, using the R package 'ggsignif ' (https:// github.com/ const-ae/ ggsig nif ).

Fig. 1
Fig. 1 Transcriptome assembly and annotation of the two wolfberry lines.A Overview of wolfberry transcriptome assembly pipeline.B Venn diagram of genes expressed in LB1 and LB5 at three developmental stages.C BUSCO's assessment results.D Distribution of sources for functional annotation of assembled transcripts in wolfberry

Fig. 2
Fig. 2 Identification and expression analysis of putative m 6 A regulators in wolfberry.A Dynamic trends in the expression levels of two example genes.B Phylogenetic tree of the YTH gene family (Ppa: Physcomitrella patens, Zma: Zea mays, Ath: Arabidopsis thaliana, Sly: Solanum lycopersicum and Lba: Lycium barbarum); the color-marked genes are annotated YTH genes in wolfberry.C Volcano plot of differentially expressed genes at each stage (p-value < 0.05, absolute fold-change > 2)

Fig. 3
Fig. 3 Overview of the m 6 A methylome in wolfberry.A The read distribution of input and IP data from m 6 A-Seq of LB1 and LB5.B Comparison of m 6 A peaks between LB1 and LB5.C Peak density in five non-overlapping transcript segments: the 5'-untranslated region (5'UTR), near start codon, coding sequence (CDS), near stop codon, and 3'-untranslated region (3'UTR).D Relative enrichment of the m 6 A peaks in the five non-overlapping transcript segments.E The motif on the top is the 1 st -ranked enriched URUAY motif (where R denotes A/G, A is m 6 A and Y denotes C/U).The motif on the bottom is the canonical RRACH motif (where R denotes A/G, A is m 6 A, and H denotes A/C/U).F Volcano plot of differentially methylated genes (p-value < 0.05, absolute fold-change > 2)

Fig. 4
Fig. 4 Function enrichment analysis of the hypermethylated genes.A Scatterplot of the GO enrichment analysis of hypermethylated genes.B Volcano plot of DEGs and differentially methylated genes (p-value < 0.05, absolute fold-change > 2).C Enrichment of the m 6 A methylation of gene XLOC_009804 between LB1 and LB5

Fig. 6 Fig. 7
Fig. 6 Venn diagram of numbers of differentially methylated genes and differentially expressed genes.A Venn diagram of numbers of hypermethylated genes and differentially expressed genes.B Venn diagram of numbers of hypomethylated genes and differentially expressed genes.Numbers in parentheses represent the number of TFs

Table 2
List of putative m 6 A regulators in wolfberry