An isofemale population of T. palmi maintained at Advanced Centre for Plant Virology, Indian Agricultural Research Institute (IARI), New Delhi, was used in the study. The identification of T. palmi was based on the morphological keys and mitochondrial cytochrome oxidase subunit I (mtCOI) sequences amplified in PCR using primer pairs, LCO 1490 and HCO 2198 (Folmer et al., 1994; Jangra et al., 2020). The population was established from a single adult female of T. palmi on eggplant (var. Navkiran, Mahyco). The population has been maintained at the thrips rearing facility under controlled environmental conditions at 28 ± 1°C temperature, 60 ± 10% relative humidity, and 16 h of light–8 h of the dark since 2018. Fresh virus-free healthy eggplants were supplied as and when required. The isofemale population has been used in replicate throughout the experiment.

The initial inoculum of GBNV was collected from the pure culture maintained at Advanced Centre for Plant Virology, IARI. Healthy cowpea plants (var. Pusa Komal) were raised from seeds in a plant growth chamber (A2000, Conviron, Canada). The healthy cowpea plants at a two-leaf stage were sap-inoculated with GBNV as described by Ghosh et al. (2021). Briefly, the sap was extracted by grinding the inoculum in a pre-sterilized mortar on ice. Ice-cold sodium phosphate buffer (0.01 M, pH 7.0) was added at tissue to buffer ratio of 1:6 wt/vol. β-mercaptoethanol (0.2%) was added while grinding the tissue. A pinch of Celite was dusted on the cowpea plants. The sap was applied to the leaves by gently rubbing with gloved hands. The inoculated plants were washed with sterile distilled water after a while. The inoculated plants were maintained at 25 ± 1°C, 60 ± 10% relative humidity, and 16 h of light–8 h of dark under insect-proof conditions. Nutrient solution containing macro, meso, and trace elements at appropriate concentrations was provided to plants and regularly monitored for symptom appearance. All the inoculated plants were tested in reverse-transcriptase PCR (RT-PCR) using GBNV-specific primers, AG109F and AG110R (Supplementary Table 1) as described below.

The total RNA was isolated from cowpea plants using RNeasy Plant Mini Kit (Qiagen, Germantown, MD, United States). Complementary DNA (cDNA) was synthesized with random primers using FIREScript RT cDNA synthesis kit (Solis Biodyne, Estonia). The reaction mixture contained 1X RT reaction buffer, 1.0 μg template RNA, 5 μM random primers, 2 mM dNTP mix, 10-unit FIREScript RT, and 20-unit RiboGrip RNase inhibitor. The reverse transcription was carried out in a T100 thermocycler (Bio-Rad, United States) by primer annealing at 25°C for 10 min, reverse transcription at 50°C for 30 min and followed by enzyme inactivation at 85°C for 5 min. PCR was carried out in 25 μl reactions containing 2 μl cDNA, 2.5 μl 10X PCR buffer (Thermo Fischer Scientific), 0.4 μM each of forward and reverse primers, 260 μM dNTP mix (Thermo Fischer Scientific), and 2 U DreamTaq polymerase (Thermo Fisher Scientific). PCR was performed in a T100 Thermal Cycler with 95°C for 5 min followed by 30 cycles of 94°C for 30 s, 59°C for 45 s, and 72°C for 1 min followed by a final extension at 72°C for 7 min. RT-PCR products were resolved on 1% agarose gel stained with GoodView (BR Biochem, India) and observed in a gel documentation system (MasteroGen Inc., Taiwan). The purified product was sequenced bidirectional. The sequences were processed by BioEdit, and BLASTn was performed to check the species homology. The consensus sequence was submitted to GenBank.

