The D. citri specimens used in this study were collected from Sun Yat-sen University, Guangzhou, China, in 2014. The strain was reared on Citrus reticulata Blanco cv. Shatangju in a greenhouse under a 14 h: 10 h light: dark cycle at 26 ± 1°C and 65–70% relative humidity for more than 25 generations. The newly emerged female adults were fed for 14 days on the citrus plants infected with HLB, and they were confirmed to carry CLas through sampling detection. Meanwhile, the female adults that was continuously fed on healthy citrus plants for 14 days as control.

Total RNA (1 μg) was reverse-transcribed into first-strand cDNA using a PrimeScript^TM RT reagent kit (Takara Bio, Inc., Otsu, Shiga, Japan). Quantitative PCR (qPCR) was performed using a 10 μL reaction containing 1 μL cDNA, 0.3 μL each of 10 μmol⋅L^–1 forward and reverse primers (Supplementary Table S1), and 5 μL SYBR FAST Universal qPCR mix (KAPA Biosystems, Woburn, MA, United States) in a LightCycler 480 system (Roche Diagnostics GmbH, Mannheim, Germany), with the following amplification conditions: 5 min at 95°C, followed by 45 cycles at 95°C for 10 s, 60°C for 20 s, and 72°C for 20 s. The qPCR experiments were performed for each sample using three biological and three technical replicates. The expression levels of selected genes were normalized to the expression levels of D. citri actin-1. The differential gene expression was calculated using the 2^{–△ △ Ct} method (Livak and Schmittgen, 2001), and the results were expressed as mean ± SE. The differences in gene expression levels between infected and uninfected D. citri were analyzed using t-tests in the SPSS 18.0 statistical software (P < 0.01).

All 20 insect individuals were pooled and ground in liquid nitrogen. Total RNA was extracted using a total RNA extract kit (Omega, Norcross, GA, United States) according to the manufacturer’s instructions. The cDNA libraries from triplicates of each insect line were generated using the Illumina TruSeq Stranded mRNA kit (Illumina, San Diego, CA, Untied States) and were size-selected to an average fragment size of 450 bp. The libraries were sequenced on the Illumina HiSeq Xten platform producing paired-end (PE) 150 bp reads.

Low quality reads from Illumina sequencing were filtered out using Cutadapt software v1.16 with parameters of—discard-trimmed-n3-e0.1 (Martin, 2011). The obtained clean reads were mapped to the D. citri genome (NCBI accession No. GCA_000475195.1) with the Hisat2 software (Kim et al., 2015). Read counts for all genes were extracted with HTSeq-count (Anders et al., 2015) and normalized using the R package DESeq2 (Love et al., 2014). The differentially expressed genes (DEGs) between the CLas-infected and non-infected lines were estimated by DESeq2 according to the threshold of | log2 ratio| > 1 and adjusted P < 0.05 (BH adjustment). The Blast2GO software was used to generate GO annotation of the DEGs. The DEGs were further annotated in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database by KEGG Automatic Annotation Server. The GO functional enrichment and KEGG pathway enrichment of the DEGs were analyzed using the Fisher’s exact test with Benjamini and Hochberg (BH) adjustment.

Twelve genes were selected for qPCR validation based on the differentially expressed level and specific functions (Supplementary Table S1). Following the user manual, the first cDNA strand fragments were synthesized from total RNA using the PrimeScript^TM RT Master Mix kit (Takara, Japan). Specific primers for selected genes are listed in Supplementary Table S1. qRT-PCR was performed on an ABI7500 real-time fluorescent quantitative PCR instrument (Applied Biosystems, Foster City, United States). The relative mRNA levels of the selected genes were estimated according to the aforementioned method.

For metabolite extraction, 20 female adults from each insect line were pooled and surface washes of sterile water for 1 min were performed three times. Then the samples were ground to a powder form with an electronic grinding machine. A total of 20 mg of ground sample was placed in an EP tube, and a 1,000 μL extract solution (acetonitrile: methanol: water = 2: 2: 1, with isotopically labelled internal standard mixture) was added. After 30 s of vortex mixing, the samples were homogenized at 35 Hz for 4 min and sonicated for 5 min in ice-water bath. The homogenization and sonication cycle were repeated three times. Then, the samples were incubated for 1 h at −40°C and centrifuged at 12,000 rpm for 15 min at 4°C. The resulting supernatant was transferred to a fresh glass vial for analysis. The quality control sample was prepared by mixing an equal aliquot of supernatant from all the samples. Six biological replications were performed for each insect line.

