In silico prediction of candidate gene targets for the management of African cassava whitefly (Bemisia tabaci, SSA1-SG1), a key vector of viruses causing cassava brown streak disease

Sample collection and RNA isolation

The cassava B. tabaci population (SSA1-SG1-Ug) was derived from a single virgin female and male collected from Namulonge, Uganda. Namulonge is one of the areas in Uganda with superabundant B. tabaci populations and a high incidence of both cassava brown streak disease and cassava mosaic disease (Mugerwa et al., 2018). The B. tabaci SSA1-SG1-Ug population (member species name ‘SSA1-SG1’) was maintained at the Natural Resources Institute (UK), on aubergine (Solanum melongena cv. Black Beauty) under quarantine insectary conditions (60% relative humidity, 28 °C and 12 h light:12 h dark regime) with light intensity of 500 µmol m⁻² s⁻¹ PAR (Photosynthetic Active Radiation). Approximately 10,000 bacteriocytes and 300 guts were dissected in phosphate-buffered saline (PBS, 1X) from 500 and 300 female adult whiteflies respectively using fine pins and a stereomicroscope (Fig. 1). Gut samples were stored in RNAlater^TM stabilization solution (Invitrogen, Waltham, MA, USA) to preserve their RNA quality. A total of 50 female adult whiteflies were also collected in three replicates, flash-frozen in liquid nitrogen and stored at −80 °C prior to analysis of whole body gene expression in these insects.

RNA extraction from dissected guts and RNAseq library preparation

Total RNA was isolated from B. tabaci gut, bacteriocytes and whole body samples using TRIzoL^TM reagent (Invitrogen), with slight modifications made in the manufacturer’s instructions to ensure purity and quantity of RNA for downstream analysis for low tissue amounts. Specifically, we included a second clean-up with 100 µL chloroform, a second phase separation steps and addition of 1 µl linear acrylamide to the aqueous phase (before the isopropanol step) to improve the precipitation of RNA. Total RNA was further treated with ezDNase^TM enzyme (Invitrogen) according to the manufacturer’s instructions, to remove contaminating genomic DNA. RNA quantity and integrity were confirmed using a Qubit fluorometer (Life Technologies, Carlsbad, CA, SA).

Total RNA was depleted of cytoplasmic and mitochondrial rRNA using the Ribo-Zero Gold rRNA removal epidemiology kit (Illumina, San Diego, CA, USA) according to manufacturer’s instructions. For the whole body samples, 1 µg total RNA was used for rRNA depletion, while 367 and 230 ng total RNA were used for the bacteriocyte and gut samples. RNAseq libraries were prepared with the NEBNext Ultra II RNA Library Prep Kit (New England Biolabs Inc, Ipswich, MA, USA), according to manufacturer’s instructions. For the whole body samples, two technical replicates of 15 ng rRNA-depleted RNA were used as input while for the bacteriocyte and gut samples, all the rRNA-depleted RNA (not quantified) was used as input. The final libraries were quantified with a Qubit HS DNA assay (Thermo Fisher Scientific, Waltham, MA, USA) while the size distribution was determined using a Fragment Analyzer instrument (Agilent Technologies, Inc., Santa Clara, CA, USA). Equimolar amounts of each library were pooled for sequencing.

RNA sequencing and RNAseq analysis

The RNAseq libraries were sequenced and >40 million single-end 75 nt reads were generated on a NEXTseq500 sequencing platform at the Cornell Genomics Facility (Biotechnology Resource Center). Raw reads were filtered to remove low-quality reads and adaptor sequences with Cutadapt (Martin, 2011). The cleaned RNAseq reads were aligned against the SSA1-SG1-Ug genome (GCA_902825415.1) using HISAT2 (Kim, Langmead & Salzberg, 2015), which runs Bowtie2 in the background. The accepted hits were used, firstly to develop gene annotation files containing gene models based on RNAseq aligned reads for guts, bacteriocytes and whole body using cufflinks and cuffmerge version 2.2.1 (Trapnell et al., 2010). The accepted hits were then assembled and aligned reads counted using featureCounts (Liao, Smyth & Shi, 2014). Differential gene expression in cassava B. tabaci (SSA1-SG1-Ug) gut or bacteriocytes relative to the whole body was performed based on feature counts using edgeR (Robinson, McCarthy & Smyth, 2009) version 3.28.0 with three considerations. These include (i) false discovery rate controls based on the Benjamin and Hochberg method (Benjamini & Hochberg, 1995), (ii) trimmed mean of M-values normalization to account for library size variation between the samples and (iii) normalization for each gene, where the counts for each gene were mean centred and scaled to unit variance. The statistical power of this experimental design as calculated in RNAseqPower (Hart et al., 2013) was one, with a minimum sample size of eight whiteflies per sample to detect a two-fold difference in gene expression. Genes with enriched or depleted expression in the gut samples were identified by the general linear model (likelihood ratio test) analysis of log2 fold difference (logFC) data.