The freshly emerged first instar larvae (L1, less than 24 h old) of T. palmi were collected using a Camel hairbrush from the isofemale population on eggplants and used for GBNV acquisition in three replicates to generate the viruliferous population. The virus acquisition setup was established as described by Ghosh et al. (2021). About 50 T. palmi L1s in each replicate were placed on a GBNV-infected leaf and allowed to feed for 24 h at 28 ± 1°C and 60 ± 10% relative humidity. The setup was continuously monitored to confirm the feeding of larvae on the leaf. After the acquisition access period, the larvae were transferred on a detached healthy cowpea leaf. The petiole of the detached leaf was inserted in a slant of 0.8% dextrose agar within an insect breeding dish (10 cm diameter, 4 cm height, SPL Life Sciences, Korea). The GBNV-exposed T. palmi were reared up to the adult stage on cowpea leaves at 28 ± 1°C, 60 ± 10% relative humidity, and 16 h of light–8 h of dark. The breeding dishes were regularly monitored and supplied with fresh leaves as and when required. To generate nonviruliferous T. palmi adults, the L1s were released on the healthy cowpea leaves placed on the agar slant within the insect breeding dishes in three replicates as described above. Adults were collected immediately after their emergence and used in the study. Three sets of viruliferous (designated as TpTrI1, TpTrI2, and TpTrI3) and nonviruliferous (TpTrH1, TpTrH2, and TpTrH3) T. palmi adult populations were developed. Each set was divided into two parts. One part was used for transcriptome sequencing and another part was preserved at −80°C for gene expression analysis in reverse transcriptase-quantitative real-time PCR (RT-qPCR).

Infection of GBNV in T. palmi was confirmed by randomly collecting five individual adults from each of the viruliferous and nonviruliferous insect populations and testing in RT-PCR. Total RNA from adults of T. palmi was isolated using NucleoSpin RNA XS (Macherey-Nagel, Germany), and RT-PCR was performed as described above.

Three sets of viruliferous (TpTrI1, TpTrI2, and TpTrI3) and nonviruliferous (TpTrH1, TpTrH2, and TpTrH3) adults of T. palmi were subjected to total RNA isolation using TRIzol (Invitrogen, United States). Each set comprised 10 adult individuals that included both females and males (1:1).

RNA quality was checked using RNA 6000 Nano Kit (Agilent Technologies, United States) on 2100 Bioanalyzer (Agilent Technologies) with a minimum RNA Integrity Number (RIN) value of 7. RNA concentrations were determined with a NanoDrop ND-8000 spectrophotometer (Thermo Fischer Scientific). RNA-Seq libraries for all samples were prepared using NEBNext UltraII RNA library preparation kit for Illumina (New England Biolabs, United States), and sequencing was done in a single HiSEQ 4,000 (Illumina, Inc. United States) lane using 150 bp paired-end chemistry. The library preparation and sequencing were done by commercial service providers (NxGen Bio Life Sciences, India). Briefly, total RNA was used to purify poly(A) messenger RNA (mRNA) using oligo-dT beads. Magnetic beads were used for two rounds of purification. During the second elution of the poly-A RNA, the RNA was also fragmented into 200–500 bp pieces in the presence of divalent cations at 94°C for 5 min using an ultrasonicator. The cleaved RNA fragments were copied into first-strand cDNA using SuperScript-II Reverse Transcriptase (Thermo Fischer Scientific) and random primers. After second-strand cDNA synthesis, fragments were end-repaired and A-tailed, and indexed adapters were ligated. The products were purified and enriched with PCR to create the final cDNA library. The tagged cDNA libraries were pooled in equal ratios and used for 2 × 150 bp paired-end sequencing on a single lane of the Illumina HiSeq4000. Illumina clusters were generated and were loaded onto Illumina Flow Cell on Illumina HiSeq 4000 instrument and sequencing was carried out. After sequencing, the samples were demultiplexed and the indexed adapter sequences were trimmed using the CASAVA v1.8.2 software (Illumina, Inc.).

