1 Introduction
The Rosaceae family is one of the most economically and ecologically important flowering plant families worldwide, encompassing diverse fruit crops, ornamental plants, and medicinal species [
1–
3]. Within this family, the
Rosa genus comprises over 200 species with global distribution, renowned for their unparalleled ornamental value, profound cultural significance, and rich genetic diversity that serves as a core resource for horticultural and medicinal breeding [
4,
5]. Evolutionary and functional genomic studies of
Rosa species have thus attracted extensive attention to unlock their breeding potential and elucidate the genetic basis of key agronomic traits.
Rosa laevigata Michx., commonly known as JīnYīngZǐ or Cherokee rose, is a perennial woody plant native to southern China, with unique biological, medicinal, and horticultural values [
6]. As a classic medicinal material recorded in the
Chinese Pharmacopoeia, its fruits have been used in traditional Chinese medicine for thousands of years [
7]. Modern phytochemical studies have revealed its abundant bioactive components, including extremely high vitamin C content, diverse flavonoids, polyphenols, and triterpenoid saponins, with proven antioxidant, anti-inflammatory, and hepatoprotective activities, highlighting its great potential in functional food and pharmaceutical development [
8–
10]. Additionally,
R. laevigata exhibits excellent tolerance to drought, cold, and pathogens, making it a valuable genetic resource for resistance breeding. Notably, it is a diploid species with a karyotype of 2
n = 2
x = 14, providing an ideal model for genomic and functional studies within the
Rosa genus.
High-quality reference genomes are fundamental for deciphering plant evolutionary history, genetic regulatory mechanisms underlying key traits, and advancing molecular breeding. Existing fragmented assemblies for related
Rosa species fail to resolve highly repetitive genomic regions, including centromeres, telomeres, and long terminal repeat retrotransposons (LTR-RTs) [
11,
12]. These limitations severely restrict in-depth exploration of
R. laevigata’s unique biological characteristics, functional verification of medicinal metabolite biosynthetic pathways, and development of efficient molecular breeding strategies.
Telomere-to-telomere (T2T) gap-free genome assembly, the current gold standard for genomic research, provides complete and accurate representation of all chromosomal regions from telomere to telomere, including previously intractable highly repetitive regions [
13–
15]. This strategy enables precise genome annotation, comprehensive structural variation detection, and high-resolution comparative genomic and evolutionary analyses with unprecedented accuracy [
16]. To overcome the limitations of existing genomic resources and unlock the genetic potential of
R. laevigata, we integrated PacBio HiFi long reads, Oxford Nanopore ultra-long reads, and high-throughput chromosome conformation capture (Hi-C) technology to construct the first chromosome-level, gap-free T2T reference genome of
R. laevigata. We further performed comprehensive genomic, transcriptomic, and metabolomic analyses to elucidate its evolutionary history and the transcriptional regulatory network of flavonoid biosynthesis. This high-quality T2T genome will serve as an invaluable resource for dissecting the molecular basis of its unique medicinal and agronomic traits, accelerating
Rosa molecular breeding, and providing critical insights into Rosaceae evolutionary dynamics.
2 Materials and Methods
2.1 Plant materials and sample collection
All wild healthy individuals of R. laevigata were collected from their native habitats in Hunan Province, China, and authenticated by Central South University of Forestry and Technology. Young, disease-free tender leaves were harvested from a single highly homozygous individual for genome sequencing. For sequencing-library construction, approximately 2 g of young leaf tissue was used for genome survey sequencing, 3 g for PacBio HiFi sequencing, 5 g for Hi-C library construction, and 15 g for ultra-long ONT sequencing. Five distinct tissue types (roots, stems, leaves, mature fruits, and seeds) were collected at the fully mature fruit stage for transcriptomic and metabolomic analyses, with approximately 0.5 g of material collected for each omics sample. For the spatiotemporal fruit development analysis, fruit samples were collected at seven consecutive key developmental stages from July to late December, and approximately 0.5 g of material was collected for each transcriptomic or metabolomic sample. The five tissue types were represented by four biological replicates for both transcriptomic and metabolomic analyses. For the seven fruit developmental stages, transcriptome sequencing used four biological replicates per stage, whereas metabolomic profiling used six biological replicates per stage. All samples were immediately snap-frozen in liquid nitrogen after collection and stored at –80 °C until DNA, RNA, and metabolite extraction.
2.2 Phenotypic observation and data recording
Digital images of various organs and fruits at different developmental stages of R. laevigata were captured and recorded for morphological documentation and comparative analysis. Fruits were continuously monitored throughout the entire developmental period, and key growth parameters including longitudinal diameter, transverse diameter, fresh weight per fruit, and dry weight per fruit were quantitatively analyzed.
2.3 Genome sequencing, assembly, and annotation
2.3.1 Genome size estimation
Genome size, heterozygosity, and repeat content of
R. laevigata were estimated via
K-mer analysis of Illumina sequencing data and flow cytometry [
17]. For
K-mer analysis, 17-mer frequency profiles were generated from high-quality Illumina reads using Jellyfish v2.3.0, and genomic characteristics were calculated using GenomeScope v2.0 [
18]. Assembly statistics were evaluated using QUAST v5.2.0 with default assembly-statistics settings. For flow cytometry analysis, intact nuclei were isolated from fresh
R. laevigata leaves using a plant nuclear extraction kit, stained with propidium iodide (PI), and analyzed on a BD FACSCalibur flow cytometer (BD Biosciences, San Jose, CA, USA) with
Solanum lycopersicum (tomato) as the internal reference standard [
19]. Genome size was calculated based on the ratio of mean fluorescence intensity of the
R. laevigata sample to that of the reference.
2.3.2 DNA extraction and sequencing
Genomic DNA was extracted from young leaves of R. laevigata. The purity and concentration of the extracted DNA were determined using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and a Qubit 4.0 fluorometer (Invitrogen, Carlsbad, CA, USA), respectively. Four types of sequencing libraries were constructed for the telomere-to-telomere (T2T) genome assembly: (1) an Illumina paired-end library with an insert size of 350 bp was constructed using the TruSeq DNA PCR-Free Library Preparation Kit (Illumina, San Diego, CA, USA) and sequenced on the Illumina NovaSeq 6000 platform with a 150 bp paired-end sequencing strategy; (2) a PacBio HiFi library with an insert size of 15–20 kb was constructed using the SMRTbell Express Template Preparation Kit 2.0 (Pacific Biosciences, Menlo Park, CA, USA) and sequenced on the PacBio Sequel II platform; (3) an Oxford Nanopore Technologies (ONT) 1D sequencing library was constructed following the manufacturer's standard protocol (Oxford Nanopore Technologies, Oxford, UK) and sequenced on the PromethION platform; (4) a Hi-C library was constructed using the Arima-HiC Kit (Arima Genomics, San Diego, CA, USA) following the manufacturer’s standard protocol. This Hi-C library was also sequenced on the Illumina NovaSeq 6000 platform with a 150 bp paired-end strategy.
2.3.3 T2T genome assembly
PacBio HiFi reads were first processed with CCS using a minimum predicted read-quality threshold of min-rq = 0.99, and
de novo assembly of primary contigs was performed using Hifiasm v0.19.8 with default parameters [
20]. For chromosome-level scaffold construction, filtered Hi-C reads were aligned to the assembled contigs using BWA-MEM v0.7.17, and chromosome anchoring, ordering, and orientation were performed using the ALLHiC workflow followed by manual inspection and correction with Juicebox [
21,
22]. ONT ultra-long reads were then incorporated during the gap-filling stage: ultra-long reads were mapped to the Hi-C assembly using minimap2, reads spanning gap-flanking regions were extracted with samtools, and the remaining gaps were filled using TGS-GapCloser v1.2.1 followed by Racon v1.4.20 correction. Telomeric repeats were identified using Tandem Repeats Finder v4.09 with the canonical plant telomeric repeat motif (AAACCCT) as the query sequence, confirming the presence of telomeres at both ends of each pseudochromosome [
23,
24]. Centromeric regions were identified based on the enrichment of centromere-specific satellite repeats and long terminal repeat (LTR) retrotransposons [
25]. Gap filling was performed using PacBio HiFi sequencing data to eliminate remaining gaps, resulting in the final gap-free telomere-to-telomere chromosome-level genome assembly [
26]. Basic assembly statistics (including contig N50, scaffold N50, and total gap number) were calculated using QUAST v5.2.0 [
27]; Assembly consistency and potential exogenous contamination were assessed using short-read mapping, coverage statistics, and GC-depth distribution; no obvious foreign-contamination signal was detected. The final gap-free telomere-to-telomere assembly was evaluated using Benchmarking Universal Single-Copy Orthologs (BUSCO) v5.4.7 analysis against embryophyta_odb10 and by Merqury v1.3 analysis [
28,
29].
2.3.4 Genome annotation
Repeat sequence annotation was performed using a combined approach of
de novo prediction and homology-based alignment: a custom repeat library was constructed using RepeatModeler v2.0.4 and LTR_FINDER_parallel v1.3 [
30,
31], which was then merged with the RepBase database for homology-based repeat masking using RepeatMasker [
32]. For protein-coding gene prediction, an integrated strategy combining
de novo prediction, homology-based prediction, and transcriptome evidence was employed. For
de novo prediction, Augustus v3.5.0, SNAP v1.0, and Genscan v1.0 were used [
33,
34]. For homology-based prediction, protein sequences from five closely related Rosaceae species (
Fragaria vesca,
Rosa chinensis,
Prunus persica,
Pyrus bretschneideri, and
Malus domestica) and
Arabidopsis thaliana were aligned to the
R. laevigata genome to predict gene structures. For transcriptome-based evidence analysis, high-quality RNA sequencing reads from all tissue and developmental stage samples in this study were aligned to the assembled genome using HISAT2 v2.2.1; transcripts were assembled using StringTie v2.2.1, and PASA v2.5.2 was used to optimize untranslated region (UTR) prediction and identify alternative splicing events [
35–
37]. All data obtained from the three strategies were integrated using EVidenceModeler v1.1.1 to generate the consensus gene set [
37]. Functional annotation of the final protein-coding gene set was achieved by aligning protein sequences against public databases including the NCBI non-redundant (Nr) database [
38], Swiss-Prot [
39], Kyoto Encyclopedia of Genes and Genomes (KEGG) [
40], Gene Ontology (GO) [
41], Pfam [
42], and Eukaryotic Orthologous Groups (KOG) [
43] using BLASTP with an
E-value threshold of 1e-5. Protein domain annotation was performed using InterProScan [
44]. Candidate genes involved in the phenylpropanoid/flavonoid biosynthetic pathway were identified by integrating genome annotation, KEGG pathway assignment, BLASTP homology searches against characterized phenylpropanoid/flavonoid genes, conserved-domain confirmation, and comparative synteny with related Rosaceae species. For multi-copy gene families, candidate copies were retained and prioritized based on conserved catalytic domains, sequence similarity, chromosomal/syntenic context.
2.4 Gene family and phylogenomic analysis
2.4.1 Genome-wide synteny analysis
Using MCScanX with default parameters [
45], we identified both intragenomic syntenic blocks within the
R. laevigata genome and intergenomic syntenic blocks between
R. laevigata,
R. rugosa and
R. idaeus, based on all-against-all BLASTP alignments of protein sequences with an
E-value cutoff of ≤ 1 × 10
−5. The identified syntenic blocks, including collinear gene pairs and chromosomal rearrangement events, were visualized using TBtools v2.476 [
46].
2.4.2 Synonymous substitution rate analysis
Synonymous substitution rate (Ks) values of paralogous gene pairs located in intra-genomic syntenic blocks of
R. laevigata and orthologous gene pairs between
R. laevigata and related species were calculated using KaKs_Calculator v3.0 with the YN model [
47]. The Ks frequency distribution was plotted to identify potential whole-genome duplication (WGD) events in the evolutionary history of
R. laevigata.
2.4.3 Analysis of genome evolution
To investigate the evolutionary position of
R. laevigata, gene family clustering was performed using OrthoFinder v2.5.5 with protein sequences from
R. laevigata and nine other plant species, including eight Rosaceae species (
Rosa rugosa,
Fragaria vesca,
Rubus idaeus, Prunus campanulata, Prunus persica,
Pyrus communis,
Malus ioensis,
Crataegus pinnatifida) and one eudicot species (
Vitis vinifera) as the outgroup [
48]. Single-copy orthologous genes shared by all 10 species were extracted, and a maximum likelihood (ML) phylogenetic tree was constructed using IQ-TREE v2.2.6 with ModelFinder-based model selection (-m MFP), automatic thread allocation (-T AUTO), and branch-support evaluation using ultrafast bootstrap and SH-aLRT tests with 1,000 replicates where applicable [
49]. The final ML inference used the best-fit amino-acid substitution model selected by ModelFinder. The divergence time between species was estimated using MCMCTree in the PAML v4.9j package, with calibration points obtained from the TimeTree database [
50,
51]. The calibration constraints were
Vitis vinifera -
Rosa rugosa, 109.8-124.4 Mya;
Malus ioensis -
Prunus persica, 34.4-67.2 Mya;
Prunus persica -
Rosa rugosa, 49.2-77.1 Mya; and
Rubus idaeus -
Rosa rugosa, 28.7-67.1 Mya. Gene family expansion and contraction analysis was performed using CAFE v5 based on the phylogenetic tree and gene family clustering results [
52].
2.5 Transcriptome sequencing and data analysis
Transcriptome sequencing was performed for all five tissue types (roots, stems, leaves, fruits, and seeds) and seven fruit developmental stages of
R. laevigata, with four biological replicates per sample. Total RNA was extracted using Trizol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions. After cDNA library construction, sequencing was carried out on the Illumina NovaSeq 6000 platform with a 150 bp paired-end sequencing strategy. Raw RNA-seq reads were preprocessed using Trimmomatic v0.39 with the following settings: ILLUMINACLIP:TruSeq3-PE.fa:2:30:10, LEADING:3, TRAILING:3, SLIDINGWINDOW:4:15, and MINLEN:36 [
53]. Clean reads were aligned to the assembled
R. laevigata genome using HISAT2 v2.2.1 with default parameters [
35], and gene expression levels were normalized to transcripts per kilobase million (TPM). Differentially expressed genes (DEGs) between different tissues and between consecutive fruit developmental stages were identified using DESeq2 v1.38.3 [
54], with the screening criteria of|log
2(fold change)| ≥ 2 and adjusted
P value < 0.05.
2.6 Metabolome profiling and data analysis
Untargeted metabolomic analysis was performed on the same set of samples used for transcriptome sequencing, including five tissue types and seven fruit developmental stages, with at least four biological replicates per sample. Metabolite separation and detection were accomplished using an ultra-performance liquid chromatography (UPLC) system (Agilent 1290 Infinity II, Agilent Technologies, Santa Clara, CA, USA) coupled with a triple quadrupole-linear ion trap mass spectrometer (AB Sciex QTRAP 6500+, AB Sciex, Framingham, MA, USA) [
55]. Metabolite identification was achieved by comparing retention times, precursor ion masses, and MS/MS fragmentation spectra against an in-house standard metabolite database, the Metlin database, and the Human Metabolome Database (HMDB) [
56]. Metabolite quantification was performed using the multiple reaction monitoring (MRM) mode. Raw metabolomic data were processed using Analyst v1.6.3 and MultiQuant v3.0.3 software [
57]. Differentially accumulated metabolites (DAMs) were screened based on the following criteria: variable importance in projection (VIP) value ≥ 1 in the orthogonal partial least squares discriminant analysis (OPLS-DA) model, |log
2 (fold change) | ≥ 1, and
P value < 0.05 from Student’s
t-test.
2.7 Co-expression and regulatory network analysis
Weighted gene co-expression network analysis (WGCNA) was performed to identify key regulatory modules and hub genes associated with the biosynthesis of key flavonoid compounds in
R. laevigata [
58]. The TPM expression matrix derived from RSEM was used for WGCNA, comprising four biological replicates per developmental stage and 28 RNA-seq samples in total. Before network construction, genes were filtered using a mean expression threshold of > 1 and a coefficient-of-variation threshold of > 0.1. A co-expression network was constructed with the following parameters: soft power = 8, minModuleSize = 30, minKMEtoStay = 0.3, and mergeCutHeight = 0.25. Module-trait relationships were calculated by Spearman correlation analysis using continuous metabolite-trait data for the contents of eight key metabolites, and modules significantly associated with the accumulation of flavonoids and flavan-3-ols were identified based on the combined evaluation of correlation coefficients and corresponding
P values [
59]. The co-expression relationships between structural genes and transcription factors were calculated using the Pearson correlation coefficient (PCC) [
60]. The co-expression regulatory network was visualized using Cytoscape v3.9.1, and hub transcription factors with high connectivity were screened as candidate regulators of flavonoid biosynthesis in
R. laevigata [
61].
2.8 Measurement of physiological indicators
Key physiological and nutritional bioactive indicators of
R. laevigata fruits at seven developmental stages were systematically determined, with three independent biological replicates set for each sample. The measured indicators included the contents of flavonoids, phenolics, saponins, proanthocyanidins, ascorbic acid (AsA), and tannins. In addition, the total antioxidant capacity of the samples was assessed using the ferric reducing antioxidant power (FRAP) assay and the 2,2-diphenyl-1-picrylhydrazyl (DPPH) radical scavenging assay. All procedures were performed following the manufacturer’s instructions of the assay kits from Comin Biotechnology Co., Ltd. (Suzhou, China) [
62,
63].
2.9 Bioinformatics analysis of promoter cis-elements
Genomic sequences of 2000 bp upstream of key structural genes involved in flavonoid and flavan-3-ol biosynthetic pathways were extracted from the genome of
R. laevigata. Cis-acting regulatory elements in these promoter sequences were predicted and annotated using the PlantCARE database [
38,
64]. Statistical analysis was performed on the functional cis-elements, and the distribution of key cis-acting elements on each promoter was visualized using TBtools.
2.10 Bioinformatics analysis of RlmMYB gene family
MYB members in
R. laevigata were comprehensively identified using the Hidden Markov Model (HMM) profile (PF00249), the SMART database, and the NCBI Conserved Domain Database (CDD) [
65,
66]. Chromosomal localization was analyzed and visualized based on the genome annotation information. Multiple sequence alignment of full-length RlmMYB protein sequences was performed using MEGA-X, and a maximum likelihood phylogenetic tree was constructed in combination with MYB protein sequences from
Arabidopsis thaliana [
67]. Expression patterns of
RlmMYB genes in different tissues and fruit developmental stages were extracted from the transcriptome data, and a hierarchical clustering heatmap of the expression profiles was constructed. Meanwhile, molecular docking simulation between the protein sequence of RlmMYB32 and the promoter binding site of RlmLAR was performed using AlphaFold3 [
68].
2.11 Statistical analysis
All quantitative analyses in this study were performed with at least three independent biological replicates, and no quantitative comparison was based on fewer than three biological replicates. Transcriptomic analyses were conducted with four biological replicates per sample. Metabolomic analyses used four biological replicates for tissue samples and six biological replicates for fruit developmental-stage samples, and physiological measurements were performed with three independent biological replicates per stage. Quantitative data are presented as mean ± standard deviation (SD). Statistical analyses were conducted using SPSS v26.0 and Microsoft Excel software. Differences between two groups were evaluated using Student’s t-test, and differences among multiple groups were assessed using one-way analysis of variance (ANOVA) followed by Duncan’s multiple range test. P < 0.05 was considered statistically significant.
3 Results
3.1 Generation and annotation of a gapless genome for R. laevigata
R. laevigata Michx. is an economically and medicinally important Rosaceae species renowned for its bioactive secondary metabolites. Here, we report a high-quality gap-free telomere-to-telomere (T2T) genome assembly of R. laevigata (Fig. 1, Table 1). First, 17-mer K-mer analysis using 25 Gb Illumina reads estimated a genome size of 514.43 Mb, with 0.87% heterozygosity and 54.11% repeat content (Fig. 1H, Table S1). De novo assembly was performed using 51.8 Gb ONT ultra-long reads, 39.1 Gb PacBio HiFi reads, and 50 Gb Hi-C reads (Tables S2 and S3). Hi-C scaffolding anchored contigs into seven pseudochromosomes, yielding a 503.66 Mb gap-free T2T assembly with contig and scaffold N50 values of 67.24 Mb (Tables 1 and S4). Rigorous quality validation confirmed high assembly integrity: BUSCO completeness reached 98.9% (Table S5); read mapping showed 99.39% alignment rate and 99.94% genome coverage (Table S6); Merqury analysis yielded a consensus quality value (QV) of 46.78 (Table S7); and extremely low homologous SNP rate (0.000048%) and heterozygous SNP rate (0.443024%) demonstrated high base accuracy (Table S8). All seven pseudochromosomes possess complete (AAACCCT)n telomeric repeats at both ends, and seven distinct centromeric regions were annotated (Figs. S9 and S10). Genome annotation identified 56.76% repetitive sequences (predominantly transposable elements) (Table S9) and 30,259 high-confidence protein-coding genes with average transcript and CDS lengths of 2816.78 bp and 1152 bp, respectively (Table 1). This high-quality reference genome provides a solid foundation for future genetic and functional studies of R. laevigata.
3.2 Phylogenetic relationships and whole-genome duplication events of R. laevigata
We performed comprehensive phylogenomic analyses using the gap-free T2T genome assembly of R. laevigata generated in this study, combined with genome sequences of eight representative Rosaceae species and Vitis vinifera. A total of 10,957 orthologous gene clusters were identified across the ten species, from which 1151 high-confidence single-copy orthologs were selected to construct the phylogenetic tree, and divergence times were estimated using multiple fossil calibration points from the TimeTree database (Fig. 2A). Phylogenetic analysis revealed that R. laevigata and R. rugosa are the closest sister species within the genus Rosa; molecular dating indicated that R. laevigata diverged from R. rugosa approximately 10.67 million years ago (Mya), the genus Rosa separated from the Fragaria lineage around 22.19 Mya, and the core Rosaceae lineage diverged from the outgroup V. vinifera about 117 Mya. We further analyzed the dynamics of gene family expansion and contraction along the phylogenetic tree, identifying 358 significantly expanded and 278 contracted gene families in R. laevigata (Fig. 2A). KEGG functional enrichment analysis of expanded gene families showed significant enrichment in key biological pathways related to plant secondary metabolism and environmental adaptation (Fig. 2B). Gene family clustering analysis further revealed that among the 17,889 genes assigned to R. laevigata orthologous clusters, 369 were species-specific (Fig. 2C). Orthologous gene sharing analysis among five Rosaceae species (Fig. 2E) identified 471 species-specific genes in R. laevigata, whose KEGG enrichment analysis demonstrated significant overrepresentation in photosynthesis, biosynthesis of secondary metabolites, and plant hormone signal transduction (Fig. 2F), suggesting that these expanded gene families and species-specific genes may contribute to the environmental adaptability and growth characteristics of R. laevigata.
Transposable elements (TEs) are major drivers of plant genome size variation, structural evolution and functional innovation, with their dynamic insertion and amplification playing pivotal roles in shaping plant genome evolutionary trajectories [
69]. Repeat annotation showed that LTR retrotransposons (LTR-RTs) are the most abundant repetitive elements in the
R. laevigata genome (Table S9). In
R. laevigata, 224,367 intact LTR-RTs were identified, including Ty3/Gypsy (13%), Ty1/Copia (23%), and unclassified LTR-RTs detected but not assigned to Copia or Gypsy in this LTR-focused classification (64%) (Fig. 3A). Ty3/G and Ty1/Copia displayed distinct chromosomal distribution patterns (Fig. 3B). Phylogenetic analysis showed that Ty3/Gypsy mainly comprised the Athila, Ogre, Retand, Tekay, Galadriel, Reina, and CRM subfamilies (Fig. 3C), whereas Ty1/Copia included the Angela, Bianca, SIRE, Ikeros, Tork, TAR, Ivana, Ale, and Alesia subfamilies, with Ale and Bianca being the most abundant (Fig. 3D). Insertion time analysis revealed a recent burst of LTR-RT amplification in
R. laevigata, peaking at 0–2 Mya and showing the lowest insertion density among the four compared species (Fig. 3E). Notably, Ty3/Gypsy peaked later than Ty1/Copia and showed a much higher insertion density, suggesting that recent genome expansion in
R. laevigata was mainly driven by gypsy amplification (Fig. 3E). Consistently, comparative TE analysis showed that
R. laevigata had the highest abundance of specific LTR-RTs among the five species (Fig. 3F).
Whole-genome duplication (WGD) events are core evolutionary forces driving plant speciation, gene functional innovation and adaptive evolution [
70]. Ks density distributions of paralogous gene pairs in
R. laevigata and other Rosaceae species showed two peaks at approximately 0.2 and 1.5. The younger peak likely reflects a recent WGD event shared by Rosaceae species (Fig. 3G), while the single peak at Ks ≈ 1.5 represents the ancient γ whole-genome triplication event [
15,
71]. Corrected Ks distributions of orthologous gene pairs further supported a divergence order consistent with the phylogenetic tree (Fig. 2A). Whole-genome synteny analysis with
Rosa rugosa,
Rubus idaeus,
Crataegus pinnatifida, and
Vitis vinifera showed that
R. laevigata had the highest level of collinearity with
R. rugosa, while collinearity decreased progressively with increasing phylogenetic distance, accompanied by more rearrangements, inversions, and translocations (Fig. 3H). These results suggest the high degree of genome conservation within the genus Rosa and the progressive structural divergence of
R. laevigata from more distantly related species.
3.3 Tissue-specific transcriptomic and metabolomic profiling of R. laevigata
The unique medicinal and biological functions of plant tissues are largely determined by their specific gene expression patterns and metabolite accumulation profiles [
72]. To systematically dissect the spatial transcriptional landscape of
R. laevigata, we performed high-throughput RNA sequencing (RNA-seq) on five representative tissues: roots, stems, leaves, seeds and fruits (Fig. 4). Differential expression analysis identified numerous DEGs among tissues (Figs. 4B and 4C). KEGG enrichment of tissue-specific upregulated genes revealed clear functional differentiation (Fig. 4D). Leaf-upregulated genes were mainly enriched in photosynthesis-related pathways, consistent with the primary role of leaves in carbon assimilation. Root-upregulated genes were enriched in pathways associated with environmental adaptation, transport, hormone signaling, and flavonoid biosynthesis, reflecting their roles in nutrient uptake, stress response, and metabolite transport. Fruit-upregulated genes were significantly enriched in secondary metabolite biosynthetic pathways, particularly flavonoid biosynthesis, phenylpropanoid biosynthesis, amino sugar and nucleotide sugar metabolism, and ascorbate and aldarate metabolism, indicating that fruits are major sites for the biosynthesis and accumulation of medicinally active compounds in
R. laevigata. Consistently, key structural genes in the flavonoid and phenylpropanoid pathways displayed strong tissue-specific expression patterns (Fig. 4E). Collectively, these results demonstrate that tissue-specific gene expression patterns in
R. laevigata are tightly linked to functional differentiation, and the transcriptional activation of secondary metabolite biosynthetic pathways in reproductive tissues, particularly fruits, provides the molecular basis for the accumulation of medicinal active ingredients in fruits.
To further investigate tissue-specific metabolic features, we performed UPLC–MS/MS-based metabolomic profiling of the same five tissues used for transcriptome analysis. Partial least squares discriminant analysis (PLS-DA) based on metabolomic data revealed clear metabolic separation among tissues and high clustering of biological replicates, indicating strong organ specificity in metabolite accumulation (Fig. 4F). Differentially accumulated metabolites (DAMs) were identified between tissues (Fig. 4G), and Venn analysis showed that 65.2% of metabolites were shared across all five tissues, whereas seven metabolites were fruit-specific (Fig. 4H). A stacked bar chart further showed clear differences in metabolite class composition among tissues, with flavonoids, lipids, and terpenoids together accounting for more than 40% of the detected metabolites. KEGG enrichment of upregulated DAMs was largely consistent with the transcriptomic results (Fig. 4J). In fruits, these metabolites were mainly enriched in flavonoid biosynthesis, phenylpropanoid biosynthesis, isoflavonoid biosynthesis, glycolysis/ gluconeogenesis, galactose metabolism, and unsaturated fatty acid biosynthesis. We further examined the accumulation of ten key bioactive metabolites with documented medicinal importance, including apigenidin, liquiritigenin, naringenin, luteolin, quercetin, kaempferol, catechin, corosolic acid, gallic acid, and rutaretin, all of which showed distinct tissue-specific accumulation patterns (Fig. 4K). Overall, integrated transcriptomic and metabolomic analyses revealed the tissue-specific regulatory and metabolic landscapes of R. laevigata, identified fruits as a major tissue for medicinal compound biosynthesis and accumulation, and provided a valuable set of candidate genes and metabolites for functional studies and genetic improvement.
3.4 Dynamic transcript-metabolome and metabolite changes during fruit ripening
Fruits are the primary medicinal organs of
R. laevigata, and their nutritional and medicinal value largely depends on the dynamic accumulation of bioactive metabolites during ripening [
10,
73,
74]. To characterize this process, fruits were collected at consecutive developmental stages from April (young ovary stage) to December (fully mature stage). Phenotypic observations showed continuous fruit development throughout ripening (Fig. 5A). Consistently, fresh weight, dry weight, transverse diameter, and longitudinal diameter all increased significantly over time (Figs. 5B–5E). Measurements of physiological and bioactive traits further showed that eight core nutritional and medicinal indicators, including total flavonoids, total phenolics, total saponins, oligomeric proanthocyanidins (OPCs), total tannins, ascorbic acid, ferric reducing antioxidant power (FRAP), and DPPH radical scavenging activity, increased progressively during fruit ripening (Figs. 5F–5M). Notably, the accumulation of these bioactive components was closely synchronized with fruit growth, with the fastest increase occurring during the color transition and maturation stages from September to December, indicating that fully mature fruits possess the highest nutritional and medicinal value, consistent with the traditional harvest period.
To investigate metabolic dynamics during fruit ripening, we performed untargeted metabolomic profiling of fruits from seven developmental stages using UPLC–MS/MS. Principal component analysis (PCA) revealed clear metabolic separation among stages (Fig. 6A). Analysis of differentially accumulated metabolites (DAMs) showed that the number of upregulated DAMs generally increased during ripening, whereas downregulated DAMs remained relatively low, suggesting that fruit ripening is mainly characterized by the continuous accumulation of functional metabolites (Fig. 6B). Metabolite classification identified flavonoids, prenol lipids, organic oxygen compounds, carboxylic acids and derivatives, and fatty acyls as the major metabolite classes. Notably, the relative abundances of flavonoids, carboxylic acids and derivatives, phenylpropanoids, and terpenoids increased steadily with ripening, in agreement with the observed physiological changes (Fig. 6E). Among metabolites with high variable importance in projection (VIP) values, the top-ranked compounds were mainly flavonoid derivatives, phenolic acids, triterpenoid saponins, and polysaccharide precursors, all of which showed marked accumulation at late ripening stages (Fig. 6F). Correlation-based network analysis further revealed strong associations among flavonoid, phenylpropanoid, and carbohydrate biosynthetic pathways (Fig. 6G). Time-series clustering divided DAMs into ten modules, of which M5, M6, and M7 showed accumulation patterns highly consistent with bioactive metabolite trends (Fig. 6H). KEGG enrichment analysis showed that these modules were significantly enriched in multiple secondary metabolic pathways related to medicinal activity, as well as primary metabolic pathways such as starch and sucrose metabolism that contribute to sugar and polysaccharide accumulation (Figs. 6I–6K).
To further uncover the transcriptional basis of these metabolic changes, we performed RNA-seq on the same fruit samples used for metabolomic analysis (Fig. 6C). Differential expression analysis showed that the number of DEGs also generally increased during fruit ripening (Fig. 6D). Time-series clustering identified modules G8, G9, and G10 as continuously upregulated during ripening, with expression patterns closely matching the accumulation of bioactive metabolites (Fig. 6L). KEGG enrichment analysis indicated that these genes were significantly enriched in flavonoid biosynthesis, phenylpropanoid biosynthesis, plant hormone signal transduction, starch and sucrose metabolism, and secondary metabolite biosynthesis. These enrichment patterns largely overlapped with those of DAMs, suggesting that transcriptional activation of these pathways drives the synthesis and accumulation of the corresponding metabolites (Figs. 6M–6O). Together, these results reveal the dynamic transcriptional and metabolic regulatory networks underlying R. laevigata fruit ripening and provide a foundation for functional studies of medicinal metabolite biosynthesis.
3.5 Transcriptional and metabolic landscape of the flavonoid biosynthetic pathway
Flavonoids are major bioactive constituents of
R. laevigata and are known to possess antioxidant, anti-inflammatory, and hepatoprotective activities [
75]. The upstream phenylpropanoid pathway provides the core precursor supply for flavonoid biosynthesis. Based on the T2T genome annotation, we systematically identified the key structural genes involved in the phenylpropanoid pathway in
R. laevigata, most of which displayed clear tissue-specific and developmental stage-specific expression patterns (Fig. 7). Among them, PAL, encoding the first rate-limiting enzyme of the phenylpropanoid pathway, showed relatively high expression in roots and mature fruits. Homologs of C4H and 4CL were continuously upregulated during the middle and late stages of fruit development, suggesting an enhanced precursor supply for downstream flavonoid biosynthesis.
(A–C) Expression heatmap of key structural genes involved in the phenylpropanoid biosynthesis pathway across five different tissues and seven fruit developmental stages. The color gradient represents the normalized expression level (Z-score) of each gene. (B) Proposed biosynthetic pathways for phenylpropanoids, terpenoids and flavonoids in R. laevigata. Three major metabolic branches are highlighted in different colored backgrounds. (D) Expression heatmap of key structural genes involved in flavonoids biosynthesis across five different tissues and seven fruit developmental stages. (E) Expression heatmap of key structural genes involved in terpenoid biosynthesis across five different tissues and seven fruit developmental stages. (F) Relative content heatmap of phenylpropanoid and flavonoid metabolites across five different tissues (left panel) and seven fruit developmental stages (right panel). Metabolites related to flavan-3-ol and proanthocyanidin biosynthesis are marked with red stars. The color gradient represents the normalized relative content (Z-score) of each metabolite.
In the core flavonoid pathway (Fig. 7B),
p-coumaroyl-CoA and malonyl-CoA are sequentially converted by CHS and CHI to produce naringenin, which is subsequently directed into different branches, including flavonol, anthocyanin, and flavan-3-ol biosynthesis, through the actions of F3H, F3'H, and DFR [
76]. Transcriptomic analysis showed that the core structural genes in this pathway exhibited coordinated expression patterns. In particular, CHS, CHI, and F3H, which function in the early steps of flavonoid biosynthesis, were significantly upregulated during the middle and late stages of fruit development, and their expression trends were consistent with the continuous increase in total flavonoid content in fruits (Fig. 7D). Notably, LAR (leucoanthocyanidin reductase), a key enzyme catalyzing the conversion of leucoanthocyanidins to flavan-3-ols, showed transcript abundance peaking at the late stage of fruit development. This pattern was highly consistent with flavonoid accumulation, suggesting that LAR plays an important role in determining the accumulation of major medicinally active flavan-3-ols in
R. laevigata fruits (Fig. 7D).
Based on correlation analysis between gene expression and metabolite accumulation, we constructed a transcriptional-metabolic co-regulatory network for the flavonoid biosynthetic pathway in R. laevigata (Fig. S12). Among the metabolites, afzelechin and (−)-epiafzelechin showed high connectivity in the network. Their accumulation increased significantly during fruit development and was closely associated with LAR expression, supporting the importance of the flavan-3-ol branch in fruit ripening (Fig. 7F).
In addition, we examined the terpenoid biosynthetic branch, which also derives from acetyl-CoA-related metabolic inputs, and identified key genes in this pathway, including 3-hydroxy-3-methylglutaryl-CoA reductase (HMGR), mevalonate kinase (MVK), phosphomevalonate kinase (PMVK), diphosphomevalonate decarboxylase (MVD), isopentenyl-diphosphate delta-isomerase (IDI), farnesyl diphosphate synthase (FDPS), squalene synthase (SQS), squalene epoxidase (SQE), lupeol synthase (LUP4) and cycloartenol synthase (CAS1). These genes showed dynamic expression changes during fruit development, indicating that the phenylpropanoid-flavonoid and terpenoid-saponin pathways are co-activated during R. laevigata fruit ripening (Fig. 7E). Together, these findings provide a molecular framework for understanding the coordinated biosynthesis of flavonoids and terpenoid-related active compounds during R. laevigata fruit ripening.
3.6 Identification and characterization of RlmMYB32 as a candidate regulator of flavan-3-ol biosynthesis in R. laevigata
The biosynthesis of flavonoid bioactive compounds in plants is a complex process regulated by hierarchical transcriptional networks. To identify key genes and core modules associated with flavan-3-ol accumulation in R. laevigata, we performed weighted gene co-expression network analysis (WGCNA) using transcriptome data from consecutive fruit developmental stages, with the accumulation levels of eight core flavonoid metabolites as phenotypic traits (Fig. 8). Module-trait relationship analysis identified three core modules that were highly positively correlated with target flavonoid metabolites (Fig. 8C). KEGG enrichment analysis showed that genes in these modules were significantly enriched in phenylpropanoid biosynthesis, flavonoid biosynthesis, amino sugar and nucleotide sugar metabolism, terpenoid backbone biosynthesis, and starch and sucrose metabolism (Fig. 8D). These results were consistent with our previous pathway analyses, further supporting the central roles of these modules in flavonoid biosynthesis and providing a reliable gene set for screening key regulators.
To further dissect the transcriptional regulation of flavan-3-ol biosynthesis in R. laevigata, we annotated, classified, and quantified transcription factors (TFs) in the three key modules (Figs. 8E and 8F). Among them, the MYB and MYB-related families were the most abundant, with 27 members in total. Based on Pearson correlation coefficients of gene expression profiles, we constructed a co-expression network linking TFs in the core modules with key structural genes in the flavonoid biosynthetic pathway (Fig. 8G). Network analysis showed that core structural genes were positively co-expressed with numerous TFs from the MYB, bHLH, ERF, NAC, and WRKY families. Notably, LAR, a key gene involved in flavan-3-ol biosynthesis, showed particularly strong co-expression with multiple MYB, WRKY, and MYB-related TFs, while CHS and CHI, the initiating genes of the flavonoid pathway, also formed major regulatory nodes. In addition, MYB32 exhibited the highest connectivity in the network and was identified as a candidate core regulatory transcription factor.
MYB transcription factors constitute one of the largest transcription factor families in plants and play central roles in the regulation of secondary metabolism [
77]. Based on the high-quality genome assembled in this study, we performed a genome-wide identification of the MYB gene family in
R. laevigata and identified 94 R2R3-MYB members (Fig. S13A, Table S16). Notably,
RlmMYB32, the candidate regulator highlighted in this study, clustered with the well-characterized flavonoid regulators
AtMYB12 (AT2G47460) and
AtMYB11 (AT3G62610) from
Arabidopsis thaliana [
78,
79], suggesting a conserved role in flavonoid biosynthesis.
We further performed comparative synteny analysis of key structural genes in the flavonoid biosynthetic pathway between R. laevigata and two related Rosaceae species, R. rugosa and R. idaeus (Fig. S13B). Core genes, including PAL, C4H, 4CL, CHS, CHI, F3H, LAR, and ANS, showed high collinearity among the three species, indicating strong conservation of the flavonoid biosynthetic pathway during Rosaceae evolution. To evaluate the regulatory potential of RlmMYB proteins, we analyzed cis-acting elements in the 2000-bp promoter regions upstream of 14 key structural genes (Fig. S13C). Most core genes contained MYB-binding sites, and the promoters of RlmLAR, RlmFLS, and RlmDFR each harbored three MBS motifs, suggesting that these genes may represent major targets of RlmMYB-mediated regulation.
Expression profiling across seven consecutive fruit developmental stages further divided RlmMYB genes into several distinct expression modules. Among them, RlmMYB32 showed low expression during the early stages of fruit development, followed by a marked increase from November onward and a peak at full maturity in December. This pattern was strongly positively correlated with that of its putative target gene RlmLAR, supporting a potential regulatory relationship between RlmMYB32 and RlmLAR (Fig. S13E). Molecular docking analysis suggested that RlmMYB32 can form a stable interaction with the MBS motif through multiple non-covalent interactions (Figs. S13D and S14). Together, these results reveal the evolutionary characteristics of the RlmMYB gene family, identify RlmMYB32 as a key candidate regulator of flavonoid accumulation during fruit ripening, and provide a basis for future functional validation and breeding of high-flavonoid R. laevigata cultivars.
4 Discussion
Previous genomic studies in Rosaceae have focused mainly on ornamental plants and fruit crops, leaving medicinal species underrepresented in genomic resources [
2,
5,
80–
82]. In this study, we generated the first T2T gap-free genome of
R. laevigata, providing a valuable foundation for future evolutionary and functional studies. Transposable elements (TEs) are major contributors to genome structural variation and functional innovation in plants. We found that 56.76% of the
R. laevigata genome is composed of repetitive sequences, with LTR-RTs representing the dominant fraction (Table 1). A recent burst of LTR-RT amplification occurred 0–2 Mya, consistent with patterns reported in many angiosperms [
83]. Notably, the Ty3/Gypsy superfamily showed a later insertion peak and a higher insertion density than Ty1/Copia, suggesting that recent Ty3/Gypsy amplification was a major contributor to genome expansion in
R. laevigata (Fig. 3). Gene family evolution analysis further supported this genomic trajectory, as significantly expanded gene families and species-specific genes were mainly enriched in pathways related to secondary metabolism and environmental adaptation.
Functional specialization of plant tissues is largely shaped by spatiotemporal gene expression and differential metabolite accumulation [
84]. Integrated transcriptomic and metabolomic analyses of five representative tissues showed that fruits are the major tissues for the biosynthesis and accumulation of medicinally active compounds in
R. laevigata, which agrees with the
Chinese Pharmacopoeia, where the mature fruit is recorded as the medicinal part [
85]. Fruit-upregulated genes were significantly enriched in flavonoid, phenylpropanoid, and ascorbate-related pathways, while metabolomic analysis identified seven fruit-specific metabolites and showed that core medicinal compounds, including catechin, quercetin, and kaempferol, accumulated most strongly in mature fruits (Fig. 4). These findings not only support the traditional medicinal use of
R. laevigata fruits at the molecular level, but also provide a basis for the standardized utilization and quality evaluation of its medicinal materials.
The dynamic accumulation of medicinal compounds during fruit development is directly related to the optimal harvest stage. By systematically characterizing fruits across seven consecutive developmental stages, we found that morphological development was closely synchronized with the accumulation of medicinally relevant metabolites (Fig. 5). Total flavonoids, total phenolics, and proanthocyanidins all increased progressively during ripening, with the most rapid accumulation occurring during the color transition and maturation stages from September to December (Fig. 5). These results indicate that fully mature fruits harvested in December likely represent the optimal stage for medicinal use. Dynamic transcriptomic and metabolomic analyses further suggested that the continuous accumulation of medicinal compounds during ripening is driven by coordinated activation of multiple biosynthetic pathways (Fig. 6). Both DEGs and DAMs were significantly enriched in flavonoid, phenylpropanoid, terpenoid, and starch and sucrose metabolism pathways, and time-series clustering identified several gene and metabolite modules whose patterns closely matched the accumulation of medicinal components. Notably, the phenylpropanoid-flavonoid pathway and the terpenoid-saponin pathway were co-upregulated during fruit development, suggesting partial coordination between these two major biosynthetic systems. This coordinated behavior may provide a useful framework for molecular breeding aimed at improving multiple classes of medicinal compounds simultaneously.
Flavan-3-ols are among the representative bioactive compounds in
R. laevigata, and their biosynthesis is regulated by both structural genes and transcription factors. Based on the T2T genome, we systematically annotated the key structural genes involved in this pathway and found that RlmLAR (leucoanthocyanidin reductase), a key gene in flavan-3-ol biosynthesis, showed high and fruit-biased expression during the middle and late stages of fruit development, with its expression peak coinciding with flavan-3-ol accumulation (Fig. 7). This result suggests that transcriptional activation of RlmLAR plays an important role in promoting flavan-3-ol accumulation in
R. laevigata fruits. In plants, flavonoid biosynthesis is commonly regulated by the MYB-bHLH-WD40 (MBW) complex, in which MYB transcription factors often confer pathway specificity [
86]. Through weighted gene co-expression network analysis (WGCNA), we identified RlmMYB32 as the hub transcription factor with the highest connectivity in the flavan-3-ol-associated network (Fig. 8). Phylogenetic analysis further showed that RlmMYB32 clustered with the well-characterized flavonoid regulators AtMYB12 and AtMYB11 from Arabidopsis thaliana, suggesting potential functional conservation [
78,
79]. Additional evidence supported a regulatory relationship between RlmMYB32 and RlmLAR. First, the promoter of RlmLAR contains three canonical MBSs, indicating potential cis-regulatory targets for MYB-mediated control. Second, the expression pattern of RlmMYB32 was strongly positively correlated with that of RlmLAR. Third, AlphaFold3-based molecular docking suggested that RlmMYB32 can form stable interactions with the MBS motifs in the RlmLAR promoter through multiple non-covalent interactions (Fig. S13). Collectively, these results support the hypothesis that RlmMYB32 is a key candidate regulator of flavan-3-ol biosynthesis in
R. laevigata fruits, likely acting through regulation of RlmLAR transcription. This regulatory module provides an important framework for future functional validation and for the molecular improvement of medicinal quality in
R. laevigata.
5 Conclusions
In conclusion, this study presents the first gap-free T2T reference genome of R. laevigata, a medicinally important Rosaceae species and a valuable bioactive forest biomass resource. By integrating PacBio HiFi reads, Oxford Nanopore ultra-long reads, and Hi-C chromatin interaction data, we assembled a 503.66 Mb genome comprising seven fully gapless pseudochromosomes, with a contig N50 of 67.24 Mb and 98.9% BUSCO completeness. This assembly resolves highly repetitive genomic regions, including centromeric regions, telomeric regions, and LTR retrotransposon-rich regions. Genome annotation showed that repetitive sequences account for 56.76% of the genome and identified 30,259 high-confidence protein-coding genes. Integrated spatiotemporal transcriptomic and metabolomic analyses across five tissues and seven fruit developmental stages confirmed that fruits are the major tissues for the accumulation of bioactive compounds and revealed a strong synchronization between fruit development and metabolite accumulation. Furthermore, we identified RlmMYB32 as a candidate key transcription factor potentially involved in flavan-3-ol biosynthesis, likely through the regulation of RlmLAR, a key gene in this pathway, as supported by promoter element analysis, co-expression patterns, and AlphaFold3-based molecular docking simulations. Overall, this T2T genome and multi-omics atlas provide a valuable resource for comparative genomics and evolutionary studies in Rosaceae, and establish a molecular foundation for functional validation, quality improvement, molecular breeding, and high-value utilization of R. laevigata as a bioactive biomass resource.
The Author(s) 2026. This article is published by Higher Education Press.
PDF
(32359KB)
