The transcriptome of the two types of peel of C. limon [smooth (C) and rough (D)] harvested 30, 60, 90, 120, and 150 days after bloom (Figure 1) was sequenced using a Novaseq (Illumina, San Diego, CA, United States) platform. The sequencing of 27 libraries resulted in 40,211,564 to 47,275,744 and 39,568,768 to 46,284,626 raw and clean reads, respectively. The Q20, Q30, and GC% was >98, >94, and >44%, respectively (Supplementary Table 1). The gene expression level (FPKM) of fruit peels harvested after 30 days was high both for C and D (Supplementary Figure 1A). The Pearson’s correlation coefficient between the samples ranged from.632 to.98 (Supplementary Figure 1B). Principal component analysis (PCA) showed that the treatment replicates grouped together, indicating the reliability of sampling (Figure 2A).

A total of 62, 16, 17, and 19 genes were exclusively expressed in C peels harvested after 60, 90, 120, and 150 days, respectively. The genes that were expressed in C2 were related to flavonoid biosynthesis pathway, flowering-related pathway, carbohydrate metabolism, fruit ripening (CYP71A1), calcium homeostasis, and heat stress. The genes expressed in C3 were related to reactive oxygen species production, auxin signaling, cell number regulation, ethylene signaling, amino acid biosynthesis, and stomata differentiation.

Similar to C2 and C3, we noticed the expression of ethylene related genes (e.g., 1-aminocyclopropane-1-carboxylate oxidase homolog 1), flavonoid biosynthesis pathway, carotenoid biosynthesis (geranylgeranyl pyrophosphate synthase 7 chloroplastic-like, higher GGP biosynthesis means higher carotenoid biosynthesis), cuticle development (probable lysophospholipase BODYGUARD 3), fruit ripening (protein PELPK1, (Rashid, 2014)), and organ abscission [TF RAX3 (Pei et al., 2016)] in C4. The exclusively expressed genes in C5 as compared with D4 were mostly related to cell-wall, flavonoid biosynthesis, stomata closure related genes, and anthocyanin biosynthesis related genes and TFs, e.g., myb-related protein 308 (Supplementary Table 3).

The DEGs in C2 vs. D1 were significantly enriched in three pathways i.e., MAPK signaling pathway-plant, DNA replication, and glyoxylate and dicarboxylate metabolism. In addition to MAPK signaling pathway-plant, the DEGs in C3 vs D2 were enriched in amino sugar and nucleotide sugar metabolism, pentose and glucuronate interconversion, phenylpropanoid biosynthesis, and ubiquinone and other terpenoid-quinone biosynthesis pathways. With the increase in the number of days after flowering, i.e., 120 days after flowering (DAF), the DEGs (in C3 vs. D4) were enriched in the terpenoid backbone biosynthesis, phenylpropanoid biosynthesis, amino sugar and nucleotide sugar metabolism, flavonoid biosynthesis, stilbenoid, diarylheptanoid and gingerol biosynthesis, and phenylalanine, tyrosine, and tryptophan biosynthesis pathways. The DEGs in C4 vs. D5 were significantly enriched in the ribosome pathway (Supplementary Figure 2). The enrichment results suggest the changes of hormonal and developmental signaling, sugar and amino acid contents, and pigment and scent compounds during lemon peel development at 60 and 90 days. In the latter stages of lemon peel development, the biosynthesis of phenylpropanoids, terpenoids, phenylalanine, tyrosine and tryptophan, and amino sugars was differentially regulated. Thus, below, we present pathway-specific transcriptomic changes in C and D to explore more about key changes in both peel types.

The quantitative expression of 10 C. limon genes in D1 and C2, and D4 and C5 was confirmed by quantitative real-time (qRT)-PCR analysis. The expression profiles were similar to those observed in the RNA-sequencing (Figure 4 and Supplementary Tables 2–4). Thus, these expression profiles confirm that the transcriptome sequencing results are valid and reliable.