The obtained low-quality raw reads were preprocessed to remove the adaptor contamination using Trim Galore v0.4.1¹ with the quality score <20%, and trimming ambiguous “N” nucleotides with a ratio of “N” > 5%. A reference genome index was established using BWAv0.7.5 (Li and Durbin, 2009). The processed clean paired-end reads were used to map the reference genome of T. palmi (Guo et al., 2020) by BWA-MEM (Li, 2013). Read numbers mapped to every gene were counted using Samtools v0.1.19 software. Differential expression between the viruliferous and nonviruliferous T. palmi populations was analyzed using the DESeq R package.² An absolute value of log₂ fold change >2 and value of p <  0.05 were used as the thresholds to judge the significance of gene expression differences.

Gene annotations and functional enrichment analysis including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) biological pathways were performed to identify the DEGs that were significantly enriched in GO terms or biological pathways post-GBNV infection. GO enrichment analysis for the identified GO terms was performed using the Fisher’s exact test available in the Blast2Go program (Conesa and Götz, 2008). Gene annotations against the Uniprot GO database³ were performed by aligning DEGs to the NR database using Blast2Go v5.2.5 software. GO terms were also assigned using the QuickGO online search tool.⁴ KEGG pathway enrichment analysis of DEGs was performed using the KEGG database resource⁵ to identify the pathways that were differentially regulated between viruliferous and nonviruliferous T. palmi with corrected value of p < 0.05. KEGG pathway enrichment was performed by calculating the enrichment factor for all the identified KEGG pathways. The top 20 highly enriched pathways were plotted based on the rich factor.

The top 20 genes each from up- and downregulated DEGs were selected for the construction of gene regulatory networks (GRNs). These DEGs were filtered based on the log fold change values. The network was analyzed and visualized using Cytoscape v3.7.2, where gene correlation was computed using the Pearson’s correlation coefficient of the normalized expression values (Shannon et al., 2003). Network Analyzer plug-in was used to estimate the network centrality and topology.

To validate the differential gene expression data by RNA-Seq analysis, a few highly up- and downregulated DEGs of T. palmi were selected to assess the gene expression in RT-qPCR. A part of the preserved viruliferous and nonviruliferous T. palmi populations was used in RT-qPCR analysis. Initially, three sets of primers were designed for each of the target genes, and finally, one pair of primers for each gene was selected after validation. β-tubulin was considered as an endogenous control. All the primer sets were validated and optimized in a gradient PCR. Following thermal cycling conditions were followed for validation and optimization of target genes and endogenous control primers: initial denaturation at 95°C for 5 min, then 35 cycles of 95°C for 40 s, 50°C–60°C (depending upon the Tm of primer pairs) for 40 s, and 72°C for 40 s followed by a final extension at 72° C for 10 min. The primer pairs that produced a single specific sharp band at the same PCR conditions for β-tubulin primers were further selected for RT-qPCR experiment. The list of primers that were validated and used in RT-qPCR assay is included in Supplementary Table 1.

The relative expressions of target genes were estimated by RT-qPCR following 2^−ΔΔCT method (Livak and Schmittgen, 2001). Firstly, total RNA was isolated from viruliferous (TpTrI1, TpTrI2, and TpTrI3) and nonviruliferous (TpTrH1, TpTrH2, and TpTrH3) thrips samples stored at −80°C using NucleoSpin RNA XS. cDNA was synthesized with oligo (dT) primer using FIREScript RT cDNA synthesis kit. The reaction mixture contained 1X RT reaction buffer, 1.0 μg template RNA, 5 μM oligo dT primer, 2 mM dNTP mix, 10-unit FIREScript RT, and 20-unit RiboGrip RNase inhibitor. The reverse transcription was carried out in a T100 thermocycler at 50°C for 30 min followed by enzyme inactivation at 85°C for 5 min. RT-qPCR was performed in an Insta Q 48M (Himedia, India) using the DyNAmo ColorFlash SYBR Green qPCR Kit (Thermo Fisher Scientific). A 20 μl of RT-qPCR reaction mixture consisted of 10 μl of 1X DyNAmo ColorFlash SYBR Green Master Mix, 0.25 μM ROX passive reference dye, 0.2 μM each forward and reverse primer, and 2 μl template cDNA. Thermal cycling was performed as initial denaturation at 95°C for 5 min, 35 cycles of 95°C for 40 s, annealing at 53–59°C (depending upon the Tm of primer pairs) for 40 s, and 72°C for 40 s. Since SYBR Green I dye binds nonspecifically to any double-stranded DNA, a dissociation or melting stage was added after every reaction to determine the specificity of the amplicons based on the melting curves. The RT-qPCR was performed with three biological and two technical replicates. Log₂ fold change value was calculated and relative expression of mRNA was determined by normalized log₂ 2^−ΔΔCT (Livak and Schmittgen, 2001) values of viruliferous (TpTrI) in comparison to nonviruliferous (TpTrH) T. palmi. The ΔC_T was calculated as C_T of targeted gene—C_T of endogenous control. ΔΔCT value was calculated as average ΔC_T of nonviruliferous replicates—ΔCT of each replicate. The average 2^-(ΔΔCT) was calculated for both the viruliferous and nonviruliferous populations and transformed to log₂ 2^-(ΔΔCT). The log₂-fold change of the target gene expression in viruliferous T. palmi population was estimated by normalizing the fold change values of the nonviruliferous population. The relative expressions of target genes post-GBNV exposure in RT-qPCR assay were compared with log₂-fold changes obtained through RNA-Seq analysis, and Pearson’s correlation coefficient was calculated using CORREL function in MS Excel.