Identification of osmoregulation genes

Nucleotide sequences for the genes enriched in B. tabaci guts were extracted from the SSA1-SG1-Ug genome (GCA_902825415.1) using Samtools faidx tool (Li et al., 2009). Blastx search against NCBI protein reference non-redundant database (nr) with a cutoff E-value of 1e-3 (Altschul et al., 1997) mapping and annotation of these sequences were done using Blast2GO suite (Conesa et al., 2005). The suite was also used to generate and characterize gene ontology terms based on BLAST output. Gene ontology annotation was characterized into three groups, namely molecular function, biological process and cellular component. Selected sequences with metabolic functions were further translated to protein sequences using the ExPASy molecular biology server of the Swiss Institute of Bioinformatics (Gasteiger et al., 2003). The translated protein sequences with Aamy domain (PF00128, IPR006047) for glycosyl hydrolases (family 13) and aquaporin domain (PF00230, IPR000425), for osmoregulation genes were selected from the blast2go output. The selected sequences were further analyzed using SignalP 4.0 server (Petersen et al., 2011), NCBI conserved domain database (Marchler-Bauer et al., 2015), HMMER V.3 (Eddy, 2011) and InterProScan V5.52-86.0 (Zdobnov & Apweiler, 2001) to confirm the functional domain for sucrase and aquaporin. The selected protein sequences were also analyzed using NCBI Batch web CD-search tool to determine if the conserved domains were complete or incomplete at N-terminus or C-terminus.

Phylogenetic analysis

Phylogenetic analyses compared selected sequences to experimentally validated osmoregulation genes in other B. tabaci species and aphids (NCBI accession number: NP_001119607.1). Protein sequences were aligned by ClustalW (Larkin et al., 2007) implemented in MEGA X using default parameters. The best-fit amino acid substitution model to describe the substitution pattern was selected based on Bayesian information criterion (BIC) scores, computed using MEGA X (Kumar et al., 2018). The best amino acid substitution model with the lowest BIC score was selected and used in the phylogenetic reconstruction implemented in Bayesian Evolutionary Analysis Sampling Trees (BEAST version 1.10.2) (Suchard et al., 2018). BEAST was run with four gamma categories and strict molecular clock model at 1.0 clock rate. Yule model with uniform birth rate (Gernhard, 2008) and inverse gamma were used to provide prior on the tree. Markov Chain Monte Carlo (MCMC) was set with a chain length of 10,000,000. The generated tree was visualized using FigTree V1.4.3. The phylogenetic trees were used to study the relationship of selected genes with experimentally validated osmoregulation genes in other B. tabaci and aphid species (Price et al., 2007; Shakesby et al., 2009) to select the most suitable gene target.

Genome-scale reconstruction of Portiera and SSA1-SG1 metabolic models

Single compartment models for Portiera and B. tabaci SSA1-SG1 were reconstructed following the procedures of Thiele & Palsson (2010). A pseudo reaction for biomass generation was developed as the objective function for the model as described by Ankrah, Luan & Douglasa (2017). Stoichiometric coefficients in this reaction describe the relative abundance of essential building blocks, focusing on amino acid synthesis (amino acids, ATP and water), that enables Portiera and B. tabaci to grow and survive (Table S4). This reaction was selected as the objective in the calculation of optimal flux distribution for the Portiera and B. tabaci SSA1-SG1 model.