The PCA based on metabolites showed grouping of treatments of the same experimental group, indicating that the results are reliable (Figure 5A). We identified a total of 1,083 metabolites in the C. limon fruit peels (Figure 5B and Supplementary Table 5). For treatment comparisons D1 vs. C2, D2 vs. C3, D3 vs. C4, and D4 vs. C5, there were 331, 314, 332, and 245 metabolites that were significantly differentially accumulated (DAMs), respectively (Figure 5C). Thirteen common metabolites were accumulated in both peels collected after all the sampling dates (Figure 5D).

The DAMs in the treatment comparisons were enriched in multiple pathways (Supplementary Figure 3). Particularly, we observed that the DAMs in D2 vs. C1 were enriched in 36 pathways such as those observed in transcriptome-based KEGG pathway enrichment. Interestingly, we observed the enrichment of DAMs in the flavonoid biosynthesis, plant hormone signal transduction, anthocyanin biosynthesis, and amino acid related pathways. The DAMs in D2 vs. C3, D3 vs. C4, and D4 vs. C5 were enriched in 33, 20, and 25 pathways, respectively. The results of the transcriptome-based KEGG pathway enrichment analysis showed that the DAMs were enriched in some of the pathways (Supplementary Figures 2, 3). These observations confirm the transcriptome-based results that the C. limon peels harvested after 60 and 90 days in spring and autumn have a differential regulation of hormone signaling, amino sugar and amino acid-related pathways, and pigment related pathways. While in the latter DAF, the C. limon peels have a differential regulation of the phenylpropanoid and terpenoid, and tyrosine and tryptophan biosynthesis pathways.

A total of 331 metabolites were differentially accumulated in D1 vs C2; 197 were downregulated and 134 were upregulated in D1 as compared with C2. The downregulated metabolites belong to fatty acids and conjugates, flavones and flavonols, flavanons, amino acids, peptides and analogs, benzoic acid and derivatives, carbohydrates, and terpenoids (Supplementary Table 5). The metabolites that were increasingly accumulated in D1 as compared with C2 mostly belonged to flavonoids and flavonols, and amino acids, peptides, and analogs. In the case of the treatment comparison D2 vs. C3, 314 metabolites accumulated differentially; the content of 181 and 133 metabolites was decreased and increased, respectively, in D2 as compared with C3. The metabolites with decreased content in D2 were classified as fatty acids and conjugates, fatty alcohols, fatty amides, fatty esters, glycerophosphocholines, flavonoids, isoprenoids, and steroids. Additionally, we found lignans, terpene lactones, amino acids, peptides and analogs, sesquiterpenoids, and monoterpenoids in reduced quantities in D2 as compared with C3. On the contrary, flavonoids, sterols, isoprenoids, and several fatty acids and conjugates were increasingly accumulated in D2 as compared with C3. The differentially accumulated metabolites in D3 vs. C4 were classified as alcohols and polyols, benzene and substitutive derivatives, carboxylic acids and derivatives, coumarins, fatty acyls, flavonoids, carbohydrates, and terpenoids and derivatives. A similar class of compounds was differentially accumulated in the treatment comparison of the C. limon peels collected 150 DAF i.e., D4 vs. C5 (Supplementary Table 5).

An important aspect of differential fruit development in D and C could be the temperature at the time of pollen and ovary formation, flowering, and early fruit formation in watermelon (Lin, 2007; Wei et al., 2015), tomato (Yang et al., 2009), and cucumber (Xiaoyan and Zhilong, 2008). Since the early growth period of the D fruits was in low temperatures (see material and methods section), therefore, the differences in fruit (particularly peel) development could be due to lower temperature. This assumption is based on earlier observation in longan where the fruits experiencing lower temperatures resulted in larger fruits due to its effects on cell division (Yang et al., 2010a,b). In this regard, earlier studies have shown the effect of cold on citrus leaves and different plant parts or plants as a whole have added much knowledge (Jian-huaa and Su-jinb, 2010; Jiang et al., 2016, 2021). However, the effect of cold on the differential peel formation in lemon is an open question for future studies. Specific experiments where lemon trees grown under controlled conditions could delineate the specific role(s) of temperature on the peel growth and variation. Not only, the cold temperature but the overall growth conditions (light, water, and nutrition, etc.) can also regulate the plant growth differently and might exert comprehensive influence on the regulation of genes and the accumulation of metabolites (Seymour et al., 2013). Thus, in further studies on the differences in lemon peels grown and harvested in different seasons, the other abiotic and biotic factors will help to further enhance our understanding of the mechanism of different peel traits.

