Introduction
The citrus greening disease, also known as Huanglongbing (HLB), caused by the phloem-limited, Gram-negative causative bacterium ‘Candidatus Liberibacter asiaticus’ (CLas), is the most destructive disease of citrus plants worldwide (Hall et al., Reference Hall, Richardson, Ammar and Halbert2013). It causes tremendous economic losses each year in most citrus-growing regions by reducing fruit yield and quality, promoting tree decline, and eventually results in plant death (Slisz et al., Reference Slisz, Breksa, Mishchuk, McCollum and Slupsky2012; Munir et al., Reference Munir, He, Wu, He, Khan, Huang, Cui, He and He2018). Thus HLB has become the most serious threat to the citrus industry for decades worldwide, and the CLas bacterium is currently unculturable, there are no effective prevention and control methods once established (da Graca et al., Reference da Graca, Douhan, Halbert, Keremane, Lee, Vidalakis and Zhao2016).
The Asian citrus psyllid (ACP), Diaphorina citri Kuwayama, a phloem-sucking pest of citrus, is considered as the major transmitting vector of CLas bacterium, the causal agent of HLB (George et al., Reference George, Ammar, Hall, Shatters and Lapointe2018). The psyllid transmits HLB pathogen in a persistent propagative manner, and both nymphs and adults can transmit CLas once they acquired CLas at their nymphal stages (Canale et al., Reference Canale, Tomaseto, Haddad, Della Coletta and Lopes2017; Andrade et al., Reference Andrade, Pang, Achor, Wang, Yao, Singer and Wang2020). The effective control of ACP is recognised as a key measure to prevent HLB spreading, and it is crucial to investigate novel targets and develop new strategies for ACP management (Vanaclocha et al., Reference Vanaclocha, Jones, Tansey, Monzo, Chen and Stansly2019).
Like other phloem-sucking insects, D. citri feed on the phloem sap of citrus plants; both nymphs and adults puncture the epidermis using a specialised piercing-sucking style, and continuously secrete saliva into plant tissues before ingesting juice from the phloem (George et al., Reference George, Ammar, Hall and Lapointe2017). Saliva is primarily derived from salivary glands (SG) and contains a variety of bioactive components with function in food ingestion, digestion, lubrication, tissue penetration, detoxification, and suppressing host plant defences (Perilla-Henao and Casteel, Reference Perilla-Henao and Casteel2016; Aljbory and Chen, Reference Aljbory and Chen2018). SG, secreting gelling and watery saliva, act as reservoirs and facilitate the transmission of plant pathogens by inoculating saliva into healthy plants via stylet during insect feeding. Many genes and proteins in the SG are associated with the essential biological processes of insect host that facilitate feeding behaviour or repress host immune responses. Several SG-secreted proteins acting as effector proteins have been reported, such as MED, SP75, Bt56, and vitellogenin (Yang et al., Reference Yang, Guo, Chu, Ding, Wei, Cheng and Wan2017; Matsumoto and Hattori, Reference Matsumoto and Hattori2018; Xu et al., Reference Xu, Qian, Wang, Shao, Hong, Liu and Wang2019; Ji et al., Reference Ji, Fu, Shi, Li, Jing, Wang, Yang, Tian, Wang, Ju, Guo, Liu, Dou, Hoffmann, Zhu-Salzman and Fang2021). Additionally, in several vector-borne hemipterans, like Acyrthosiphon pisum Harris, Bemisia tabaci Gennadius, Nilaparvata lugens Stal, and D. citri (Carolan et al., Reference Carolan, Caragea, Reardon, Mutti, Dittmer, Pappan, Cui, Castaneto, Poulain, Dossat, Tagu, Reese, Reeck, Wilkinson and Edwards2011; Su et al., Reference Su, Li, Li, Luan, Ye, Wang and Liu2012; Ji et al., Reference Ji, Yu, Fu, Chen, Ye, Li and Lou2013; Yu and Killiny, Reference Yu and Killiny2018), SG, and secreted saliva play critical roles mediating interactions among insects, insect-borne pathogens, and their host plants (Perilla-Henao and Casteel, Reference Perilla-Henao and Casteel2016; van Bell and Will, Reference van Bell and Will2016).
The transmission of HLB pathogen by D. citri is associated with the properties of SG during insect feeding. Both D. citri nymphs and adults can transmit CLas by acquiring the pathogen, and the CLas-infected psyllids have the ability to transmit HLB pathogen throughout their lives. CLas is transmitted by D. citri more efficiently when it is acquired by nymphs rather than adults (Mann et al., Reference Mann, Fattah-Hosseini, Ammar, Stange, Warrick, Sturgeon, Shatters and Heck2018). Prior studies have reveal that the newly produced pathogen can be detected in various internal organs of CLas-infected D. citri, including the alimentary canal, haemolymph, filter chamber, midgut, SG, Malpighian tubules, testis, ovaries, muscle, and fat tissues (Ammar et al., Reference Ammar, Shatters and Hall2011a; Coletta et al., Reference Coletta, Daugherty, Ferreira and Lopes2014). Among these tissues, the proportion of infected SG was significantly lower than that of other tissues in the insect body, suggesting that the SG may serve as an important barrier to CLas transmission in the psyllid vector (Ammar et al., Reference Ammar, Shatters, Lynch and Hall2011b). After the acquisition, CLas is initially reproduced in the alimentary canal, and subsequently migrates from midgut to haemolymph and other tissues (Chen et al., Reference Chen, Li, Liu, Chi, Jia and Wei2022). It eventually moves into SG to proliferate, from where newly produced pathogens will be injected into healthy plant tissues along with saliva when sucking juice from the phloem (Hall et al., Reference Hall, Richardson, Ammar and Halbert2013; Wu et al., Reference Wu, Luo, Xu, Wu, Qureshi and Cen2016). In this process, the SG may act as the last, and most important, barrier that the pathogen must breach before entry into new plants, suggesting that the SG could play an important role in the transmission of HLB pathogen.
In the last decade, high-throughput RNA-Seq technology has been widely employed to elucidate the molecular mechanisms of insect–pathogen interactions, especially in the research of host cell responses to exogenous pathogenic infection. By comparative transcriptomic analysis, Kaur et al. identified a large number of differentially expressed genes (DEGs) associated with metabolic pathways, signal transduction, transport and catabolism, receptors, glucose transporters, α-glucosidases and the uric acid pathway between B. tabaci fed on Tomato chlorosis virus (ToCV) infected and uninfected tomato plants (Kaur et al., Reference Kaur, Chen, Zheng, Hasegawa, Ling, Fei and Wintermantel2017). By RNA-Seq analysis, Kruse et al. investigated the interactions between CLas and the ACP gut, and CLas exposure resulted in changes in pathways involving the tricarboxylic acid (TCA) cycle, iron metabolism, insecticide resistance, and the insect immune system (Kruse et al., Reference Kruse, Fattah-Hosseini, Saha, Johnson, Warwick, Sturgeon, Mueller, MacCoss, Shatters and Heck2017). Yu et al. also performed a comparative transcriptome-wide analysis of D. citri midgut upon CLas infection, and a total of 778 DEGs involved in ubiquitination, the immune response, the ribosome, endocytosis, the cytoskeleton, and insecticide resistance were identified (Yu et al., Reference Yu, Li, Zeng, Song, Yu, Su, Chen, Yi and Lu2020).
Additionally, the functions and potential roles of the SG in transmission of plant pathogens have been studied in several other hemipteran insects by RNA-Seq and transcriptomic analysis, such as A. pisum (Carolan et al., Reference Carolan, Caragea, Reardon, Mutti, Dittmer, Pappan, Cui, Castaneto, Poulain, Dossat, Tagu, Reese, Reeck, Wilkinson and Edwards2011), B. tabaci (Su et al., Reference Su, Li, Li, Luan, Ye, Wang and Liu2012), Sitobion avenae, Empoasca fabae Harris and N. lugens (DeLay et al., Reference DeLay, Mamidala, Wijeratne, Wijeratne, Mittapalli, Wang and Lamp2012; Zhang et al., Reference Zhang, Fan, Sun, Francis and Chen2017; Rao et al., Reference Rao, Zheng, Liu, Guo, Guo, Wu, Shangguan, Wang, Wu, Wang, Hu, Xu, Jiang, Huang, Shi and He2019). For ACP insect, SG play an important role in HLB pathogen colonisation and transmission. Although a lot of differential genes and proteins in D. citri saliva and SG have been identified (Yu and Killiny, Reference Yu and Killiny2018; Liu et al., Reference Liu, He, Guan, Zhong, Pang and Han2021), the functions and properties of these differentially expressed SG components in CLas-free and CLas-infected psyllid are still limited. The mechanism underlying CLas infection in the SG of D. citri remains to be addressed.
In the current study, a comparative transcriptomic analysis was performed between the SG in CLas-infected and uninfected ACP using an Illumina RNA-Seq technology. An abundance of DEGs in the SG of D. citri were identified, revealing their potential roles in response to CLas invasion. This study provides valuable transcriptomic data for further studies of functional genes in SG, which will contribute to our understanding of the molecular mechanisms between CLas and D. citri.
Materials and methods
Citrus psyllid culture
The colony of ACP used in this study was originally collected from a field population on the orange jessamine Murraya exotica L. Mant at the Engineering Research Center of Ministry of Education, South China Agricultural University, Guangzhou, China, in September 2020. A subcolony of ACP was reared in mesh cages (60 cm × 60 cm × 90 cm) with environmental conditions of 26 ± 1°C, 65 ± 5% relative humidity, and a 14:10 (light:dark) photoperiod. Healthy D. citri adults were inoculated onto CLas-infected Citrus reticulata Blanco cv. Shatangju and maintained under the same rearing conditions, and they were infected with HLB pathogen by sucking juice from the phloem during feeding.
To calculate the infection rate of CLas in ACP population, genomic DNA from the head of individual ACP was extracted using a TIANamp Micro DNA Kit (TIANGEN Biotech, Beijing, China), and a 16S ribosomal DNA fragment of CLas was amplified by using specific primer set OI1/OI2c (Shokrollah et al., Reference Shokrollah, Abdullah, Sijam and Abdullah2010). Afterwards, agarose gel electrophoresis was conducted to determine the CLas infection rate. When the proportion of CLas-infected psyllid reached ≥80%, the infected ACPs were collected for tissue dissection.
Tissue collection and total RNA extraction
For SG collection, the uninfected and CLas-infected ACP adults (male/female = 1:1) were anaesthetised on ice, and their SG were dissected in ice-cold sterile 1× phosphate-buffered saline (PBS) solution under a stereomicroscope (ZEISS, Germany). After the dissection, SG tissues were washed with precooled diethypyrocarbonate (DEPC) treated water to remove other remaining tissues, then immediately transferred into RNAlater ® stabilisation solution (Invitrogen, USA) and frozen at −80°C. For transcriptome sequencing, each sample consisted of approximately 1000 pairs of SG, and three biological replicates of the CLas-infected and uninfected ACP subcolonies were processed. The total RNA of citrus psyllid SG was extracted using a miRNeasy® Micro Kit (QIAGEN, Germany) based on the manufacturer's directions. RNA concentration and purity were analysed using a NanoPhotometer Spectrophotometer (Nanodrop Technologies, USA). The integrity of total RNA was detected by RNase-free agarose gel electrophoresis, and assessed on an Agilent 2100 Bioanalyzer system (Agilent Technologies, Palo Alto, CA, USA). The flowchart of transcriptome experiments used in this study is shown in fig. 1.
Figure 1. Flowchart illustrating the procedure used in the current study.
cDNA library construction and Illumina sequencing
After quality inspection, the mRNA of citrus psyllid SG was separated from total RNA with oligo (dT) beads (Epicentre, Madison, WI, USA), and the enriched mRNA was then fragmented into approximately 200–700 nt RNA short fragments by using RNA Fragmentation Reagents (Invitrogen, USA). The cDNA was synthesised using the NEBNext® Ultra™ II RNA Library Prep Kit for Illumina® (New England Biolabs, Ipswich, MA, USA) following the instruction manual. The bead-purified double-stranded cDNA fragments were subjected to end repair, ‘A’ base added, and then ligated into Illumina sequencing adapters. The ligation reaction was purified with the AMPure XP Beads (Beckman Coulter, Brea, CA, USA), and subsequently was size-selected by agarose gel electrophoresis and polymerase chain reaction (PCR) amplification. The quality of library was monitored on the Agilent 2100 Bioanalyzer system, and the concentration was assessed by Qubit® 2.0 Fluorometer (Invitrogen, USA). After the index cluster generation, the resulting cDNA libraries were then deep-sequenced using an Illumina Hiseq™ 4000 platform, and 2 × 150 bp paired-end reads were generated.
Data filtering and read mapping
The raw reads of FASTQ format were produced from the sequencing data using CASAVA software. Raw sequencing reads were monitored for quality using fastp software (v0.18.0), and Q20, Q30, fuzzy base ratio, error rate, and the GC content of raw data were calculated simultaneously. The clean reads were obtained by removing reads containing adapters, poly A/T tail, reads containing more than 10% of unknown nucleotides (poly-N), and low-quality reads containing more than 50% of low-quality bases (Q-value ≤20). To build the index of the reference genome, the paired-end clean reads were mapped to the D. citri reference genome (Genome assembly version 1.1, NCBI RefSeq sequence: GCF_000475195.1) using HISAT2 (v2.2.4) with ‘-rna-strandness RF’ and other parameters set as default. The mapped reads were then assembled by using StringTie (v1.3.1) in a reference-based approach.
Relationship analysis of SG transcriptome samples
To evaluate the reliability and operational stability of the SG transcriptome samples from CLas-free and CLas-infected ACP, correlation analysis was conducted by R Correlation. The Pearson's correlation coefficient R 2 between two SG transcriptomes was calculated based on fragments per kilo bases per million (FPKM) values. The closer R 2 gets to 1, the better the repeatability between two parallel experiments. Furthermore, a principal component analysis (PCA) from the SG transcriptomes was performed with R package gmodels according to the gene expression of correlated variables. Samples from different effective treatments tend to exhibit separate aggregated distributions, and the more similar of the sample compositions in each treatment, the closer the distance reflected in the PCA diagram.
Identification of DEGs
For each transcription region, a FPKM value was calculated to estimate its expression abundance and variations by RSEM. DESeq2 R software (v1.16.1) was used to perform differential expression analysis between CLas-free and CLas-infected subcolonies. The transcripts with the parameter of a false-discovery rate (FDR)-adjusted P-value <0.05 and an absolute value log 2(fold-change) of 0 were defined as the threshold for significantly differential expressed.
Gene ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) enrichment analysis of DEGs
GO functional annotation of DEGs between the SG of CLas-infected and uninfected ACP was conducted by the topGO package. All DEGs were mapped to GO terms in the GO database, and GO terms with a FDR-corrected P-value <0.05 were defined as significantly enriched. KEGG pathway enrichment analysis for DEGs was implemented by using the clusterProfiler R package, and a corrected P-value of 0.05 was set as the threshold for significantly enriched pathways in DEGs.
Reverse-transcription quantitative PCR (RT-qPCR) analysis of DEGs
The transcript abundance of 16 selective DEGs was analysed by RT-qPCR. Primer-BLAST program was used to design specific RT-qPCR primers, and all of the primers used are listed in Supplementary table S1. The first strand cDNA was synthesised with approximately 1 μg total RNA of each sample by using the PrimerScript™ RT reagnet Kit with gDNA Eraser (Perfect Real Time) (Takara, Dalian, China). RT-qPCR reactions were carried out using TB Green® Premix Ex Taq™ II (Tli RNaseH Plus) (Takara, Dalian, China), and the reaction mixture contained 5 μl of Premix Ex TaqII (2×), 0.5 μl of forward primer (10 μM), 0.5 μl of reverse primer (10 μM), 1 μl of cDNA template, and 3 μl of double-distilled water (ddH2O) in a final reaction volume of 10 μl. qPCR was performed with the CFX96 Real-Time PCR Detection System (Bio-Rad, USA), and the standard thermal cycle parameters were as follows: 95°C for 2 min; 40 cycles of 95°C for 15 s, 60°C for 20 s, 72°C for 30 s. The reference gene D. citri Actin (GenBank accession number: XM_008470467) was used as an internal control for normalisation, and the relative gene expression was calculated using the 2−ΔΔCt method. A total of three biological replicates were performed for each treatment, and each biological replicate consisted of three technical replicates.
Statistical analysis
Experimental data from quantitative detection of CLas in the SG assay and RT-qPCR validation of DEGs were subjected to independent sample t-test using SPSS 25.0 software. Statistically significant differences were established as *P < 0.05 or **P < 0.01.
Results
CLas infection in the SG of citrus psyllid
The SG samples were dissected in sterile PBS and then washed with DEPC-treated water to remove excess tissues. Thus, the detached SG showed relatively intact morphology (fig. 2). To verify the presence of CLas in the SG, PCR was performed using CLas-specific primer set OI1/OI2c with genomic DNA isolated from pooled CLas-infected SG samples. The results of gel electrophoresis showed that the CLas-infected SG exhibited clear OI1/OI2c target fragments of 1167 bp in length, while that were not present in the control (fig. 3a). Meanwhile, the RT-qPCR test further validated the CLas infection in the SG of exposed ACP subcolony, which was consistent with the gel imaging results (fig. 3b). Overall, these results indicated that the selected SG samples could be used for the subsequent transcriptome sequencing.
Figure 2. Morphology of the ACP SG that observed under a Discovery V20 stereomicroscope. (a) The position of the SG in the head of an ACP adult; (b) paired SG of an adult ACP; and (c) the structure of the whole excised ACP SG. pg, principal gland; ag, accessory gland.
Figure 3. Qualitative and quantitative detection of CLas in the SG of CLas-free and CLas-infected citrus psyllid subcolonies. (a) Qualitative detection PCR. M, 2000 bp DNA marker; PC, positive control; NC, negative control. The CLas-exposed SG and ddH2O were used as the DNA templates of positive control and negative control, respectively. α-tubulin: the reference gene for the presence of DNA template in each lane. (b) Quantitative detection RT-qPCR. The relative titre of CLas was calculated using the 2−ΔΔCt method. The data represent mean ± SEM (n = 3), **P < 0.01 (independent t-test).
Read assembly and Illumina sequencing
For quality assessment, fastp software was used to trim adaptor sequences and low-quality reads. In total, an average 6,686,101,700 and 6,680,688,100 raw reads from the SG of CLas-infected and uninfected ACP were generated, respectively. After data filtering, an average 6,568,553,270 and 6,565,561,078 high-quality clean reads were obtained with an average GC content of 44.18 and 44.35% (table 1). The mean Q20 (Phred score >20) and Q30 (Phred score >30) values were more than 96.8 and 91.4%, respectively. The paired-end clean reads were subsequently mapped to the ACP reference genome, which resulted in an average mapping rate of 80.26 and 80.12% from the two subcolonies (table 2). The results showed that the accuracy of the sequencing data was adequate for further analysis.
Table 1. Summary of data filtering statistics from two different treatments of citrus psyllid SG
Table 2. Statistics of the number of reads and their ratios mapped to citrus psyllid reference genome in the two SG transcriptomes
Quality control analyses of SG transcriptome samples
A correlation test was conducted to evaluate the reliability of two SG transcriptome samples from CLas-free and CLas-infected ACP as well as operational stability. The average Pearson's correlation coefficient R 2 values between transcriptome samples from CLas-free and CLas-infected ACP reached more than 99.69 and 99.68% (figs 4a and b), respectively, indicating that the SG transcriptome samples exhibited a high stability and repeatability.
Figure 4. Relationship analyses of the SG transcriptome samples from CLas-free and CLas-infected ACP. The repeatability and diversity of the SG transcriptome samples were calculated according to correlation heat map (a), scatter plots, (b) and PCA (c). The oblique line represents the fitting line of the scatter distribution. The Pearson's correlation coefficient was calculated based on the FPKM value of each sample.
In addition, the uncorrelated variable values converted from all transcriptome samples were subjected to PCA, and the total cumulative contribution rates of the first and second principal components to sample variance were 95.8% (fig. 4c), indicating that the quality of SG transcriptomes was adequately reliable. Accordingly, these results suggested that the correlation test and PCA were consistent with respect to the RNA-Seq datasets of CLas-infected and uninfected APC subcolonies.
Identification of DEGs
Significant differences in the SG of CLas-infected D. citri were founded compared with the uninfected D. citri. As a result, the comparative expression analysis resulted in a total of 25,624 genes in the SG transcriptomes of the two ACP subcolonies (table S2). Of these, 861 DEGs were identified in the CLas-infected ACP SG compared with the uninfected ACP, including 202 upregulated and 659 downregulated genes (fig. 5a, table S3). In addition, based on the transcript abundance of each sample, the hierarchical clustering analysis of all the DEGs was performed to evaluate the expression patterns of the identified genes (fig. 5b). These results suggested that CLas infection significantly altered the transcriptional profiles of the DEGs in the ACP SG.
Figure 5. Identification and hierarchical cluster analysis of DEGs between CLas-free and CLas-infected ACP. (a) Volcano plot for the expression pattern of each gene. The abscissa represents the fold change of gene expression in different samples, and the ordinate indicates the significance level of difference in DEGs. The red, blue, and grey points represent upregulated, downregulated, and no significant difference genes, respectively. (b) Hierarchical clustering of DEGs. Rows represent different DEGs, and columns indicate different samples. Red bands and blue bands indicate high and low gene expression level, respectively.
GO and KEGG enrichment analysis of DEGs
GO enrichment analysis was performed to investigate the functions of the DEGs, and a corrected P-value <0.05 was defined as significant enrichment. In total, the GO assignments generated 7516 GO terms, which were functionally grouped into three main categories: biological processes (5544 unigenes), cellular component (775 unigenes), and molecular function (1197 unigenes) (table S4). For each category, the top 20 most significantly enriched GO terms were selected for further analysis (fig. 6a, table S5). In biological processes, the majority of DEGs were associated with cellular processes, single-organism processes, and metabolic processes. A total of 52 upregulated DEGs were related to organic substance metabolic process, 44 upregulated DEGs were mainly involved in cellular metabolic process, and 40 upregulated DEGs were associated with nitrogen compound metabolic processes. In addition, a total of 46 DEGs were mainly involved in the immune system process, including 9 upregulated and 37 downregulated genes, indicating that CLas may suppress the host's immune responses by reducing the expression of immune-related genes. Among molecular function terms, binding, catalytic activity, and transporter activity were the most abundant of 13 categories, while in cellular component, the top three common categories were cell, cell part, and membrane. Among the top 20 assigned GO terms based on GO annotation, the DEGs were mostly enriched in GO functional categories of plasma membrane and cell periphery (fig. 6b, table S4). These results suggested that the cellular activities in D. citri SG were metabolically activated after CLas infection.
Figure 6. GO classification and enrichment analysis of DEGs between CLas-free and CLas-infected ACP. (a) Statistics of the enriched GO terms and genes in GO database according to hierarchical classification. The abscissa represents three main functional categories: BP, biological process; CC, cellular component; and MF, molecular function. For each category, the top 20 GO terms were selected for functional analysis. The ordinate indicates the number of DEGs. (b) GO enrichment analysis of DEGs. The outer ring represents the enriched top 20 GO terms, and different colours indicate different ontologies. Ring 2 represents the gene counts and Q-value in the corresponding GO terms background. The bar length represents the number of DEGs, and the colour changes represent the various Q-values. Ring 3 represents the counts for upregulated and downregulated genes in the corresponding GO terms. The inner ring represents the rich factor values (the number of DEGs divided by the number of background genes in the corresponding GO terms). Each cell in the background grid lines represents 0.1. The scatter diagram of top 20 enriched GO categories. The abscissa indicates the gene ratio. The ordinate indicates different categories.
By assigning to the reference pathway in KEGG, the biological pathways affected by the DEGs were further analysed. According to the KEGG pathway classification, the DEGs were functionally grouped into 42 KEGG categories, including 276 pathways (fig. 7a, table S6). The top 20 KEGG-enriched pathways are shown in fig. 7b, and the majority of DEGs were predicted to be involved in mannose type O-glycan biosynthesis, oxytocin signalling pathway, peroxisome proliferator-activated receptor signalling pathway, and salivary secretion (table S7). Among these pathways, a total of eight, seven, three, and two DEGs were significantly enriched in endocytosis, lysosome, peroxisome, and phagosome, respectively. In addition, the enrichment analysis revealed the assignment of DEGs in multiple carbohydrate metabolism pathways, including pyruvate metabolism, glycolysis/gluconeogenesis, citrate cycle (TCA cycle), starch and sucrose metabolism, pentose and glucuronate interconversions, ascorbate and aldarate metabolism, amino sugar, and nucleotide sugar metabolism. Among amino acid metabolism pathways, arginine and proline metabolism, alanine, aspartate and glutamate metabolism, and cysteine and methionine metabolism were significantly altered. These results suggested that the SG were active in carbohydrate and amino acid metabolism, and the utilisation of these organics may be altered by CLas. Furthermore, several immune system-related pathways, like platelet activation, chemokine signalling pathway, nucleotide oligomerisation domain (NOD)-like receptor signalling pathway, and Toll and immune deficiency (Imd) signalling pathway were found to be significantly affected by the DEGs, indicating that CLas infection may modulate the immune defence in citrus psyllid SG.
Figure 7. KEGG classification and enrichment analysis of DEGs. (a) Statistics of the enriched pathways and genes in KEGG database based on hierarchical classification. The X-axis indicates the number of DEGs. The Y-axis indicates different pathways. (b) KEGG enrichment analysis of DEGs. The outer ring indicates the enriched top 20 pathways, and different colours indicate different clusters. Ring 2 indicates the pathway counts and Q-value in the corresponding DEG background. The bar length represents the number of DEGs, and the colour changes represent the various Q-values. Ring 3 indicates the counts for upregulated and downregulated DEGs in the corresponding pathway. The inner ring indicates the rich factor values (the number of DEGs divided by the number of background genes in the corresponding pathway). Each cell in the background grid lines represents 0.1. The scatter diagram of top 20 enriched KEGG pathways. The X-axis indicates the gene ratio and the Y-axis indicates different categories.
RT-qPCR validation of DEGs
In order to verify the accuracy of the transcriptome sequencing data, the relative expression abundance of 16 DEGs (eight upregulated and eight downregulated genes) were randomly selected for RT-qPCR analysis (fig. 8a). All of the selective target genes are shown in table S8. A comparative analysis suggested that the RT-qPCR results were consistent with our RNA-Seq data (fig. 8b). For instance, the genes of vitellogenin-2-like, acetylcholine receptor subunit alpha-like 1, and E3 ubiquitin-protein ligase MYCBP2-like in CLas-infected subcolony exhibited significant downregulation in both transcriptome data and the RT-qPCR analysis, with a similar fold change of expression level. Moreover, eight genes were upregulated in the SG after CLas infection, whereas only one gene (Adenosine triphosphate [ATP] dependent RNA helicase ddx17) had no significant change in the RT-qPCR analysis. In addition, we performed a linear regression analysis to determine the correlation between RT-qPCR and transcriptome data. The analysis resulted in a correlation coefficient R 2 value of 0.9408 and a corresponding slope of 0.8227 (fig. 8c), indicating a significant correlation between RT-qPCR and the transcriptome data. In summary, these results confirmed the reliability of the RNA-Seq data.
Figure 8. Correlation analysis of the RNA-Seq data and RT-qPCR validation. (a) The relative transcription abundance of 16 DEGs from the SG transcriptomes between CLas-free and CLas-infected ACP. The relative expression levels were calculated using the 2−ΔΔCt method, and the significant differences are indicated by *P < 0.05 or **P < 0.01 using independent t-test. (b) Comparison of the transcriptome data and RT-qPCR results. The data represent mean ± SEM (n = 3). (c) The linear regression analysis between the transcriptome and RT-qPCR data. The abscissa indicates the ratios obtained by RT-qPCR, and the ordinate indicates the ratios obtained by RNA-Seq.
Analysis of DEGs associated with ubiquitination, cytoskeleton, and salivary secretion
Based on functional annotation, the DEGs involved in ubiquitination, cytoskeleton, and salivary secretion showed significant differences in the SG transcriptome between CLas-free and CLas-infected ACP subcolonies (fig. 9, table 3). In the present study, a total of 12 ubiquitination-related DEGs were obtained from the SG transcriptome. Among them, ten DEGs were downregulated in the SG after CLas infection, while two DEGs associated with ubiquitin mediated proteolysis were upregulated, including RING-box protein HRT1-like and cullin-2. For cytoskeleton analysis, a total of 15 DEGs were identified, and 13 DEGs were significantly downregulated, except for troponin T-like and troponin T isoform X11. In addition, we also observed decreased levels of DEGs involved in salivary secretion in the CLas-infected D. citri SG transcriptome compared with that in the control. These results revealed that infection with CLas led to lower levels of DEGs involved in ubiquitination, the cytoskeleton, and salivary secretion in ACP SG, which may be important in the process of CLas invasion.
Figure 9. Hierarchical cluster analysis of DEGs involved in ubiquitination, cytoskeleton, and salivary secretion between CLas-free and CLas-infected ACP SG. Red colour indicates high gene expression level, and green colour indicates low gene expression level.
Table 3. Presence of the DEGs upon CLas infection involved in ubiquitination, cytoskeleton, and salivary secretion between CLas-free and CLas-infected D. citri SG transcriptome
Analysis of DEGs associated with endocytosis, lysosome, and xenobiotics biodegradation and metabolism
According to the GO and KEGG enrichment analyses, many DEGs were classified and significantly enriched in endocytosis, lysosome, and xenobiotics biodegradation and metabolism (fig. 10, table 4). In this instance, a total of 21 DEGs related to endocytosis were found to be significantly downregulated in CLas-infected ACP subcolony. For the lysosome analysis, seven DEGs exhibited significant downregulation in the SG under CLas infection, including four upregulated DEGs and three downregulated DEGs. For xenobiotics biodegradation and metabolism, a total of 12 DEGs showed significant changes; of these DEGs, five genes were significantly upregulated and seven genes were significantly downregulated in CLas-infected D. citri compared with the uninfected control. These results indicated that the endocytosis, lysosome, and xenobiotics biodegradation and metabolism may play crucial roles in the SG response to CLas infection.
Figure 10. Hierarchical cluster analysis of DEGs associated with endocytosis, lysosome, and xenobiotics biodegradation and metabolism between CLas-free and CLas-infected ACP SG. Red colour represents high gene expression level, and green colour represents low gene expression level.
Table 4. Presence of the DEGs upon CLas infection involved in endocytosis, lysosome, and xenobiotics biodegradation and metabolism between CLas-free and CLas-infected D. citri SG transcriptome
Discussion
D. citri is currently the principal transmission vector of CLas bacterium that causes the destructive citrus greening disease (HLB) worldwide. The SG is an important secretory tissue in D. citri with functions in insect feeding, resisting against pathogens, and overcoming plant defences, and it also acts as a barrier to the transmission of CLas (Giron et al., Reference Giron, Huguet, Stone and Body2016). However, the molecular mechanism of SG defence against CLas infection remains unclear. Recent findings have uncovered the transcriptomic and proteomic differences between CLas-infected and uninfected D. citri involving the cytoskeleton, ubiquitination, metabolism, insecticide resistance, and immune system (Kruse et al., Reference Kruse, Fattah-Hosseini, Saha, Johnson, Warwick, Sturgeon, Mueller, MacCoss, Shatters and Heck2017, Reference Kruse, Ramsey, Johnson, Hall, MacCoss and Heck2018; Yu et al., Reference Yu, Li, Zeng, Song, Yu, Su, Chen, Yi and Lu2020; Liu et al., Reference Liu, He, Guan, Zhong, Pang and Han2021). In the current study, we focused on the transcription level changes of the SG conferred by CLas infection in D. citri by using a comparative transcriptomic approach. The transcriptome dataset revealed the functional annotation of DEGs involved in the lysosome, endocytosis, cellular processes, metabolic processes, the immune system, and digestion system. These results suggested that the HLB pathogen may manipulate multiple essential biological processes in host vector to promote its transmission and proliferation.
In host–pathogen interactions, many pathogens interact with highly conserved host components for extracellular recognition and attachment, such as cytoskeleton. The cytoskeleton is a major host structural component, which is essential for controlling cell shape, adhesion, migration, and intracellular transport (Wen et al., Reference Wen, Zhang, Lin, Shi and Jiu2020; Benoit et al., Reference Benoit, Baillet and Poüs2021). The bacterial pathogens modulate the host cytoskeleton to thrive and promote numerous infection events, including invasion, intracellular replication, and dissemination (Colonne et al., Reference Colonne, Winchell and Voth2016). Sarkar et al. revealed that the vector-borne ‘Candidatus Liberibacter solanacearum’ employed the actin cytoskeleton of gut cells for successful movement and transmission in the host vector (Sarkar et al., Reference Sarkar, Kontsedalov, Lebedev and Ghanim2021). Moreover, several Liberibacter species inducing cytoskeletal rearrangements of host cell for invasion and systemic infection in psyllids have also been reported (Ghosh et al., Reference Ghosh, Jassar, Kontsedalov, Lebedev, Wang, Turner, Levy and Ghanim2019; Tang and Tamborindeguy, Reference Tang and Tamborindeguy2020). In the current study, many DEGs were mostly enriched in the pathway of regulation of actin cytoskeleton. In addition, a total of 15 DEGs associated with the cytoskeleton were identified, and most were significantly downregulated in the CLas-infected groups compared with uninfected groups. For example, actin and tubulin, the major protein constituents of cytoskeletal filaments involved in cell motility and signal transduction (Gouin et al., Reference Gouin, Welch and Cossart2005), were significantly downregulated after CLas infection. Rap proteins belong to small GTPases involved in cellular signal transduction, adhesion, and the regulation of cytoskeletal dynamics. Interestingly, the Rap 1 signalling pathway, which involves the metabolism of talin, integrin, actin, and other cytoskeletal proteins (Bromberger et al., Reference Bromberger, Zhu, Klapproth, Qin and Moser2019; Maywald et al., Reference Maywald, Picciotto, Lepa, Bertgen, Yousaf, Ricker, Klingauf, Krahn, Pavenstadt and George2022), was significantly assigned to the top 20 pathways in KEGG (fig. 7). Furthermore, a previous study had revealed the disrupted architecture of actin cytoskeleton within the D. citri midgut following exposure to CLas (Ghanim et al., Reference Ghanim, Fattah-Hosseini, Levy and Cilia2016). This may result from the reduced expression of cytoskeleton-related proteins induced by the CLas infection. Remarkably, vitellogenin transcription was also responsive to CLas-infection in the SG because of its involvement in host recognition (table S3). Vitellogenin is one of the major known pattern recognition receptors with functions in insect fecundity and eliciting host immunity (Han et al., Reference Han, Wang, Song, Zhang and He2022; Sarkar and Ghanim, Reference Sarkar and Ghanim2022). It often acts as a pattern recognition molecule to recognise pathogens, enhance macrophage phagocytosis and autophagy (Brumin et al., Reference Brumin, Lebedev, Kontsedalov and Ghanim2020; Sarkar and Ghanim, Reference Sarkar and Ghanim2022). Previous reports have showed that CLas infection causes the upregulation of vitellogenin genes in the haemolymph, gut, and whole body of D. citri (Ramsey et al., Reference Ramsey, Chavez, Johnson, Hosseinzadeh, Mahoney, Mohr, Robison, Zhong, Hall, MacCoss, Bruce and Cilia2017; Kruse et al., Reference Kruse, Ramsey, Johnson, Hall, MacCoss and Heck2018; Liu et al., Reference Liu, He, Guan, Zhong, Pang and Han2021). However, the majority of vitellogenin genes were unexpectedly downregulated in the SG conferred by CLas, indicating that CLas may induce low expression of vitellogenin genes for successful evasion and proliferation. This result is consistent with a previous study that showed the reduced vitellogenin in CLas-infected D. citri (Vyas et al., Reference Vyas, Fisher, He, Nelson, Yin, Cicero, Willer, Kim, Kramer, May, Crow, Soderlund, Gang and Brown2015), and our hypothesis contributes to explaining the potential upregulation of the CLas titre by silencing of a vitellogenin-like protein in D. citri (Peng et al., Reference Peng, Yuan, Huang, He, Fu, Duan, Yi, Yuan, Yuan, Wang and Zhou2023). In general, our findings suggested that CLas might inhibit the expression of cytoskeleton-associated proteins and vitellogenin genes in the SG to avoid the host immune recognition and facilitate its transmission.
The replication and circulation of the HLB pathogen in SG depend on the suppression of host immune defence and overcoming the SG barrier. The bacterium possibly manipulates the vector host ACP directly and/or indirectly to enhance its spread and transmission (Galdeano et al., Reference Galdeano, Pacheco, Alves, Granato, Rashidi, Turner, Levy and Machado2020), and the immune defence mechanism of SG in response to CLas invasion remains unclear. In this instance, many DEGs are significantly enriched in pathways involved in immune defences, such as the NOD-like receptor signalling pathway, ubiquitin-mediated proteolysis, and Toll and Imd signalling pathway. Among these pathways, numerous immune-related genes, like heat shock proteins, serine/threonine-protein kinases, cathepsins, ubiquitin-ligases, and so forth, were significantly altered in the CLas-infected SG. Cathepsins belong to the papain-like family of cysteine proteases involved in various physiological processes, including apoptosis, intracellular protein degradation, and the immune response (Yang et al., Reference Yang, Zhang, Zhang, Liu, Li, Gong and Hou2020; Chen et al., Reference Chen, Guo, Li, He, Chen, Shen, Zhang, Mao and Zhang2023). For example, Sun et al. revealed that cathepsins O and L may be involved in the innate immune system of Antheraea pernyi Guer (Sun et al., Reference Sun, Tang, Wang, Abbas, Tian, Zhu and Liu2018, Reference Sun, Zhu, Tang, Sun, Chen, Abbas, Wang, Qian, Wei and Liu2017); Guo et al. revealed a novel function for a cathepsin-type protease cathepsin B3 in aphid saliva that elicited effective host plant defences in phloem (Guo et al., Reference Guo, Zhang, Tong, Ge, Wang, Zhao, Zhu-Salzman, Hogenhout, Ge and Sun2020); Yu et al. identified two D. citri cathepsin genes induced by CLas, indicating their potential roles in the innate immune response (Yu et al., Reference Yu, Huang, Li, Xie, Zhou and Lu2019). In this study, a total of 37 cathepsin genes were identified in the SG, and one cathepsin L-like gene was significantly downregulated in CLas-infected ACP. Our result suggested that cathepsins may be important in the SG defence against CLas invasion. For the immune response analysis, we also observed that many DEGs associated with ubiquitination pathway were altered in the obtained SG transcriptome. Ubiquitination is an enzymatic cascade reaction catalysed by several classes of enzymes, including ub-activating enzyme E1, ub-conjugating enzymes E2, and ub-protein ligase E3 (Swatek and Komander, Reference Swatek and Komander2016). Ubiquitination involves membrane trafficking and endocytosis, and it plays critical roles in bacterial recognition, targeting, and clearance (Grumati and Dikic, Reference Grumati and Dikic2018; Mukherjee and Dikic, Reference Mukherjee and Dikic2022). Furthermore, the ubiquitin proteasome system also acts as a key signalling pathway in the host immune response to infection (Ligeon et al., Reference Ligeon, Temime-Smaali and Lafont2011). In this study, we identified 12 ubiquitination-related DEGs in the SG transcriptome. Of those genes, ten DEGs were significantly downregulated after CLas infection, including five E3 ubiquitin-protein ligase genes. Yu et al. also identified 13 DEGs involved in ubiquitination in the D. citri midgut after exposure to CLas (Yu et al., Reference Yu, Li, Zeng, Song, Yu, Su, Chen, Yi and Lu2020). Thus, ubiquitination-correlated genes might play an important role in the process of CLas infection. Moreover, previous studies have revealed that many bacterial effectors enable them to interfere with the ubiquitination pathway to avoid triggering the immune responses (Munro et al., Reference Munro, Flatau and Lemichez2007; Spallek et al., Reference Spallek, Robatzek and Gohre2009). Overall, we considered that CLas may suppress cathepsins and host ubiquitination process to reduce immune-related proteins in the SG of D. citri.
The CLas bacterium is an obligate intracellular pathogen with a relatively small genome (Duan et al., Reference Duan, Zhou, Hall, Li, Doddapaneni, Lin, Liu, Vahling, Gabriel, Williams, Dickerman, Sun and Gottwald2009). Comparative genome sequencing of CLas with the closest cultured relative revealed that CLas is not able to synthesise several essential elements, such as the aromatic amino acids, fatty acids, vitamins, and lipopolysaccharides (Lai et al., Reference Lai, Davis-Richardson, Dias and Triplett2016). Consequently, CLas needs to rely on the insect vector for these essential nutrients. Based on KEGG enrichment analysis, the significant changes conferred by CLas infection were with respect to the pathways involved in nutrient and detoxification metabolism. Trehalose is the major non-reducing circulatory disaccharide in the haemolymph of insects, and provides instant energy substance and defence against abiotic stresses (Iordachescu and Imai, Reference Iordachescu and Imai2008; Shukla et al., Reference Shukla, Thorat, Nath and Gaikwad2015). Trehalase, a hydrolase that catalyses a trehalose molecule into two glucose molecules, is synthesised in the fat body and then transported to the haemolymph via trehalose transporter (Wang et al., Reference Wang, Gou, Guo, Zhou and Liu2021). It has been well documented that insect trehalase can regulate numerous biological processes and metabolic pathways involved in insect homoeostasis, flight, growth, development, nutrition, energy metabolism, and chitin synthesis during moulting (Zhao et al., Reference Zhao, Yang, Shen, Liu, Shi, Wang and Tang2016; Yasugi et al., Reference Yasugi, Yamada and Nishimura2017). In our study, seven trehalase genes were identified in SG transcriptome (table 5), and the induced expression of trehalases post-CLas exposure indicated that CLas may regulate trehalose metabolism pathway by stimulating trehalase activity for nutrition availability. The transcriptomic results also suggested that CLas modulated the detoxification metabolism in D. citri SG. Insect detoxification enzymes, such as cytochrome P450 monooxygenases (CYPs), UDP-glycosyltransferases (UGTs), and glutathione S-transferases (GSTs), are well known to be involved in metabolic detoxification of xenobiotics (Lu et al., Reference Lu, Song and Zeng2021; Gao et al., Reference Gao, Qiao, Wei, Moussian and Wang2022). CYPs and GSTs have been well documented in the SG of several hemipteran insects. For instance, 41 CYPs and 4 GSTs were predicted in the SG of E. fabae (DeLay et al., Reference DeLay, Mamidala, Wijeratne, Wijeratne, Mittapalli, Wang and Lamp2012), and 59 CYPs and 20 GSTs genes were found with various expression levels in Nephotettix cincticeps Uhler (Matsumoto et al., Reference Matsumoto, Suetsugu, Nakamura and Hattori2014). Su et al. identified eight CYPs and five GSTs genes in SG transcriptome of B. tabaci, and a number of genes involved in metabolism and transport were significantly enriched (Su et al., Reference Su, Li, Li, Luan, Ye, Wang and Liu2012). Ji et al. also performed transcriptomic analyses of two N. lugens populations with different virulence; in total, nine and ten GST genes were obtained from the SG, respectively (Ji et al., Reference Ji, Yu, Fu, Chen, Ye, Li and Lou2013). Moreover, 14 CYPs genes (including nine genes upregulated and five genes downregulated) were found to be differentially expressed between the brown planthopper (BPH) populations, and genes related to metabolism, binding, and transport were significantly active. In our DEGs, a total of 21 GST genes were identified in the SG transcriptome. Two CYP genes and two UGT genes were significantly upregulated after CLas infection. Other detoxification enzymes, like esterases (ESTs), were also considered to be responsible for the metabolism of xenobiotics in insects (Xu et al., Reference Xu, Chang, Lu, Tie, Dong, Chen, Ma, Liu and Li2021). Five EST genes were significantly downregulated in the SG of CLas-infected D. citri. These results are in accordance with previous reports that showed the altered detoxification-related genes in CLas-positive D. citri (Tiwari et al., Reference Tiwari, Pelz-Stelinski, Mann and Stelinski2011; Kruse et al., Reference Kruse, Fattah-Hosseini, Saha, Johnson, Warwick, Sturgeon, Mueller, MacCoss, Shatters and Heck2017). In conclusion, our observations suggested that CLas may manipulate trehalase activities and the transcription of detoxification enzymes during infection to facilitate its proliferation in the SG.
Table 5. Candidate genes identified from the SG transcriptomes of ACP
The transmission of CLas involves the propagation and circulation of the bacterium within host vectors, and the pathogen may alter the fitness, survival, and behaviour of its insect host. Previous study reported that the CLas-infected psyllids showed a significant reduction in total non-probing time, salivation time, and phloem ingestion time, indicating that D. citri were at a higher hunger level and they tended to forage more often after CLas infection (Killiny et al., Reference Killiny, Hijaz, Ebert and Rogers2017). During psyllids feeding, salivary secretions play crucial roles in tissue penetration, lubrication, digestion, and nutrient ingestion. Mucin-like proteins (Mul) are a group of glycoproteins that cover the epithelial cell layer of many cuticle-free tissues in insects, such as the SG, midgut and Malpighian tubules (Lou et al., Reference Lou, Shen, Li, Huang, Lu and Zhang2019; Zhao et al., Reference Zhao, Zhang, Yang, Niu, Zhang and Yang2020; Ahmed et al., Reference Ahmed, Seo and Kim2021). Mul are abundantly secreted as the principal components of saliva during insect feeding, and specifically expressed in the SG. In the transcriptome herein obtained, we identified 59 mucin-like proteins from dissected ACP SG, and the most abundant transcript encodes a mucin-5AC-like protein with an average FPKM value of 8597.05 (table S6). In addition, one digestive system-related pathway, salivary secretion, was found to be significantly enriched in the top 20 pathways in KEGG (fig. 7, table S6). Among these pathways, ten DEGs were significantly downregulated after CLas infection. These results indicated that the CLas pathogen may manipulate the secretory activity of saliva in the SG to alter the feeding behaviour of its psyllid vector. Interestingly, recent studies have revealed the potential roles of Mul in salivary sheath formation and modulation of plant defence responses (Shangguan et al., Reference Shangguan, Zhang, Liu, Zhao, Wang, Wang, Guo, Rao, Jing, Guan, Ma, Wu, Hu, Chen, Du, Zhu, Yu and He2018; Huo et al., Reference Huo, Zhao, Meng, Yang, Zhang, Liu, Fang and Zhang2022; Liu et al., Reference Liu, Yi, Jia, Miao and Hou2022). Thus, these mucin-like proteins may also play dual roles in the formation of salivary sheath and response to plant immune defence in the SG of D. citri.
In summary, by using comparative transcriptomic analysis, we identified 861 DEGs between the SG of CLas-free and CLas-infected ACP colonies. The results of functional analysis revealed that CLas made noticeable effects on the SG during infection involved in host recognition, the immune response, nutrient and detoxification metabolism, and salivary secretion. This study provides scientific insights into the further research on the response mechanisms of the citrus psyllid to the invasion of CLas, which is much valuable for the formulation of consistent control strategy of citrus psyllid.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0007485324000038.
Data availability
All information of transcriptome sequencing has been deposited in the NCBI database with a bioproject accession no. PRJNA1013779. The raw sequence reads from six salivary gland transcriptome samples are available in the NCBI Sequence Read Archive (SRA) database under the accession numbers SRR25937708–SRR25937713.
Acknowledgements
The authors thank the editor and reviewers for the modifications and suggestions for the manuscript.
Author contributions
Conceptualisation: S. Z. and B. Q.; methodology: S. Z. and B. Q.; software: S. Z. and B. Q.; validation, S. Z., X. R., Y. H., W. S., B. E. D. and B. Q.; formal analysis: S. Z. and B. Q.; investigation: S. Z., X. R., Y. H., W. S., B. E. D. and B. Q.; resources: S. Z., W. S., and B. Q.; data curation: B. Q.; writing – original draft preparation: S. Z.; writing – review and editing: S. Z., B. E. D. and B. Q.; visualisation: B. Q.; supervision: B. Q.; project administration: W. S., B. E. D. and B. Q.; funding acquisition: W. S. and B. Q. All authors have read and agreed to the published version of the manuscript.
Financial support
This study was supported by the project of Guangdong Laboratory for Lingnan Modern Agriculture (NT2021003), the National Joint Fund for Regional Innovation and Development (U22A20481), and the Program of Agricultural Science and Technology Innovation for the 14th Five-year Plan of Guangdong Province (2023SDZG06) to B.-L. Q.
Competing interests
None.
