B. cockerelli colonies free of CLso [CLso(-)] or carrying CLso (haplotype B) [CLso(+)] were reared on potato plants and maintained at the Texas A&M AgriLife Research & Extension Center, Weslaco, Texas. The absence [CLso(-)] or presence [CLso(+)] of the CLso bacterium in the psyllid salivary glands and ovaries was determined by employing molecular diagnostics; polymerase chain reaction (PCR) using 16S rDNA primers specific to the CLso (Avila et al., 2019). Further, salivary glands and ovaries of teneral adult psyllids were carefully dissected under the microscope. Briefly, total RNA from three biological replicates representing the salivary glands and ovaries of CLso(-) and CLso(+) psyllids were extracted using a Direct-zol RNA Miniprep kit (Zymo Research, Irvine, CA). The quantity and quality of the RNA samples were estimated by NanoDrop-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA) and Bioanalyzer (Agilent Technologies, Santa Clara, CA), respectively. The four RNA samples representing CLso(-) and CLso(+) conditions of the salivary glands and ovaries were further processed to prepare libraries and sequenced in paired-end mode to obtain 150 bp long reads per Illumina’s recommendations (San Diego, CA) at Texas A&M AgriLife Genomics and Bioinformatics Services (College Station, TX). The raw files of RNA-Seq data of the four samples, each with three biological replicates, have been deposited at NCBI’s Gene Expression Omnibus (GEO) under the PRJNA1076531 BioProject accession number.

The quality of the paired-end reads was analyzed using the NGS QC Toolkit using the default parameters (Patel and Jain, 2012). The high-quality reads were mapped to the reference genome of B. cockerelli (Kwak et al., 2023) using the STAR aligner at default settings (Dobin et al., 2013). The number of reads in each gene was estimated using the reference annotations (Kwak et al., 2023) and the feature Counts tool at the default settings (Liao et al., 2014). Differentially expressed genes (DEGs) under CLso(+)/CLso(-) [CLso(-) as reference] in each organ and/or between the two organs (ovaries/salivary glands; salivary glands as reference) under CLso(+) and CLso(-) conditions were estimated with the default normalization algorithms implements in DESeq2 tool (Love et al., 2014). In addition, each gene’s expression level was estimated using the normalized read count data by considering read depth and gene size (fragments per kilobase per million mapped fragments; FPKM).

Principal component analysis (PCA) was performed using the expression values (FPKM) with the Factoextra package in the R program to examine the relationship among the four samples. The Pearson correlation coefficient (R) among the four samples was estimated and shown via heatmap using Python scripts.

To examine the chromosome-wide expression patterns, we first classified all genes into expressed (≥0.5 FPKM) or not-expressed (<0.5 FPKM) categories. Subsequently, the frequency of expressed genes per 100 KB was estimated and plotted via Circos-plot (Krzywinski et al., 2009). Likewise, the chromosome-wide expression level was estimated and visualized using box plots.

To examine organ-specific expression under CLso(+) and CLso(-) conditions, we estimated the tissue specificity index (TSI) (Kryuchkova-Mostacci and Robinson-Rechavi, 2017) among the four samples. Only the expressed (≥0.5 FPKM) genes in any of the four samples were selected and used as input to determine TSI in the R program. Those genes with ≥0.9 TSI in any of the four samples were considered specifically expressed, as described previously (Jain et al., 2022a).

We performed a weighted gene co-expression network (WGCNA; v1.61) (Langfelder and Horvath, 2008) analysis among the four distinct samples to identify sets of co-expressed genes exhibiting significant correlation. Genes with ≥0.1 FPKM in any of the four samples were selected and transformed into a log₂ scale. Those genes with ≥0.1 variances (based on their log₂ transformed FPKM values) were further used to determine significantly correlated modules. A cut-off of 9 β-value (soft-threshold power) and scale-free topology index (R²) was used to estimate the adjacency matrix and sets of co-expressed modules with high correlation (R) and significance (P-value) as described earlier (Garg et al., 2017). The relationship of the co-expressed modules was determined using the default parameters implemented in the WGCNA program and displayed as a dendrogram.

