Integrated Transcriptomics and Metabolomics Analyses Provide Insights Into the Response of Chongyi Wild Mandarin to Candidatus Liberibacter Asiaticus Infection

Candidatus Liberibacter asiaticus (CLas) is the causative agent of Huanglongbing (HLB), which has caused great economic losses to the citrus industry. The molecular mechanism of the host response to CLas in wild citrus germplasm has been reported less. Eighteen weeks after inoculation via grafting, all the CLas-inoculated Chongyi wild mandarin (Citrus reticulata) were positive and showed severe anatomical aberrations, suggesting its susceptibility to HLB. Transcriptomics and metabolomics analyses of leaves, barks, and roots from mock-inoculated (control) and CLas-inoculated seedlings were performed. Comparative transcriptomics identified 3,628, 3,770, and 1,716 differentially expressed genes (DEGs) between CLas-infected and healthy tissues in the leaves, barks, and roots, respectively. The CLas-infected tissues had higher transcripts per kilobase per million values and more genes that reached their maximal expression, suggesting that HLB might cause an overall increase in transcript accumulation. However, HLB-triggered transcriptional alteration showed tissue specificity. In the CLas-infected leaves, many DEGs encoding immune receptors were downregulated. In the CLas-infected barks, nearly all the DEGs involved in signaling and plant-pathogen interaction were upregulated. In the CLas-infected roots, DEGs encoding enzymes or transporters involved in carotenoid biosynthesis and nitrogen metabolism were downregulated. Metabolomics identified 71, 62, and 50 differentially accumulated metabolites (DAMs) in the CLas-infected leaves, barks and roots, respectively. By associating DEGs with DAMs, nitrogen metabolism was the only pathway shared by the three infected tissues and was depressed in the CLas-infected roots. In addition, 26 genes were determined as putative markers of CLas infection, and a hypothesized model for the HLB susceptibility mechanism in Chongyi was proposed. Our study may shed light on investigating the molecular mechanism of the host response to CLas infection in wild citrus germplasm.


INTRODUCTION
Also known as citrus greening disease, Huanglongbing (HLB) is one of the most destructive diseases of citrus and has resulted in huge economic losses to the citrus industry. It is caused by the phloem-limited, gram-negative, alpha-proteobacteria Candidatus Liberibacter spp. [i.e., C. L. asiaticus (CLas), C. L. africanus (CLaf), and C. L. americanus (CLam)] (Albrecht and Bowman, 2008). CLas is transmitted via grafting with CLas-infected scions and phloem-feeding psyllids [Asian citrus psyllid (ACP) Diaphorina citri] (Bové, 2006). Once infected, the anatomical aberration in newly growing flushes was observed at the early invasion of the pathogen (Folimonova and Achor, 2010). As HLB develops, infected citrus plants will manifest a series of typical symptoms, such as leaves with an asymmetrical pattern of blotchy yellowing or mottling and malformed fruits with aborted seeds (Fu et al., 2016). Eventually, CLas can reduce resource availability for roots and then leads to starvation, which results in widespread tree death due to the decreased photoassimilate transport (Wang and Trivedi, 2013).
To investigate the underlying pathogenetic mechanism of HLB, studies have scanned changes in gene expression levels in CLas-infected citrus (Wang et al., 2016;Rawat et al., 2017;Zhao et al., 2019). In CLas-infected leaves, many differentially expressed genes (DEGs) involved in carbohydrate metabolism, particularly starch pathways, were found to be significantly upregulated (Albrecht and Bowman, 2008). A study performed on citrus phloem by microarray revealed that HLB could disrupt various biological functions, such as sugar metabolism, plant defense, phytohormones, and cell wall metabolism (Kim et al., 2009). The transcriptome of the CLas-positive citrus roots showed that HLB could alter some biological processes, such as cell wall modifications, protease-involved protein degradation, carbohydrate metabolism, hormone synthesis and signaling, transcription activities, and stress responses (Zhong et al., 2015). Nevertheless, there are rare reports on the comparison of transcript differences among various CLas-positive tissues, especially in wild citrus germplasms.
Moreover, bioinformatic analysis of four RNA-Seq datasets on HLB response and tolerance in leaves revealed that sucrose and starch metabolism was highly linked with disease symptoms, but that DEGs in each data set varied because of the conditions of leaves, infection stages, and environments (Balan et al., 2018). Therefore, analysis on transcriptional level alone might be difficult to screen the critical molecular signals in response to CLas infection. Metabolomics could identify low-molecularweight metabolites (<1 kDa) that alter their accumulation in response to physiological challenges (Padhi et al., 2019;Huang et al., 2020). Integrated transcriptomics and metabolomics analysis may identify key HLB-related pathways, in which DEGs and differentially accumulated metabolites (DAMs) are both enriched.
To evaluate the influences of HLB on various tissues, RNAseq and metabolome libraries were constructed with leaves, barks, and roots from healthy and CLas-infected Chongyi wild mandarin (Citrus reticulata), wild citrus germplasm in Jiangxi province, China (Peng et al., 2021). Based on the number of reads per gene, the number of DEGs between CLas-infected vs. healthy leaves/barks/roots was estimated, and their putative functions were investigated. Parallel gene expression changes among tissues were compared by scaling transcripts per million. Meanwhile, relationships between DEGs and DAMs in typical HLB-affected pathways were analyzed.