LC-MS/MS analyses were performed using an UHPLC system (Vanquish, Thermo Fisher Scientific) with a UPLC BEH Amide column (2.1 mm × 100 mm, 1.7 μm) coupled to a Q Exactive HFX mass spectrometer (Orbitrap MS, Thermo Fisher Scientific). The mobile phase consisted of 25 mmol/L ammonium acetate and 25 ammonia hydroxide in water (pH = 9.75) (A) and acetonitrile (B). The analysis was carried with elution gradient as follows: 0–0.5 min, 95% B; 0.5–7.0 min, 95–65% B; 7.0–8.0 min, 65–40% B; 8.0–9.0 min, 40% B; 9.0–9.1 min, 40–95% B; 9.1–12.0 min, 95% B. Column temperature was 30°C. Auto-sampler temperature was 4°C, and the injection volume was 2 μL. The QE HFX mass spectrometer was used for its ability to acquire MS/MS spectra on information-dependent acquisition mode in the control of the acquisition software (Xcalibur, Thermo Fisher Scientific). In this mode, the acquisition software continuously evaluates the full scan MS spectrum. The electrospray ionization (ESI) source conditions were set as follows: sheath gas flow rate 50 Arb, Aux gas flow rate 10 Arb, capillary temperature 320°C, full MS resolution 60,000, MS/MS resolution 7,500, collision energy 10/30/60 in NCE mode, spray voltage 3.5 kV (positive) or −3.2 kV (negative).

For peak detection, extraction, alignment, and integration, the raw MS/MS spectra data were converted to the mzXML format using ProteoWizard and processed with an in-house program, developed using R and based on XCMS. Then, an in-house MS2 database (BiotreeDB) was used for metabolite annotation, with a cutoff value of 0.3.

The difference in metabolite contents between the two groups was calculated using analysis of variance (ANOVA) and the Mann-Whitney U-test in SPSS 13.0. Statistical significance was set at P < 0.05. MetaboAnalyst¹ was further used for the hierarchical cluster analysis of all samples and overlapped metabolites. PCA and partial least squares-discriminant analysis (PLS-DA) were conducted to investigate the relationships among the test samples. For the pathway enrichment analysis, the differential metabolites between tested groups were assigned to metabolic pathways using the tool MetPA in MetaboAnalyst. The −log(P) and impact values of each metabolic pathway were calculated using the hypergeometric test.

The transcriptome and metabolome pathways that overlapped were considered, and the corresponding KEGG xml files (KGML) were downloaded from the KEGG pathway database. The interactions among these pathways were estimated using a custom Perl script and the KGML files as input files (Chen et al., 2016; Pang et al., 2018). Finally, the pathway–pathway network was visualized using the Cytoscape software (V.3.6.0).

The CLas-specific fragments had 1,167 bp, and the gel imaging results showed that D. citri carried CLas after feeding on CLas-infected plants for 14 days (Figure 1A). The qPCR test also confirmed the CLas infection in the infected insect line, which showed a much higher level of CLas mRNA than the control non-infected line (Figure 1B). Thus, the target insect lines could be used for the subsequent transcriptomic and metabolomic analyses.

The RNA-seq experiments yielded an average 42,029,244 raw paired-end reads (range from 39,283,042 to 44,690,980) for each sample with an average Q30 value of 93.94% (Supplementary Table S2). After quality filtering, we obtained an average 39,538,101 high-quality clean reads for each sample. The average mapping rate of the clean reads to the D. citri genome was 81.91% (range from 81.10 to 82.66%), and the mapped reads were used for subsequent analysis.

The correlation test and principal component analysis (PCA) results were consistent regarding the triplicate RNA sequencing (RNA-seq) datasets of CLas-infected and non-infected insect lines (Supplementary Figure S1). The comparative analysis revealed that CLas infection resulted in 1,194 DEGs (4.97% of the total detected genes in the D. citri genome), including 662 up- and 532 downregulated genes (Figure 2A and Supplementary Table S3). Twelve randomly selected DEGs were validated by qPCR (Figure 2B), suggesting the reliability of the RNA-seq results. Among the upregulated genes, we found eight vitellogenin genes and two hemocyanin genes (Supplementary Table S2); the upregulation of these proteins has been observed in previous proteomic studies (Ramsey et al., 2017; Kruse et al., 2018), and our findings further confirmed this at a transcriptomic level.