To examine the functional relevance of the DEGs and co-expressed genes, we performed gene ontology (GO) enrichment analysis using the BINGO plug-in implemented in Cytoscape (v3.9.1) (Shannon et al., 2003). The enriched GO terms with <0.05 P-value were considered significant. The number of genes and significance level (P-value) in each enriched GO term were presented via bubble plots generated in the R program. Moreover, comparative GO analyses among the different sets of genes were examined using the EnrichmentMap tool implemented in Cytoscape (v3.9.1) (Shannon et al., 2003).

To assign transcription factors (TFs), we performed a local blastP between amino acid sequences of the psyllid proteins with amino acid sequences of Drosophila melanogaster TFs available in AnimalTFDB (v3.0) (Shen et al., 2023). The most significant match with a cutoff of <0.005 P-value was assigned for each gene of D. melanogaster to determine TFs in B. cockerelli.

Psyllid tissue-specific transcriptome rendered in a total of 8.44-11.78 million paired-end reads for each sample. Of these, 93-95% were high-quality reads, as determined by the NGS QC Toolkit (Patel and Jain, 2012). The reads were next mapped to the B. cockerelli reference genome (Kwak et al., 2023). The current B. cockerelli genome draft is the first release and consists of ~4409 scaffolds of at least >1 kb size, with 13 biggest scaffolds representing 13 distinct chromosomes, covering ~92% of the psyllid genome (Kwak et al., 2023). Due to noncontiguous sequences in the smaller scaffolds, we only mapped the RNA-Seq reads to the 13 biggest scaffolds. In total, ~4.71-7.48 million read pairs were mapped uniquely (~53-66%) to the 13 biggest scaffolds of the reference genome (Kwak et al., 2023) ( Supplementary Table 1 ). The mapping percentage could be marginally higher if the reads were mapped to all other scaffolds. Approximately 19,699 genes were annotated on the 13 chromosomes (Kwak et al., 2023). All the downstream transcriptome analyses were performed with these reference gene annotations.

First, we examined the broader relationships of the four distinct samples representing the salivary glands and ovaries of CLso(-) and CLso(+) psyllids using the FPKM expression metric. The transcriptome profiles of the salivary glands and ovaries were quite distinct, with Pearson’s correlation of 0.34 and 0.38 under CLso(-) and CLso(+) conditions, respectively. In response to CLso [CLso(-) vs CLso(+)] within the specific organ(s), a correlation of 0.97 and 0.74 were detected in the salivary glands and ovaries, respectively ( Figures 1A, B ). These correlations were as expected since transcriptomes of distinct tissue types/organs would be more disparate than the overall gene expression profiles influenced by CLso’s presence. Recently, a tissue-specific transcriptome analysis of another closely related psyllid, Diaphorina citri (Asian citrus psyllid), vectoring Candidatus Liberibacter asiaticus (CLas), the presumptive causal agent of Huanglongbing disease was reported (Mann et al., 2022). The results showed unique transcriptome profiles in the distinct organs responding to CLas (Mann et al., 2022). Our results in B. cockerelli in response to CLso were broadly similar to those observed for the tissue-specific transcriptomes of the Asian citrus psyllid vectoring CLas in the virtue of organ-specific transcriptome signatures under the influence of fastidious bacteria (Mann et al., 2022).

The 13 chromosomes of B. cockerelli were assigned into 12 autosomes and one sex chromosome (Kwak et al., 2023). The size of the autosomes ranged from 23.61 65.45 Mb, while the sex chromosome spanned 22 Mb (Kwak et al., 2023). The genomic attributes of the ~19,699 B. cockerelli annotated genes and their expression dynamics remain unexplored. Here, we analyzed the length distribution of the ~19,699 B. cockerelli genes. A substantial fraction (26.3%) of the genes were ∼1-3 kb long ( Supplementary Figure 1A ). Next, we examined the chromosome-wide expression profiles based on the FPKM metric. Genes with ≥0.5 FPKM were considered expressed (Jain et al., 2022b). Among the 12 autosomes, ~980-2888 genes were located per chromosome. Of them, ~535-1,598 (~37.96-61.2%) genes were found to be expressed in any of the four distinct samples analyzed. A similar analysis was performed with the single sex chromosome, too. The percentage of expressed genes on the sex chromosome was greater (58.59-68.98%; 531 out of 693 annotated genes) than the autosomes ( Supplementary Figure 1B ; Supplementary Table 2 ). Further, we estimated the frequency of expressed genes in a window size of 100 kb, as shown via Circos-plot. In general, the expressed genes were found to be evenly distributed across the length of most of the chromosomes ( Supplementary Figure 2A ). Interestingly, ChrX harbored the highest frequency (58.59-68.98%) of expressed genes (≥0.5 FPKM) and the least frequency in Chr12 (37.96-44.59%) and Chr1 (39.44-43.39%) in any of the four samples analyzed ( Supplementary Figures 1B , 2A ). Moreover, the expression level of genes located in ChrX (1.18-2.69; log₂ transformed FPKM median expression) also exhibited higher expression than autosomes (0.16-0.63; log₂ transformed FPKM median expression) in all four samples analyzed ( Supplementary Figure 2B ). The high-level expression of genes on the sex chromosome may be critical for regulating biological processes related to sex determination and governing reproductive development processes in B. cockerelli.