Plant Materials and Inoculation
Seeds of Chongyi wild mandarin (C. reticulata) were sampled in Chongyi County, Jiangxi province, China. One-year-old seedlings were grafted with healthy (control) and CLas-infected buds from sweet orange trees. All the plants were maintained in the same greenhouse (at 25-28 • C). The leaves were monitored biweekly by nested PCR to detect CLas. After 18 weeks, leaves, barks, and roots from CLas-positive and healthy seedlings were collected for paraffin sectioning, RNA-Seq, and ultra-high performance liquid tandem chromatography quadrupole time of flight mass spectrometry (UPLC-QTOFMS).

Paraffin Sectioning
Tissues were fixed in formaldehyde-acetic acid-ethanol fixative, dehydrated in ethanol, and then embedded in paraffin. After being deparaffinized with xylene and stained with a fast green, the sliced sections were analyzed with a DM6 B (Leica, Wetzlar, Germany) microsystem combined with a DFC 7000T (Leica, Wetzlar, Germany) camera.

Transcriptome Profiling and Data Processing
The RNA-Seq libraries were constructed, and high-throughput sequencing was carried out by Novogene (Beijing, China) using HiSeq 2000 (Illumina, San Diego, CA, United States). All sequencing adaptors, unknown nucleotides, and low-quality reads (quality score < 20) were removed from the raw reads, and the resulting clean reads were used in subsequent analyses. HISAT2 (Kim et al., 2015) was performed to map the clean reads to the sweet orange (Clonorchis sinensis) genome (Xu et al., 2013). StringTie was used to estimate the abundance of transcripts in each sample (Pertea et al., 2016), and then the information of all the transcripts was assembled into a single gene transfer format (GTF) file using the "merge" option. Using the resulting GTF file as reference, StringTie was again used to obtain total gene expression information. Subsequently, the R subread package was applied to compute the number of uniquely mapping reads unambiguously attributed to each gene (Liao et al., 2019), which was used as input in DESeq2 (Love et al., 2014). A matrix of TPM per gene across all samples was generated to compare gene expression differences among the tissues.