“Yunning No.1” lemon is mainly cultivated in tropic areas of Yunnan province, with sporadic cultivation in Sichuan, Guangxi, Guangdong, and Hainan provinces of China. Six-year-old grafted plants (disease- and insect-free) of Citrus limon (L.) Burm. f. “Yunning No.1” on Poncirus trifoliata L., growing at the Lujiangba Experimental base of the Institute of Tropical and Subtropical Economic Crops, Yunnan Academy of Agricultural Sciences (21°59′N, 98°53′E) were selected for fruit peel sampling. The experimental station lies in a subtropical dry-hot climate zone with an average annual temperature of 21.3°C (10-year average), average annual sunshine of 2,329.7 h, average radiation of 138,449 cal/(year.cm²), and an average rainfall of 755.3 mm. The quality of soil in the location is sandy loam with medium fertility. The experimental group consisted of lemon fruits that flowered in winter and ripened in spring. The samples were collected 60, 90, 120, and 150 days after flowering. The control group consisted of lemon fruits that flowered in spring and ripened in autumn. The samples were collected 30, 60, 90, 120, and 150 days after flowering (Figure 1). Since the temperatures in winter are low (on an average of 13.9°C), the peel and flesh are not very distinguishable, and it is not feasible to separate the peel in the case of the experimental group. However, the peel can be easily separated only after 60 days after flowering. Therefore, the sample at 30 days after flowering was not processed for the experimental group. The temperatures in both growing seasons are given in Supplementary Table 8. Three plants were selected for each time point, and triplicate samples were harvested for each time point in each season. The harvested fruits were surface-cleaned with distilled water and wiped to dry. Fruit peels were removed, wrapped in aluminum foil, and immersed and stored in liquid nitrogen.

Ribonucleic acid extraction, cDNA library construction, and sequencing were performed at Novogene (Beijing NuoheZhiyuan Technology Co., Ltd). Briefly, total RNAs were extracted using the Spin Column Plant total RNA Purification Kit (Sangon Biotech, Shanghai, China); their purity was assessed on 1% agarose gels and quantified using an Agilent Bioanalyzer 2100 (Agilent Technologies, Sta. Clara, CA, United States) as reported earlier (Chen et al., 2019). Sequencing libraries were created using the NEB Next Ultra RNA Library Prep Kit following the instructions of the manufacturer. Twenty-four paired-end cDNA libraries with an insert size of 300 bp were constructed. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego, CA, United States) according to the instructions of the manufacturer. After cluster generation, the library preparations were sequenced on a Novaseq (Illumina, San Diego, CA, United States) platform, and 150 bp paired-end reads were generated.

Raw data (raw reads) of fastq format were first processed through in-house Perl scripts, and clean data (clean reads) were obtained. The checking of sequencing error rate distribution was performed according to Q_phred = −10log10 (e), where e = sequencing error rate and Q_phred = is the based quality value of Illumina. The Q20, Q30, and GC content of the clean data was calculated as reported earlier (Goldstein et al., 2016). All the downstream analyses were based on clean data with high quality.

An index of the reference genome was built using Hisat2 v2.0.5, paired-end clean reads were aligned to the reference genome (Xu et al., 2013; Kim et al., 2019; Li et al., 2020), and read numbers mapped to each gene were counted using the feature Counts v1.5.0-p3 (Liao et al., 2014). We then calculated the FPKM (fragments per kilobase of transcript sequence per millions base pairs sequenced) of each gene, and the differential expression analysis was performed using the DESeq2 R package (1.20.0). The P-values were adjusted using Benjamini and Hochberg’s approach (Benjamini and Hochberg, 1995) in order to control the FDR (false discovery rate). The threshold for differentially expressed genes (DEGs) was P-value < 0.05. The transcripts were then annotated in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa and Goto, 2000). We used the cluster Profiler R package to test the statistical enrichment of differential expression genes in KEGG pathways.