Transcription factors (TFs) play an important role in regulating the expression of downstream genes via binding in their cis-regulatory elements (CREs) (Priest et al., 2009; Wittkopp and Kalay, 2011). However, the genes encoding TFs in B. cockerelli are not well annotated or described in the released genome (Kwak et al., 2023). Therefore, we annotated and characterized the TFs in B. cockerelli based on the most significant match with drosophila’s TFs available in AnimalTFDB (v3.0) (Shen et al., 2023). A total of 801 genes encoding TFs were detected in B. cockerelli. Of these, 528-575 TFs were found to be expressed with ≥0.5 FPKM ( Figure 2A ). The percentage of expressed TFs (65.92-71.79%) was substantially higher than non-TF genes (46.71-52.39%) across the four different samples ( Figure 2B ; Supplementary Table 3 ). Among the expressed TFs, zf-H2C2 (30.37%), zf-C2H2 (22.43%), BTB (9.97%), homeodomain (5.3%), and HLH (4.21%) were the five topmost represented classes of TFs ( Figure 2C ; Supplementary Table 3 ).

Further, to examine genes expressed in an organ-specific manner in the salivary glands and ovaries in the CLso(-) or CLso(+) conditions, we estimated the tissue specificity index (TSI) of the four distinct samples using their FPKM values (Kryuchkova-Mostacci and Robinson-Rechavi, 2017). The relevance and estimation methods are quite different from differential expression analysis. In TSI, a set of uniquely expressed genes is determined in a tissue/organ/condition-specific manner with reference to all the other remaining samples. In contrast, differential expression analysis is performed between two samples, representing only one reference. TSI analysis showed that 123 and 20 genes exhibited tissue-specific expression under the CLso(-) condition in the salivary glands and ovaries, respectively. Likewise, 119 and 41 genes showed a tissue-specific expression under the CLso(+) condition in the salivary glands and ovaries, respectively ( Figure 3A ). Gene ontology (GO) analysis revealed significant enrichment of processes involved in localization/transport, signaling, and endocytosis preferentially/specifically under CLso(+) condition in the salivary glands. Likewise, specific enrichment of GO terms involved in the primary metabolic process, macromolecule metabolic process, and protein metabolic processes were observed under the CLso(+) condition in the ovaries. In contrast, significantly enriched processes under the CLso(-) condition in both organs were represented with only a few genes, rendering only a trivial correlation in the absence of CLso ( Figure 3B ). In addition, other attributes of GO enrichment analysis, such as molecular function and cellular component, were interrogated. Molecular function GO terms involved in DNA binding, signal transducer activity, ATP binding, and membrane transporter activity were enriched in specifically expressed genes in salivary glands under CLso(+) condition, and localized in the membrane, plasma membrane, membrane-bound vesicles ( Supplementary Figures 3A, B ). Likewise, GO terms related to catalytic activity, peptidase/endopeptidase activity, and localization in intracellular spaces were enriched under the CLso(+) condition in ovaries. Conversely, the enrichment of the molecular function and cellular component GO terms were underrepresented in CLso(-) condition in both the salivary glands and ovaries ( Supplementary Figures 3A, B ). We suggest that the influence of CLso might be implicated in determining unique signatures of the two distinct organs via regulating stress-responsive biological processes and functions.

