Transcriptome profiling and identifying differentially expressed genes (DEGs)

The transcriptomes of AZ-C from Dd, Rd, Dh, and Rh (Fig. 1a) were examined using the RNA-seq. Before RNA-sequencing analysis, CLas titer in AZ-C tissue was measured and the results verified that the fruits from HLB-diseased trees (Dd and Rd) were infected with CLas, while the fruits from healthy trees (Dh and Rh) were free of CLas (Fig. S1). Although the CLas titer in Dd was little bit higher than that in Rd, they were not statistically different (p > 0.05) (Fig. S1). Eight cDNA libraries were constructed and sequenced, two libraries for each of Dd, Rd, Dh, and Rh. The raw data have been deposited in NCBI Sequence Read Archive (SRA) through Gene Expression Omnibus (GEO) (access number: GSE101381). After filtering the low-quality and adapter sequences, a sequencing depth of 58–67 million 100 bp paired-end reliable reads per library was reached (Table S1). The clean reads were aligned to the reference genome sequences of Citrus sinensis v1.129. The RNA-Seq reads were mapped to 21,781–22,452 Citrus sinensis transcripts. The mapped transcripts and the corresponding Arabidopsis orthologs as well as their annotations are listed in Table S2.

Fig. 1: Images/Diagrams illustrating the fruit calyx abscission zone (AZ-C) tissue, correlations/distances among transcriptomes and the number of differentially expressed genes in AZ-C of different types of fruit.


a The fruit types and abscission zone tissue used in the experiment. b Correlation matrix and cluster dendrogram of the whole dataset of the mapped reads. The analysis was performed by comparing the values of the entire transcriptome (23310 transcripts) in all eight samples with two biological replicates. Correlation analysis (coefficients R2) and hierarchical cluster analysis were performed using R software. Red color indicates a stronger correlation and green weaker (R2). c The number of differentially expressed genes (DEGs) resulting from pairwise comparisons of transcriptomes between Dd and Rd (Dd/Rd), Dd and Dh (Dd/Dh), Dh and Rh (Dh/Rh), Rd and Rh (Rd/Rh)

Hierarchical cluster analysis of the reads was conducted to evaluate the distances among samples based on correlations regarding gene expression, and the results indicate that gene expression correlations between the biological replicates had high coefficients (R2 = 0.982–0.989, n = 23,310) (Fig. 1b), demonstrating the reliability of the data produced and illustrating the consistency of the transcriptional changes within each type of fruit. Dd groups were farthest away from all other groups, indicating that the transcriptome profile in the AZ-C of dropped fruit from HLB-diseased trees (Dd) was quite different from all other types of fruit. Although both from HLB-diseased trees, the similarity between retained (Rd) and dropped fruit (Dd) was less than between HLB and healthy fruit (Rd vs. Rh or Dh). Meanwhile, the Dh groups were closer to Rh than to Dd or Rd, indicating that the dropped fruit from healthy trees (Dh) was more similar to retained fruit on healthy trees (Rh) than to HLB-affected fruit (Rd or Dd).

Pairwise comparisons between Dd and Rd, Dh and Rh, Dd and Dh, and Rd and Rh resulted in four sets of regulated genes, respectively, for Dd vs. Rd, Dh vs. Rh, Dd vs. Dh, and Rd vs. Rh. Those genes with a fold change of at least 2 (|log2FC|≥1), and with a p-value less than 0.05 were considered differentially expressed genes (DEGs). There were 1047 (403 up- and 644 downregulated), 1599 (777 up- and 822 downregulated), 813 (553 up- and 260 downregulated), and 764 (521 up- and 243 downregulated) DEGs, respectively, for Dd vs. Rd, Dd vs. Dh, Dh vs. Rh, and Rd vs. Rh (Fig. 1c).

Significantly regulated functional gene groups

To gain a holistic understanding of the functional groups that were involved in the abscission of HLB-affected fruit, the DEGs were subjected to Wilcoxon rank-sum test analysis by PageMan30 integrated in the MapMan software31. The functionally related groups of genes were identified, which showed significant patterns of regulation compared with the complete collection of genes under analysis. A graphical summary, illustrating differentially regulated functional groups, is shown in Fig. 2. DEGs in each group and the results of Wilcoxon test are listed in Table S3.

