Transcriptome differential expression analysis of defoliation of two different lemon varieties

Experimental materials

Previously, our group conducted observation records and research on two lemon varieties (‘Allen Eureka’ and ‘Yunning No. 1’) for three consecutive years (Figs.1A and 1B). As shown in Fig. 1, both lemon varieties began defoliation from November 11–November 20, with peak defoliation from November 21-December 10, and defoliation ending between December 11-December 20. Five rootstock types were selected for comparison of rootstock combinations with ‘Allen Eureka’ scions, and the five rootstocks were selected from C.junos (Sieb.) Tanaka, C.grandis (L.) Osbeck, C.limonia, Valencia orange (nondeciduous rootstock in winter), and Poncirus trifoliate (L.) Raf (deciduous rootstock in winter). Nine plants were grafted with each rootstock type, and an average plant height of 60 cm was transplanted. A total of 1 year after transplanting, it was found that different rootstocks had a great influence on the growth of ‘Allen Eureka’. The ‘Allen Eureka’-scion/C.junos (Sieb.) Tanaka-rootstock combination had the weakest growth, with an average plant height of 90 cm; the ‘Allen Eureka’/Valencia orange had the strongest growth, with an average plant height of 180 cm; and the ‘Allen Eureka’-scion/Poncirus trifoliate (L.) Raf-rootstock combination had medium growth, with an average plant height of 120 cm. Abnormal winter defoliation was also observed in all five rootstocks grafted to ‘Allen Eureka’, and the heaviest defoliation was in the ‘Allen Eureka’-scion/Poncirus trifoliate (L.) Raf-rootstock combination, with cumulative winter defoliation of 180 leaves/plant. The lowest defoliation was in the ‘Allen Eureka’-scion/C.grandis (L.) Osbeck-rootstock combination with a cumulative winter defoliation of 117 leaves/plant, which indicated that this defoliation state of ‘Allen Eureka’ was stable and that the defoliation state of ‘Allen Eureka’ could not be solved by changing the rootstocks. Second, the defoliation of the ‘Yunning No. 1’-scion/Poncirus trifoliate (L.) Raf-rootstock was also compared with that of the ‘Allen Eureka’-scion/Poncirus trifoliate (L.) Raf-rootstock, and the cumulative winter defoliation of ‘Yunning No. 1’-scion/Poncirus trifoliate (L.) Raf-rootstock, at the same time and under the same cultivation conditions, was 23 leaves/plant, which is a normal physiological defoliation of old leaves. The results of the period study showed that ‘Allen Eureka’ and ‘Yunning No. 1’ exhibited significant varietal differences.

For experimental materials, the scions of two lemon varieties, ‘Allen Eureka’ (Allen Eureka) and ‘Yun Lemon No. 1’ (Yuning No. 1), were selected, and ‘Yun Lemon No. 1’ of the same defoliation period was used as a control and grafted onto the annual Poncirus trifoliate (L.) Raf rootstocks in November 2019. When the grafted seedlings of the two lemon varieties reached an average plant height of 60 cm, the seedlings of the two lemon varieties with good growth conditions and basically the same size were selected for transplanting and cultivation, respectively. The culture container was a plastic bucket (diameter of the mouth of the bucket 60 cm, diameter of the bottom of the bucket 55 cm, height of the bucket 60 cm), and the culture substrate was cocopeat: red loam soil: organic fertilizer = 6:3:1. One hundred and eighty plants of each variety were planted (a total of 360 plants) and placed in the Institute of Tropical Subtropical Cash Crops, Yunnan Academy of Agricultural Sciences in a rain-sheltered arched shed (97° 87′ E, 24° 2′ N). The height of the shed was 3.2 ~ 3.5 m. General management, the average height of the plants was 200 cm and the average crown width was 150 cm, for the experimental sampling.

Sampling