A few previous studies analyzed the differential expression of genes in response to CLso in the whole psyllid (B. cockerelli) (Nachappa et al., 2012; Huot et al., 2018). The psyllid salivary glands and ovaries are vital organs for pathogenesis and transmission to the host plants. Therefore, a tissue-specific differential expression analysis in these organs in response to CLso may provide blueprints for insect-pathogen interactions. We identified DEGs between CLso(-) and CLso(+) conditions [CLso(-) as reference] within the distinct organ(s). A total of 110 and 93 genes exhibited up- and downregulation under the CLso(+)/CLso(-) condition in the ovaries ( Figure 4A ; Supplementary Table 4 ). Likewise, 514 and 479 genes showed up- and downregulation in response to CLso in the salivary glands, respectively ( Figure 4B ; Supplementary Table 4 ).

Next, we identified common and unique sets of genes that were up- and downregulated in response to CLso in the two psyllid organs. A total of 54 and 461 genes were upregulated explicitly in the ovaries and salivary glands, respectively. Likewise, 66 and 449 genes exhibited specific downregulation in the ovaries and salivary glands, respectively ( Figure 4C ). To examine the functional relevance of these uniquely up- and downregulations in response to CLso, we analyzed the enrichment of GO terms. Our results showed that specifically upregulated genes in CLso-infected ovaries were related to GO terms involved in primary/macromolecule/protein metabolic processes, proteolysis, and cuticle development. For instance, the GO term ‘proteolysis’ contained the genes encoding trypsin 1, carboxypeptidase N, and cathepsin L1. These are critical proteolytic enzymes involved in diverse biological processes, such as reproduction, embryo development, and immunity (Law et al., 1977; Krautz et al., 2014; Jagdale et al., 2017; Yu et al., 2019; Ferrara et al., 2020). The upregulation of cathepsin L1 (BcCathL1) in psyllid ovaries might suggest that the host is mounting a tissue-specific immune response against CLso. Also, it could be associated with ovarian diapause, reducing the number of eggs oviposited, as was demonstrated in Coccinella septempunctata (Linnaeus) (Chen et al., 2022). The transcriptional modulation of proteolytic enzymes (upregulation) determined in our study might also be related to a previous study by Albuquerque Tomilhero Frias et al., 2020, wherein CLso-infected females (haplotype B, same as used in this study) showed reduced oocyte development and eggs oviposited compared to uninfected psyllids. We hypothesize that reduced progeny could be associated with the survival and reproduction of the CLso-infected psyllids in the host plants. A body of evidence revealed the preference–performance hypothesis (PPH), which determines the choice of host plant selection by the phytophagous insect to ensure female oviposition and successive transmission of the offspring (Mayhew, 1997; Gripenberg et al., 2010). The phenomenon of PPH is supposedly implicated in the CLso-infected phytophagous psyllids for their successive transmission. Likewise, GO terms related to cellular/metabolic processes, response to stimuli/stress, reproductive cellular process, oxidation-reduction, cell cycle phase, membrane organization, vesicle-mediated transport, RNA splicing, phagocytosis, and programmed cell death were enriched in specifically upregulated genes in the salivary glands ( Figure 4D ). For instance, transcripts of Carboxypeptidases (E and M) from the GO term ‘vesicle-mediated transport’ were up-regulated explicitly in CLso-infected salivary glands. Carboxypeptidases belong to a large family of enzymes governing the cleavage of C-terminal residues (Gomis-Rüth, 2008) and they play an important role in several biological processes, such as control of peptide activity at the cell surface and post-translational processing of extracellular proteins and peptides (Reznik and Fricker, 2001). The function of Carboxypeptidases is still unclear in tomato/potato-psyllid-CLso interactions. However, these proteins could be plausibly involved in the defense response of the host plants. We hypothesize that the colonization of salivary glands by CLso may modulate the expression of host vesicle-mediated transport-related genes to evade plant-immune defenses towards the pathogen, degrading critical proteins and peptides necessary to reduce insect feeding behavior and the probability of CLso transmission. An electrical penetration graph study on CLso-infected psyllids showed a significant increase in salivation, phloem ingestion, and several probes in comparison to the uninfected psyllids (Valenzuela et al., 2020), supporting our hypothesis; further investigations will unveil the function of vesicle-mediated transport-related genes in the potato psyllid-CLso interaction. Conversely, GO terms involved in the metabolic process, response to stimuli, small molecule metabolic process, and phagocytosis were enriched in the specifically downregulated genes in the ovaries. Likewise, GO terms involved in RNA splicing/processing, locomotion, and homeostatic processes were enriched in specifically downregulated genes in the salivary glands ( Figure 4D ).