Ten genes were randomly selected for verification of the transcriptome results through quantitative real-time PCR (qRT-PCR) reactions. We used gene-specific primers for qRT-PCR (Supplementary Table 7). The primers were designed in Primer3Plus (Untergasser et al., 2007). A housekeeping gene i.e., Actin, was used as an internal control. All the reactions were carried out in a Rotor-Gene 6000 (Qiagen, Shanghai, China) machine. The thermal cycles consisted of 50 (2 min) and 95°C (2 min) followed by 40 cycles of 95 (3 s) and 60°C (30 s). This was followed by the verification of single production amplification by melting curve analysis with temperatures ranging from 55 to 95°C by increasing 1°C every step. The total reaction volume was 10 μl (30 ng cDNA, 5 μl 1 × SYBR^® Select Master Mix (Applied Biosystem, Carlsbad, CA, United States), and.2 μl (20 μM) of each primer). We analyzed three biological replicates for each time point. The rest of the conditions for reactions and the analysis of the expression data were performed as reported earlier (Sendín et al., 2012; Chin et al., 2021).

Tissues (100 mg) were separately ground with liquid nitrogen and mixed with 80% methanol and.1% formic acid. The mixture was incubated on ice for 5 min and then centrifuged at 15,000 g for 20 min at 4°C. The supernatant was diluted to a final concentration containing 53% methanol by adding LC-MS grade water. The solution was again centrifuged at 15,000 g for 20 min at 4°C, and the resulting supernatant was used for liquid chromatography with tandem mass spectrometry (LC-MS/MS) (Want et al., 2006).

The ultra-high performance (UHP) LC-MS/MS analysis was conducted at Novogene Co., Ltd. (Beijing Novogene Technology Co., Ltd)¹ following their standard procedures. The raw data files obtained by UHPLC-MS/MS analyses were processed using the Compound Discoverer 3.1 (CD3.1; Thermo Fisher Scientific, Waltham, MA, United States) to perform peak alignment, peak picking, and quantitation for each metabolite. Briefly, the data were screened by retention time and mass-to-charge ratio. To make the metabolite identification accurate, peak alignment was performed according to retention time deviation and mass deviation (part per million, ppm). Following this, signal-to-noise ratio, adduct ion, and peak area were quantified in ppm. We then identified the metabolites by comparing the quantified data with mzCloud, mzVault, and the MassList primary database search library. Metabolites with a coefficient of variation less than 30% (Dai et al., 2017) in the QC sample were retained as the final identification result for subsequent analyses (Fernie et al., 2011).

Statistical analyses were performed using the statistical software R (R version R-3.4.3), Python (Python 2.7.6 version), and CentOS (CentOS release 6.6); when data were not normally distributed, normal transformations were attempted using area normalization method.

Based on the total putative identified metabolites as background in each pathway, we performed a hypergeometric test to find the significantly enriched KEGG pathways for the differentially accumulated metabolites with the R software. The p-values were adjusted by FDR with a threshold of less than or equal to.05.

The identified metabolites were annotated using the KEGG² (Kanehisa and Goto, 2000), HMDB³ (Wishart et al., 2018), and LIPIDMaps databases⁴. The principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA) were performed in metaX (Boulesteix and Strimmer, 2007; Wen et al., 2017). We performed a univariate analysis (t-test) to calculate the statistical significance (P-value). The metabolites with VIP > 1 and P < 0.05 and fold change ≥ 2 or FC ≤ 0.5 were considered to be differentially accumulated metabolites (DAMs). Volcano plots were used to filter metabolites of interest, which were based on log2 FC and −log10 (p-value) of the metabolites.

For clustering heatmaps, the data were normalized using z-scores of the intensity areas of differential metabolites and were plotted with the Pheatmap package in the R language. The Pearson’s correlation between differential metabolites was analyzed by cor in the R language. P-value < 0.05 was considered as statistically significant, and correlation plots were plotted with the corrplot package in the R language. The functions of these metabolites and metabolic pathways were studied using the KEGG database. The metabolic pathway enrichment of differential metabolites was performed; when ratio was satisfied by x/n > y/N, a metabolic pathway was considered as enrichment; when the P-value of a metabolic pathway < 0.05, the metabolic pathway was considered as statistically significant enrichment.