Based on the observation and recording of the phenological period of two lemon varieties, ‘Allen Eureka’ and ‘Yunning No. 1’, sampling was carried out according to the time of defoliation and the amount of defoliation in three periods, namely, the pre-defoliation stage (the beginning of defoliation, CQ), the mid-defoliation stage (a large amount of defoliation, CZ), and the post-defoliation stage (near the end of defoliation, CH). The sampling times were November 15, 2021, December 1, 2021, and December 15, 2021 (Table S1). Each variety was divided into nine groups (three groups for enzyme activity determination, three groups for hormone determination, and three groups for transcriptome determination). First, 1.0 g of the middle and lower mature petioles of barrel-planted plants from different treatment periods were collected from the AZ (0.3–0.5 cm at the base of the petiole) (Fig. S1), and the samples were transferred to liquid nitrogen quickly frozen after sample collection and placed at −80 °C for freezing and preservation for future use.

Measurement of enzyme activities

A kit was used to determine superoxide dismutase (SOD) activity, peroxidase (POD) activity, catalase (CAT) activity, malondialdehyde (MDA) content, cellulase activity, and pectinase activity, and was biologically repeated three times. The primary test instrument was an MD SpectraMax 190 full-wavelength enzyme labeler.

Hormone determination

In this experiment, the 1-aminocyclopropane-1-carboxylic acid (ACC), abscisic acid (ABA), auxin (IAA), and gibberellin (GA3) contents of the two lemon varieties at three periods were measured using high-performance liquid chromatography-tandem mass spectrometry (HPLC‒MS/MS), repeated three times. The main testing instruments were an Agilent 1290 high-performance liquid chromatograph (Agilent, Santa Clara, CA, USA) and a SCIEX-6500 Qtrap (MSMS) (AB SCIEX, Framingham, MA, USA).

RNA extraction, RNA-seq library preparation, and transcriptome sequencing

Total RNA was extracted from 18 lemon samples using a plant RNA extraction kit (QIAGEN, Hilden, Germany), and three periods of two lemon varieties, three biological replicates (named ACQ-1, ACQ-2, ACQ-3, ACZ-1, ACZ-2, ACZ-3, ACH-1, ACH-2, ACH-3, YCQ-1, YCQ- 2, YCQ-3, YCZ-1, YCZ-2, YCZ-3, YCH-1, YCH-2, and YCH-3), and three technical replicates were performed for each sample. After that, mRNA was purified using magnetic beads and Oligo (dT), and the mRNA was broken up into small pieces by adding fragmentation buffer. Using six-base random hexamers as a template and mRNA as the starting material, the primary cDNA strand was created. The second cDNA strand was subsequently produced by adding buffer, dNTPs, RNase H, and DNA polymerase I. After purification with the QiaQuick PCR kit and elution with EB buffer, the second cDNA strand was end-repaired, poly (A) was added, the sequencing junctions were connected, the fragment sizes were selected by agarose gel electrophoresis, and then PCR was performed. The sequencing library was constructed and then sequenced with Illumina HiSeq 6000.

Transcriptome data analysis