Next, we performed a comparative GO among the four sets of specifically up- and downregulated genes in response to CLso. The comparative GO enrichment analysis enables us to identify relatively more significant GO terms among the multiple sets of genes with greater stringent criteria than the conventional GO enrichment analysis from a single set of genes. Amino acid/histamine/guanosine metabolic processes, cuticle development, and protein repair terms were exclusively enriched in specifically upregulated genes in the ovaries. Likewise, in response to stress, DNA metabolic and DNA replication-related GO terms were enriched in specifically upregulated genes in the salivary glands ( Figure 4E ). Conversely, no GO terms were enriched in specifically downregulated genes in the ovaries. However, metabolic processes related to nucleic acid/nitrogen compound/small molecule were enriched in specifically downregulated genes in the salivary glands ( Figure 4E ). Further analysis of the molecular function and cellular component GO terms revealed the enrichment of metal ion binding, hydrolase activity, peptidase activity, and localization in intracellular spaces in the upregulated genes in response to the presence of CLso bacterium in ovaries. Likewise, GO terms involved in metal ion binding, hydrolase activity, oxidoreductase activity, ATPase binding, GTP binding, and helicase activity were preferentially enriched and localized in intracellular spaces, cytoplasm, cytoskeleton, and ribonucleoprotein complexes in the presence of CLso in the salivary glands ( Supplementary Figures 4A, B ). Conversely, enriched GO terms in downregulated genes in the presence of CLso in the ovaries included metal ion binding and oxidoreductase activity. Likewise, cofactor binding, transferase activity, lyase activity, nuclease activity, UDP-glucosyltransferase activity, and localization into the cytoskeleton, nucleus, microtubule-associated complex, and lipid particle were identified in the set of downregulated genes in the presence of CLso in the salivary glands ( Supplementary Figures 4A, B ). Overall, the results suggest that distinct transcriptional changes in multiple biological processes and functions in salivary glands and ovaries are crucial for the CLso-potato/tomato psyllid interactions.

After initial acquisition by B. cockerelli, CLso multiplies and circulates in the body of the psyllid. About two weeks after the acquisition, the CLso reaches the salivary glands and acquires potency to infect the host plants (Sengoda et al., 2013; Rashed et al., 2014; Sengoda et al., 2014). In contrast, replication and vertical transmission of CLso in successive generations are mostly confined to the ovaries and eggs of the psyllid (Hansen et al., 2008; Casteel et al., 2012; Kean et al., 2019). Therefore, we next performed a comparative analysis between the two organs (ovaries/salivary glands) in the presence or absence of the bacterium. In the control psyllids [CLso(-)], a total of 2,481 and 3,113 genes exhibited up- and downregulations, respectively ( Figure 5A ; Supplementary Table 5 ). Likewise, 2,526 and 3,043 genes showed up- and downregulations, respectively, under the CLso(+) condition ( Figure 5B ; Supplementary Table 5 ). Further, these four sets of up- and downregulated genes were used to identify uniquely up- and downregulated genes. In total, 493 and 689 genes exhibited specific up- and downregulation between the ovaries and salivary glands under the CLso(-) condition. Likewise, 535 and 622 genes showed exclusively up- and downregulation, respectively, under the CLso(+) condition ( Figure 5C ).

We performed an enrichment analysis of GO terms to gain insights into the functional relevance of these specifically up- and downregulated genes. Under the CLso(-) condition and upregulation in the ovaries, GO terms involved in macromolecule/nitrogen compound/nucleic acid metabolic processes, cellular catabolic process, and proteolysis were enriched. Conversely, GO terms involved in localization, signaling, neurological system process, and polysaccharide biosynthetic process were enriched under CLso(-) condition and downregulated in the ovaries ( Figure 5D ). Likewise, in response to CLso(+) condition and upregulation in the ovaries, GO terms involved in nucleic acid/cellular macromolecule/DNA metabolic processes were enriched. However, GO terms related to nucleic acid/nitrogen compound metabolic processes, membrane invagination, and endocytosis were enriched in the set of downregulated genes under CLso(+) condition ( Figure 5D ).