The isofemale T. palmi population generated from a single adult female was confirmed on the basis of morphological characters. The adults of T. palmi were yellow in color measuring about 0.8 mm–1 mm in length. Three brick red ocelli were visible in a triangular fashion between the compound eyes. A pair of interocellar setae originated outside the ocellar triangle. The head was quadrangular in shape with seven segmented antennae. Female adults were comparatively larger than males. Female adults had sharp ovipositor at the apex of the abdomen while males were with rounded apex. Further, the nucleotide sequence of mtCOI confirmed the identity as T. palmi. A 657 bp fragment of mtCOI was amplified using primer pairs LCO 1490 and HCO 2198 (Folmer et al., 1994). The nucleotide sequence of the amplified product showed 100% homology to other T. palmi isolates (MN594549, NC_039437, and KX622379) available in GenBank. The sequence can be retrieved from the GenBank with the Accession number MW020349.

The initial inoculum of GBNV was collected from the pure culture (tomato isolate) maintained at Advanced Centre for Plant Virology, IARI. Healthy cowpea plants (Vigna unguiculata var. Pusa Komal) were sap-inoculated with GBNV. Characteristic symptoms of GBNV were recorded on sap-inoculated cowpea plants 10–14 days post-inoculation (dpi). Local lesions appeared on inoculated leaves and yellow necrotic spots with multiple rings-like structures were visible on systemic leaves. RT-PCR with GBNV-specific primers produced an amplicon of 1.8 kb for all symptomatic samples that confirmed the GBNV infection in inoculated plants. No amplification was observed in RT-PCR for healthy plant samples. Further, the nucleotide sequence of the amplified product showed more than 97% identity to GBNV. The sequence can be retrieved from the GenBank with Accession No. MN566913.

The isofemale T. palmi adults in three replicates (TpTrI1, TpTrI2, and TpTrI3) exposed to GBNV-infected cowpea plants during the L1 stage were tested in RT-PCR. A 1.8 kb product was visualized on agarose gel that confirmed the infection of GBNV in T. palmi adults. T. palmi populations (TpTrH1, TpTrH2, and TpTrH3) that were exposed to healthy cowpea plants did not produce any GBNV-specific amplification in RT-PCR.

T. palmi transcriptome data generated a total of 156,127,441 raw reads of 150 bp from three replicated cDNA libraries of each viruliferous and nonviruliferous T. palmi population. Around 46.84 Gb of the sequence was generated. The average reads per library were about 26,021,240 (Supplementary Table 2). After preprocessing, a total of 156,074,069 clean reads of 20–150 bp were filtered from all six libraries. The preprocessed reads ranged between 99.95 and 99.98%. After removing the low-quality reads from all the sets, the high-quality reads were used for further analyses. These were mapped with the T. palmi genome (Guo et al., 2020) yielding an average mapping percent of 97.31%.