DEGs were mostly enriched in gene ontology (GO) functional categories of carboxylic acid metabolic process, organic acid metabolic process, and cofactor binding (Figure 3A and Supplementary Table S3). However, the functional enrichment of upregulated and downregulated genes was largely different. For example, the downregulated genes were mostly related to functions of photoreceptor activity, protein-chromophore linkage, and hydroxymethyl-, formyl-, and related transferase activity, whereas the upregulated genes were preferentially enriched in the functions of lipid transporter activity and lipid transporter activity (Supplementary Table S4).

We further analyzed the pathways affected by the DEGs. The analysis showed the enrichment of DEGs in multiple carbohydrate metabolism pathways, including the citrate cycle (TCA cycle), ascorbate and aldarate metabolism, pentose, and glucuronate interconversions, and pyruvate metabolism (Figure 3B and Supplementary Table S5). This indicated that CLas infection may affect the carbon utilization of D. citri. Moreover, two immune system-related pathways, antigen processing, and presentation and NOD-like receptor signaling pathway, were found to be significantly enriched with DEGs, suggesting that CLas inevitably modulated D. citri immunity. Additionally, we noticed that several cofactor and vitamin metabolism pathways were significantly affected by the DEGs. This is in accordance with the GO enrichment in the functional category of cofactor binding. The affected cofactors and vitamins included folate, riboflavin, and retinol.

The metabolome of insects from CLas-infected and non-infected lines was measured using LC-MS/MS, which identified a total of 315 metabolites across all samples (Supplementary Table S6). The PCA showed that almost all samples were located within the 95% confidence interval, and these samples were distinctly separated into two groups according to their grouping situation (Figure 4A). The substantial difference between CLas-infected and non-infected insect lines was further confirmed by OPLS-DA score plots (Supplementary Figure S2). Of all detected metabolites, 105 showed significant difference in abundance between CLas-infected and non-infected lines (Mann-Whitney U-test, P < 0.05, Figure 4B). This result revealed that infection with CLas led to higher level of steroids (ergosterol), flavonoids (castavinol), and nucleosides such as cytosine, inosine, 5′-methylthioadenosine and 2-methylguanosine. In addition, we observed the upregulation of fatty acid molecules (linoleic acid, polyoxyethylene dioleate, 2-methylbutyroylcarnitine, L-acetylcarnitine, butyrylcarnitine, and polyoxyethylene dioleate) but the downregulation of fatty amides (linoleamide, palmitic amide, and oleamide) in CLas-infected insects (Figure 4B).

The metabolic pathway analysis revealed that the pathways of vitamin B6 metabolism, glycerophospholipid metabolism, and purine metabolism were largely affected by the differentially expressed metabolites under CLas infection (Figure 4C and Supplementary Table S7). Notably, most of the detected metabolites in these pathways were significantly upregulated, indicating the potentially positive effect of CLas in the metabolism of these compounds.

We further determined the gene regulatory network associated with CLas infection using an integrated analysis of the transcriptome and metabolome. Nine pathways were found to be affected by both differentially expressed genes and metabolites (Supplementary Table S8), and these pathways could be linked to each other through differential factors except the tryptophan metabolism pathway (Figure 5). The glycine, serine, and threonine metabolism contained the most differential factors, and this pathway interacted with those of the cysteine and methionine metabolism, glycerophospholipid metabolism, sphingolipid metabolism, and purine metabolism. The purine metabolism pathway also harbored multiple differential factors, most of which were significantly upregulated under CLas infection.

Among all related differential factors, the downregulated gene LOC103515380 (branched-chain-amino-acid aminotransferase) was involved in both pantothenate and CoA biosynthesis and cysteine and methionine metabolism pathways, indicating its reduced function in the regulatory network under CLas infection. In addition, we observed increased levels of the metabolite phosphatidylcholine, an intermediate product of the glycerophospholipid metabolism and the origin of the arachidonic acid metabolism (Supplementary Figure S3). Overall, the network analysis showed that CLas infection mainly impacted the metabolism of amino acids, lipids, and cofactors of D. citri.