Fig. 2: PageMan display of results of Bin-wise Wilcoxon rank sum test for significant MapMan gene categories that were regulated in calyx abscission zones of D compared to R (Dd vs Rd and Dh vs Rh) and counterparts between HLB-diseased and healthy (Dd vs Dh and Rd vs Rh).


Colored boxes indicate statistically significant groups (p-value < 0.05). The color scale represents regulation of gene expression, with red indicating a trend within the group for upregulation of expression, and green, downregulation. The arrows shown with the same color scheme, and the text alongside each row, provide MapMan annotation of differentially regulated gene classes

Although the total number of upregulated genes were less than that of downregulated genes in DEGs of Dd vs. Rd and Dd vs. Dh (Fig. 1c), more significantly upregulated than downregulated gene groups were identified by the Wilcoxon test (Fig. 2). The most noteworthy upregulated gene groups in AZ-C of Dd fruit include: lipid metabolism, amino acid and secondary metabolism, hormone metabolism–ET, and –jasmonate, signaling-receptor kinases, etc; while the most noteworthy downregulated gene groups in Dd were “transport” related and a group of non-annotated (unknown) genes (Fig. 2, lines Dd/Rd and Dd/Dh). Quite different from Dd, few gene groups regulated in Dh were identified by the test, with a mainly upregulated gene group related to amino acid metabolism (Fig. 2, line Dh/Rh). The comparison of Rd with Rh showed that gene groups related to major carbohydrate metabolism, cell wall degradation, and cell cycle were upregulated, while photosynthesis, abiotic stress, transport, and a group of functionally unknown genes were downregulated in Rd compared with Rh (Fig. 2, line Rd/Rh).

Since the goal of the study was to determine the mechanism of HLB-related fruit abscission, the analyses and discussions afterward will be focused on the comparison between dropped and non-dropped fruit from HLB-trees (Dd vs. Rd), and compared with those from healthy trees (Dh vs. Rh).

Enrichment analyses of DEGs for the significant biological processes and pathways

To understand the significant biological processes and pathways that are involved in the HLB-related fruit drop, the upregulated and downregulated DEGs from the comparisons between dropped and retained fruit (Dd vs. Rd and Dh vs. Rh) were subjected to the gene ontology (GO) and KEGG pathway enrichment analyses, which are illustrated in Figs. 3 and 4, respectively for Dd vs. Rd and Dh vs. Rh. All the GO terms and KEGG pathways are listed in Table S4. For the GO enrichment results illustrated in Fig. 3a, b, and Fig. 4a, only the GO terms with multiple-test adjusted p-value below 0.001 (highlighted by colors) are shown due to the big numbers of GO terms with multiple-test adjusted p-value below 0.05.

Fig. 3: Diagrams showing significant terms by gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of differentially expressed genes (DEGs) between dropped diseased (Dd) vs retained diseased (Rd).


a, b Significant GO terms under the category of “Biological Process” for upregulated (a) and downregulated (b) genes. Statistical test method: hypergeometric, and multi-test adjustment method: Yekutieli (FDR under dependency). c, d The significant KEGG pathways related to upregulated (c) and downregulated genes (d). Statistical test method: two-sided hypergeometric test

Fig. 4: Diagrams showing significant terms by the gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of differentially expressed genes (DEGs) between dropped healthy (Dh) vs. retained healthy (Rh).


a, b Significant GO terms under the category of “Biological Process” for upregulated (a) and downregulated (b) genes. Statistical test method: hypergeometric, and multi-test adjustment method: Yekutieli (FDR under dependency). c, d The significant KEGG pathways related to upregulated (c) and downregulated genes (d). Statistical test method: two-sided hypergeometric test

