Abstract ‘Fumei’ apple is characterized by high anthocyanin content and thick wax layer. Long non-coding RNAs (lncRNAs) play essential roles in the growth and development of various plants via regulation of gene expression. This study explored the potential mechanism underlying anthocyanin accumulation and cuticular wax formation during the development of ‘Fumei’ apple fruit. The results demonstrated that anthocyanin accumulation correlates with fruit coloration, while wax content drives wax layer formation. A total of 6039 and 3410 differentially expressed genes (DEGs), as well as 230 and 131 differentially expressed lncRNAs (DELs) were identified in the M1/M2 and M2/M3 pairs, respectively, by using RNA-seq. In the M1/M2 pair, the DEGs were mainly enriched in the ‘photosynthesis’ and ‘flavonoid biosynthesis’ pathways; in the M2/M3 pair, the DEGs were significantly enriched in the ‘photosynthesis’ and ‘cutin, suberine and wax biosynthesis’ pathways. Furthermore, the structural and regulatory genes involved in anthocyanin and cuticular wax biosynthesis were investigated, and the potential lncRNAs and genes that may control the anthocyanin and cuticular wax biosynthesis were identified. This study provides candidate lncRNAs and potential regulatory genes associated with both the regulation of anthocyanins and wax during apple development. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-025-06545-3. Keywords: Apple, Anthocyanin, Wax, lncRNAs Background Apple is one of the most widely cultivated fruits in the world, particularly in China. The cultivation area and yield of apple in China exceed more than 50% of the global area and yield. Fruit quality is an important economic trait of apple, which directly affects the market value of apple. The fruits of most apple varieties are red, smooth and waxy as those traits are attractive to consumers [[40]1, [41]2]. Additionally, red apples are rich in anthocyanins, the major compounds in the flavonoid pathway, while waxy apples are covered with cuticular wax [[42]3–[43]5]. It has been demonstrated that the content of anthocyanins determines the red pigmentation of apples. Anthocyanins are biosynthesized on the cytosolic surface of the endoplasmic reticulum (ER) by anthocyanin biosynthetic enzymes, which are encoded by structural genes [[44]6]. They are transported into the vacuoles by the anthocyanin transporter proteins [[45]7–[46]9]. These structural genes can be divided into early biosynthetic genes (EBGs, such as CHS, CHI, and F3H) and late biosynthetic genes (LBGs, such as DFR, ANS, and UFGT) [[47]6, [48]10]. To date, three anthocyanin transporters (ATP-binding cassette (ABC), glutathione S-transferase (GST), and multidrug and toxic compound extrusion (MATE)) have been identified in plants [[49]9, [50]11, [51]12]. Anthocyanin biosynthesis is mainly controlled by regulatory genes, including MYBs, bHLHs and WD40 proteins, via generation of the MBW ternary complexes [[52]13]. For instance, MdMYB1 is a key regulator in apple as it dominates anthocyanin biosynthesis on the fruit skin and interacts with MdbHLH3/33 [[53]14, [54]15]. On the other hand, some negative regulators have been identified in apple (such as MdMYB6 and MdMYB16), strawberry (such as FvMYB1) and ginkgo (such as GbMYBF2) [[55]16–[56]19]. Furthermore, other transcription factors (TFs) can directly or indirectly regulate anthocyanin biosynthesis. For instance, PpBBX16 regulates light-induced anthocyanin biosynthesis by binding to the PpMYB10 promoter in pear [[57]20]; and MdEIL1 promotes anthocyanin biosynthesis by binding to the MdMYB1 promoter in apple [[58]21]. Covering the fruit surface, cuticular wax can effectively reduce post-harvest water loss and prevent pathogen invasion, thereby improving the appearance of fruits. In plants, cuticular wax biosynthesis generally comprises three steps: (1) de novo synthesis of C16 and C18 fatty acids [[59]22]; (2) extension of very long chain fat acids (VLCFAs) [[60]23, [61]24]; (3) formation of VLCFA derivatives, including alkanes, aldehydes, alcohols, ketones, and esters [[62]25]. Different enzyme-coding genes have been isolated for each step. For instance, fatty acyl-ACP thioesterases (FATA and FATB) trigger the hydrolysis of C16–C18 acyl-acyl carrier proteins (ACP) into C16–C18 fatty acids in plastids, and the products are then activated by long-chain acyl-CoA synthetases (LACS) and exported to the ER [[63]26]. The fatty acid elongase (FAE) complex is responsible for catalyzing the sequential reactions, including dehydration, reduction, condensation and second reduction, wherein KCS, KCR, ECR and HCD are the functioning enzymes [[64]23, [65]24]. Finally, the wax components synthesized on the ER are transported to the plasma membrane (PM) by the ABC transporters and LTPGs [[66]27, [67]28]. The wax biosynthesis is essentially regulated by transcriptional regulation in plants. For instance, the AP2-type TF WIN1/SHN1 promotes cuticular wax accumulation by regulating KCS1, CER1 and CER2 in Arabidopsis [[68]29, [69]30]. Three homologous genes (AtMYB96, MdMYB94 and CsMYB96) are involved in drought stress response and wax biosynthesis [[70]31–[71]33]. In apple, the MYB TF MdMYB106 can improve the wax content by interacting with the MdCER2 and MdCER2L1 promoters and undergo MdBT2-mediated degradation under nitrogen signals [[72]34]. The molecular mechanisms underlying anthocyanin accumulation and cuticular wax biosynthesis have been extensively studied in various plant species. However, the interplay between these two pathways remains largely unexplored. In this study, coloration and cuticular wax biosynthesis of apples were investigated based on the ‘Fumei’ apple. Specifically, transcriptome analysis was employed to identify DEGs and lncRNAs at different developmental stages of ‘Fumei’ apple. Functional analysis revealed that the DEGs were mainly enriched in the ‘Flavonoid’ and ‘Cutin, suberin and wax biosynthesis’ pathways. The cis-regulated target genes of DELs were significantly enriched in the ‘response to UV-B’ pathway and other stress-related pathways. Also, MdMYB94 was found to participate in both the biosynthesis of anthocyanins and wax. This study provides candidate genes and lncRNAs for the improvement on coloration and wax formation of apple fruits. Materials and methods Materials ‘Fumei’ apple fruit were obtained from a local orchard, which is a breeding base of Qingdao Agricultural University in Qingdao, China. A total of twenty ‘Fumei’ apple trees with a row spacing of 1.2 m × 3.5 m (4-year-old) were selected for fruit collection. The fruit were bagged on the 25th May 2022 (at 30 days after full bloom; DAFB) and the bags were removed at 160 DAFB. The fruits were collected at five stages including 0 days after bags removal (DABR), 6 DABR, 12 DABR, 18 DABR and 24 DABR. At each time 20 fruits of uniform size without mechanical damages were sampled. The peels were collected and stored at − 80 °C. Measurement of total content of anthocyanins The content of anthocyanins was measured by a method reported previously [[73]35]. First, 0.5 g of the fruit peel was collected and ground into a fine powder, followed by extraction in 10 mL methanol with 1% HCl for 24 h. Then, 1 mL of the extract was added into 4 mL of KCl buffer (0.025 M; pH = 1.0) and 4 mL of NaAc buffer (0.4 M; pH = 4.5), respectively. The two solutions were mixed and kept untouched at 4 °C for 15 min in the dark. Finally, the absorbance of the solution was measured at 510 and 700 nm, respectively. The content of anthocyanin was calculated by OD = (A[510]– A[700]) − 0.1 × (A[510]– A[700]). Characterization by scanning electron microscopy The fruit skin was cut into a 0.5 × 0.5 cm pieces and fixed in the FAA solution. Then, the pieces were dehydrated using acetone. After the critical-point dried and gold spraying, the samples were checked by SEM (JSM-7500 F). Extraction and GC-MS analysis of cuticular wax The fruit samples were immersed in chloroform to extract the cuticular wax. Then, the solution was concentrated using a rotary evaporator at 40 °C. After drying under N[2] flow, the waxes were isolated and stored at − 80 °C until use. The wax samples were dissolved in the chloroform and methanol solution (v: v = 9: 1) and mixed with 1 mg of n-tetracosane (Sigma Aldrich, China) as internal standard. After drying under N[2] flow, the samples were mixed with bis-N, O-(trimethylsilyl) trifluoroacetamide (BSTFA, Sigma Aldrich, China) at 70 °C for 40 min. After that, the samples were injected into the GC-MS (QP-2020, Shimadzu, Tokyo, Japan) for further analysis. The NIST 17 MS Library was employed to identify the detected compounds. Additionally, the content of cuticular wax was calculated by using the internal standard method. RNA isolation, library construction, and sequencing RNAprep Pure Plant Kit (Tiangen, Beijing, China) was used to isolate total RNA of the samples. Firstly, DNA and rRNA were removed from the total RNA samples. Then, the RNA was broken into 300–350 bp fragments and used as templates for synthesis of the first-strand cDNA by the random primer method. After removing the RNAs using the RNase H, the second-strand cDNA was synthesized. The double stranded cDNA was purified and end-repaired, followed by addition of the adapters. Finally, PCR was employed to construct a specific library. The Illumina NovaSeq6000 platform was employed for sequencing. Read mapping and transcriptome assembly The raw data were generated by the sequencing, and adapters and low-quality reads were eliminated to obtain clean reads. The clean data were aligned to the GDDH13 apple genome [[74]36] using HISAT2 software [[75]37]. The StringTie software was employed to reconstruct the single sample transcript, while the Fragments Per Kilobase of transcript per Million fragments mapped (FPKM) was calculated by using the StringTie [[76]38]. Additionally, the StringTie-Merge tool was used to combine transcripts from different samples. The raw sequencing data are publicly accessible under GSA No. CRA016451 ([77]https://ngdc.cncb.ac.cn/gsa/browse/CRA016451) at [78]https://ngdc.cncb.ac.cn/gsa. Identification of LncRNAs The FEELnc tool was used to identify lncRNAs, and the protein-encoding capability of the lncRNAs was predicted by using computational tools including CPC2, PLEK, CNCI, Pfam and HMMER. The lncRNAs candidates were screened based on the following criteria: transcript length ≥ 200 bp, with two or more exons, and FPKM ≥ 0.1. Differential expression analysis of LncRNAs and encoding genes The expression of lncRNAs and encoding genes was analyzed by the DESeq2 R package [[79]39]. DEGs were analyzed with the following parameters:|log[2]FC| ≥ 1 and p value ≤ 0.05. Then, the DEGs were investigated by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis with the enrichGO and clusterProfiler packages, respectively. Prediction of target genes for LncRNAs The predicted target genes were divided into cis- and trans-regulated genes. Herein, genes within 100 Kb upstream and downstream of lncRNAs were taken as potential cis-regulated target genes, while the genes whose expression is correlated with the lncRNAs were regarded as trans-regulated target genes. RT-PCR analysis Total RNA of the apple peels was extracted by using the RNAprep pure Plant Kit (Tiangen, Beijing, China) and cDNA was synthesized by using the PrimeScript™ 1st strand cDNA Synthesis Kit (Takara, Beijing). All primers were provided by Shanghai Sangon Biotech Co. in China (Table [80]S1). RT-PCR was performed with three biological and three technical replicates for each sample, and the results were analyzed by using the CFX96 system (BioRad, Hercules, CA, USA). The program was set as follows: 95 °C for 30 s, followed by 38 cycles of 95 °C for 10 s, and finally 58 °C for 30 s. MdACTIN was used as internal controls for gene expression. The transcription levels of genes were determined by a method reported previously [[81]40]. Statistical analysis Data are presented as mean ± standard deviation (SD) from three independent experiments. Significant differences between groups were determined by one-way ANOVA followed by Tukey’s post-hoc test (p < 0.05) using the SPSS 22 software (IBM Corp., Armonk, NY, USA). Results Content of anthocyanins during fruit coloration ‘Fumei’ apple was selected as the study subject due to its characteristic red appearance and thick wax layer when harvested late (Fig. [82]S1). In this study, we collected ‘Fumei’ fruit samples at five developmental stages (0, 6, 12, 18, 24 DABR). The results showed that after bag removal, the fruit color gradually turned red and the fruit appearance became shiny (Fig. [83]1A). Then, we measured the contents of anthocyanins in the fruit at the five developmental stages. The results showed that the content of anthocyanins gradually increased steadily until 18 DABR and decreased with time from 24 DABR (Fig. [84]1b). Fig. 1. [85]Fig. 1 [86]Open in a new tab Development pattern of ‘Fumei’ apple fruit. A, Fruit phenotypes of ‘Fumei’ at different development stages. B, Anthocyanin contents in fruit peels of ‘Fumei’. C, Total wax contents in fruit peels of ‘Fumei’. Error bar means ± SD, n ≥ 3 and different letters indicate significant differences as assessed by one-way ANOVA (with Tukey’s test) (p < 0.01) Changes in cuticular wax structure of ‘fumei’ fruit The cuticular wax structure on ‘Fumei’ apple was investigated after bag removal. The SEM images (250×) showed deposit of cuticular wax on the fruit surface after bag removal. Also, crystals and small plates were observed at 6 DABR and the plates grew larger (Fig. [87]2A). Moreover, we also checked the structure of cuticular wax by high-power microscopy field (2500×). The wax layer on the fruit surface grew larger and thicker with time, especially at 6 and 12 DABR, but no wax films were generated at this stage. From 18 DABR, wax films were generated, which were stacked into layers covering the fruit surface (Fig. [88]2A). Fig. 2. [89]Fig. 2 [90]Open in a new tab SEM images and the wax components of fruit peels of ‘Fumei’. A, SEM images of fruit peels of ‘Fumei’. The contents of alkane (B), alcohol (C) and wax ester (D) in fruit peels of ‘Fumei’. Error bar means ± SD, n ≥ 3 and different letters indicate significant differences as assessed by one-way ANOVA (with Tukey’s test) (p < 0.01) Cuticular wax content and composition on ‘fumei’ fruit The total cuticular wax content was measured. The results indicated that the content of cuticular wax increased from 18 DABR and reached the peak at 24 DABR (Fig. [91]1C). Interestingly, the total content of cuticular wax was high at the beginning. Furthermore, the composition of cuticular wax on ‘Fumei’ fruit was determined by gas chromatography-mass spectrometry (GC-MS). The results showed that alkane was the predominant wax component in ‘Fumei’ fruit. Specifically, the content of C14 alkane was the highest at 18 DABR; the contents of C21 and C24 alkane increased at the last two stages (Fig. [92]2B); and those of C11, C19 and C21 wax ester increased drastically at the final two stages (Fig. [93]2C). However, the alcohol content decreased with over time during the development of ‘Fumei’ fruit (Fig. [94]2D). Overall, the increased content of alkane contributed to the accumulation of cuticular wax in ‘Fumei’ fruit. RNA-seq analysis of different developmental stages of ‘fumei’ fruit To explore the DEGs of the anthocyanin and cuticular wax pathways in ‘Fumei’ fruit peels, the samples were performed by RNA sequencing at different developmental stages. Herein, the content of anthocyanins in the fruit increased drastically at 12 DAFR, and that of cuticular wax was maximized at 24 DAFR (Fig. [95]1A). Thus, 3 stages (0, 12,24 DAFR marked as M1, M2, M3) with three replicates for each stage were used for RNA-Seq analysis. Transcript abundance with at least two-fold changes and p < 0.05 denoted DEGs. A total of 6039 (2857 upregulated and 1389 downregulated) and 3410 (1389 upregulated and 2021 downregulated) DEGs were identified in the M1/M2 and M2/M3 pairs, respectively (Fig. [96]S2 A, B, C). GO analysis classified DEGs into three major categories: biological process, cellular component and molecular function, respectively. The results indicated that most genes related to the biological process were involved in the ‘response to chitin’, ‘photosynthesis’ and ‘photosynthesis, light reaction’ pathways. Many genes corresponding to the cellular component were related to the ‘plastid thylakoid membrane’, ‘chloroplast thylakoid membrane’ and ‘photosynthetic membrane’ pathways. Genes corresponding to molecular function were related to ‘DNA-binding transcription factor activity, RNA polymerase II-specific’ and ‘monooxygenase activity’ pathways (Fig. [97]S3 A, B). The KEGG pathway enrichment analysis was also carried out. The results demonstrated that the DEGs were mainly enriched in ‘photosynthesis’ and ‘flavonoid biosynthesis’ pathways in the M1/M2 pair, while in the ‘photosynthesis’ and ‘cutin, suberine and wax biosynthesis’ pathways in the M2/M3 pair (Fig. [98]S3 C, D). Additionally, we investigated the co-expressed genes in the M1/M2 and M2/M3 pairs. A total of 2122 DEGs were co-expressed in both M1/M2 and M2/M2 pairs (Fig. [99]3A). GO analysis showed that these genes were mainly involved in biological processes and cellular component (Fig. [100]3B). For instance, most genes corresponding to biological processes were related to photosynthesis and responsive to light, suggesting that photosynthesis is involved in the development of ‘Fumei’ peels. Besides, some genes were responsive to plant hormones (such as jasmonic acid and ethylene). DEGs were enriched in the ‘response to the fatty acid’ pathway, indicating that the formation of cuticular wax on the ‘Fumei’ fruit is related to these genes (Fig. [101]3B). The KEGG pathway enrichment analysis showed that DEGs were mainly enriched in the ‘circadian rhythm’, ‘carbon fixation in photosynthetic organism’, ‘photosynthesis’ and ‘carotenoid biosynthesis’ pathways (Fig. [102]3C). Fig. 3. [103]Fig. 3 [104]Open in a new tab Differentially expressed genes (DEGs) and GO and KEGG analyses. A, Venn diagram indicating the number of DEGs in the M1/M2 and M2/M3 pairs. B, GO enrichment analysis of DEGs in the M1/M2 and M2/M3 pairs. C, KEGG pathway analysis of DEGs in the M1/M2 and M2/M3 pairs DEGs related to the anthocyanin pathway in ‘fumei’ fruit To further investigate anthocyanin biosynthesis, DEGs related to biosynthesis and transport pathways were analyzed (Fig. [105]4A). As demonstrated, eight CHS genes, seven CHI genes, three DFR genes, one ANS gene and five UFGT genes were differentially expressed in the M1/M2 pair, while one CHS gene, two CHI genes and two UFGT genes were differentially expressed in the M2/M3 pair (Fig. [106]4B). Besides the structural genes, the transporter coding genes were also analyzed. The results suggested that 25 GST gens, 11 MATE genes and six ABC were differentially expressed in the M1/M2 pair, while 15 GST genes, seven MATE genes and four ABC genes were differentially expressed in the M2/M3 pair (Fig. [107]4B). In the M1/M2 pair, four CHS genes (MD04G1003300, MD13G1285100, MD04G1003400, MD04G1003000), three CHI genes (MD01G1167300, MD07G1186300, MD07G1233400), two DFR genes (MD08G1028200, MD15G1024100), two UGFT genes (MD07G1306900, MD01G1234400), one GST gene (MD17G1272100), one MATE gene (MD07G1277300) and one ABC gene (MD11G1313800) were highly expressed at Stage 3, indicating that the expression of these genes contributes to anthocyanin accumulation. In other words, the anthocyanin biosynthetic genes were up-regulated after bag removal. Fig. 4. [108]Fig. 4 [109]Open in a new tab Expression patterns of genes related to biosynthesis and transport in the anthocyanin pathway in the M1/M2 and M2/M3 pairs. A, A model of anthocyanin biosynthetic pathway. B, Expression patterns of genes involved in the anthocyanin pathway in the M1/M2 and M2/M3 pairs. Value of log[2] (fold change) was used in the heatmap. Red colour represents up-regulation, and blue colour represents down-regulation DEGs related to cuticular wax in ‘fumei’ fruit DEGs involved in cuticular wax biosynthesis, including the synthesis and transport path, were also explored (Fig. [110]5A). The results showed that two LACS genes, ten KCS genes, one HCD gene, one ECR gene, two CER4 genes, three WSD1 genes, four CER1 genes, seven CER3 genes, and two CYPB5 genes were differentially expressed in the M1/M2 pair, while two FATB genes, three KCS genes, one HCD gene, one ECR gene, two CER4 genes, and three CER3 genes were differentially expressed in the M2/M3 pair (Fig. [111]5B). Also, two transporters, including four LTPG genes and 27 ABCG genes, were differentially expressed in the M1/M2 pair, while two LTPG genes and 24 ABCG genes showed differential expression in the M2/M3 pair. Notably, one ECR gene (MD02G1261000) and two CER4 genes (MD03G1141400, MD11G1221200) were up-regulated at 12 DABR but down-regulated at 24 DBAR (Fig. [112]5B). Fig. 5. [113]Fig. 5 [114]Open in a new tab Expression patterns of genes related to cuticular wax biosynthesis in the M1/M2 and M2/M3 pairs. A, A model of cuticular wax biosynthetic pathway. B, Expression patterns of genes involved in the cuticular wax pathway in the M1/M2 and M2/M3 pairs. Value of log[2] (fold change) was used in the heatmap. Red colour represents up-regulation, and blue colour represents down-regulation Differentially expressed TFs involved in anthocyanin and cuticular wax accumulation Anthocyanin biosynthesis is typically controlled by regulatory genes in plants. Hence, the regulatory genes involved in the anthocyanin biosynthetic pathway were investigated (Table [115]1). As indicated, MYB1 (MD09G1278600), a key regulatory factor, was differentially expressed in both M1/M2 and M2/M3 pairs, and was up-regulated at 12 DABR and down-regulated at 24 DABR. Moreover, most MYB genes were differentially expressed in the M1/M2 pair and down-regulated in the development process. The bHLH3 gene (MD11G1286900) exhibited a similar expression pattern to the MYB1 gene, and interacted with MYB1 to generate a MBW complex in apple. MYB9 was differentially expressed in the M2/M3 pair, indicating that MYB9 might regulate anthocyanin biosynthesis at the late stage of fruit development. Two MYB_SH[AL]QKY[RF] transcription factors, MdLUX (MD03G1086800) and MdPCL-like (MD11G1095900), exhibited opposite expression patterns and were down-regulated and then up-regulated (Table [116]1). Table 1. Transcription factors related DEGs that were involved in anthocyanin pathway Name ID log[2] foldchange in M1/M2 log[2] foldchange in M2/M3 MYB1 [[117]2] MD09G1278600 4.82 -1.52 MYB3 [[118]42] MD15G1411200 2.00 -4.04 MYB9 [[119]43] MD08G1070700 nd 3.38 MYB11 [[120]43] MD09G1184000 2.14 nd MYB90-like [[121]44] MD09G1278400 3.93 nd MYB114 [[122]35] MD17G1261100 5.89 nd MYB306-like [[123]41] MD17G1050900 1.15 -1.39 bHLH3 [[124]45] MD11G1286900 1.33 -1.14 LUX [[125]46] MD03G1086800 -2.87 2.43 PCL-like [[126]46] MD11G1095900 -1.89 1.59 [127]Open in a new tab Note: ‘nd’ means ‘not detected’ and representatives that the gene is not differently expressed in the comparison Regulatory genes involved in the cuticular wax biosynthetic pathway were also investigated (Table [128]2). As indicated, four MYB genes and one AP2 TF were differentially expressed. The MYB16 gene (MD17G1051700) was differentially expressed in the M2/M3 pair, while the MYB41 gene (MD09G1097700) was differentially expressed in the M1/M2 pair. MYB94 (MD17G1050900), MYB96 (MD14G1200100) and WRI4 (MD17G1226700) exhibited consistent expression patterns, all of which were differentially expressed in both M1/M2 and M2/M3 pairs; moreover, they were up-regulated at 12 DABR and down-regulated at 0 DABR and 24 DABR. Notably, MYB306-like genes and MYB94 genes belonged to the same ID (MD17G1050900), suggesting that they may simultaneously regulate the anthocyanin and cuticular wax biosynthesis in apples [[129]32, [130]41]. Table 2. Transcription factors related DEGs that were involved in cuticular wax pathway Name ID log[2] foldchange in M1/M2 log[2] foldchange in M2/M3 MYB16 [[131]17] MD17G1051700 nd -1.31 MYB41 [[132]47] MD09G1097700 2.80 nd MYB94 [[133]32] MD17G1050900 1.15 -1.39 MYB96 [[134]48] MD14G1200100 2.17 -1.39 WRI4 [[135]49] MD17G1226700 2.40 -1.50 [136]Open in a new tab Note: ‘nd’ means ‘not detected’ and representatives that the gene is not differently expressed in the comparison Identification of LncRNAs and DEGs in ‘fumei’ fruit To identify all potential lncRNAs in ‘Fumei’ apple, 4080 lncRNAs were included in this study (Table [137]S2). Compared with mRNAs, lncRNAs had shorter sequences, and most of them had sequence lengths less than 2500 bp (Fig. [138]S3A). Moreover, most lncRNA transcripts had at least two exons, which were fewer than those of protein-coding mRNAs (Fig. [139]S3B). Notably, the overall expression level of lncRNAs was lower than that of mRNAs (Fig. [140]S3C). According to the transcriptome data, 230 and 131 lncRNAs were differentially expressed in the M1/M2 and M2/M3 pairs, respectively (Fig. [141]6A). Specifically, 121 lncRNAs were up-regulated and 99 lncRNAs were down-regulated at 12 DABR compared with 0 DABR (Fig. [142]S3A, C). In addition, 66 lncRNAs were up-regulated and 65 lncRNAs were down-regulated at 24 DABR compared with 12 DABR (Fig. [143]S3B, C). Meanwhile, a total of 80 lncRNAs were differentially expressed in both pairs (Fig. [144]6A). It has been demonstrated that LncRNAs can regulate the expression of coding genes via cis and trans patterns. In this study, we also investigated the cis target genes. The results showed 178 and 80 lncRNA-targeted DEGs in the M1/M2 and M2/M3 pairs, respectively (Table [145]S3 and S4). GO analysis showed that most of these target genes were involved in the ‘regulation of meristem development’, ‘response to UV-B’, ‘microtubule cytoskeleton’, and ‘calmodulin binding’ pathways in the M1/M2 pair (Fig. [146]6B). KEGG pathway enrichment analysis showed that these target genes were mainly enriched in the ‘Trypophan metabolism’, ‘NOD-like receptor signaling pathway’ and ‘Oxidative phosphorylation’ pathways (Fig. [147]6C). For the M2/M3 pair, GO analysis indicated that most target genes were involved in different processes related to cell death and extra-membrane component (Fig. [148]6D); KEGG pathway enrichment analysis showed that most target genes were enriched in the ‘RNA polymerase’, ‘Pyrimidine metabolism’, ‘Purine metabolism’ and ‘Nucleocytoplasmic transport’ pathways (Fig. [149]6E). Overall, most target genes were involved in amino acid synthesis and disease-resistant processes. Fig. 6. [150]Fig. 6 [151]Open in a new tab Numbers of DELs and GO and KEGG analyses of their target genes in the M1/M2 and M2/M3 pairs. A, Venn diagram indicating the number of DELs in the M1/M2 and M2/M3 pairs. B, GO enrichment analysis of target genes of DELs in the M1/M2 pair. C, KEGG pathway analysis of target genes of DELs in the M1/M2 pair. D, GO enrichment analysis of target genes of DELs in the M2/M3 pair. E, KEGG pathway analysis of target genes of DELs in the M2/M3 pair To clarify whether the DELs target the anthocyanin and cuticular wax biosynthetic genes, we investigated all targeted genes in the M1/M2 and M2/M3 pairs. The results showed that no genes were present in both anthocyanin and cuticular wax biosynthesis pathways. However, MdMYB1 (MD09G1278600), a key regulatory gene in red apple, was targeted by the lncRNA MSTRG.27523.1 (Table [152]S3), indicating that MSTRG.27523.1 may regulate anthocyanin biosynthesis via affecting MdMYB1 in red apples. In addition, some target genes were involved in biotic and abiotic stress responses in apple. For instance, MdERF114 (MD16G1140800) interacts with MdMYB8 to improve resistance to F. solani, resulting in enhanced binding activity to the MdPRX63 promoter in apple [[153]50]. Two NLR coding genes (MD02G1112400 and MD02G1112800) were also present in Table [154]S3, suggesting that the biosynthesis of anthocyanin and cuticular wax at the late stage helps to improve plant resistance. Validation of gene expression To validate the gene expression results obtained from the RNA-Seq, 15 mRNAs and 5 lncRNAs were selected for qRT-PCR verification. MD02G1112400 and MD02G1112500 were downregulated at 12 DABR and upregulated at 24 DABR, whereas their corresponding lncRNAs exhibited the opposite pattern. MD09G1278600 and MD16G1140800 showed expression patterns that closely matched those of their associated lncRNAs (Fig. [155]7). Despite some quantitative differences in expression levels, the trends of the expression levels of genes were similar in both the RNA-Seq and qRT-PCR data. These results confirm the reliability of the RNA-Seq data. Fig. 7. [156]Fig. 7 [157]Open in a new tab Validation of the expression of genes and their lncRNAs at different developmental stages. Error bar means ± SD, n ≥ 3 Discussion As a key fruit quality, the color of fruit significantly affects its economic value of apple. ‘Fumei’ is an apple variety with a fully-red phenotype and high wax content. In this study, we investigated the key regulator factors related to the high contents of anthocyanin and wax, and identified the lncRNA-mRNA regulator modules involved in these pathways. Firstly, the contents of anthocyanin and cuticular wax in ‘Fumei’ were monitored during fruit development. Similar to most red apple varieties [[158]2, [159]51], ‘Fumei’ showed high anthocyanin content at the late developmental stage (Fig. [160]1B). Anthocyanin is an important flavonoid determining the phenotype of red apples and can be readily detected. Indeed, high content of anthocyanin can enhance the appearance quality of fruits such as apple. Cuticular wax can improve the fruit glossiness and affect the sensory quality, conferring the fruit with higher economic value. ‘Fumei’ exhibited an attractive appearance due to their high anthocyanin and wax accumulation (Fig. [161]S1). The composition of wax was different at different stages (Fig. [162]2B, C, D) due to its tissue or organ specificity. For instance, the leaves and stems of Camelina sativa have high contents of wax esters and triterpenoids/sterols, respectively, while its seed coats have high levels of primary alcohols [[163]52]. In ‘Golden Delicious’ and ‘Red Delicious’ apples, alkane is the dominant wax component, which is consistent with our findings in ‘Fumei’ [[164]53]. Anthocyanin is biosynthesized by structural genes on ER and transported to the vacuoles for storage. In this study, we analyzed the RNA-seq data at different developmental stages (0, 12, 24 DABR). The results showed that most DEGs were enriched in pathways related to ‘photosynthesis and response to light’ (Fig. [165]S3A, B), indicating that light induces anthocyanin biosynthesis after bag removal. Previous studies have demonstrated that blue light can induce the expression of anthocyanin biosynthetic genes and anthocyanin accumulation in pear and strawberry [[166]54, [167]55]. Moreover, in non-red apple cultivars, light-induced anthocyanin accumulation has been associated with epigenetic modifications, such as DNA methylation and histone modification [[168]56, [169]57]. Overall, these findings collectively highlight the close relationship between light and anthocyanin biosynthesis in plants. To further investigate the regulatory mechanisms underlying anthocyanin biosynthesis, we identified seven transcription factors (TFs), including MYBs and bHLHs, that may play key regulatory roles in this pathway (Table [170]1). Unlike anthocyanins, cuticular wax was present from 0 DAFB, and its biosynthesis appears to be independent of light (Fig. [171]1C). However, previous studies have suggested that specific light wavelengths influence the accumulation and composition of cuticular wax in faba bean [[172]58]. The cuticular wax in plants can prevent water loss and damages caused by biotic and abiotic stresses, including pathogens, pest, and UV radiation [[173]59–[174]61]. lncRNAs are ncRNAs with sequence lengths longer than 200 nt, and play key roles in flavonoid accumulation, fruit ripening, reproductive tissue development and responses to biotic and abiotic stresses in plants. For instance, COLDWRAP represses FLC expression by regulating H3K27m3 of the FLC promoter and develops an intragenic chromatin loop upstream of the FLC gene [[175]62]. MdLNC499 and MdLNC610 were identified as a cis-regulator in light-induced anthocyanin accumulation [[176]63, [177]64]. BplncSIR1 affects salt resistance of Betula platyphylla by regulating NpNAC2 to induce ROS scavenging and stomatal movement [[178]65]. In this study, we analyzed DELs and their target genes at different developmental stages. A total of 230 and 131 DELs were identified in the M1/M2 and M2/M3 comparisons, respectively (Fig. [179]6A). Notably, MdMYB1, a key regulator of anthocyanin biosynthesis, was predicted to be regulated by MSTRG.27523.1 (Table [180]S3). Additionally, a number of transcription factors were identified as putative targets of lncRNAs (Table [181]S3); however, their specific roles in anthocyanin and cuticular wax biosynthesis remain unclear. Interestingly, no structural genes associated with either anthocyanin or cuticular wax biosynthetic pathways were identified among the predicted lncRNA targets. Overall, this study provides potential candidate lncRNAs and TFs that may regulate anthocyanin and cuticular wax biosynthesis during apple fruit development. Conclusions In this study, we identified key lncRNAs and genes involved in apple coloration and wax formation through transcriptome analysis. The findings revealed that multiple protein-coding genes contribute to anthocyanin and wax biosynthesis, and that MdMYB1, a central regulator of anthocyanin biosynthesis, may be regulated by MSTRG.27523.1. Furthermore, we identified several potential lncRNA-mRNA regulatory interactions that may provide new insights into the molecular mechanisms governing anthocyanin and cuticular wax biosynthesis during apple fruit development. Electronic supplementary material Below is the link to the electronic supplementary material. [182]12870_2025_6545_MOESM1_ESM.jpg^ (162.5KB, jpg) Supplementary Material 1: Figure S1. Genealogy and phenotype of ‘Fumei’ apple [183]12870_2025_6545_MOESM2_ESM.tif^ (18.2MB, tif) Supplementary Material 2: Figure S2. Numbers of DEGs identified in the M1/M2 and M2/M3 pairs. Volcano maps showing the upregulated and downregulated DEGs in the M1/M2 and M2/M3 pairs (A, B). Column charts showing numbers of upregulated and downregulated DEGs in the M1/M2 and M2/M3 pairs [184]12870_2025_6545_MOESM3_ESM.tif^ (19.9MB, tif) Supplementary Material 3: Figure S3. GO enrichment and KEGG analysis of DEGs in the M1/M2 and M2/M3 pairs. A, GO enrichment analysis of DEGs in the M1/M2 pair. B, KEGG pathway analysis of DEGs in the M1/M2 pair. C, GO enrichment analysis of DEGs in the M2/M3 pair. D, KEGG pathway analysis of DEGs in the M2/M3 pair [185]12870_2025_6545_MOESM4_ESM.tif^ (21.3MB, tif) Supplementary Material 4: Figure S4. Characteristics of lncRNAs identified in ‘Fumei’ apple. A, Length statistics of lncRNA and mRNA. B, Numbers of exons in lncRNAs and mRNAs. C, Average expression level of lncRNAs and mRNAs [186]12870_2025_6545_MOESM5_ESM.tif^ (187.3MB, tif) Supplementary Material 5: Figure S5. Numbers of DELs identified in the M1/M2 and M2/M3 pairs. Volcano maps showing the upregulated and downregulated DELs in the M1/M2 and M2/M3 pairs (A, B). Column charts showing numbers of upregulated and downregulated DELs in the M1/M2 and M2/M3 pairs [187]12870_2025_6545_MOESM6_ESM.xlsx^ (99.7KB, xlsx) Supplementary Material 6: Table S1. Primers used in this study. Table S2. IDs of lncRNAs identified in this study. Table S3. Details of lncRNA-targeted DEGs in M1/M2 pair. Table S4. Details of lncRNA-targeted DEGs in M2/M3 pair Acknowledgements