The two single-compartment models generated (for Portiera and B. tabaci) were combined to form a two-compartment model of B. tabaci SSA1-SG1 and Portiera, which was used for analysis and identification of essential symbiosis genes. The reactions in single models were renamed with a prefix SSA1_Por for Portiera model and SSA1_Bt for the B. tabaci SSA1-SG1 model. Metabolites in Portiera and B. tabaci models were localized in (p) and (c) compartments, respectively. Transport reactions were assigned to ensure that dead-end metabolites from each compartment are exchanged between the two compartments and import essential metabolites into the compartment. Exchange reactions were also included so that the metabolites or products produced by the two compartments are exported to the external environment (B. tabaci hemolymph).

Analysis and evaluation of models

Several parameters of the model were assessed. These included the optimal steady-state flux profile and essential reaction/genes in the Portiera and Portiera/B. tabaci model that would best predict the amino acid requirement for growth of Portiera and B. tabaci respectively. The model was analysed using flux balance analysis (Varma & Palsson, 1994) to determine the steady-state flux profile that optimizes the objective function and to simulate gene essentiality of the two-compartment model.

The amino acid requirement for growth of Portiera and B. tabaci (SSA1-SG1) were determined using flux balance analysis (FBA), with simulations performed under aerobic conditions with a maximum oxygen uptake rate of 20 mmol gDW⁻¹ h⁻¹. The objective function for the two-compartment model was developed (Table S4) to optimize the amino acid requirement for both Portiera and B. tabaci. FBA was also used to determine the effect of gene/reaction deletion on the metabolic model. Each reaction was perturbed, one at a time, by constraining both the upper and lower bound to zero, thereby restricting it from carrying flux. Gene or reaction knock out were implemented and growth simulated using SingleGeneDeletion or SingleReactionDelection function of COBRA Toolbox respectively. A reaction was classified as essential or critical if when removed from the model, the resultant value of biomass flux was zero (the network does not support growth) or non-essential if perturbation caused no effect on the network (network supports growth). Genes mediating reactions in essential amino acid biosynthesis were selected from the set of identified essential genes as potential gene targets for controlling cassava B. tabaci through RNAi. Minimization of metabolic adjustment (MOMA) (Segrè, Vitkup & Church, 2002) and regulatory on/off minimization (ROOM) (Shlomi, Berkman & Ruppin, 2005) were also used to confirm the essential genes selected.

Validation of selected genes using real-time quantitative PCR

The expression of the selected genes in the whitefly gut or bacteriocytes was validated using real-time quantitative PCR (qPCR). In brief, three technical samples of guts from 300 whiteflies, bacteriocytes from 500 whiteflies and whole body samples from 50 whiteflies were collected and preserved immediately at −80 °C. Total RNA was extracted as described above and its quantity and integrity were confirmed using a Qubit fluorometer (Life Technologies). Only samples with RIN values between 6–8 were considered for further analysis. A total RNA for whitefly gut, bacteriocyte, and whole body samples were normalized to 120 ng, to ensure an equal amount of starting RNA template. Normalized RNAs were then treated with 1 µl ezDNase^TM enzyme (Invitrogen) for 2 min at 37 °C to remove genomic DNA. The RNA sample was further incubated at 55 °C for 5 min with 10 mM DTT to inactivate the enzyme. First-strand cDNA synthesis was carried out using SuperScript™ III Reverse Transcriptase (Invitrogen) according to the manufacturer’s instructions. A three-step thermocycling protocol was performed on C1000 Thermocycler with CFX96 Real-time detection system (Bio-Rad, Hercules, CA, USA) with CFX Manager Version 3.1.1517.0823. The three steps included a polymerase activation and DNA denaturation at 95 °C for 3 min followed by 40 cycles of both denaturation at 95 °C and annealing/extension at 55.2 °C for 45 s. Both 60S ribosomal protein L13a (RPL13) and β-tubulin were used as the internal controls (Fig. S5, Table S5). No template controls (NTC) were also added to evaluate the quality of analysis. A total of 20 µl reaction mix containing 10 µl iQ™ SYBR^® Green Supermix (2x), 1 µl forward and reverse primer (Table S1), 1 µl cDNA and 7 µl water was used. Melting curve analysis at the end of the reaction was done by increasing temperature from 55–95 °C in increments of 0.5 °C every 5 s to assess the dissociation characteristics of double-stranded DNA during heating (Table S5).