For Dd (Fig. 3), the significantly upregulated biological process (Fig. 3a) were those related to defense response (GO:0006950, GO:0006952, GO:0050896, GO:0042221, GO:0010033, GO:0009743, and GO:0010200), secondary metabolism (GO:0019748, GO:0006575, GO:0006725, and GO:0009698), and hormone signaling (GO:0009719, GO:0009725). The significantly upregulated pathways (Fig. 3c) include plant hormone jasmonic acid biosynthesis (“alpha-linolenic acid metabolism” pathway (KEGG:00592, p = 4.6E-5)) and the pathways related to secondary metabolism and defense response to the pathogen (“phenylpropanoid biosynthesis” (KEGG:00940, p = 9.9E-3), “phenylalanine metabolism” (KEGG:00360, p = 4.3E-2), “ubiquinone and other terpenoid-quinone biosynthesis” (KEGG:00130, p = 2.7E-2), “tyrosine metabolism” (KEGG:00350, p = 6.3E-3), and “plant–pathogen interaction” (KEGG:04626, p = 3.5E-2)). Among these biological processes or pathways, “response to chitin” was the most significant biological process (multi-test adjusted p = 9.95E-13), and “alpha-linolenic acid metabolism” (jasmonic acid biosynthesis) was the most significant pathway (p = 4.6E-5) that were upregulated in Dd compared with Rd.

The significantly downregulated biological processes in Dd (Fig. 3b) were related to plant-type cell wall modification (GO:0042545, GO:0009664, and GO:0009827), cell surface receptor-linked signaling (GO:0007166, GO:0007167, and GO:0007169), organ development (GO:0032502, GO:0007275, GO:0048856, GO:0048731, and GO:0048513), and lipid localization (GO:0051179, GO:0010876); while the downregulated pathways in Dd (Fig. 3d) included “starch and sucrose metabolism” (KEGG:00500, p = 4.4E-2) and “zeatin biosynthesis” (KEGG:00908, p = 2.0E-2).

For Dh (Fig. 4), the significantly upregulated biological processes were related to cell death (GO:0016265, GO:0009987, GO:0008219, GO:0012501, and GO:0006915) and stress/defense response (GO:0050896, GO:0006950, and GO:0006952) (Fig. 4a); while the upregulated pathways (Fig. 4c) were related to protein degradation, including “ubiquitin-mediated proteolysis” (KEGG:04120, p = 6.4E-3) and transcription processing, including “spliceosome” (KEGG:03040, p = 4.7E-2). The downregulated biological processes for Dh (Fig. 4b) include “response to chitin” (GO:0010200) and GO terms related to lipid biosynthesis (GO:0008610, GO:0032787, GO:0006631, and GO:0006633). The downregulated pathways in Dh (Fig. 4d) include the “alpha-linolenic acid metabolism” pathway (KEGG:00592, p = 1.2E-3) and “plant–pathogen interaction” (KEGG:04626, p = 1.4E-3), both of which were contrarily upregulated in Dd (Fig. 3c).

Taken together, the enrichment analysis results indicated that defense response (especially “response to chitin”), secondary metabolism, and hormone signaling (in particular jasmonic acid) were upregulated (Fig. 3a, c), while organ development and cell wall organization were downregulated in Dd (Fig. 3b, d), which is quite different from its healthy tree counterpart Dh. In Dh (Fig. 4), “response to chitin” and jasmonic acid (JA) biosynthesis were downregulated, while “cell death” and protein degradation were upregulated. Cell death in Dh may be linked to wounding caused by the strong force applied to the healthy trees (in order to get enough numbers of “dropped fruit”). The cell wall degradation in Dd may be directly linked to the abscission process activated by plant hormones (such as JA) resulting from defense response.

The regulated genes under GO terms of “response to chitin” (GO:0010200) and “secondary metabolic process” (GO:0019748) are listed, respectively, in Table S5 and Table S6. Hormone-related genes will be further analyzed separately. Twenty chitin-responsive genes and 21 secondary metabolism-related genes were upregulated in Dd, but none of them were upregulated in Dh; and instead, some of them (six chitin-responsive genes and four secondary metabolism-related genes) were downregulated in Dh (Tables S5, S6). More than half (12 out of 21) of the secondary metabolism-related genes that were upregulated in Dd are from the phenylpropanoid pathway (Table S6, highlighted in red).