The differential expression analysis identified a total of 28,165 transcripts in the viruliferous vs. nonviruliferous T. palmi populations. Out of these, a total of 2,363 transcripts were significantly expressed in viruliferous T. palmi at stringent filtering parameters (log₂ fold change > ± 2 and adjusted value of p < 0.05). A total of 1,383 differentially expressed transcripts were upregulated, while 980 transcripts were downregulated in viruliferous (pooled TpTrI) as compared to nonviruliferous (pooled TpTrH) populations (Figure 1). Table 1 represents the top 20 up- and downregulated genes of T. palmi under the stringent filtering criteria in response to GBNV-infection. The putative genes of T. palmi such as cytochrome P450 (CYP), gamma-aminobutyric acid (GABA) receptor, heat shock protein (hsp70), nephrin, sperm-associated antigen 6 (spag6)-like, and ubiquitin-like containing PHD and RING finger domains 1 (UHRF-1)-binding protein 1, etc. were found highly upregulated after exposure to GBNV. Significant downregulations of ABC transporter, arrestin domain-containing 3 (ARRDC3), elongation of very-long-chain fatty acid (ELOVL), dynamin-1-like, serpin H1-like, tyrosine–protein kinase receptor Tie-1 (TIE1), larval cuticle protein A2B-like, pupal cuticle protein C1B-like, and endocuticle structural glycoprotein SgAbd-2-like, etc., were also observed.

The 1,383 (58.52%) up- and 980 (41.48%) downregulated DEGs were categorized under the three major classes, namely, cellular components, biological processes, and molecular functions. The highest number of GO terms were categorized into molecular functions (70%) followed by cellular components (24%), and biological processes (6%; Figure 2A). In the molecular functions category, there was an abundance of genes associated with catalytic activity followed by binding, molecular transduction activity, antioxidant activity, and transporter activity. Similarly, genes associated with the membrane, membrane part, cell, cell part, protein-containing complex, organelle, synapse, synapse part were differentially enriched under the cellular components category. In the biological processes, genes involved in the cellular process, response to stimuli, signaling, regulation of the biological process, localization, metabolic process, biological regulation, and signaling were enriched. Within each of the functional groups, the top enrichment terms for both up- and downregulated DEGs identified as per Fisher’s exact test value of p are listed in Figures 2B,C, respectively.

Kyoto Encyclopedia of Genes and Genomes pathway analysis of DEGs showed that the genes were involved in the functions like metabolic pathways, biosynthesis of secondary metabolites, cellular processes, endocytosis, organismal systems, and signaling pathways. Majorly affected metabolic pathways were glycerolipid, glycerophospholipid, purine, fructose and mannose, galactose, glutathione, and carbon metabolism. The differentially expressed signaling pathways were calcium, cyclic adenosine monophosphate (cAMP), glucagon, hypoxia-inducible factor 1 (HIF-1), mitogen-activated protein kinase (MAPK), neurotrophin signaling, phospholipase D, Ras-proximate-1 (Rap1) signaling, and Ras signaling pathways. The most affected cellular and organismal systems were dopaminergic synapses, endocytosis, glutamatergic synapse, lysine degradation, neuroactive ligand-receptor interaction, oocyte meiosis, and synaptic vesicle cycle. In KEGG pathway enrichment analysis, the highest number of DEGs was annotated in drug metabolism (map00982, 00983) under the class xenobiotics biodegradation and metabolism followed by glycerolipid metabolism (map00561), methane metabolism (map00680), glycerophospholipid metabolism (map00564), and fructose and mannose metabolism (map00051; Figure 3). The pathway enrichment of DEGs showed the DEGs were significantly enriched in ascorbate and aldarate metabolism (map00053), steroid hormone biosynthesis (map00140), and metabolism of xenobiotics by cytochrome P450 (map00980) pathways. The pathways identified in GO enrichment analysis were consistent with the findings of the KEGG pathway study.