Expression patterns of the gut and bacteriocytes of Bemisia tabaci SSA1-SG1

RNAseq analysis was conducted to identify candidate osmoregulation and symbiosis genes with enriched expression in the whitefly gut and bacteriocytes respectively. A total of 49,426,054, 47,305,213 and 47,646,265 raw reads was generated for the whitefly gut, bacteriocyte and whole body samples respectively. These samples had closely similar percentage mapping against the B. tabaci SSA1-SG1-Ug genome of 78.1, 83.84 and 78.0% for the gut, bacteriocyte and whole body samples, respectively. RNAseq analysis revealed a total of 10,027 transcripts (64% of what has been reported in B. tabaci species Middle East Asia Minor 1 (MEAM1)) (Chen et al., 2016). Of these transcripts, 1,316 and 1,909 had enriched expression in the whitefly gut and bacteriocytes relative to the whole body. This gene set was then assessed for candidate genes with osmoregulatory function, specifically genes encoding sucrase enzymes (SUC genes) and water-specific aquaporin channels (AQP genes) and symbiosis genes (amino acid synthesis, transport and horizontally transferred genes).

Genes encoding sugar metabolism with enriched expression in the SSA1-SG1 gut

A total of 24 transcripts encoding sugar processing enzymes with an AmyAc_Maltase (cd11328) catalytic domain and belonging to glycosyl hydrolase family 13 (EC 3.2.1.20) had significantly enriched expression in whitefly gut relative to the whole body (Table 1). Of these, seven transcripts had the highest expression with log₂CPM ranging from 10.18–11.86, equivalent to a 17–30-fold enrichment in the B. tabaci SSA1-SG1 gut relative to whole body. Further analysis revealed that 14 of the 24 α-glucosidase genes with enriched expression in the gut had both a catalytic site residue, required for enzymatic activity, and a predicted signal peptide, required for the enzyme to function as an extracellular enzyme in the gut lumen (Table 1). Phylogenetic analysis of the 14 selected protein sequences identified three clusters (Fig. 2), with some protein sequences aligning most closely with the experimentally validated pea aphid sucrase protein (SUC1). The protein of SSA1-SG1-Ug that aligned most closely with the pea aphid SUC1 is ENSSSA1UGG000718 (SSA-ECA-Ssa05164), which has a predicted signal peptide and intact active site. The transcript was enriched 26-fold in the gut relative to the whole body. We predict that, as for the pea aphid SUC1, this protein has a key osmoregulatory function. We name the protein SSA1_SUC1 and its cognate gene SSA1_SUC1. Another protein that aligned closely to the experimentally validated pea aphid sucrase protein and its corresponding gene and had a high gene expression in the SSA1-SG1 gut was ENSSSA1UGG009974 (Ssa05154). We name this protein SSA1_SUC2 and its cognate gene SSA1_SUC2.

Water cycling within cassava Bemisia tabaci SSA1-SG1 gut

