Differential Transcriptional Responses in Two Old World Bemisia tabaci Cryptic Species Post Acquisition of Old and New World Begomoviruses

Begomoviruses are transmitted by several cryptic species of the sweetpotato whitefly, Bemisia tabaci (Gennadius), in a persistent and circulative manner. Upon virus acquisition and circulative translocation within the whitefly, a multitude of molecular interactions occur. This study investigated the differentially expressed transcript profiles associated with the acquisition of the Old World monopartite begomovirus, tomato yellow leaf curl virus (TYLCV), and two New World bipartite begomoviruses, sida golden mosaic virus (SiGMV) and cucurbit leaf crumple virus (CuLCrV), in two invasive B. tabaci cryptic species, Middle East-Asia Minor 1 (MEAM1) and Mediterranean (MED). A total of 881 and 559 genes were differentially expressed in viruliferous MEAM1 and MED whiteflies, respectively, compared with their non-viruliferous counterparts, of which 146 genes were common between the two cryptic species. For both cryptic species, the number of differentially expressed genes (DEGs) associated with TYLCV and SiGMV acquisition were higher compared with DEGs associated with CuLCrV acquisition. Pathway analysis indicated that the acquisition of begomoviruses induced differential changes in pathways associated with metabolism and organismal systems. Contrasting expression patterns of major genes associated with virus infection and immune systems were observed. These genes were generally overexpressed and underexpressed in B. tabaci MEAM1 and MED adults, respectively. Further, no specific expression pattern was observed among genes associated with fitness (egg production, spermatogenesis, and aging) in viruliferous whiteflies. The weighted gene correlation network analysis of viruliferous B. tabaci MEAM1 and MED adults identified different hub genes potentially implicated in the vector competence and circulative tropism of viruses. Taken together, the results indicate that both vector cryptic species and the acquired virus species could differentially affect gene expression.


Introduction
Bemisia tabaci (Gennadius) (Hemiptera: Aleyrodidae), commonly known as the sweetpotato whitefly, is widely distributed in the subtropics/tropics, where it can be an agricultural pest [1][2][3]. Feeding by this pest induces physiological symptoms and results in the transmission of a devastating group of plant viruses belonging to the genus Begomovirus [4][5][6][7][8]. Bemisia tabaci comprises several cryptic species with distinct phylogeographical distributions [9][10][11][12][13]. These cryptic species and their haplotypes are morphologically In this study, differential gene expression in B. tabaci MEAM1 and MED was assessed by transcriptome analysis post acquisition of three begomoviruses, the Old World monopartite TYLCV, and the New World bipartite SiGMV and CuLCrV. Libraries were prepared for viruliferous and non-viruliferous whiteflies given a 72 h acquisition access period (AAP) on the respective begomovirus-infected or non-infected host plants, followed by a 72 h feeding access on cotton, a non-host of the three begomoviruses. The objectives of this study were to (i) assess the differences in gene expression between B. tabaci MEAM1 and MED upon the acquisition of Old and New World begomoviruses, and (ii) locate putative hub genes and co-expressed genes (modules) responsive to the acquisition of Old and New World begomoviruses using weighted gene correlation network analysis (WGCNA).
The B. tabaci MEAM1 colony (Genbank accession number: MN970031) was established from individuals collected in 2009 from cotton plants in Tifton, Georgia. This colony was reared on cotton plants in 10 cm diameter × 8 cm tall pots in whitefly-proof cages in the greenhouse under the conditions described above. The B. tabaci MED colony (GenBank accession number: MZ469725) was established on cotton plants in 2017 from individuals collected from poinsettia (Euphorbia pulcherrima Wild. Ex Klotzsch) plants from a nursery located in North Georgia, USA and maintained by Prof. Ronald D. Oetting, College of Agriculture and Environmental Sciences, University of Georgia, Griffin, GA. The B. tabaci MED colony was reared under the same conditions as those reported for the B. tabaci MEAM1 colony. The identities of both colonies were examined by PCR amplification and by partially sequencing the mitochondrial cytochrome oxidase I (COI) gene in the 3 barcode region [55,56]. The B. tabaci MEAM1 and MED colonies were maintained in separate greenhouses to avoid population admixtures, and their identities were periodically screened (every alternate month) to ensure no contamination had occurred.
The virus isolates used in this study are described extensively in Gautam et al. [36]. Briefly, the TYLCV isolate was collected in 2009 from a commercial tomato farm located in Montezuma in Macon County in GA, USA. TYLCV has since been maintained in TYLCVsusceptible tomato cultivar (Florida 47) through B. tabaci MEAM1-mediated inoculation. The SiGMV isolate was collected in 2018 from prickly sida plants on field edges at the UGA Horticulture Hill Farm in Tifton (Tift County, GA, USA). SiGMV has since been maintained in prickly sida plants through B. tabaci MEAM1-mediated inoculation. The CuLCrV isolate used in this study was collected in 2016 from infected yellow summer squash plants from the UGA Horticulture Hill Farm in Tifton (Tift County, GA, USA). CuLCrV has since been maintained in yellow summer squash plants through B. tabaci MEAM1-mediated inoculation. Virus-infected plants were obtained by providing B. tabaci MEAM1 a 72 h acquisition access period (AAP) on TYLCV/SiGMV/CuLCrV-infected plants. Subsequently, viruliferous B. tabaci MEAM1 (~100/plant) were attached to the non-infected~10 cm tall seedlings (tomato, squash, or prickly sida) using clip cages. Inoculated individual plants were placed in separate insect-proof cages under the conditions described above. A week post inoculation, all of the leaves were removed. The infection status of plants in insect-free newly-emerged leaves was evaluated at~3 weeks post inoculation by endpoint PCR using species-specific primers [36,39].