Further, we performed comparative GO analysis among the specifically up- and downregulated genes between the ovaries and salivary glands under the CLso(-) and CLso(+) conditions. Metabolic processes related to nucleic acid/RNA/nitrogen compound/macromolecule were enriched in the upregulated genes in the ovaries under both the CLso(-) and CLso(+) conditions. However, metabolic processes involved in amino acid/carboxylic acid/ketone bodies/organic acid were exclusively enriched in response to CLso(+) condition in the ovaries. Conversely, processes involved in membrane invagination and endocytosis were specifically enriched in the set of downregulated genes under the CLso(+) condition ( Figure 5E ). Previous studies showed that the reproductive fitness of Bactericera cockerelli is adversely affected by CLso’s presence (Nachappa et al., 2014; Albuquerque Tomilhero Frias et al., 2020). Such antagonistic relationships have been studied earlier (Schwenke et al., 2016). The CLso replication and psyllid reproduction processes are physiologically and energetically demanding and governed by insect hormones, 20-hydroxyecdysone, and insulin/insulin-like growth factor via endocrine-regulated energy metabolism pathway (Lorenz and Gerd Gäde, 2009; Schwenke et al., 2016). Our results suggest that the perturbation of various metabolic processes in response to CLso in the ovaries could be corroborated by a trade-off between lower reproductive fitness and enhanced immunity of the psyllids. Further analysis of the molecular function and cellular component GO terms revealed enrichment of hydrolase activity, nucleotide binding, and transferase/methyltransferase activity and localization into membrane-bound organelle, nucleus, nuclear lumen, and nucleoplasm in the set of upregulated genes in ovaries in the presence of CLso bacterium. Likewise, nucleotide binding, transferase activity, ATP binding, and hydrolase/helicase/ligase activity and localization into the membrane and membrane-bound organelle were identified among the upregulated genes in the ovaries in the absence of CLso bacterium ( Supplementary Figures 5A, B ). Conversely, the enriched GO terms in downregulated genes were minimal in the presence of CLso while a few GO terms, such as cation binding, transcription regulator activity, and mono-oxygenase activity and, localization into the membrane-bound organelle, and membrane were enriched in the absence of the CLso bacterium ( Supplementary Figures 5A, B ).

Next, we analyzed the differential expression of genes encoding TFs between the ovaries and salivary glands under CLso(-) and CLso(+) conditions. In total, 31 and 25 TFs under the CLso(-) condition showed up- and downregulation in the ovaries, respectively. Likewise, a total of 19 and 23 TFs were specifically up- and downregulated, respectively, under CLso(+) condition in the ovaries ( Supplementary Figure 6A ). Among them, TFs encoding zinc-finger (zf-H2C2 and zf-C2H2) types, BTB, homeodomain, HMG box, HLH, and bZIP were predominant ( Supplementary Figure 6B ). Zinc finger and bZIP proteins are involved in diverse processes, including growth, development, and immune response in different organisms (Alves et al., 2013; Fu and Blackshear, 2017; Rakhra and Rakhra, 2021; Zheng et al., 2021). Likewise, HLH and HMG TFs were implicated in defense responses (Malarkey and Churchill, 2012; Lee et al., 2014; Chen et al., 2017; Cook et al., 2020). The TFs perturbed in the psyllid organs and associated with CLso could also mediate signal transduction and contribute to psyllid growth, development, and stress responses. TFs play a crucial role in regulating the expression of several associated downstream genes and their dynamic expression between the ovaries and salivary glands of the psyllids in the absence or presence of CLso bacterium may provide blueprints and mechanistic insights into the pathogenesis and transmission of the diseases to the host plants.