In order to remove junctions and low-quality data from the transcriptome data, the raw data were filtered using fastp software (https://github.com/OpenGene/fastp) to remove reads containing junctions, reads with an N-ratio of 10% or more, and reads with a low-quality (quality value less than 20) base ratio of more than 50%. Finally, FastQC software was used to check the quality of the clean data, and subsequent analysis was carried out after passing the quality control. All downstream analyses were based on sequence data. We downloaded lemon reference genome and gene model annotation files directly from the Genome website (https://www.citrusgenomedb.org/jbrowse/index.html?data=data/Climon_Alt_v1.0f) (Bao et al., 2023). The index of the reference genome was built using HISAT2 v2.0.5, and paired-end clean reads were aligned to the reference genome using HISAT2 v2.0.5 (Kim, Langmead & Salzberg, 2015; Shirasawa et al., 2017). The mapped reads of each sample were assembled by StringTie (v1.3.3b) (Pertea et al., 2015) in a reference-based approach. With |log₂Fold Change| ≥ 1 and padj ≤ 0.05 as the thresholds, differential expression analysis of two conditions/groups (two biological replicates per condition) was performed using the DESeq2 R package (1.20.0). Gene Ontology (GO) and KEGG pathway enrichment analysis of differentially expressed genes (DEGs) was implemented using the clusterProfiler R package.

RT-qPCR analysis

The expression levels of five DEGs selected from the RNA-seq data were verified by RT-qPCR using the cDNA samples used for RNA-Seq library construction. Primer 5.0 software was used to design primer pairs. The primer sequence was synthesized by Beijing Tian Yi Hui Yuan Biotechnology Co., Ltd. (Beijing, China), and the primer information is shown in Table S2. RT‒qPCR was performed using a Roche LightCycler 480 with the following procedure: 95 °C for 5 min, followed by 40 cycles at 95 °C for 10 s, 60 °C for 10 s, and 60°C for 30 s. The housekeeping gene Actin was used as an internal control. The relative expression level was calculated by the 2^−ΔΔCt method.

Statistics and analysis of experimental data

Gene expression data were expressed as the mean ± standard deviation. The significance between means was statistically analyzed using a t-test (p < 0.05). Microsoft Office Excel 2016 was utilized for statistical processing of data as well as graphing, and the results were expressed as the mean ± standard deviation. Gene coexpression networks were constructed using the R software WGCNAv1.69 package.

Changes in enzyme activities during different defoliation periods in two lemon varieties

The SOD, MDA, cellulase, and pectinase activities of the two lemon varieties, ‘Allen Eureka’ and ‘Yunning No. 1’, all gradually increased with prolonged defoliation time. The content of SOD in ‘Allen Eureka’ was significantly higher than that of ‘Yunning No. 1’ in the mid- and post-defoliation stages of defoliation. CAT activity tended to increase and decrease with prolonged defoliation, but the magnitude of change was small in both varieties. POD activities all showed a significant decreasing trend, with the prolongation of defoliation time being significantly lower in the mid-defoliation stage (CZ) and post-defoliation stage (CH) than in the pre-defoliation stage (CQ), and the differences between the two varieties were distinct (Fig. 2). The protective enzymes SOD, POD, and CAT showed an upward trend in a short period of time when the plants faced adversity, but as the stress intensified or the stress time prolonged, the activities of the protective enzymes eventually declined, which may be attributed to the fact that the duration and intensity of the stress exceeded the upper limit of the plant’s tolerance. Although both SOD and POD were effective in removing excess reactive oxygen species and superoxide anion radicals from the cells, the MDA content of both lemon varieties kept increasing, indicating that the activities of the two enzymes were not sufficient to remove reactive oxygen species and superoxide radicals produced in the cells throughout the defoliation process. Increased cellulase and pectinase activities can both lead to faster cell wall degradation and accelerated leaf abscission.

Changes in endogenous hormone contents during different defoliation stages in two lemon varieties

The ACC, ABA, and GA3 content of the two lemon varieties, ‘Allen Eureka’ and ‘Yunning No. 1’, gradually increased with the prolongation of defoliation time, and the ACC, ABA, and GA3 content of ‘Yunning No. 1’ was significantly higher than that of ‘Allen Eureka’ in both the mid- and post-defoliation stages. The IAA content of the two lemon varieties showed a decreasing trend from the pre-defoliation stage to the mid-defoliation stage. At the mid-defoliation stage, the IAA content of ‘Allen Eureka’ was significantly higher than that of ‘Yunning No. 1’, and the IAA content of the two lemon varieties increased again from the mid-defoliation stage to the post-defoliation stage (Fig. 3). Increased levels of ACC and ABA promote organ abscission. IAA both inhibits and promotes abscission, and its effect on organ abscission is related to the concentration, timing, and application site of the growth hormone used. The role played by GA in plant organ abscission is still controversial.

Statistics of off-region transcriptome sequencing data

After sequencing quality analysis, 18 libraries were constructed from three replicates across three periods of ‘Allen Eureka’ and ‘Yunning No. 1’. A total of 116.35 Gb of clean data was obtained for 18 samples with Q30 > 89.15%. The number of raw reads ranged from 39,849,322 to 46,685,856. Both clean and uniquely mapped reads made up more than 70% (Table 1; Table S3), indicating that the reference genome was selected appropriately. The Pearson correlation coefficients between the biological replicates for each sample ranged from 0.89 to 0.97. The 18 samples were subjected to principal component analysis (PCA). After dimensionality reduction, it was found that the samples of the same variety and period were clustered together (Fig. 4), and the results of the analysis showed good reproducibility of all the samples used for this experiment.

Number of DEGs at different defoliation stages in two lemon varieties

Two-by-two comparisons of the ‘Allen Eureka’ and ‘Yunning No. 1’ genotypes were performed at the pre-defoliation, mid-defoliation, and post-defoliation stages of the two lemon varieties to identify lemon defoliation-related genes. When compared to ‘Yunning No. 1’, 898 genes in ‘Allen Eureka’ were significantly differentially expressed, 320 genes were significantly upregulated, and 578 genes were significantly downregulated in the pre-defoliation stage (CQ). At the mid-defoliation stage (CZ), 4,856 genes were significantly differentially expressed, 2,513 genes were significantly upregulated, and 2,343 genes were significantly downregulated. At the post-defoliation stage (CH), 3,126 genes were significantly differentially expressed, with 1,336 genes significantly upregulated and 1,790 genes significantly downregulated (Fig. 5A), and the most DEGs in the mid-defoliation stage (CZ) out of the three periods (Table S4). A Venn diagram (Fig. 5B) was plotted based on all the DEGs across the three periods, and a total of 143 DEGs in the three periods compared with 400 (CQ), 3,621 (CZ), and 1,990 (CH) genes were expressed in only one period.

GO enrichment analysis of DEGs at different defoliation stages in two lemon varieties

Additionally, GO enrichment analysis showed that 99 (pre-defoliation stage CQ), 338 (mid-defoliation stage CZ), and 378 (post-defoliation stage CH) DEGs were annotated to the GO database for the ‘Allen Eureka’ and ‘Yunning No. 1’ lemons in the three periods, respectively. According to the significance level of GO term enrichment, the differential genes of the two varieties were mainly enriched in the molecular functional classifications of xyloglucosyl transferase activity, glucosyltransferase activity, transferring glycosyl groups, iron ion binding, oxidoreductase activity, hydrolase activity, DNA-binding transcription factor activity, transcriptional regulatory activity, and classification of cellular components of the apoplast, extracellular region, cell wall, external encapsulating structure, and cell periphery. In addition, biological processes included cellular glucan metabolic, glucan metabolic, cellular polysaccharide metabolic, polysaccharide metabolic, and cellular carbohydrate metabolic processes (Fig. 6). The related differential genes in molecular functions, cellular components, and biological process classifications were significantly different in all three periods, with most of them enriched to significantly higher differential genes in the mid-defoliation stage (CZ) GO than in the pre-defoliation stage (CQ) and post-defoliation stage (CH).

KEGG enrichment analysis

The results of KEGG metabolic pathway enrichment of DEGs between two lemon varieties, ‘Allen Eureka’ and ‘Yunning No. 1’, at three periods (Fig. 7) showed that the number of DEGs and related metabolic pathways enriched at the mid-defoliation stage (CZ) was the highest, which mainly focused on phenylpropanoid biosynthesis, plant hormone signal transduction, protein processing in the endoplasmic reticulum, amino sugar and nucleotide sugar metabolism, and glutathione metabolism, cyanoamino acid metabolism. The number of differential genes and related metabolic pathways enriched in the post-defoliation stage (CH) was high, which focused on phenylpropanoid biosynthesis, plant hormone signal transduction, protein processing in the endoplasmic reticulum, amino sugar and nucleotide sugar metabolism, pentose and glucuronate interconversions, and alpha-linolenic acid metabolism pathways. Relatively few metabolic pathways were enriched in the pre-defoliation stage (CQ), which mainly focused on protein processing in the endoplasmic reticulum, plant hormone signal transduction, amino sugar and nucleotide sugar metabolism, and phenylpropanoid biosynthesis pathways.

DEGs of the plant hormone signal transduction pathway

The KEGG enrichment analysis revealed that the DEGs between ‘Allen Eureka’ and ‘Yunning No. 1’ lemons were more abundant in all three periods of the plant hormone signal transduction pathway, and comparing the three periods, the mid-defoliation stage (CZ) had the most DEGs in the plant hormone signal transduction pathway. Among them, Aux/IAA (CL9G066930012_alt, CL3G046631012_alt, and CL3G046634012_alt), SAUR (CL0G069923012_alt), GH3 (CL2G041491012_alt), and the two-component response regulator ARR (CL9G067813012_alt) were downregulated in the mid-defoliation stage compared to the pre-defoliation stage, and the relative expression of ‘Allen Eureka’ was significantly higher than that of ‘Yunning No. 1’ (Fig. 8A). The expression of growth hormone-responsive genes Aux/IAA, SAUR, and GH3-related genes were all down-regulated, suggesting that the IAA signaling pathway was inhibited, and that the massive leaf abscission in the middle stage of defoliation might be related to the down-regulation of these growth hormone-responsive genes.

Construction of gene coexpression networks

To gain insight into the coexpressed genes of the two lemon varieties at three defoliation stages (CQ, CZ, and CH) and to mine genes related to defoliation height, weighted gene co-expression network analysis (WGCNA) was performed. By screening the weight values, β = 14 was finally chosen to construct the network, the dynamic shear tree method was used to merge the modules with similar expressions, and a total of 11 coexpression modules were obtained (Fig. 9). The turquoise module contained the highest number of genes (2,926 genes), followed by the blue module containing 1,923 genes, and the lowest number of genes was within the gray module, which contained 69 genes. The average number of genes contained in each module was 737.

Four of the 11 modules were highly specific to the sample. One was present in each of the pre-mid-defoliation stages, and two were present in the post-defoliation stage (Fig. S2). Five reported defoliation-associated genes (CL2G044356012_alt, CL6G058882012_alt, CL9G068078012_alt, CL1G038862012_alt, and CL3G046093012_alt) were found in the blue module by enrichment analysis. Therefore, the present study focused attention on the blue module, screened the genes with high connectivity in the module, and visualized the network using Cytoscape software (Fig. S3). Six genes of ribosomal proteins (CL6G060070012_alt, CL1G038121012_alt, CL6G059277012_alt, CL7G060527012_alt, and CL9G068716012_alt) and gibberellin-regulated proteins (CL5G054544012_alt) were found to have the connectivity of the module ranked in the top 1% of the blue module, which can be used as the hub genes in the module. In addition, WRKY genes (CL6G058882012_alt and CL1G038862012_alt) ranked low in connectivity in the module, but these two genes were presumed to be key genes in the leaf abscission regulatory pathway. Organ shedding involves two important processes, cell division and hydrolase induction, both of which are based on active metabolism. The protein and RNA content of the AZ significantly increase during abscission, and various metabolic inhibitors and protein synthesis inhibitors inhibit the formation of the AZ. Gibberellin regulates organ shedding and has a range of effects that are more pronounced at higher concentrations. WRKY transcription factors are mainly expressed during organ abscission in response to abscisic acid and ethylene, which in turn regulate organ abscission.

RT-qPCR validation of DEGs

In this study, we also used RT‒qPCR to validate five genes involved in plant hormone synthesis, signaling, and other processes related to plant organ abscission: peroxidase (CL2G044126012_alt), Aux/IAA gene family (CL3G046634012_alt, CL9G066930012_alt), weak ethylene-insensitive protein (CL2G044587012_alt), and protein phosphatase (CL8G065380012_alt). The results (Fig. 10; Table S5) showed that the changes in the fluorescence quantitative expression levels of the five DEGs converged with the changes in transcriptome gene abundance.