Feeding Assay, DNA & RNA Isolation, and RNA Sequencing
Approximately 1000 B. tabaci MEAM1 or MED newly emerged (1-3 days old) whiteflies were collected from the above established colonies and then introduced to either noninfected or virus-infected tomato, prickly sida, and squash plants for an acquisition access period (AAP) of 72 h. After 72 h on either non-infected or virus-infected tomato, prickly sida, or squash plants, the whiteflies (MEAM1 and MED) were transferred to cotton plants for another 72 h to minimize the host effect from either non-infected or infected plants. After 72 h on cotton,~100 live whiteflies in three to five biological replicates (per treatment) were collected and immediately stored at −20 • C for 20 min and thereafter used for RNA extraction on the same day. Another set of 20 individual whiteflies (MEAM1 or MED) per treatment were collected and used to confirm the presence/absence of the virus. The results indicated that 90-100%, 80-90%, or 70-90% of the B. tabaci MEAM1 and 95-100%, 65-85%, or 60-80% of the MED whiteflies were positive for TYLCV, SiGMV, or CuLCrV, respectively. As expected, individual B. tabaci MEAM1 or MED that were given a 72 h feeding access on non-infected tomato, prickly sida, or squash plants (and for 72 h on cotton plants) tested negative for TYLCV, SiGMV, and CuLCrV, respectively.
The total whitefly DNA was extracted from individuals using an InstaGene Matrix containing 6% Chelex resin (Bio-Rad, Hercules, CA, USA), as previously described by Gautam et al. [30]. The extracted DNA was stored at −20 • C until it was used. PCR reactions containing 2X GoTaq Green Master Mix (Promega, Madison, WI, USA) and primers specific to TYLCV, SiGMV, or CuLCrV (Table S1) were performed using a T-100 thermocycler (Bio-Rad, Hercules, CA, USA). The total RNA was extracted using the Ambion TRIzol Reagent (Thermo Fisher, Waltham, MA, USA), according to the manufacturer's instructions, and purified using the PureLink TM RNA Mini Kit (Thermo Fisher, Waltham, MA, USA), according to the manufacturer's instructions. Approximately 50 µL of RNA was shipped to Novogene Corporation Inc. (Sacramento, CA, USA), where RNA sample quality control (QC), mRNA library preparation, and RNA sequencing were carried out.

RNA Sequencing, Transcriptome Assembly, and Analysis
RNA sample quality control (QC) was carried out using a Nanodrop for preliminary quantitation. Agarose gel electrophoresis was used to test RNA degradation and potential contamination, while the Agilent 2100 was used to assess RNA integrity (RIN) and quantitation. The samples that passed QC (RIN value > 6.8 and concentration of >20 ng/uL) were used for library construction. In brief, the library construction involved enriching mRNA using oligo(dT) beads and the removal of rRNA using the Ribo-Zero kit. Subsequently, the mRNA was fragmented and followed by first-and second-strand cDNA synthesis. Finally, cDNA libraries were generated through adaptor ligation and PCR enrichment. The library QC consisted of three steps: testing library concentration using Qubit 2.0, testing the insert size using the Agilent 2100, and quantifying the library concentration using qPCR. The libraries that passed QC were sequenced on a NovaSeq 6000 Sequencing System (Illumina, San Diego, CA, USA) using the NovaSeq paired end 150 sequencing setting.
The weighted gene co-expression network (WGCNA) software v1.70-3 was run using the R software v4.1.0 [65,67]. Prior to the construction of the co-expression networks, filtering was performed to select the top 95% quantile of DEGs from the 25 and 27 samples of B. tabaci MEAM1 and MED, respectively. The soft-thresholding power modules (1 to 40) was tested by the gradient independent method, and the scale-independent condition of the signed R2 was set to 0.90. Gene modules were constructed with a power of 28, with a minimum module size set to 20. The interaction relationships across the co-expression modules were used to construct a topological overlap matrix (TOM) by using correlation expression values. The TOM was represented with a dendrogram by setting the parameters mergeCutHeight = 0.15 and detectCutHeight = 0.995. The module color was randomly assigned, and the first principal component from each module was used to calculate the module eigengene. A module-trait relationship heatmap was generated using the labeled-Heatmap package available in WGCNA software to show the topological overlap of the co-expression modules based on the module eigengenes. The labeledHeatmap parameters included matrix = moduleTraitCor and main = paste ("Module-trait relationships"). The network analysis of the top 30 genes from the B. tabaci MEAM1 and MED turquoise modules was visualized with Cytoscape v3.9.0 [68].

Validation of RNA Sequencing Data
To validate the differential expression analysis, 10 genes were randomly selected for each virus treatment (n = 56), and their expression levels were compared between viruliferous and non-viruliferous whiteflies (B. tabaci MEAM1 or MED) for the three viruses (TYLCV, SiGMV, and CuLCrV). For each sample, an aliquot (10 µL) of the RNA was used for validation by reverse transcription quantitative PCR (RT-qPCR), and the remainder were used for RNA Sequencing. The primer sequences used for the selected genes are listed in Table S1. The total RNA for each sample (500 ng) was reverse transcribed into cDNA using the GoScript TM Reverse Transcription System (Promega, Madison, WI, USA), according to the manufacturer's instructions. Quantitative PCR reactions were assembled in 12.5 µL duplicate reactions for three biological replicates and carried out using 2X GoTaq qPCR Master Mix (Promega, Madison, WI, USA) in a Mastercycler ® ep realplex thermal cycler (Eppendorf, Hauppauge, NY, USA). An initial denaturation step (2 min at 94 • C) was followed by 40 cycles of denaturation (20 s at 94 • C), primer annealing (15 s at 58-60 • C, depending on the primer set), and extension (15 s at 68 • C). Finally, the melting curve analysis was conducted to evaluate the specificity of the fluorescence signal. Relative expression levels were calculated using the 2 −∆∆Ct method, and the expression level of each gene was normalized to the expression level of β-actin, a whitefly reference gene [69,70].

Summary of RNA Sequencing
Three to five biological replicates were included per virus treatment, resulting in a total of 25 and 27 libraries constructed for B. tabaci MEAM1 and MED adults, respectively. Raw read pairs for the generated libraries ranged from~19 to 33 million. After trimming and removing the reads that aligned with the ribosomal RNA, the mitochondrion genome, and the three bacterial endosymbiont genomes (Candidatus Portiera aleyrodidarum, Hamiltonella, and Rickettsia),~17 to 31 million read pairs were retained in the libraries (Table 1). Bemisia tabaci MEAM1 and MED cleaned read pairs for the different libraries mapped onto the B. tabaci MEAM1 genome by~90 to 95% and 88 to 91%, respectively ( Table 1). The Pearson's correlation coefficients for all of the biological replicates ranged from 0.9 to 1 ( Table S2), suggesting that the data were highly reproducible.        The RNA sequencing-based differential gene expression results were validated via RT-qPCR on 10 randomly selected differentially expressed genes per whitefly species-virus treatment (n = 56), using the β-actin gene as an internal standard (Table S3). The results showed that the expression trends obtained from RNA sequencing and RT-qPCR were highly consistent (Table S3).