Seven aquaporin genes were identified in the transcriptome of B. tabaci SSA1-SG1-Ug (Table 2). The transcript abundance of three genes (ENSSSA1UGG010352 (SSA-ECA-Ssa03484), ENSSSA1UGG011821 (SSA-ECA-Ssa02238) and ENSSSA1UGG001892 (SSA-ECA-Ssa08000) were significantly enriched in the gut with log₂CPM of 8.49, 5.17 and 3.31 and fold difference of 32.7, 2.2 and 2.1-fold respectively. Phylogenetic analysis (Fig. 3) included the sequences of experimentally-validated AQP proteins in B. tabaci studied by Van Ekert et al. (2016), revealing that ENSSSA1UGG010352 has homology to the water channel BtDRIP. Inspection of the predicted protein sequence of ENSSSA1UGG010352 confirmed that it had the motifs of a DRIP water channel: the two NPA boxes, the ar/R constriction region with phenylalanine at position 71, histidine at position 198, alanine at position 207 and arginine at position 213, and a mercury-sensitive cysteine residue located two residues after the first NPA box.

Genes with enriched expression in the cassava Bemisia tabaci bacteriocytes

Interaction between B. tabaci and Portiera for essential amino acid provisioning takes place within the bacteriocyte. RNASeq analysis to identify genes with enriched expression in SSA1-SG1 bacteriocytes relative to whole body revealed a total of 1,909 transcripts with enriched expression log fold change (log FC) >1 in the bacteriocyte. Of these, 20 transcripts/genes encoding for proteins involved in amino acid biosynthesis and transport were selected (Table 3). Another category of genes with enriched expression in the bacteriocyte was the horizontally transferred genes (HTGs). HTGs contributing to lysine, arginine and proline biosynthesis within bacteriocytes had relatively high gene expression compared to other HGTs. Proteins encoded by these genes included; 2OG-Fe, oxygenase, allophanate hydrolase, 4-hydroxy-tetrahydrodipicolinate reductase (dapB), diaminopimelate decarboxylase (LysA), arginosuccinate lyase (ArgH) and arginosuccinate synthase (ArgG) (Table 4). The role and essentiality of these genes in the genome-scale metabolic network was analysed using constraint-based metabolic modelling techniques.

Metabolic models for Portiera and its host, SSA1-SG1

Single compartment models of both the whitefly and the primary endosymbiont were reconstructed. The reconstructed single compartment metabolic model for Portiera (SSA1-SG1) comprised of 76 intercellular reactions, 150 metabolites and a total of 90 genes hence supporting a metabolic network, iKT90 while that of B. tabaci SSA1-SG1 (iKT330) comprised 251 metabolites and 233 intracellular reactions, mediated by 330 metabolic genes. The shared metabolic interaction between Portiera SSA1-SG1 and B. tabaci SSA1-SG1 involved a total of 420 metabolic genes mediating 310 intercellular reactions and an optimal flux balance solution of 0.39 (h⁻¹), hence a two-compartment model iKT420.

Shared metabolic interaction between Portiera and Bemisia tabaci SSA1-SG1

Genome-wide reconstruction of a metabolic pathway for Portiera from SSA1-SG1 revealed that Portiera had genes that mediated the terminal reaction of three essential amino acids. These were (i) threonine, (ii) methionine and (iii) tryptophan, with reactions producing 0.101, 0.046 and 0.012 mmol gDW⁻¹ h⁻¹ flux, respectively. Furthermore, Portiera also contributes to the synthesis of metabolites for intermediate reactions of other essential amino acids. The bulk of reactions of the biosynthesis pathways of seven essential amino acids (arginine, histidine, lysine, phenylalanine, isoleucine, valine and leucine) are mediated by genes in Portiera. Metabolites produced in these reactions are exported across a symbiosome membrane into the bacteriocyte for amino acid biosynthesis. Metabolites (intermediate precursors) exported across the symbiosome membrane include, (i) 0.371 mmol gDW⁻¹ h⁻¹ N(Omega)-(L-Arginino) succinate for arginine biosynthesis, (ii) 0.049 mmol gDW⁻¹ h⁻¹ L-Histidinol phosphate for histidine biosynthesis, (iii) 0.1599 mmol gDW⁻¹ h⁻¹ LL-2,6-Diaminoheptanedioate for lysine biosynthesis, (iv) 0.2097 mmol gDW⁻¹ h⁻¹ phenylpyruvate for phenylalanine biosynthesis, (v) 0.1133 mmol gDW⁻¹ h⁻¹ (5)-3-methyl-2-oxopentanoate for isoleucine biosynthesis, (vi) 0.3107 mmol gDW⁻¹ h⁻¹ 3-methyl-2-oxobutanoate for valine biosynthesis and (vii) 0.1653 mmol gDW⁻¹ h⁻¹ 3-carboxy-4-methyl-2-oxopentanoate for leucine biosynthesis.

In addition, the single-compartment model of Portiera revealed that Portiera lacks genes that encode for metabolites required in the biosynthesis of non-essential amino acids (alanine, asparagine, aspartic acid, cysteine, glutamic acid, glutamine, glycine, proline, serine and tyrosine). It also requires eight host-derived metabolites for the biosynthesis of essential amino acids. These include, (i) 2, 3, 4, 5-tetrahydrodipicolinate, (ii) L-aspartate, (iii) alpha-D-ribose 5-phosphate, (iv) 5-methyltetrahydrofolate, (v) L-glutamine, and (vi) L-serine. In addition, Bemisia tabaci SSA1-SG1 imports D-fructose, homocysteine, thiamine diphosphate, riboflavin, pyridoxine 5-phosphate and nicotinate D-ribonucleotide from the insect hemolymph to the bacteriocyte. On the other hand, it exports 39 and 54 metabolites to hemolymph and Portiera respectively.

Terminal reactions in the essential amino acid biosynthesis pathway are important because they influence the provisioning of amino acid. The generated metabolic reconstructions revealed that the terminal reaction of isoleucine, leucine, valine, histidine and phenylalanine are all mediated by genes of intrinsic origin in the host bacteriocyte of SSA1-SG1. For example, for branched-chain amino acids, one enzyme, branched-chain amino acid transaminases (EC 2.6.1.42) encoded by BCAT, mediate the conversion of 3-methyl-2-oxobutanoate to valine, 4-methyl-2-oxopentanoate to leucine and (S)-3-Methyl-2-oxopentanoate to isoleucine in SSA1-SG1 whitefly.

Essential gene targets in different essential amino acid biosynthesis pathways in SSA1-SG1

A total of 270 reactions were characterized as essential based on the growth rate. Of these, 129 reactions were mediated by genes in Portiera while 141 reactions were mediated by genes in B. tabaci SSA1-SG1 (Fig. S2). Single gene deletion also revealed that seven of the ten biosynthesis pathways had essential genes mediating their terminal reactions. Therefore, these genes were characterized as essential genes in the provisioning of essential amino acids in the cassava B. tabaci SSA1-SG1. These include, (i) argH (EC 4.3.2.1) for arginine biosynthesis, (ii) lysA (EC 4.1.1.20) for lysine biosynthesis (Fig. 4), (iii) branched-chain amino acid transaminase (EC 2.6.1.42) for isoleucine, valine and leucine biosynthesis, (iv) hisD (EC 1.1.1.23) for histidine and (v) aspC (EC 2.6.1.58) for phenylalanine biosynthesis (Fig. 5). Of these, two genes (argH and lysA) are horizontally transferred genes of bacterial origin while three genes (BCAT, hisD and aspC) are of intrinsic origin (Fig. 6). All the selected reactions except PHETA1 for phenylalanine biosynthesis were mediated by the protein encoded by a single gene. Therefore, three factors were considered in the selection of the best candidate symbiosis gene targets for the management of cassava whitefly SSA1-SG1. These include, (i) indispensability of the gene, (ii) reaction mediated by a single gene and (iii) gene mediating a terminal reaction.

Validation of expression of selected genes in Bemisia tabaci guts and bacteriocytes using real-time quantitative PCR

In addition to the selected osmoregulation genes (AQP1, SUC1 and SUC2), genes mediating terminal reactions in the different amino acid pathways were selected for validation of gene expression in the bacteriocyte. A total of five genes were selected as critical gene targets with enriched expression in the bacteriocyte compared to the whole body. These were (i) ENSSSA1UGG003633: diaminopimelate decarboxylase (LysA), essential in lysine biosynthesis, (ii) ENSSSA1UGG001672: arginosuccinate lyase (ArgH), essential in arginine biosynthesis, (iii) ENSSSA1UGG010036: branched-chain-amino-acid aminotransferase (BCAT), critical in isoleucine, leucine and valine biosynthesis, (vi) ENSSSA1UGG004484: aspartate aminotransferase (AAT), essential in phenylalanine biosynthesis and (viii) ENSSSA1UGG004185: 4-hydroxy-tetrahydrodipicolinate reductase (dapB), essential in lysine biosynthesis. Results from RNAseq analysis and RT-qPCR were similar with fold differences in the gut ranging from 22.30 to 26.54-folds for RT-qPCR compared to 21.71 to 32.67-folds from RNAseq analysis (Fig. 7). In the bacteriocyte, gene expression for the selected gene ranged from 0.57 to 10.07-folds for RT-qPCR compared to 1.79 to 13.89-folds for RNAseq analysis.