A total of 40 nodes comprised the top up- and downregulated DEGs were considered in network analysis. The degree and betweenness centrality in the interaction network of the top DEGs were calculated. Interaction networks corresponding with the highly regulated DEGs demonstrated a total of 438 edges. The degree of certain nodes was markedly higher compared with the average degree in the networks (Figure 4). The upregulated DEGs with significantly higher degrees were calphotin-like isoform X2, proline-rich protein 36-like isoform X1, dipeptidyl-aminopeptidase B, sodium channel protein para isoform X8, nephrin-like, cAMP-specific 3′,5′-cyclic phosphodiesterase isoform X1, UHRF1-binding protein 1-like isoform X3, spag6-like, laccase-5, CYP 4C1-like, phosphoribosylformylglycinamidine synthase (PFAS), and hsp70 II-like. Among the downregulated nodes, a group of DEGs associated with cuticle development such as larval cuticle protein A2B-like, endocuticle structural glycoprotein SgAbd-2-like, and ABC transporter F family member 4-like were with higher degrees of interactions. Besides, proline-rich protein HaeIII subfamily 1-like, ARRDC3, serpin H1-like, TIE-1-like, myrosinase 1-like, alpha-tocopherol transfer protein (TTPA)-like, and dynamin-1-like protein showed significant interactions with higher degrees.

To validate the differential gene expression data, thirteen highly expressed genes of T. palmi were selected and mRNA expression levels were quantified in RT-qPCR in response to GBNV infection. A part of the samples used for RNA-Seq was preserved and used in RT-qPCR assay. Primer pairs designed for this study were initially optimized in a gradient PCR (Supplementary Table 1). Primer pairs of target genes that produced a single sharp amplicon at the same PCR conditions for endogenous control (β-tubulin) primers were selected for RT-qPCR assay. The primer pair specific to β-tubulin produced sharp bands at all annealing temperatures between 53 and 59°C. The annealing temperature of primer pairs for each target gene was standardized within the same temperature range. One best primer pair for each of the target genes was optimized for RT-qPCR assay. The annealing temperature, amplicon size, and melting temperature of each amplicon in RT-qPCR are included in Supplementary Table 1. The primer pairs for target and endogenous control genes produced single specific peaks without any secondary amplification in the RT-qPCR melting curve analysis that indicated the specificity of the reactions (Supplementary Figure 1).

The relative expression of target genes was estimated by a 2^−ΔΔCT method. The variation among the biological replicates was normalized with respect to the C_T value of endogenous control, β-tubulin. Among the selected highly expressed genes, the expression of UHRF1-binding protein 1 of T. palmi was upregulated by log₂ 5.5-fold in RT-qPCR in response to GBNV infection (Figure 5). Likewise, mRNA expression of nephrin was upregulated by log₂ 4.8-fold in viruliferous T. palmi adults in comparison to nonviruliferous adults. mRNA levels of spag6-like, hsp70, and GABA receptor were upregulated by log₂ 3.37, 0.98, and 0.80-fold, respectively, post-GBNV exposure. In viruliferous T. palmi, the expression of serpin was downregulated by log₂ 8.5-fold. Similarly, log₂ 4.3-fold downregulation of tyrosin kinase was recorded in RT-qPCR in response to GBNV infection. Expression of ABC transporter F family member 4-like gene was also downregulated by log₂ 3.57-fold, while it was log₂ 1.65-fold downregulations for ARRDC3 in viruliferous T. palmi. Some other genes such as elongation of fatty acid chain, dynamin-1-like, myrosinase 1-like showed log₂ 1.36, 0.66, and 2.57-fold, respectively, lower expression in viruliferous T. palmi as compared to nonviruliferous.

The fold change values estimated from RNA-Seq were plotted against fold change values obtained in RT-qPCR assays of the selected genes (Figure 5). The Pearson’s correlation coefficient of 0.81 indicated the acceptability of DEGs through RNA-Seq (Supplementary Table 3).