Common DEGs among B. tabaci MEAM1 and MED Adults
Of the 881 and 559 DEGs found in B. tabaci MEAM1 and MED adults, respectively, 146 were found to be differentially expressed in common following the acquisition of the three viruses ( Figure 1a). Seventy genes-almost half of the common DEGs between B. tabaci MEAM1 and MED adults-were unknown, whereas the other 76 genes were annotated (http://www.whiteflygenomics.org/cgi-bin/bta/index.cgi accessed on 15 December 2021). All of the common DEGs (146) between B. tabaci MEAM1 and MED adults were analyzed using functional analysis tools (Blast, InterProScan, Blast2GO Mapping, and Blast2GO Annotation) available in the OmicsBox version 1.4.11. Only 72 genes (69 annotated and 3 unknown) were assigned functional groups under three classification systems: biological process (44 genes), molecular function (60 genes), and cellular component (33 genes). Fourteen Gene Ontology (GO) terms were assigned under the biological process category, of which the cellular process and the establishment of localization were significant ( Figure 3a). Eight GO terms were assigned under the molecular function category, of which catalytic activity, binding, heterocyclic compound binding, and transport activity were significant ( Figure 3b). In the cellular component category, two GO terms were significant, and they were membrane and intrinsic component of membrane ( Figure 3c). The RNA sequencing-based differential gene expression results were validated via RT-qPCR on 10 randomly selected differentially expressed genes per whitefly species-virus treatment (n = 56), using the β-actin gene as an internal standard (Table S3). The results showed that the expression trends obtained from RNA sequencing and RT-qPCR were highly consistent (Table S3).

Common DEGs among B. tabaci MEAM1 and MED Adults
Of the 881 and 559 DEGs found in B. tabaci MEAM1 and MED adults, respectively, 146 were found to be differentially expressed in common following the acquisition of the three viruses ( Figure 1a). Seventy genes-almost half of the common DEGs between B. tabaci MEAM1 and MED adults-were unknown, whereas the other 76 genes were annotated (http://www.whiteflygenomics.org/cgi-bin/bta/index.cgi accessed on 15 December 2021). All of the common DEGs (146) between B. tabaci MEAM1 and MED adults were analyzed using functional analysis tools (Blast, InterProScan, Blast2GO Mapping, and Blast2GO Annotation) available in the OmicsBox version 1.4.11. Only 72 genes (69 annotated and 3 unknown) were assigned functional groups under three classification systems: biological process (44 genes), molecular function (60 genes), and cellular component (33 genes). Fourteen Gene Ontology (GO) terms were assigned under the biological process category, of which the cellular process and the establishment of localization were significant ( Figure 3a). Eight GO terms were assigned under the molecular function category, of which catalytic activity, binding, heterocyclic compound binding, and transport activity were significant (Figure 3b). In the cellular component category, two GO terms were significant, and they were membrane and intrinsic component of membrane ( Figure 3c).  following the acquisition of viruses. Cluster representatives in a two-dimensional space were derived by applying multidimensional scaling to a matrix of the gene ontology terms' semantic similarities. Bubble color indicates the p-value, and size indicates the frequency of the GO term in the underlying GOA database.
Only 28.8% (42 of 146) of the total number of common DEGs between B. tabaci MEAM1 and MED adults were annotated using the KEGG analysis [63] and mapped to 74 biochemical pathways (Table S4a,b). The most represented category pathways were (i) metabolism, (ii) organismal systems, (iii) human diseases/pathogens, and (iv) cellular processes.

Unique DEGs in B. tabaci MEAM1 and MED Adults
Of the 881 DEGs found in the B. tabaci MEAM1 adults, 735 were found to be uniquely and differentially expressed following the virus acquisition ( Figure 1a). Using the functional analysis tools available in the OmicsBox version 1.4.11, 526 of the 735 genes were assigned functional groups under three classification systems: biological process (290 genes), molecular function (387 genes), and cellular component (296 genes). Twenty-two GO terms were assigned under the biological process category, of which 20 terms were significant ( Figure 4a). Ten GO terms were assigned under the molecular function category, five (binding, protein binding, hydrolase activity, heterocyclic compound binding, and organic cyclic compound binding) of which were significant ( Figure 4b). In the cellular component category, sixteen GO terms were identified, and 12 of those were significant ( Figure 4c).  (Table S4a-b). The most represented category pathways were (i) metabolism, (ii) organismal systems, (iii) human diseases/pathogens, and (iv) cellular processes.

Unique DEGs in B. tabaci MEAM1 and MED Adults
Of the 881 DEGs found in the B. tabaci MEAM1 adults, 735 were found to be uniquely and differentially expressed following the virus acquisition (Figure 1a). Using the functional analysis tools available in the OmicsBox version 1.4.11, 526 of the 735 genes were assigned functional groups under three classification systems: biological process (290 genes), molecular function (387 genes), and cellular component (296 genes). Twenty-two GO terms were assigned under the biological process category, of which 20 terms were significant ( Figure 4a). Ten GO terms were assigned under the molecular function category, five (binding, protein binding, hydrolase activity, heterocyclic compound binding, and organic cyclic compound binding) of which were significant (Figure 4b). In the cellular component category, sixteen GO terms were identified, and 12 of those were significant ( Figure 4c). Of the 559 DEGs found in B. tabaci MED adults, 413 were found to be uniquely and differentially expressed following the virus acquisition (Figure 1a). Using functional analysis tools, 223 genes of the 413 genes were assigned functional groups under three classification systems: biological process (120 genes), molecular function (179 genes), and cellular component (107 genes). Sixteen GO terms were assigned under the biological process category, of which four (biogenesis, cellular process, cellular component organization, and localization) were significant (Figure 5a). Eight GO terms were assigned under the molecular function category, of which seven were significant (Figure 5b). In the cellular component category, eleven GO terms were identified, of which three (membrane, intrinsic component of membrane, and protein-containing complex) were significant (Figure 5c).  Only 44.8% (329 of 735) of the total number of DEGs in the B. tabaci MEAM1 adults were annotated using the KEGG analysis [63] and were mapped to 301 biochemical pathways (Table S4c-d). The most represented pathway categories were (i) metabolism, (ii) human diseases, and (iii) organismal systems. For the B. tabaci MED adults, 42.4% (175 of 413) of the total number of DEGs were annotated using the KEGG analysis [63] and were mapped to 162 biochemical pathways (Table S4e,f). Similar to the MEAM1 whiteflies, metabolism, human diseases/pathogens, and organismal systems were the most represented pathway categories.

Co-Expression Networks from B. tabaci MEAM1 and MED Adults
The co-expression of the genes from the B. tabaci MEAM1 and MED adults that acquired TYLCV, SiGMV, or CuLCrV were identified in relation to the expression data. Weighted gene co-expression network analysis (WGCNA), which clusters genes into modules based on weighted gene-gene interactions, was used to study the co-expression. For B. tabaci MEAM1, eight modules were identified (Figure 6a), and a Pearson's correlation coefficient analysis showed the connections between the three virus species ( Figure  6b). The B. tabaci MEAM1 interactions for the largest module, MEturqoiuse, were checked for the top interacting genes among the 30 identified genes (Figure 6c). The top four most highly connected genes were: Bta03319 (Tat-linked quality control protein TatD), Bta00139 (Riboflavin kinase), Bta15228 (Ubiquitin-conjugating enzyme E2 L3, putative), and Bta06777 (GRIP and coiled-coil domain-containing protein, putative). In the B. tabaci MED analysis, nine modules were identified (Figure 7a), and a Pearson's correlation coefficient analysis showed the connections between the three virus species (Figure 7b). The B. tabaci MED interactions for the largest module, MEturqoiuse, were checked for the top Cluster representatives in a two-dimensional space were derived by applying multidimensional scaling to a matrix of the gene ontology terms' semantic similarities. Bubble color indicates the p-value, and size indicates the frequency of the GO term in the underlying GOA database.
Only 44.8% (329 of 735) of the total number of DEGs in the B. tabaci MEAM1 adults were annotated using the KEGG analysis [63] and were mapped to 301 biochemical pathways (Table S4c,d). The most represented pathway categories were (i) metabolism, (ii) human diseases, and (iii) organismal systems. For the B. tabaci MED adults, 42.4% (175 of 413) of the total number of DEGs were annotated using the KEGG analysis [63] and were mapped to 162 biochemical pathways (Table S4e,f). Similar to the MEAM1 whiteflies, metabolism, human diseases/pathogens, and organismal systems were the most represented pathway categories.

Co-Expression Networks from B. tabaci MEAM1 and MED Adults
The co-expression of the genes from the B. tabaci MEAM1 and MED adults that acquired TYLCV, SiGMV, or CuLCrV were identified in relation to the expression data. Weighted gene co-expression network analysis (WGCNA), which clusters genes into modules based on weighted gene-gene interactions, was used to study the co-expression. For B. tabaci MEAM1, eight modules were identified (Figure 6a), and a Pearson's correlation coefficient analysis showed the connections between the three virus species (Figure 6b). The B. tabaci MEAM1 interactions for the largest module, MEturqoiuse, were checked for the top interacting genes among the 30 identified genes (Figure 6c). The top four most highly connected genes were: Bta03319 (Tat-linked quality control protein TatD), Bta00139 (Riboflavin kinase), Bta15228 (Ubiquitin-conjugating enzyme E2 L3, putative), and Bta06777 (GRIP and coiled-coil domain-containing protein, putative). In the B. tabaci MED analysis, nine modules were identified (Figure 7a), and a Pearson's correlation coefficient analysis showed the connections between the three virus species (Figure 7b). The B. tabaci MED interactions for the largest module, MEturqoiuse, were checked for the top interacting genes among the 30 identified genes (Figure 7c). The top four most highly connected genes were: Bta03156 (Serine protease 7, isoform), Bta03155 (Serine protease), Bta03265 (Serine protease 7, isoform A), and Bta06572 (unknown protein).
connected genes were all unknown proteins (Bta14024, Bta04891, and Bta08994). For SiGMV, six modules were identified (Figure S2a), and a Pearson's correlation coefficient analysis showed the connections between the two whitefly cryptic species ( Figure S2b). SiGMV interactions for the largest module, MEturqoiuse, were checked for the top interacting genes among the 30 identified genes ( Figure S2c). The top four most highly connected genes were all unknown proteins (Bta10403, Bta01241, and Bta04118), except for Bta13010 (Phosphatidylethanolamine-binding protein). For CuLCrV, five modules were identified (Figure S3a), and a Pearson's correlation coefficient analysis showed the connections between B. tabaci MEAM1 and MED ( Figure S3b). CuLCrV interactions for the largest module, MEturqoiuse, were checked for the top interacting genes among the 30 identified genes ( Figure S3c). The top three most highly connected genes were: Bta02748 (AP-3 complex subunit mu-1), Bta06159 (Ribosomal protein L19), and Bta14433 (unknown protein).    For TYLCV, six modules were identified (Figure S1a), and a Pearson's correlation coefficient analysis showed the connections between B. tabaci MEAM1 and MED ( Figure S1b). TYLCV interactions for the largest module, MEturqoiuse, were checked for the top interacting genes among the 30 identified genes ( Figure S1c). The top three most highly connected genes were all unknown proteins (Bta14024, Bta04891, and Bta08994). For SiGMV, six modules were identified (Figure S2a), and a Pearson's correlation coefficient analysis showed the connections between the two whitefly cryptic species ( Figure S2b). SiGMV interactions for the largest module, MEturqoiuse, were checked for the top interacting genes among the 30 identified genes ( Figure S2c). The top four most highly connected genes were all unknown proteins (Bta10403, Bta01241, and Bta04118), except for Bta13010 (Phosphatidylethanolamine-binding protein). For CuLCrV, five modules were identified (Figure S3a), and a Pearson's correlation coefficient analysis showed the connections between B. tabaci MEAM1 and MED ( Figure S3b). CuLCrV interactions for the largest module, MEturqoiuse, were checked for the top interacting genes among the 30 identified genes ( Figure S3c). The top three most highly connected genes were: Bta02748 (AP-3 complex subunit mu-1), Bta06159 (Ribosomal protein L19), and Bta14433 (unknown protein).
Following the overview of the functional annotation and pathway analyses of the B. tabaci MEAM1 and MED genes, genes that could influence virus-vector interaction were examined.
3.6. DEGs among B. tabaci MEAM1 and MED Adults Associated with Virus-Vector Interactions 3.6.1. Virus Infection Among the common DEGs between B. tabaci MEAM1 and MED adults, KEGG annotations revealed three generic biotic stress-induced genes ( Table 2) under the "human disease" category. These genes, all underexpressed except in one instance (Table 2), were associated with the infection of two viruses: Human T-cell leukemia virus 1 (Bta03437-Zinc finger protein) and Measles (Bta02903-Heat shock protein 70 and Bta09867-70 kDa heat shock protein) ( Table S4b).  In B. tabaci MEAM1 adults specifically, 27 genes were associated with different viruses including Human T-cell leukemia virus 1 (four genes), Human immunodeficiency virus 1 (6 genes), Hepatitis B virus (four genes), Hepatitis C virus (one gene), Coronavirus (one gene), Influenza A virus (two genes), Measles virus (three genes), Herpes simplex virus 1 (two genes), Human cytomegalovirus (eight genes), Kaposi sarcoma-associated herpesvirus (three genes), Epstein-Barr virus (three genes), and Human papillomavirus (nine genes) ( Table S4d). The majority of the genes (21/27) associated with virus infection were overexpressed, with the Ras-like/related (four genes) and heat shock protein (three genes) genes being the most predominant (Table 2).
Similarly, seven genes, four overexpressed and three underexpressed, were associated with the acquisition of three begomoviruses in B. tabaci MED ( Table 2). The differentially expressed genes were associated with the infection of Human T-cell leukemia virus 1 (one gene), Human immunodeficiency virus 1 (three genes), Coronavirus (one gene), Measles virus (1 gene), Herpes simplex virus 1 (four genes), Human cytomegalovirus (two genes), and Epstein-Barr virus (two genes) ( Table S4f).

Signal Transduction
Among the shared DEGs identified in B. tabaci MEAM1 and MED adults, KEGG annotation identified four genes. All but one (Bta02903) were underexpressed (Table 3) and were associated with three signal transduction pathways, namely, AMPK (Bta07645), Apelin (Bta03437), and MAPK (Bta02903 and Bta09867) ( Table S4b). Table 3. Differential expression of genes associated with signal transduction in viruliferous compared with non-viruliferous Bemisia tabaci Middle East-Asia Minor 1 and Mediterranean (MED) adults.

Gene ID Annotation FC in Whitefly & Virus Treatment
DEGs shared by MEAM1 and MED  Fifty-one genes (35 overexpressed and 16 underexpressed) were annotated using the KEGG analysis (Table 3) and associated with 29 signal transduction pathways in MEAM1 whiteflies (Table S4d). The pathways and their associated genes are included below (Figure 8). Some of the prevalent genes in the MEAM1 whiteflies signaling pathways were Ras-like/related genes, protein kinase, acyl-CoS desaturase, heat shock, and guanine nucleotide-binding (Table 3). In the B. tabaci MED adults, 11 genes (4 overexpressed and 7 underexpressed) were associated with 14 signaling pathways using KEGG annotation (Table 3) (Table S4f). These pathways and the associated genes are included below (Figure 9).  In the B. tabaci MED adults, 11 genes (4 overexpressed and 7 underexpressed) were associated with 14 signaling pathways using KEGG annotation (Table 3) (Table S4f). These pathways and the associated genes are included below (Figure 9).

Signaling Molecules and Virus Interaction
Using KEGG annotation, five putative receptor genes (two overexpressed and three underexpressed) were identified in B. tabaci MEAM1 adults. They include: Bta03110, Bta01738-Corticotropin-releasing factor receptor 1, Bta15373-Neuropeptide FF receptor 2 under the neuroactive ligand-receptor interaction, Bta04818 under the cytokine-cytokine receptor interaction, and Bta07388-Syndecan under the cell adhesion molecules (Table S4d). The genes Bta03110 and Bta04818 were also associated with both signal transduction and/or virus infection (Tables 2 and 3). Putative receptor genes were not identified in B. tabaci MED adults' transcriptome data.

Immune Systems
Among the DEGs shared by both B. tabaci MEAM1 and MED adults, 10 genes (6 upand underexpressed depending on whitefly/virus treatment and 4 underexpressed) under the organismal systems were associated with three immune system pathways by KEGG annotation ( Table 4). Seven of the ten genes were cathepsin B, while the other three were protein spaetzle and heat shock proteins (Table 4). These genes belonged to the antigen processing and presentation (seven cathepsin B genes) pathway, NOD-like receptor signaling pathway (seven cathepsin B and two heat shock protein genes), and Toll and Imd signaling pathways (protein spaetzle gene) ( Table S4b). Table 4. Differential expression of genes associated with immune system categories in viruliferous compared with non-viruliferous Bemisia tabaci Middle East-Asia Minor 1 (MEAM1) and Mediterranean (MED) adults.

Gene ID Annotation FC in Whitefly & Virus Treatment
DEGs shared by MEAM1 and MED  A total of 33 genes (21 overexpressed and 12 underexpressed) were annotated using the KEGG analysis (Table 4) and associated with 17 immune system pathways in B. tabaci MEAM1 adults (Table S4d). Cathepsin B (seven genes) and Ras-like/related (five genes) genes were the most predominant in this category (Table 4).

Cellular Processes (Apoptosis, Lysosome, and Phagosome)
Among the DEGs shared by B. tabaci MEAM1 and MED adults, eight cathepsin genes were associated with the apoptosis and lysosome pathways by KEGG annotation ( Table S4b). All of those genes were also associated with immune systems, except for Bta20005-cathepsin F ( Table 4). The Bta05445-protein transport protein Sec61 subunit β was the only gene associated with phagosome among the shared DEGs identified in B. tabaci MEAM1 (fold change: −0.93, treatment: SiGMV) and MED (fold change: 0.89, treatment: TYLCV) adults.
A total of 34 genes (17 overexpressed and 17 underexpressed) were annotated by KEGG analysis (Table 5) and associated with cellular processes in B. tabaci MEAM1 adults (Table S4d). Cathepsin F/F-like (nine genes) and cathepsin B (seven genes) genes were the most predominant in this category (Table 5), whereas, among the apoptosis, lysosome, and phagosome pathways, 22, 23, and 7 genes were identified, respectively. Table 5. Differential expression of genes associated with cellular processes (apoptosis, lysosome, and phagosome) in viruliferous compared with non-viruliferous Bemisia tabaci Middle East-Asia Minor 1 (MEAM1) and Mediterranean (MED) adults.

Gene ID Annotation FC in Whitefly & Virus Treatment
DEGs specific to MEAM1  In B. tabaci MED adults, 19 genes (6 overexpressed, 12 underexpressed, one in both directions) were associated with cellular processes by KEGG annotation (Table 5). Except for six genes, all were cathepsins (13 genes). The genes in the lysosome (17 genes) and apoptosis (13 genes) pathways were more represented than the genes associated with the phagosome category (four genes) ( Table S4f).

Longevity
Based on the KEGG pathway annotation, three genes (two overexpressed/underexpressed and one underexpressed) were associated with two aging pathways among the DEGs shared between B. tabaci MEAM1 and MED adults. These pathways included longevity regulating pathway-worm (Bta14177-fatty acyl-CoA reductase 1) and longevity regulating pathway-multiple species (Bta02903 and Bta09867) ( Table S4b). The Bta02903 and Bta09867 genes were associated with several pathways (Tables 2-4). Bta14177 was overexpressed (fold change: 1.59) in B. tabaci MEAM1 adults that acquired SiGMV, while it was underexpressed (fold change: −1.93) in B. tabaci MED adults that acquired TYLCV.
Fifteen genes (five overexpressed and ten underexpressed) were annotated using the KEGG analysis (Table 6) and associated with aging in B. tabaci MEAM1 adults (Table S4d). Acyl-CoA related (five genes) and heat shock protein (four genes) genes were the most predominant in this category (Table 6). Three, nine, and four genes were classified under longevity regulating pathway, longevity regulating pathway-worm, and longevity regulating pathway-multiple species, respectively. Table 6. Differential expression of genes associated with longevity in viruliferous compared with non-viruliferous Bemisia tabaci Middle East-Asia Minor 1 (MEAM1) and Mediterranean (MED) adults.

Gene ID Annotation FC in Whitefly & Virus Treatment
DEGs specific to MEAM1

Bta00981
Cyclic AMP-responsive element-binding protein 3-like protein 2 1.60 for SiGMV In B. tabaci MED adults, five genes (two overexpressed and three underexpressed) were identified in the aging pathways based on KEGG annotation (Table 6). Four acyl-CoA-related genes and the heat shock protein gene were grouped under longevity regulating pathway-worm and longevity regulating pathway-multiple species, respectively (Table S4f).

Discussion
This study assessed the differences in gene expression in B. tabaci MEAM1 and MED whiteflies in response to the acquisition of Old and New World begomoviruses. Both the Old and New World viruses examined in this study were transmitted by the B. tabaci MEAM1 cryptic species, and only the Old World virus (TYLCV) was transmitted by the B. tabaci MED cryptic species [36]. Gene expression differences were compared between the two cryptic species following the acquisition of TYLCV, SiGMV, or CuLCrV. Across all three virus treatments, more transcriptional changes were recorded with B. tabaci MEAM1 adults (881 genes) than with B. tabaci MED adults (559 genes), which represented only 5.6% and 3.6%, respectively, of the overall number of genes (15,668) [71]. Two other studies using the same whitefly cryptic species and virus species reported different numbers of differentially expressed genes in B. tabaci MEAM1 that acquired TYLCCNV when compared with their non-viruliferous counterparts. There were 457 genes in the study of Luan et al. [50], whereas there were 1606 genes in the study of Luan et al. [44]. Similarly, in the TYLCV pathosystem, two independent studies reported a contrasting number of differentially expressed genes (79 genes versus 1347 genes) in viruliferous B. tabaci MEAM1 compared with non-viruliferous, while another study identified 78 differentially expressed genes in viruliferous B. tabaci MED, compared with the non-viruliferous treatment [47][48][49]. Further, contrasting transcriptional changes in the whitefly have been observed in semi-persistently transmitted criniviruses. For example, for the tomato chlorosis virus (ToCV) pathosystem, one study reported 221 differentially expressed genes in viruliferous B. tabaci MED adults compared with non-viruliferous adults, whereas another study reported 1155 differentially expressed genes in viruliferous B. tabaci MEAM1 adults compared with non-viruliferous adults [45,49]. In another crinivirus pathosystem, cucurbit yellow stunting disorder virus, 262 differentially expressed genes were reported in viruliferous B. tabaci MEAM1 compared with non-viruliferous counterparts [46]. The differences in the experimental design and analysis, such as: (i) AAP from 1-7 days, (ii) viruliferous whiteflies with no gut clearing to minimize potential indirect host plant effects, (iii) low throughput compared with high throughput sequencing platforms, (iv) the different number of libraries sequenced, and (v) bioinformatics analyses, could have partly contributed to the observed variations in the number of transcriptional changes observed in whiteflies in different studies.
This study investigated the transcriptional changes in B. tabaci MEAM1 and MED induced by three plant viruses belonging to the same genus, Begomovirus, and differential gene expression patterns varied between them. More transcriptional changes were associated with TYLCV acquisition (MEAM1-462 genes and MED-413 genes), followed by SiGMV acquisition (MEAM1-459 genes and MED-165 genes), whereas CuLCrV acquisition (MEAM1-69 genes and MED-44 gene) induced the least number of transcriptional changes. Both TYLCV and SiGMV accumulated at significantly higher levels (copies per ng DNA) in both their host plants and vector's tissues (midgut, hemolymph, and salivary glands), while CuLCrV accumulated at significantly lower amounts in its host plant and vector's tissues [30,36]. Bemisia tabaci MEAM1 transmits SiGMV and CuLCrV, and their circulative tropism in whiteflies could have facilitated binding with different whitefly proteins and putative receptors, leading to increased transcriptional changes [22,72,73]. Bemisia tabaci MED did not transmit the two New World viruses and accumulated at reduced amounts in its tissues; hence, the reduced circulative tropism of New World viruses could have influenced the lower transcriptional responses compared to TYLCV. The differences in B. tabaci MEAM1 and MED in terms of acquisition and the competency to inoculate different begomoviruses observed here and in previous studies [74,75] support previous hypotheses surrounding B. tabaci-begomovirus transmission specificity.
The acquisition of three begomoviruses by B. tabaci MEAM1 and MED resulted in transcriptional changes, of which 146 genes were differentially expressed in common between the two invasive B. tabaci species. KEGG analysis annotated only 42 of the 146 genes, and the most represented pathways were classified under metabolic pathways involving carbohydrate and lipid metabolism. Virus (hepatitis C virus and HIV) infection and replication have been reported to increase carbohydrate and lipid metabolism in infected CD4 + T cells [76,77]. The majority of the genes associated with lipid metabolism in both B. tabaci MEAM1 and MED were underexpressed, and those associated with carbohydrate metabolism were underexpressed in B. tabaci MEAM1 and partly overexpressed in B. tabaci MED. The results herein are consistent with previous reports of the downregulation of genes associated with lipid metabolism in B. tabaci MEAM1 or MED post acquisition of ToCV, TYLCCNV, and TYLCV, or of TYLCV and ToCV from co-infected plants [44,49].
In contrast, some studies have reported the upregulation of genes associated with lipid metabolism in a circulative and propagative virus study system (tomato spotted wilt orthotospovirus-thrips) [78,79]. Evidence for reduced gene expression associated with carbohydrate and lipid metabolism in viruliferous B. tabaci MEAM1 and MED may offer another line of evidence to counter the hypothesis that begomoviruses replicate in their whitefly vector [80,81]. Nevertheless, the temporary replication of at least one begomovirus (TYLCV) in B. tabaci MEAM1 salivary glands has been reported in a recent study [82]. The replication affected TYLCV accumulation in the salivary glands minimally given the overall virus accumulation within whiteflies.
The weighted gene correlation network analysis (WGCNA) has been applied to analyze various biological processes for different organisms [83,84]. In this study, the gene co-expression network analysis of B. tabaci MEAM1 and MED adults that acquired TYLCV, SiGMV, or CuLCrV identified several candidate genes for studying the vector competence of two invasive whiteflies. WGCNA provided scale-free networks between genes within modules, and these modules aided in identifying genes that make B. tabaci MEAM1 a more competent vector of New World begomoviruses (SiGMV and CuLCrV) than its sister cryptic species B. tabaci MED. For instance, among the top four interacting genes identified in the largest module for B. tabaci MEAM1 was the gene Bta15228 (ubiquitin-conjugating enzyme E2 L3, putative). The ubiquitin-conjugating enzyme was required for Notch signaling activation during Drosophila wing development. In addition, this gene was involved in the endocytic trafficking of the Notch protein [85]. This gene could also be involved in the endocytic trafficking of virus particles in B. tabaci MEAM1. In B. tabaci MED, three serine protease genes (Bta03155, Bta03156, and Bta03265) were identified in the largest module among the top four interacting genes. Serine proteases are proteolytic enzymes responsible for digestion, larval-pupal molting, signal transduction, and immune responses in insects [86][87][88][89]. The identification of serine protease genes among the top interacting genes in B. tabaci MED adults could be signatures of heightened innate immune responses such as melanization and the production of antimicrobial peptides triggered by the acquisition of begomoviruses. The ubiquitin-conjugating enzyme and serine protease genes identified in B. tabaci MEAM1 and MED adults, respectively, could be important targets for future investigation.
Similarly, the gene co-expression network analysis for TYLCV, SiGMV, or CuLCrV acquired by B. tabaci MEAM1 and MED adults was examined with the purpose of identifying candidate genes involved in the acquisition of Old and New World begomoviruses. None of the top three interacting genes identified in the largest module for TYLCV, an Old World begomovirus, were functionally annotated. In contrast, a majority of the top three/four interacting genes identified in the largest module for New World begomoviruses-namely, SiGMV and CuLCrV-were functionally annotated. For instance, among the top four interacting genes identified in the largest module for SiGMV was the gene Bta13010phosphatidylethanolamine-binding protein (PEBP). The downregulation of PEBP genes in the Bombyx mori strain resistant to Bombyx mori nucleopolyhedrovirus (BmNPV) upon infection induced enhanced apoptosis, thereby repressing the ability of BmNPV to infect other cells [90]. In B. tabaci MED adults that transmitted TYLCV, the relative expression of PEBP4 was increased significantly and was shown to interact with the TYLCV coat protein, putatively influencing apoptosis and autophagy mechanisms in the whitefly [91]. The potential for the direct role of PEBP4 in B. tabaci MEAM1 capable of transmitting SiGMV remains to be investigated; however, based on the increased relative expression observed in this study and the prior inference, PEBP4 may be conferring a similar function across multiple begomovirus-vector systems. Among the three interacting genes identified in the largest module for CuLCrV was the gene Bta02748 (AP-3 complex subunit mu-1). Adaptin proteins (AP) are membrane-bound heterotetrameric complexes localized in cellular buds and vesicles, and their main role is intracellular trafficking in the trans-Golgi network [92]. The interaction of AP with the human respiratory syncytial virus matrix protein was essential for the latter's trafficking in the host cells [93]. The regulation of the AP-3 complex subunit mu-1 in the whiteflies that acquired a New World begomovirus could play a role in virus tropism, although this requires functional validation.
KEGG analysis identified many genes implicated in human virus infection that were differentially regulated in viruliferous B. tabaci MEAM1 (21 overexpressed and 9 underexpressed genes) and MED (5 overexpressed and 5 underexpressed genes) adults compared with non-viruliferous counterparts. For example, in B. tabaci MEAM1 adults that acquired SiGMV, Bta14253-protein kinase C was overexpressed compared with non-viruliferous adults. This gene was implicated in the infection of human immunodeficiency virus 1 (HIV), hepatitis B virus, coronavirus, influenza A virus, and human cytomegalovirus by KEGG analysis. Protein kinase C agonists are reported to induce latent HIV expression from viral reservoirs and protect primary CD4 + T cells from HIV infection through the down-modulation of HIV coreceptor expression [94]. Another class of genes, zinc finger proteins (Bta01535 and Bta03437) was underexpressed in B. tabaci MEAM1 adults that acquired TYLCV compared with their non-viruliferous counterparts. In addition, these genes (zinc finger proteins) were also modulated in both directions (up and down) in B. tabaci MED adults that acquired SiGMV compared with their non-viruliferous counterparts. Zinc finger proteins were also implicated in herpes simplex virus 1 infection (HSV-1). The upregulation of zinc finger proteins was reported to play a role in HSV-1 replication and binding to promoter proteins [95]. Five heat shock proteins (Bta00008, Bta02903, Bta03000, Bta06076, and Bta14532) were all, except in one instance, underexpressed in viruliferous whiteflies compared with non-viruliferous counterparts. The increased expression of heat shock proteins was reported after measles virus infection, resulting in increased cytopathic effects [96]. The role of the above-selected genes in both B. tabaci MEAM1 and MED adults following the acquisition of begomoviruses is unknown and warrants further investigation.
For the successful circulative translocation of begomoviruses in whiteflies, the viral coat protein must interact with putative receptors at the midgut, in the hemolymph, and at the primary salivary glands [22]. A few putative receptors such as the GroEL chaperone protein, heat shock proteins, midgut proteins, peptidyl-prolyl isomerase protein, and peptidoglycan recognition protein gene were identified in previous studies [23][24][25][26][27]. In addition to those, this study identified other possible putative receptors in the cytokinecytokine receptor interaction and cell adhesion molecule pathways that may play a role in the circulative movement of begomoviruses in their vectors, and they include Bta04818type I serine/threonine kinase receptor and Bta07388-syndecan. Serine/threonine kinase, an anchored protein-associated kinase in mammalian cells, was reported as an important cellular component in regulating the entry or clathrin-mediated endocytosis of rabies virus (RABV) [97]. Syndecan is a cell surface heparan sulphate proteoglycan (HSPGs), which plays several roles including virus (hepatitis E virus and herpes simplex virus) attachment and entry [98]. Hepatitis E virus, human papillomavirus, and HIV are some of the viruses reported to bind to cell surface HSPGs for their initial attachment to host cells [99][100][101]. The upregulation of the type I serine/threonine kinase receptor and syndecan genes in B. tabaci MEAM1 adults that acquired SiGMV could be associated with their binding to viral DNA, thereby enabling circulative translocation. Aminopeptidase N is another gene that was differentially expressed (underexpressed) in both B. tabaci MEAM1 and MED adults that acquired TYLCV. This gene, a plant virus receptor, was overexpressed in aphids and thrips that acquired pea enation mosaic virus and tomato spotted wilt virus (TSWV), respectively [78,102]. Whether the aminopeptidase N gene plays any role in begomovirus (TYLCV) reception in either B. tabaci MEAM1 or MED adults is unknown at this juncture.
Whiteflies have an innate immune system that has been documented against pathogens [103,104]. The RNA interference (RNAi) pathway is the major mechanism insects use against virus infection [43]. In addition, other innate antimicrobial pathways such as Imd, Toll, JAK-STAT, phagocytosis, apoptosis, and proteolysis were reported to play crucial roles in insects' (Diptera and Lepidoptera) antiviral responses [105][106][107]. A majority of the genes associated with innate immune pathways in the viruliferous B. tabaci MEAM1 adults were overexpressed, as opposed to the viruliferous B. tabaci MED adults, where they were underexpressed (Tables 4 and 5). For instance, several cathepsin B genes were identified in both B. tabaci MEAM1 and MED adults. These genes were associated with apoptosis and immune system pathways such as the NOD-like receptor signaling pathway and antigen processing and presentation pathway. Further, cathepsins are implicated to play a role in virus transmission [108,109]. All of the nine genes identified in B. tabaci MED were underexpressed, while four genes were overexpressed and three were underexpressed in B. tabaci MEAM1 adults (Table 4). Similar to our findings, several studies have identified the differential expression (both up-and downregulation) of genes such as cathepsins associated with innate immune pathways in whiteflies following begomovirus or crinivirus acquisition [44][45][46][47]49,51]. The contrasting differential expression of genes under the innate immune pathways such as cathepsins in B. tabaci MEAM1 and MED identified in this study possibly highlights differences in the immune responses to virus acquisition or other vector-virus interactions between the two B. tabaci cryptic species. Ras-like GTP-binding proteins were other immune response genes identified only in B. tabaci MEAM1 adults that acquired TYLCV or SIGMV. Like B. tabaci MEAM1 adults that acquired TYLCV or SiGMV, ras-like GTP-binding proteins were overexpressed in thrips that were exposed to TSWV [78]. The upregulation of these genes in viruliferous B. tabaci MEAM1 adults could be another mechanism evolved by whiteflies to counter begomoviruses.
The exposure to begomoviruses, especially TYLCV and CuLCrV, affected the expression of several fitness-related genes in B. tabaci MEAM1 and MED. This study identified the up-and downregulation of genes associated with egg production, spermatogenesis, and longevity in whiteflies that acquired begomoviruses. For example, Bta06648gametocyte-specific factor 1 (GTSF1) gene was overexpressed in B. tabaci MEAM1 adults that acquired TYLCV. In Drosophila, this gene was essential for P-element-induced wimpy testis-interacting RNA-mediated transcriptional repression, the histone mediated repression of transposons, and their neighboring genes in the ovary [110]. Similarly, in mice, this gene was essential for spermatogenesis and transposon suppression in mouse testes [111]. The upregulation of the GTSF1 gene in B. tabaci MEAM1 adults that acquired TYLCV may increase the development of sperm cells in their male reproductive organs, thereby influencing reproduction. Two egg production associated genes-vitellogenin (Bta07852 and Bta11903) were underexpressed in B. tabaci MEAM1 and MED adults that acquired TYLCV, whereas one vitellogenin gene (Bta11903) was overexpressed in B. tabaci MEAM1 adults that acquired CuLCrV. Further, most genes associated with aging pathways in viruliferous B. tabaci MEAM1 and MED adults were underexpressed. Positive, neutral, and negative effects of begomovirus infection on whitefly fecundity and longevity have been reported [42,[112][113][114]. Evidence at the transcriptome level for the underexpression of the egg production gene in B. tabaci MEAM1 and MED adults that acquired TYLCV, as indicated in this study, is consistent with the Pan et al. [114] study. That study reported that B. tabaci MEAM1 and MED that acquired TYLCV produced fewer eggs on cotton than their non-viruliferous counterparts [114]. In contrast, the overexpression of oviposition-related genes in B. tabaci MEAM1 that acquired CuLCrV may indicate higher fecundity, but this assumption was not supported by Gautam et al. [42]. There were no significant differences in the fecundity of B. tabaci MEAM1 females that acquired CuLCrV and their non-viruliferous counterparts [42]. Gautam et al. [36] demonstrated that the fitness effects of both B. tabaci MEAM1 and MED did not vary within each host species, and in this study, host effects were eliminated by gut clearing for 72 h. Consequently, this study predominantly provides evidence at the transcriptome level for some of the begomoviruses-induced (predominantly direct) macro-level fitness effects on their vectors-B. tabaci MEAM1 and MED.

Conclusions
The acquisition of begomoviruses by two invasive whiteflies induced variable transcriptional responses. Increased transcriptional changes were observed in B. tabaci MEAM1 adults compared to B. tabaci MED adults. Although the three plant viruses belong to the same genus, Begomovirus, the gene expression patterns varied between them. The gene co-expression network analysis for B. tabaci MEAM1 and MED that acquired begomoviruses revealed some candidate genes such as the ubiquitin-conjugating enzyme, which could be involved in the endocytic trafficking of virus particles. In addition, other possible putative receptors such as type I serine/threonine kinase receptor and syndecan were identified, and these could be involved in the circulative movement of begomoviruses in their vectors. Several genes implicated in the infection of human viruses were differentially regulated in viruliferous whiteflies. Furthermore, under the innate immune pathways, the contrasting differential expression of genes such as cathepsins occurred between B. tabaci MEAM1 and MED adults that acquired the begomoviruses. The differences in the immune responses of B. tabaci MEAM1 and MED adults upon virus acquisition possibly highlight their unique defense mechanisms against begomoviruses. Multiple genes were identified that may play a role in the New World begomovirus transmission. These could be important targets to investigate in the future. This study also provides evidence at the transcriptome level for some of the observed macro-level fitness effects of begomoviruses on their whitefly vectors.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/cells11132060/s1, Figure S1: TYLCV weighted gene co-expression network analysis. Figure S2: SiGMV weighted gene co-expression network analysis. Figure S3: CuLCrV weighted gene co-expression network analysis. Table S1: Endpoint PCR and qPCR primers used in this study. Table S2: Pearson's correlation coefficients for multiple biological replicates. Table S3: RT-qPCR concordant validation for selected differentially expressed genes from RNA-Seq data. Table S4: KEGG annotations and pathway analysis for differential expressed genes in viruliferous versus non-viruliferous whiteflies.  Data Availability Statement: The raw read sequences were submitted to NCBI with the BioProject accession number PRJNA796977.