Plant hormone-related gene expression

Plant hormones are known to mediate abscission signals15. The enrichment analyses showed “response to hormone stimulus” was the significant (p = 6.5E-4) biological process (Fig. 3a) and “jasmonic acid biosynthesis” (or “alpha-linolenic acid metabolism”) (p = 4.6E-5) was the significant pathway that were upregulated in Dd fruit (Fig. 3c). Therefore, the DEGs related to hormone biosynthesis and signaling were analyzed in detail, which are listed in Table S7. The profile of hormone gene expression was characterized as upregulation of ET and JA; and general downregulation of other hormones, including abscisic acid (ABA), auxin (AUX), brassinosteroid (BR), cytokinin (CK), and gibberellin (GA) in Dd, compared with Rd or Dh (Table S7). Different from ET and JA, salicylic acid (SA) was upregulated in dropped fruit from HLB-diseased trees (Dd) only when compared with dropped fruit from healthy trees (Dh), but not obviously regulated if compared with retained fruit on HLB-diseased trees (Rd) (Table S7).

The upregulation of ET in Dd included several 1-aminocyclopropane-1-carboxylate (ACC) oxidase genes (ACO1, EFE (or ACO4), and MJM20.4 (or ACO11)) for ET biosynthesis and eight ERF genes (ERF1, ERF4, ERF9, ERF13, ERF17, ERF20, ERF25, and ERF109) for ET signaling. ACC synthase gene (ACS) expression was not regulated; however, gene expression of an ACS inhibitor, named “ethylene-overproduction protein” (ETO1), was downregulated. ETO1 is known to promote ACS degradation and inhibit ACS enzyme activity32. The upregulation of JA included genes for JA biosynthesis (LOX3 (lipoxygenase 3 gene), OPR1 and OPR2 (12-oxophytodienoate reductase genes)), gene for methylation of jasmonate into methyljasmonate (JMT, jasmonate methyltransferase gene), and the JA-amino acid synthetase gene (JAR1). Moreover, gene expression of the key repressor of JA signaling, JAZ3 (jasmonate ZIM-domain protein 3)33, was downregulated in Dd fruit. The downregulation of ABA in Dd included the key genes for ABA biosynthesis (NCED3 and NCED4) and two ABA-responsive protein genes (HVA22E and GEM) (Table S7).

Figure 5a shows gene expression levels (mapped reads) viewed by Integrative Genomics Viewer (IGV), for some example DEGs involved in pathways of ET, JA, or ABA.

Fig. 5: Gene expression level and protein-protein interaction network for gene products related to ethylene, jasmonic acid and abscisic acid pathways.


a Gene expression level viewed by Integrative Genomics Viewer. b, c Protein–protein interaction network showing regulation of gene expression in the comparisons of Dd vs. Rd (b) or Dh vs. Rh (c). Red indicates the upregulated and green indicates the downregulated

Protein–protein interaction (PPI) network for the JA/ET pathway-related proteins was predicted based on the knowledge of PPI network of Arabidopsis34. The gene expression data (log2FC) for Dd vs. Rd and Dh vs. Rh were mapped to the PPI network, in which the regulated proteins/gene expression are illustrated by red or green color, indicating up- or downregulation, respectively; white indicates no change (Fig. 5b, c). The network includes the JA/ET pathway proteins that have been known to have interactions with each other, and proteins of other pathways/functional groups known to have connections with the JA/ET pathway proteins. For Dd vs. Rd (Fig. 5b), numerous genes of ET or JA synthesis or signaling pathways were upregulated, while genes for negative regulators of ET (ETO1) or JA (JAZ3) pathways were downregulated (Fig. 5b). Dh vs. Rh (Fig. 5c) was quite different from Dd vs. Rd, with only two genes (IAR3 and JAR1) that were similarly regulated to Dd vs. Rd, which were upregulated in both Dd and Dh (Fig. 5b, c). IAR3 is an enzyme that hydrolyzes amino acid conjugates of IAA or JA, and is induced in response to wounding35. JAR1 is known to play a role in JA-mediated defense response36, and as well, is induced in response to wounding37. Also, shown in the network are several gene products from secondary metabolism pathways that have direct connections with ET/JA, among them, are PAL1, CCoAOMT1, FAH1, and ATR2, gene expressions of which were upregulated in Dd vs. Rd (Fig. 5b), but were not downregulated in Dh vs. Rh (Fig. 5c).