Gene co-expression and regulatory networks allow a better understanding of biological functions and processes at a modular level and can uncover new insights (Iyer-Pascuzzi et al., 2011; Macneil and Walhout, 2011). We employed a weighted gene co-expression network (WGCNA) to identify modules of co-expressed genes. Eight modules (M1-M8) were identified, with the highest number of genes in the M5 module (5517). In addition, only the M5 module exhibited a significant (P-value ranging from 0.008 to 8e-04) correlation under CLso(-) and CLso(+) conditions in the salivary glands and ovaries ( Figure 6A ; Supplementary Figure 7 ). The pattern of modules was contrasting between the salivary glands and ovaries but largely similar between the CLso(-) and CLso(+) conditions within the individual organ(s). However, a subtle difference in the significance level (P-value) in the M5 module was observed between the CLso(+) and CLso(-) conditions. In the salivary glands, the significance level under CLso(+) (P-value = 0.005) was slightly higher than the CLso(-) condition (P-value = 0.008). Likewise, a lower significance level was observed under CLso(+) (P-value = 8e-04) than the CLso(-) condition (P-value = 0.004) in the ovaries. The correlation and significance level of the co-expressed genes between the CLso(-) and CLso(+) conditions of the M5 module were similar in both the salivary glands and ovaries. Conversely, a contrasting correlation and significance level of the co-expressed genes between the salivary glands (R=-0.99) and the ovaries (R=1.0) was observed ( Figures 6A, B ).

Next, we performed a comparative GO enrichment analysis to gain deeper insights into the role of preferentially co-expressed genes in the salivary glands and ovaries. Interestingly, transport, neuronal transmission, cell adhesion, light stimulus, and respiration processes were enriched in the set of genes preferentially expressed in the salivary glands ( Figure 6C ). Many of these processes are relevant to transmitting the bacterium to the vasculature of host plants. For instance, neuronal transmission in Drosophila controls immunity and governs behaviors to coexist with the pathogens (Montanari and Royet, 2021). Cell adhesion molecules facilitate the invasion and colonization of the pathogens in the host organism (Bisht and Meena, 2019). In addition, transport and respiration are physiologically and energetically demanding processes during the infection of the host organisms (Porcheron et al., 2013; Killiny et al., 2018). In agreement, previous studies in B. cockerelli have shown that salivary glands are critical for CLso transmission to host plants (Sengoda et al., 2013; Rashed et al., 2014; Sengoda et al., 2014). Likewise, GO terms involved in DNA replication, transcriptional regulation, translation, cell division, and response to stress were enriched in the set of preferentially expressed genes in the ovaries ( Figure 6C ). The enrichment of processes related to bacterial replication, specifically in the ovaries, is concurrent with the previously reported vertical transmission of the bacterium in the eggs and ovaries of the psyllid (Hansen et al., 2008; Casteel et al., 2012; Kean et al., 2019). Further, we examined the enrichment of molecular function and cellular component GO terms in the distinct sets of genes preferentially expressed in the salivary glands and ovaries. In the salivary glands, GO terms involved in cation binding, transporter, oxidoreductase, substrate-specific channel, ATPase, cofactor binding, NADH dehydrogenase, voltage-gated channel, and calmodulin activity were enriched and localized in the cytoplasm, plasma membrane, organelle membrane, mitochondria, and respiratory chain ( Supplementary Figures 8A, B ). In contrast, the enrichment of GO terms in the set of ovaries’ preferential genes included hydrolase, nucleotide binding, ATP binding, transferase activity, pyrophosphatase, ligase, helicase, RNA polymerase, and chromatin binding function, and localization into the cytoplasm, cytoskeleton, nucleus, chromosome, ribosome, and proteasome complex were identified ( Supplementary Figures 8A, B ). Overall, co-expression network analysis identified salivary gland and ovary-specific genes and uncovered biological processes and functions relevant to CLso acquisition and transmission into the host psyllid organs.

Whether the interaction between psyllid and the bacterium is of a beneficial nature or antagonistic can be debatable. In one scenario, the positive interactions between the psyllid and the bacterium could be crucial for their coexistence and confer an adaptive advantage over their natural predators. However, we and others have reported detrimental effects on the psyllid due to the presence of CLso, such as a substantial reduction in reproductive fitness (Qiu and Scholthof, 2000; Montllor et al., 2002; Nachappa et al., 2014; Heyworth and Ferrari, 2016; Doremus and Oliver, 2017; Albuquerque Tomilhero Frias et al., 2020). The trade-off is strikingly not lethal, and the association has hence prevailed.