DEG Detection and Functional Enrichment
To check for similarities among these samples, a principal component analysis (PCA) in DESeq2 was performed based on per-gene unique reads counts. Genes with a total read number across all samples ≥18 were kept. Based on the same dataset, a clustering analysis was also performed using hclust in the "the stats" package in R (http://www.R-project.org) to detect clusters among all the 18 samples. Samples in unreasonable clusters were removed.
The DESeq2 package was used for DEGs detection between CLas-infected barks (B), leaves (L), and roots (R), and their corresponding controls (CY), namely, HLB_B vs. CY_B, HLB_L vs. CY_L, and HLB_R vs. CY_R. Before analysis, genes with several reads less than the total number of samples were filtered. Genes with a false discovery rate (FDR) <0.05 and log2 (Fold Change) [|log2(FoldChange)| >1] were considered as putative DEGs. Sequence similarity comparisons of the C. sinensis genome were performed with Arabidopsis thaliana using Blastx. Using C. sinensis as background, KOBAS 3.0 (Ai and Kong, 2018) was employed to annotate the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of DEGs.

Expression Normalization and Expression Pattern
Genes with zero TPM were filtered. The absolute gene expression, using normalized TPM value, was calculated according to Brawand et al. (2011). In detail, genes with a TPM value in the interquartile range (25-75%) of each sample were selected. Variance per gene across all samples was calculated to yield 1,000 least-varying genes and obtain median TPM values in each sample. Then, these median TPM values were used to derive the scaling factor of each sample, which was then utilized to normalize all the TPM values. The normalized TPM value matrix of all the genes was then used to estimate gene expression differences among all the tissues.
To investigate gene expression variations among the tissues, common DEGs identified in all three comparisons were filtered to estimate gene expression patterns among the tissues. Based on the normalized TPM values of these DEGs, the average TPM of each gene across biological replicates was presented as the expression level of this gene in the corresponding tissue, and the average was estimated by dividing the maximal value in all the tissues to generate relative expression profiles (Darbellay and Necsulea, 2020). Finally, the relative expression profiles were clustered with the K-means algorithm in R to investigate gene expression patterns among the barks, leaves, and roots.

Metabolome Profiling and Data Processing
The UPLC-QTOFMS analysis was carried out by Biotree (Shanghai, China). Quality control (QC) samples were prepared by mixing the same volume from all the samples. The UPLC-QTOFMS analysis was performed on 1290 Infinity LC System (Agilent, Sta. Clara, CA, United States) coupled to an Agilent 6538 Accurate-Mass Quadrupole Time-of-Flight (Q-TOF) (Agilent, Sta. Clara, CA, United States) mass spectrometer. The metabolomic data were collected in a centroid mode and the mass range was set at 50-1,100 m/z using an extended dynamic range.
The raw metabolomics data were converted into a normal file (.mzXML) using MSConvert (ProteoWizard 3.0). MZmine 2.1 (Pluskal et al., 2010) was applied for identifying features, deisotopes, feature alignments, and gap-filling to avoid missing some features in the first alignment algorithm. After removing the complexes, all adducts were searched against an internal retention time metabolite library. Data on both positive and negative ions were subjected to statistical analysis. Tables with metabolite peaks (mz/rt) of all the three tissues under both healthy and symptomatic conditions were formatted as commaseparated values (.csv). The resulting three-dimensional data of the peak number, sample name, and normalized peak area were fed to the SIMCA14.1 software package (V14.1; MKS Data Analytics Solutions, Umea, Sweden) for orthogonal projections to latent structures discriminant analysis (PLS-DA).
The DAMs were selected based on a statistically significant threshold of variable influence on projection (VIP) values acquired from the PLS-DA (VIP ≥1 and p ≤ 0.05). MetaboAnalyst was employed for KEGG enrichment analysis (Xia and Wishart, 2016). The DAVID Functional Annotation tool was applied to investigate the interactions of DEGs and DAMs in pathways (Huang et al., 2009).

Putative Marker Gene Identification for Symptomatic Tissues
Gene expression specificity across tissues was estimated using the equation: tau = n i=1 (1 − ri)/(n−1), where r i represents the ratio between the expression level (normed TPM value) of a gene in sample i and the maximal expression level across samples (Liao et al., 2006). The average tau value of all the biological replicates was calculated and used as the specific expression level for a given tissue. The tau value ranged from 0 to 1. A higher tau value indicated greater variation in expressional level across tissues, namely, higher tissue specificity.
To identify the putative markers for HLB infection, genes with expression levels (TPM) ≥2 in at least one sample of each tissue were selected. The putative marker genes should have a high expression specificity index (tau ≥ 0.85) in a specific healthy tissue and maximal expression in the corresponding symptomatic tissue (Liao et al., 2006). Due to the shortage of Chongyi samples, some putative marker genes were verified by qRT-PCR with the CLas-infected and healthy leaves from C. sinensis.

Confirmation of HLB and Anatomical Aberrations in CLas-Infected Tissues
After 18 weeks of grafting, all the CLas-inoculated Chongyi seedlings were positive by Nested PCR amplification (Supplementary Figure 1).
Subsequently, anatomical differences between the CLas-positive and healthy tissues were compared. More enlarged vascular bundles of lateral veinlets and air cavities were observed in the CLas-infected leaves (Figures 1A,B). Cells of phloem, xylem, parenchyma, and cortex arranged irregularly after CLas infection (Figure 1). Particularly, more sclereids and phloem fibers accumulated in the CLas-infected barks (Figures 1E,F). In the CLasinfected roots, the cortical tissue collapsed and separated from the epidermis (Figures 1G,H). Approximately 1 year after inoculation, the remaining CLas-positive Chongyi seedlings died out. Therefore, Chongyi was considered to be susceptible to HLB.

DEGs and DAMs Induced by CLas Inoculation
Around 83.34-86.06% of the clean reads could be mapped to the C. sinensis genome (Supplementary Table 1). Genes with a total read number over 18 were used for PCA analysis. The PCA results indicated that principal component 1 (PC1) could explain about 75% variances, and PC2 could explain about 20% variances (Figure 2A). All of the samples could be divided into three main clusters in PC1 according to the tissue type, and PC2 separated these samples depending on whether a sample was healthy or CLas-infected. An Hclust analysis with all of the expressed genes yielded similar results, i.e., there were three main clusters (leaves, roots, and barks), and the root samples were clustered with the bark samples ( Figure 2B). However, the two results both indicated that HLB_L1 was clustered with the healthy leaves and that HLB_B3 was clustered with the healthy barks. Therefore, HLB_L1 and HLB_B3 were excluded from the subsequent analysis.
For the metabolites, there were 2,541 and 2,003 peaks identified in the negative and positive ion modes. Of these peaks, 130 and 193 could be mapped as known metabolites, which included amino acids, organic acids, sugars, polyamines, nitrogenous compounds, and polyphenols. Three tissue-specific groups were clustered by the sPLS-DA, whereas PC1and PC2 explained 88 and 8%variances among all of the metabolomic samples, respectively ( Figure 2C). From the heatmap of all the detected metabolites (Supplementary Figure 2), more noticeable differences were observed among the CLas-infected and healthy tissues.
There were 3,628 DEGs in the leaves, 3,770 in the barks, and 1,716 in the roots, among which 1,888 in HLB_L, 2,778 in HLB_B, and 624 in HLB_R were upregulated, and 1,740 in CY_L, 992 in CY_B, and 792 in CY_R were upregulated ( Figure 2D). Metabolite intensities of the healthy and CLaspositive samples among the three tissues were compared. Subsequently, 62 DAMs in the barks (Table 1), 71 in the leaves ( Table 2), and 50 in the roots ( Table 3) could be matched to a corresponding KEGG compound ID. The majority of DAMs in the CLas-infected barks and leaves were increased, but half in the CLas-infected roots were decreased, suggesting that HLB affected roots more greatly at the metabolic level.
Similarly, more DAMs were enriched in "metabolic pathways" and "biosynthesis of secondary metabolites" in the three tissues ( Figure 3B; Supplementary Tables 8-10). In addition, 56 pathways in the barks (Supplementary

Common DEGs in the Three Tissues Exhibited Eight Expression Patterns
To compare expression variation among the three tissues, the expression patterns of common DEGs were analyzed. There were 304 common DEGs shared by the three tissues, while 1,884, 1,930, and 6,30 were specific to the leaves, barks, and roots, respectively (Supplementary Figure 3a). The clustering result from these 304 DEGs revealed that there was an increase in the sum of squares when the number of groups (K) was reduced from nine to eight (Supplementary Figure 3b). Thus, it might be proper to divide the 304 common DEGs into eight different expression patterns. As shown in Figure 4, DEGs in clusters 2 (N = 63), 3 (N = 43), 4 (N = 33), 6 (N = 4 2), and 7 (N = 56) were upregulated in the CLas-infected tissues, but the DEGs in clusters 1 (N = 29), 5 (N = 19), and 8 (N = 18) were downregulated.

CLas-Infected Tissues Had Higher Transcripts Per Kilobase Per Million (TPM) Values and More Genes That Reached the Maximal Expression (GME)
Expression levels of a total of 19,558 genes were normalized by scaling their TPM values, and the average of log2(normed TPM) of gene expression in biological replicates was represented as the expression level of this gene in the corresponding tissue. Median TPM values in the CLas-infected leaves, barks, and roots were 0.722, 0.657, and 0.485, respectively, and those in the corresponding control tissues were 0.661, 0.467, and 0.468. However, the TPM distribution of all the expressed genes was nearly the same among all the tissue samples ( Figure 5A). As for the percentage of GMEs in each tissue, 24, 18, and 15% of the genes reached their maximal expressions in the CLas-infected leaves, barks, and roots, respectively, which were higher than  those in the controls ( Figure 5B). Higher TPM values and more GMEs in the CLas-infected tissues suggested that HLB infection might cause an overall increase in transcript accumulation. However, GMEs in eight individual pathways varied ( Figure 5C). In the pathways related to "starch and sucrose metabolism, " "galactose metabolism, " "plant-pathogen interaction, " and "plant hormone signal transduction, " the CLas-infected tissues had more genes that reached their maximal expression levels than the healthy tissues, suggesting that HLB significantly affected these pathways, especially in the CLasinfected leaves and barks ( Figure 5C). Nevertheless, "carbon fixation in photosynthetic organism, " "photosynthesis-antenna protein, " "carotenoid biosynthesis" and "nitrogen metabolism" had more GMEs in the healthy tissues than in the CLas-infected tissues, indicating that HLB suppressed related physiological processes, particularly photosynthesis in leaves and carotenoid and nitrogen metabolism in roots ( Figure 5C).

DEGs Related to Early Signal Perception and Transduction
Plants possess an enormous number of cell-surface localized and intracellular receptors to perceive early immunogenic signals, containing pattern recognition receptors (PRRs, namely, leucine-rich repeat receptor-like proteins and kinases) and nucleotide-binding-leucine-rich-repeat (NLR) receptors consisting of toll/interleukin-1 receptor/resistance (TIR) and coiled-coil (CC) types (Monteiro and Nishimura, 2018;Wan et al., 2019). There were 41, 60, and 38 DEGs related to PRRs in HLB_B/L/R, respectively (Supplementary Table 14). The majority of PRRs were upregulated in HLB_B/R but were downregulated in HLB_L (Figure 6). Surprisingly, two well-known PRRs, flagellin-sensing 2 (FLS2) and EF-Tu receptor (EFR), were identified as DEGs in HLB_L/R but not in HLB_B (Supplementary Table 14). As for NLRs, only one CC-type NLR was found in HLB_B, but 30, 23, and 12 TIR-type NLRs were identified in HLBs_B/L/R, respectively (Supplementary Table 15). All the NLRs were annotated as disease resistance proteins whose transcriptional expressions were all nearly upregulated in HLB_B, but half of them were downregulated in HLB_L (Figure 6). The transcriptional expressions of other upstream immunity-related signaling protein kinases, such ascysteine-rich receptor-like kinases (CRKs), glutamate-like receptors (GLRs), mitogen-activated protein kinase kinase kinases (MAPKKKs), and mitogenactivated protein kinases (MAPKs), were also significantly changed in HLB_L/B/R but varied among the tissues (Figure 6). Calcium-and phytohormone-mediated signal transduction was mentioned in the following parts.
In the "plant hormone signal transduction" pathway, the number of GMEs increased from 7 to 21% in the barks and 15 to 26% in the leaves, but it decreased from 19 to 12% in the roots after CLas infection ( Figure 5C). There were 31, 42, and 10 DEGs in the barks, leaves and roots, respectively (Supplementary Table 17), including four DEGs mainly related to abscisic acid (ABA), 28 related to auxin, 11 related to cytokinin, five related to ethylene, three related to gibberellin (GA), eight related to jasmonic acid (JA), and seven related to salicylic acid (SA). Among these DEGs, there were genes encoding kinases, such as Arabidopsis histidine kinase 4 and ethylene receptor 2. Some DEGs encoded transcription factors, such as ethylene insensitive 3 and TGACG-binding factors (TGAs). Several non-expresser of pathogenesis-related genes (NPR) were also identified, such as NPR3 and NPR5. Interestingly, all the five DEGs related to ethylene were universally upregulated and exclusively found in the CLas-infected leaves. Although the DEGs related to SA were mainly detected in the CLas-infected leaves, five out of the six were downregulated. Moreover, all the DEGs related to JA were upregulated regardless of tissue type. Expressions of some DEGs involved in plant-pathogen interaction and plant hormone signal transduction are illustrated in Figure 6.

HLB Affected the Accumulation of Transcripts and Metabolites Related to Carbohydrate Metabolism
At the transcriptional level, a total of 74 DEGs in "starch and sucrose metabolism" and 23 in "galactose metabolism" were detected (Supplementary Table 18), and more GMEs in HLB_L, as well as HLB_B, were observed in these two pathways ( Figure 5C). Some of the DEGs were suggested to be involved in cell wall modification, and the majority of their expressions were upregulated in HLB_L/B (Figure 6). Integrated transcriptomics and metabolomics results revealed that these two pathways were enriched in both the leaves and the barks with the detection of related DEGs and DAMs, but that no related DAMs were detected in the roots (Figure 7). To be specific, significantly increased sucrose and maltose contents in HLB_L were accompanied by the downregulation of two cell wall invertase (CWINV), β-fructofuranosidase 1 (ATBFRUCT1) and CWINV4, and by the up-regulation of vacuolar invertase, At-β-Fructosidase4 (ATBETAFRUCT4), β-amylase1 (BAM1), and BAM3in "starch and sucrose metabolism". In the "galactose metabolism" of HLB_L, the contents of sucrose, D-mannose, D-myo-inositol, α-D-galactose, and raffinose were all elevated, and so were the transcripts of α-galactosidase1 (agaL1), hexokinase-like 1 (HKL1), seed imbibition protein 1/2 (SIP1/2, genes encoding raffinose synthase. In the barks, the upregulation of sucrose synthase2 (SUS2), SUS6, and agaL2 promoted sucrose production. The upregulation ofHKL1 and hexokinase 1 (HXK1) elevated the accumulation of α-D-galactose (Figure 7).

HLB Altered the Accumulation of Photosynthesis-Related Transcripts and Metabolites
Two photosynthesis-related KEGG pathways were enriched: "photosynthesis-antenna protein" and "carbon fixation in photosynthetic organism". The genes related to photosynthesisantenna protein all belonged to the light-harvesting chlorophyll a/b-binding protein (LHCP) family (Supplementary Table 19).
These 12 LHCPs only reached their maximum expression in the healthy leaves ( Figure 5C) and were universally downregulated in the CLas-infected barks (Figure 6). Because their fragments per kilobase of exon model per million mapped fragments in the CLas-infected leaves/barks/roots were 788, 211, and 1, respectively, these results indicated that photosynthesis mainly occurred in the leaves and was hampered by HLB. There were 18 DEGs enriched in "carbon fixation in photosynthetic organism" (Supplementary Table 19) and more GMEs in the healthy leaves and barks ( Figure 5C). Their expressions were downregulated in the CLas-infected tissues (Figure 6). Although the malate content in HLB_L was significantly higher and two NADP-dependent malic enzymes (NADP-ME3/4) were upregulated, the downregulation of phosphoenolpyruvate carboxylase 4(PPC4) and Pyruvate orthophcosphate dikinase (PPDK) may significantly reduce CO 2 fixation in the leaves (Figure 7). Therefore, photosynthesis was severely impaired in the CLas-infected leaves.

HLB Regulated DEGs and DAMs Involved in Carotenoid Biosynthesis and Nitrogen Metabolism
Carotenoids are essential for physiological processes in plants, such as photosynthesis and the production of phytohormones (Cazzonelli and Pogson, 2010). The number of GMEs in the healthy and CLas-infected barks was the same, but it decreased from 40 to 7% in the roots after infection ( Figure 5C). Moreover, DEGs encoding carotenoid cleavage dioxygenases (CCDs), phytoene synthase (PSY), Zeaxanthin epoxidase (ZEP), and the iron-binding protein D27 were all downregulated in HLB_R (Supplementary Table 20). Thus, carotenoid biosynthesis was negatively influenced by HLB in the roots (Figure 6).
Nitrogen metabolism was the only pathway with detection of both DEGs and DAMs in the three tissues (Figure 7). Even though there were no DEGs directly related to it, Lglutamine was over-accumulated in all the three CLas-infected tissues. Among the DEGs related to nitrogen metabolism (Supplementary Table 20), five high-affinity nitrate transporters (NRTs) were all downregulated in HLB_R (Figure 6), which might influence the absorption of nitrogen from the external environment to the plant.

Putative Marker Genes in Tissues for Diagnosis of HLB Infection
By scanning tau values among all the tissues, a total of 26 genes that had tau > 0.85 in the healthy tissues and reached their maximal expression in the corresponding CLas-infected tissues were considered as putative markers of HLB infection, including 13 genes in the leaves, 11 in the barks, and two in the roots (Supplementary Table 21).
Based on the expression patterns of these 26 putative markers, the samples were vertically divided into two main clusters: one containing CLas-infected roots and all the controls, and the other containing CLas-infected barks and leaves ( Figure 8A). Horizontally, cluster 1 contained two markers in Root; cluster 2 included eight markers in Bark 1; cluster 3 consisted of seven in Leaf 1, three in Bark 2, and six in Leaf 2. Therefore, genes in Root, Bark 1, and Leaf 1 could be considered as specific markers in the roots, barks and leaves, respectively. In addition, genes in Bark 2 and Leaf 2 could be putative gene makers in both CLas-infected barks and leaves due to their similar expression patterns.
Marker genes in Leaf 1, Bark 2, and Leaf 2 were verified by qRT-PCR using CLas-infected and healthy leaf samples from C. sinensis. Only eight putative marker genes could be consistently amplified, and their expression all increased in the CLas-infected samples ( Figure 8B).

Signal Perception and Propagation
Plants innately evolve a two-layered immune system, consisting of pattern-triggered immunity (PTI) and effector-triggered immunity (ETI) (Zhou and Zhang, 2020). The activation of the first layer of plant immunity PTI is initiated by intercellular receptors named PRRs, and the second layer, ETI, is conditioned by intracellular receptors named NLRs (Wan et al., 2019). NLRs contain TIR-type and CC-type NLRs depending on N-terminal singaling domains (Monteiro and Nishimura, 2018). In our study, many PRRs and NLRs were differentially expressed in the CLas-infected tissues (Supplementary Table 14; Figure 6), such as the two best characterized PRRs: flagellin-sensing 2 (FLS2) and EF-Tu receptor (EFR). Consistent with a recent report that PTI and ETI mutually potentiate to activate strong defenses against pathogens (Ngou et al., 2021), the majority of PRRs and NLRs in HLB_B/R were upregulated, but those in HLB_L were downregulated (Figure 6). A critical early signaling event connecting PTI and ETI is the production of reactive oxygen species (ROS) by NADPH oxidase respiratory burst oxidase homologs (RBOHs), such as RBOHB/D (Adachi et al., 2015;Yuan et al., 2021). However, the phosphorylation of WRKYs, essential components of both PTI and ETI, by MAPK is required for RBOHs-mediated ROS bursts (Adachi et al., 2015;Ngou et al., 2021). In addition, RBOHD could be directly activated by cysteine-rich receptor-like kinase 2 (CRK2) by phosphorylation, and the crk2 mutant was impaired in the full elicitor-induced ROS bursts as well as defense against pathogen attack (Kimura et al., 2020). In the barks, nearly all the above-mentioned genes were upregulated after CLas inoculation (Figure 6). Although FIGURE 4 | Expression patterns of common DEGs among all tissues. CY_B, CY_L, and CY_R represent healthy barks, leaves, and roots, respectively. HLB_B, HLB_Y, and HLB_R represent CLas-infected barks, leaves, and roots, respectively. Normalized TPM values were averaged across replicates for each tissue. Dots represent the average profiles of genes belonging to each cluster. Gray lines represent the profiles of individual genes from a cluster. the majority of CRKs (especially CRK2) were downregulated in HLB_L/R, the up-regulation of WRKYs, MAPKKKs, and MAPKs may compensate for the activation of RBOHB/C/D, whose transcripts were all significantly increased in these two tissues (Figure 6; Supplementary Table 16).
Most of the DEGs encoding other upstream immunity-related signaling protein kinases, glutamate-like receptors (GLRs), were upregulated in all the HLB tissues (Figure 6), which were suggested to be involved in various physiological pathways, such as Ca 2+ signaling (Grenzi et al., 2021). Ca 2+ is the most prominent second messenger that functions in a wide variety of environmental responses and developmental processes in plants (Lee and Seo, 2021). EF-hand calcium-binding proteins play a principal role in calcium signaling in plants, such as calmodulin (CaM), cyclic nucleotide-gated channel (CNGC), calmodulin-binding protein (CML), ATPase E1-E2 protein, and calcium-dependent protein kinase (CPK) (Mohanta et al., 2019). One CaM and three CNGCs were identified as DEGs in our study (Supplementary Table 16), which were reported to be critical in linking pathogen patterns to plant immunity (Tian et al., 2019). In total, 10 out of the 13 CMLs were upregulated (Supplementary Table 16), further proving their close relationship with plant defense (Lv et al., 2019). CPKs are calcium sensors and play important roles in ETI. For instance, CPK5 overexpressing Arabidopsis exhibited TN2dependent autoimmunity and enhanced disease resistance (Liu et al., 2017). Meanwhile, CPK5 directly phosphorylated RBOHD in vivo, which was critical for rapid defense signal transduction and ROS-mediated cell-to-cell communication (Dubiella et al., 2013), consistent with the upregulation of CPK5 and RBOHD in HLB_L/B (Supplementary Table 16). In addition, ATPase E1-E2 protein, whose transcript significantly accumulated in HLB_L/B, was involved in Ca 2+ -translocation (Hayter and Peterson, 2004) (Supplementary Tables 2, 3). Since more than half of the DEGs identified in the plant-pathogen interaction pathway were calcium-related (Supplementary Table 16), calcium-based defense signal propagation is probably a pivotal event in CLas infection.

Multi-Immune Responses
To prevent biotrophic pathogen invasion, plant cells would initiate immune responses, such as autophagy, the ubiquitinproteasome system, and programmed cell death (Greenberg and Yao, 2004;Üstün et al., 2018). There was one upregulated DEG annotated as autophagy-related protein (ATG) in each CLas-infected tissue, namely, ATG11 in the barks, ATG18H in the leaves and roots (Supplementary Tables 2-4). The ubiquitin-proteasome system contains various components, such as plant U-box (PUB) E3 ligases, FLS2, RBOHD (Trujillo, 2021), whose transcripts nearly all significantly accumulated in the CLas-infected tissues. These results suggested that HLB may trigger autophagy and the ubiquitin-proteasome system in Chongyi. In addition, transcripts of hypersensitiveinduced reaction (HIR)1 and HIR2 significantly accumulated (Supplementary Table 2), which were reported to participate in the hypersensitive response (a programmed cell death phenomenon) and positively contributed to basal resistance against pathogen attack in the enhanced disease susceptibility 1 (EDS1)-and salicylic acid (SA)-dependent pathways (Qi et al., 2011;Li et al., 2019). Meanwhile, HIRs could contribute to ETI by forming complexes with a membrane-associated disease resistance protein RPS2 (Qi et al., 2011), whose transcriptional expression was also upregulated in the CLasinfected barks (Supplementary Table 2). Thus, multi-immune responses were triggered. FIGURE 6 | Summary of defense-related DEGs in CLas-infected leaves, barks, and roots. Red, upregulation; blue, downregulation. PRRs, pattern recognition receptors; RBOHs, plant respiratory burst oxidase homologs; ROS, reactive oxygen species; MAPKKKs, mitogen-activated protein kinase kinase kinases; MAPKs, mitogen-activated protein kinases; CMLs, calmodulin-like proteins; CPKs, calcium-dependent protein kinases; CNGCs, cyclic nucleotide-gated channels; CRKs, cysteine-rich receptor-like kinases; (Continued) FIGURE 6 | GLRs, glutamate-like receptors; NLRs, nucleotide-binding-leucine-rich-repeat receptors; ETI, effector-triggered immunity. MAPKKKs, MAPKs, CPKs, CRKs, and GLRs are immunity-related signaling protein kinases. These DEGs are listed in Supplementary Tables 16-20. It is to be noted that the phytohormone-related DEGs presented in this Figure are involved in signaling instead of synthesis or degradation.
The EDS1 family of immunity regulators contains three proteins: EDS1, phytoalexin deficient 4 (PAD4), and senescenceassociated gene 101 (SAG101), which could form mutually exclusive EDS1-SAG101 and EDS1-PAD4 heterodimers (Ngou et al., 2021). EDS1-SAG101 dimers have a helper NLR to mediate cell death and pathogen resistance (Lapin et al., 2019). EDS1-PAD4 functions as an integrator of cell surface and intracellular receptor signaling, which delivers signals from TIR-NLRs to transcriptional defense and host cell death (Dongus and Parker, 2021). Meanwhile, the EDS1-PAD4 module participates in phytohormone-mediated pathogen containment by inhibiting MYC2, a negative regulator of SA biosynthesis, and enhances the SA defense sector (Cui et al., 2018). Therefore, EDS1 is considered an immune regulatory hub (Lapin et al., 2020). In our study, EDS1 and PAD4 were both upregulated in the CLasinfected leaves and barks. Although SAG101 was not found, its two homologs were identified as DEGs, and the upregulation of SAG20 and SAG21 was consistent with their induced expression upon pathogen attack (Keates et al., 2003;Salleh et al., 2012).

Phytohormone Signal Pathways and Defense Response
Phytohormones, such as SA, jasmonic acid (JA), ethylene (ET), abscisic acid (ABA), auxin, and gibberellins (GAs), participate in defense responses (Berens et al., 2017). Among them, SA and JA are major defense-related phytohormones, while the other phytohormones participate in hormone signaling mainly through interactions with SA and JA (Pieterse et al., 2009). The transcriptional levels of DEGs involved in plant hormone signal transduction varied among the CLas-infected tissues, especially those related to SA (Supplementary Table 17; Figure 6). Nonexpressor of pathogenesis-related genes 1 (NPR1), the main protein required for SA perception, interacts with TGACGbinding (TGA) transcription factors, which modulate SA biosynthesis positively (Budimir et al., 2021). Although NPR1 was not found in our study, Blade-on-petiole 2 (BOP2, a paralog of NPR1) and TGA1/7/9 were all downregulated in the CLasinfected leaves. Meanwhile, Cs2g10790 (NPR3), whose protein transcriptionally repressed genes involved in promoting SAmediated basal immunity, such as systemic acquired resistance deficient 1 (SARD1) and WRKY70, was upregulated in the CLasinfected leaves and barks (Supplementary Table 17) (Ding et al., 2018). However, SARD1 and over 80% of the WRKYs were upregulated in HLB_L/B, probably because of the compensatory mechanism(s) to circumvent weakened SA signaling, such as MAPK signaling (Tsuda et al., 2013). In addition, three DEGs annotated as phenylalanine ammonium-lyase 1 (PAL1), FIGURE 7 | KEGG pathways with DEGs and DAMs in leaves, bark, and roots. In the rectangular box, pathways are represented in bold black, and black numbers represent the KEGG ID of enzymes. Upregulated genes are in red bold italicized words, while downregulated genes are in blue bold italicized words. The red bold represents upregulated metabolites, while the blue bold represents downregulated metabolites. Gray represents metabolites that were not detected.
encoding the key enzyme of the PAL pathway contributing to SA biosynthesis in immunity (Berens et al., 2017), were upregulated in HLB_B (Supplementary Table 2). Therefore, the above information indicated that SA-mediated defense response might be partially depressed in the leaves but not in the barks after CLas inoculation.
The pathogen of citrus canker appeared to promote JA accumulation through hijacking host ABA production, which could suppress SA-mediated defenses and benefit pathogen development (Long et al., 2019). Consistently, DEGs involved in JA, ET, and ABA signaling were nearly all upregulated in the CLas-infected leaves and barks (Figure 6), and so were those participating in JA biosynthesis, such as allene oxide cyclase 4 (AOC4), lipoxygenase 3(LOX3), 12-oxophytodienoate reductase 2(OPR2), and peroxisomal acyl-coenzyme A oxidase 1 (ACX1, Supplementary Tables 2, 3). Hence, JA signaling was probably dominated in the CLas-infected leaves and barks, which may suppress SA-mediated pathogen defense.

Alteration of Carbohydrate Metabolism and Cell Wall Composition
Huanglongbing severely influences functions related to carbohydrate biosynthesis and metabolism . The swelling of the middle lamella among cell walls surrounding sieve elements would drive starch accumulation in the leaves and phloem, and decrease starch transport to roots (Folimonova and Achor, 2010). In this study, the "starch and sucrose metabolism" and "galactose metabolism" pathways were upregulated in the CLas-infected leaves and barks but were down-regulated in the roots (Figure 5C). Consistently, the metabolomics results revealed that sucrose accumulation occurred in both the CLasinfected leaves and barks, but that no DAMs were detected in the roots (Figure 7), which might result in low carbohydrate storage in the roots.
Interestingly, many carbohydrate-related DEGs were annotated to participate in cell wall modification and remodeling. Cellulose, various hemicelluloses, and pectins are the three major types of polysaccharides in the cell wall of Arabidopsis leaves (Liepman et al., 2010). CEL1, an endo-1,4-β-glucanase, was required for cellulose formation of the cell wall during cell growth and expansion (Tsabary et al., 2003), and its transcripts significantly accumulated in the CLas-infected barks and roots (Supplementary Table 18). UDP-D-glucuronate 4-epimerases (GAEs) provide important precursors for the synthesis of pectins, which are the major component of the plant primary cell wall and critical for its integrity as well as immunity (Mølhøj et al., 2004;Bethke et al., 2016). GAE1/6 were depressed upon Botrytis cinerea and Pma ES4326 challenge in Arabidopsis (Bethke et al., 2016), but they were induced in our study (Supplementary Table 18). Meanwhile, several DEGs encoding pectin lyase-like superfamily proteins, which were reported to participate in plant cell wall modification (Hocq et al., 2020), were also identified FIGURE 9 | Hypothesized model for the whole plant and tissue-specific HLB susceptible response of Chongyi. Red arrows indicate induction, and blue ones indicate repression. Table 18). Pectinesterases (PMEs) contribute to cell wall integrity (Bethke et al., 2014), which may explain the upregulation of 11 PMEs (Supplementary Table 18).

(Supplementary
UDP-glucose 4-epimerase (UGE) provides UDP-galactose for cell wall biosynthesis (Rösti et al., 2007). UDP-glucose 6dehydrogenases 2 (UGD2) and UGD3 were required for the formation of the cell wall in growth (Reboul et al., 2011). Inferred from electronic annotation on the website of UniProtKB (https:// www.uniprot.org/uniprot/), hexosyltransferases (GAUTs) may be involved in pectin and/or xylans biosynthesis in cell walls, and alpha-galactosidase 1/2 (AGAL1/2) may function in cell wall loosening and expansion. The expression of UGE1, UGD2, UGD3, GAUTs, and AGAL1/2 was all significantly upregulated (Supplementary Table 18). Meanwhile, trehalose was suggested to play a role in cell wall modification (Bae et al., 2005). Its biosynthesis in plants needs two indispensable enzymes: trehalose-phosphate synthase (TPS) and trehalose-6-phosphate phosphatase (TPP). The tps1 mutant in Arabidopsis showed an altered cell wall structure and starch accumulation (Gómez et al., 2006), suggesting that the downregulation of two TPS1 genes in the CLas-infected leaves may also cause starch accumulation. Therefore, cell wall was highly affected in Chongyi after CLas infection, which was considered as one of the typical characteristics of HLB-susceptible genotypes (Curtolo et al., 2020).

HLB Negatively Affects Carotenoid and Nitrogen Metabolism in Root
Carotenoids play essential roles in plant physiological processes, and phytoene synthase (PSY) is considered the most important regulatory enzyme in carotenogenesis (Cazzonelli and Pogson, 2010). PSY was suppressed in the CLas-infected roots (Supplementary Table 20), and so was zeaxanthin epoxidase (ZEP), encoding the other enzyme of carotenoid biosynthesis. Carotenoid cleavage dioxygenases (CCDs) can catalyze carotenoids into various carotenoid cleavage products and provide substrates for phytohormone biosynthesis (Gao et al., 2021). CCD7 and CCD8 were all downregulated in the roots after CLas infection (Supplementary Table 20), probably because the suppression of carotenoid biosynthesis resulted in the reduction of degradation. Interestingly, strigolactones are carotenoid-derived phytohormones, whose biosynthesis requires CCD7 and CCD8 as well as the iron-binding protein D27 (Alder et al., 2012;Wang et al., 2020). D27 was downregulated in CLas-infected roots, too. Consequently, HLB may repress strigolactone production in roots. However, root secrets strigolactones into the rhizosphere as a pre-symbiotic signal for arbuscular mycorrhizal symbiosis, which is critical for Citrus to obtain essential nutrients from the soil (Akiyama et al., 2005;An et al., 2018). Therefore, the suppression of carotenoid metabolism may cause a nutrient deficiency in roots, such as nitrogen (N) deficiency. Consistently, high-affinity nitrate transports (NRTs) 2.4, NRT2.5, NRT2.6, and NRT2.7, which could reduce the absorption of N from soil to roots, were all downregulated (Figure 7; Supplementary Table 20). Thus, the suppression of carotenoid and nitrogen metabolism may accelerate the deterioration of the CLas-infected roots and weaken the defense response of the whole plant.

A Hypothesized Model for HLB Susceptibility Mechanism in Chongyi Wild Mandarin
The comparison of transcriptomic profiles from the CLassusceptible citrus genotypes and CLas-tolerant/resistant genotypes revealed the transcriptional alterations of susceptible plants, such as downregulation of receptors and WKRYs, induction of GA synthesis, and suppression of GA degradation (Curtolo et al., 2020). As mentioned above, the majority of receptors were depressed in the CLas-infected leaves (Figure 6). Although most WRKYs were upregulated in our study (Figure 6), WRKY70, which was suggested to play an important role in Poncirus trifoliata tolerance to HLB, was excluded (Peng et al., 2020). Two gibberellin 2-beta-dioxygenase 8 genes in the CLasinfected leaves and one, which was involved in GA degradation, in the barks were downregulated (Supplementary Tables 2, 3) . No GA synthesis-related DEGs were identified, but several Gibberellin stimulated Arabidopsis genes (GASAs) were upregulated in all the CLas-infected tissues (Supplementary Tables 2-4). GASAs were GA-induced and could promote GA response (Rubinovich and Weiss, 2010), suggesting that GA synthesis was probably induced in CLasinfected Chongyi. Thus, these results were consistent with the characteristics of susceptible Citrus genotypes, as summarized by Curtolo et al. (2020). Furthermore, CLas-enhanced JA synthesis and signaling may hamper defense response by antagonizing SA, making Chongyi more vulnerable to HLB. Based on the above arguments, along with severe anatomical aberrations in the CLas-infected tissues, a hypothesized model of the whole plant and tissue-specific HLB susceptible response of Chongyi is illustrated in Figure 9.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
TP and B-LZ conceived the study. X-TX, F-TC, and X-JZ conducted the management of grafting and seedling. TP and J-LK contributed to data analysis and writing of the manuscript. W-SD, MW, Z-YL, and H-NS significantly contributed to the writing of the manuscript. All authors discussed the results and commented on the manuscript.