Consistent with the upregulation of ET and JA, substantial numbers of downstream JA/ET-responsive defense-related genes were upregulated in Dd compared with both Rd and Dh fruit (Table S8). These defense genes include plant defensin type 1 gene (PDF1.4)38, an antifungal chitin-binding protein gene called pathogenesis-related 4 (PR4, or HEL), polygalacturonase-inhibiting protein gene (PGIP1), beta-1,3-glucanase gene (BG1), sulfotransferase 2A gene (ST2A), alpha-dioxygenase 1 gene (DOX1), cysteine proteinase gene (RD21A), several chitinase genes (CHIB1, CHIC, and CHIV), and plant U-box protein genes (PUB21, PUB22, PUB23, PUB24, and PUB29). None of these defense response genes were upregulated in Dh compared with Rh. However, several (5 out of 21) of them were higher in HLB-affected fruit (Rd) compared with healthy fruit (Rh) (Table S8).

Quantitative RT-PCR validation of DEGs

To verify the RNA-Seq data by qRT-PCR, 30 genes were selected as representatives covering each of the DEG categories identified by RNA-Seq analysis, including genes related to the seven categories of hormones (with an emphasis on ET and JA), secondary metabolism, and JA/ET-mediated defense responses (Fig. S2a). The qRT-PCR results for the 30 genes were comparable with the RNA-Seq results, with a squared correlation-coefficient value of 0.87 (Fig. S2b).

Phytohormone production measurements

To verify the gene expression data at final product level, phytohormone ET, JA, abscisic acid (ABA), indole-3-acetic acid (IAA), and salicylic acid (SA) were measured. Gaseous ET was directly measured for individual fruit. Being a gaseous compound ET affects fruit and plant parts that are remote from site of synthesis, and fruit ET production has been shown to be a predictor of abscission39,40. Fruit ET production was detected in 39 out of 60 Dd fruit, but none of the Rd, Dh, or Rh fruit exhibited detectable ET production under the measurement conditions used (Fig. 6a). The average ET production of the 60 Dd fruit was 0.188 (µL fruit-1 h-1), and statistical analysis indicated that ET production from Dd fruit was significantly higher than the other three types of fruit (p = 3.84E-11), which had no ET detected (Fig. 6a).

Fig. 6


Barplot showing phytohormone production by Dd, Rd, Dh, or Rh

Non-volatile acidic phytohormones JA, ABA, SA, and IAA were extracted from AZ-C tissue and derivatized to volatile methyl esters. The GC–MS results showed significantly higher level of JA (Fig. 6b) and lower level of ABA (Fig. 6c) in Dd compared with Rd, Dh, or Rh AZ-C tissue. IAA (Fig. 6d) and SA (Fig. 6e) showed no significant differences among Dd, Rd, Dh, and Rh; however, there was a trend where IAA was lower in Dd fruit and SA higher compared with Rd, Dh, or Rh.

Detection of fungus Diplodia infection in AZ-C

Since the RNA-seq analysis revealed a gene expression profile of antifungal defense (upregulation of “response to chitin”, “secondary metabolism” and “JA/ET signaling”, and “JA/ET-activated defense response”) in AZ-C of Dd fruit, and a higher incidence of fungus Diplodia had been reported in the AZ-C of HLB-affected citrus fruit compared with healthy controls27, we tested Diplodia infection in AZ-C of Dd, Rd, Dh, and Rh fruit. The qPCR results (Fig. S3) showed that Diplodia level in Dd was significantly higher than in Rd, Dh, or Rh (p < 0.001); while there was no difference among Rd, Dh, and Rh (p > 0.05) (Fig. S3).