Genome-Wide Profiling of Cardinium-Responsive MicroRNAs in the Exotic Whitefly, Bemisia tabaci (Gennadius) Biotype Q

Although the bacterial symbiont Cardinium has profound effects on the ecological adaptation of its host, the whitefly Bemisia tabaci (Gennadius) biotype Q (hereafter referred to as B. tabaci Q), the molecular mechanism underlying interactions between these two organisms is not yet fully understood. In this study, sRNA libraries were constructed, amplified, and sequenced for Cardinium-infected (C⁺) and uninfected (C^∗−) B. tabaci Q with identical genetic backgrounds. Subsequently, the genes targeted by the differentially expressed miRNAs were predicted by integrating the B. tabaci Q genome data. A total of 125 known and 100 novel miRNAs were identified, among which 23 significant differentially expressed miRNAs were identified in both libraries. It is noteworthy that an analysis of target genes showed that Cardinium-responsive miRNA-regulated genes were functional in apoptosis, reproduction, development, immune response, thermotolerance and insecticide resistance. GO and KEGG pathway analysis revealed that some miRNA-target genes are closely associated with energy metabolism. A major finding of this study was the identification of several miRNAs that may be involved in physiological processes in response to environmental stress, i.e., insecticides and high temperatures. This information will provide a foundation to help further elucidate the functions of Cardinium in B. tabaci Q.

Introduction

Most phloem-feeding insects are infected with maternally inherited intracellular bacteria that have been categorized as either obligate symbionts or facultative symbionts (Gaelen et al., 2009). One of these facultative symbionts, Candidatus Cardinium hertigii (hereafter referred as Cardinium), was initially characterized in Encarsia wasps (Zchori-Fein et al., 2004). The symbiont has been shown to be a reproductive manipulator in several arthropod taxa through various means including cytoplasmic incompatibility (CI), feminization, and induction of parthenogenesis (White et al., 2011; Penz et al., 2012). Further, more recent research has demonstrated that Cardinium can also affect the fitness of its host (Harris et al., 2010; Fang et al., 2014).

Our long-term field surveys showed that the infection rate by Cardinium remains at a low level (12.2%) in the exotic sweet potato whitefly, Bemisia tabaci (Gennadius) Q (hereafter referred to as B. tabaci Q) in the Shandong Province, China (Chu et al., 2011). Further research revealed that competitive ability and fitness are diminished in Cardinium-infected B. tabaci populations compared to the Cardinium-uninfected population (Fang et al., 2014). Fitness (survival and reproductive success) of a field collected strain of Cardinium-infected B. tabaci Q, however, was shown to have unexpectedly improved after being reared under constant laboratory conditions for > 3 years (Li Hong-Ran, unpublished data). We speculated that the potential negative effect of Cardinium on its host whitefly may be the result of environmental stress under field conditions, i.e., in response to insecticides and high temperatures. To date, little information is available regarding the molecular mechanism underlying the Cardinium-mediated phenotypes of their whitefly hosts (Li et al., 2018). A more thorough understanding of this mechanism will help understand the ecological adaptation and population dynamics of this introduced pest.

MicroRNAs (abbreviated as miRNAs) are numerous non-coding small RNAs (sRNAs) (20–24 nt) involved in the regulation of gene expression post-transcriptionally, through repression of mRNA translation or degradation of mRNA in the cytoplasm (Bühler et al., 2006; Bartel, 2009; Iliopoulos et al., 2009). As a result of ongoing research, the significant roles of miRNAs in the manipulation of metabolism, development, and epigenetic inheritance are becoming better understood (Esau et al., 2006; Carlsbecker et al., 2010; Wang et al., 2018). Accumulating evidence has indicated that miRNAs are involved in regulating diverse biological processes such as apoptosis, hematopoiesis, and patterning of the nervous and immune systems in various groups of insects (Bushati and Cohen, 2007; Asgari, 2013; Xu et al., 2017). More recently, dramatic changes in the expression levels of miRNAs have also been documented in response to symbiont and host interactions. For instance, the microRNA, aae-miR-2940-5p, which is highly enriched in Wolbachia-infected mosquitoes, may be involved in mediating the regulation of pelo in Aedes aegypti (Asad et al., 2018). Also, Wolbachia is able to use host miRNAs to regulate host gene expression, and manipulation of the mosquito metalloprotease gene via aae-miR-2940 is crucial for efficient maintenance of the endosymbiont (Hussain et al., 2011). An analysis of target genes in the two-spotted spider mite, Tetranychus urticae, showed that Wolbachia-responsive miRNAs regulate genes associated with lysosome function, apoptosis, sphingolipid metabolism, and lipid transport in both sexes, as well as reproduction in females (Rong et al., 2014). Because miRNAs are involved in infections by bacterial symbionts, such as Wolbachia, we postulated that the host miRNAs might also be involved in Cardinium-host interactions. A test of this hypothesis would provide insights into Cardinium-mediated phenotypes in host whiteflies.

Based on our current knowledge, sufficient data are available to characterize the profiles of miRNAs in whiteflies. Our previous results showed different expression levels of miRNAs in B. tabaci B vs. Q (Guo et al., 2013). Furthermore, Wang et al. (2016) investigated the expression profiles of miRNAs in B. tabaci B in response to viral infection. Still, insufficient research has been directed at exploring the potential role of miRNAs in B. tabaci following symbiont infection. In this study, we first established a Cardinium-infected B. tabaci Q strain (abbreviated as C⁺) and an uninfected strain (abbreviated as C^∗−) with identical genetic backgrounds. Deep sequencing of four sRNA libraries was then conducted for both C⁺ and C^∗−. Subsequently, bioinformatics analyses were carried out to categorize the sRNAs, assess differential expression of known and novel miRNAs, predict their targets and perform their functional annotation. These data will provide novel insight into the mechanism of insect-symbiont interactions, especially for Cardinium in B. tabaci Q.

Materials and Methods

Whitefly Colony

The B. tabaci Q colony used in this study was established from individuals originally collected in the Shandong Province, China, in July 2012, and maintained in separate cultures for 3 years on potted cotton plants (Lu-Mian-Yan 21 cultivar). Whiteflies were kept in isolated screen cages under constant conditions of 27 ± 1°C, with a 16/8 h light/dark photoperiod, and 70 ± 5% RH. An introgressive backcrossing scheme was used to homogenize nuclear genetic backgrounds of infected and uninfected whiteflies for 6 generations, using the method described by Turelli and Hoffmann (1991). After completion of the introgression series, > 98% of nuclear alleles were expected to be shared between the C⁺ and C^∗− whiteflies; the two strains were then regarded as having a nearly identical genetic background (Figure 1). The one day-old adult females that were prepared for RNA isolation were selected from these Cardinium-infected and uninfected B. tabaci Q strains. To obtain newly emerged adult females, all adult whiteflies were removed from cotton leaves containing whitefly pupae each evening (at 6:00 pm). The following morning (at 8:00 am), newly emerged adults were individually collected and transferred into plastic tubes and sexed using a stereomicroscope (Luan et al., 2008; López and Andorno, 2009). All samples were stored at −80°C until needed for RNA isolation.

FIGURE 1

FIGURE 1. Experimental design and schematic diagram of the workflow used in this study.

Cardinium Detection in Populations

Presence of the symbiont Cardinium in each population was verified every 30 days from a sample of 20 adults using PCR primers to amplify a portion of the 16S rRNA gene (Figure 2). The primers used were CLO-F (5′- GCGGTGTAAAATGAGCGTG -3′) and CLO-R1 (5′- ACCTMTTCTTAACTCAAGCCT -3′) (Weeks et al., 2003). All PCRs were performed using 13 μL samples of a solution containing 1 × buffer, 0.16 mM of each dNTP, 0.5 mM of each primer, 0.5 unit Taq DNA polymerase (Takara, Dalian, China), and 2 μL template DNA. Cycling conditions consisted of initial denaturation at 95°C for five min, followed by 35 cycles of 1 min at 94°C for denaturation, 1 min at 58°C for annealing and 1 min at 72°C for elongation, plus a final extension step at 72°C for 7 min. All amplicons were electrophoresed, along with a negative (sterile water instead of DNA) and positive controls (DNA from previous sequencing) of the symbiont on a 1.0% agarose gel, and visualized using Gelview staining.

FIGURE 2

FIGURE 2. Amplification product of Cardinium-infected and uninfected Bemisia tabaci Q using specific primers. 1–10: Cardinium-infected B. tabaci Q; 11–20: Cardinium-uninfected B. tabaci Q; + : positive control; -: negative control; M: DL2000 DNA marker.

Small RNA Library Construction for Illumina Sequencing

Total RNA was isolated from the C⁺ and C^∗− strains using the Trizol reagent (Invitrogen Catalog No. 15596-026) with a slight modification of the procedure recommended by the manufacturer. In this procedure, RNA was purified several times with isopropanol and 75% ethyl alcohol. Total RNA (>3 μg) of good quality was used to construct an sRNA library for each sample by using a TruSeq small RNA Sample Pre Kit (Illumina). Briefly, total RNA was ligated with 5′ and 3′ adapters followed by reverse transcription using RT primers. After PCR amplification of the cDNAs, amplified PCR products within the 130–160 bp size range were purified on an 8% polyacrylamide gel (100 V, 80 min). The libraries were constructed with purified RNAs and were sequenced on an Illumina HiSeq 2500/2000 platform at the Novogene Company, Beijing, China.

Bioinformatics Analysis of Sequencing Data

After Illumina sequencing, raw data were processed using Novogene’ s Perl and Python scripts. Clean data were screened to remove reads containing more than three N (undetermined bases), reads with 5′ adapter contaminants, reads without 3′ adapter or the insert tags, those containing poly A, T, G, or C and low quality reads obtained from the raw data. Then, sRNA sequences of 18–35 nt were selected to conduct all downstream analyses. To prevent every unique sRNA mapping to multiple non-coding RNA (ncRNAs), we used the following priority rule: known miRNA > rRNA > tRNA > snRNA > snoRNA > repeat > gene > novel miRNA so that every unique sRNA mapped to only one annotation. The Bowtie software was used to map the sRNA tags to the B. tabaci Q genome¹) with less than 2 bp mismatch to analyze their distribution and expression (Langmead et al., 2009).

Next, the mappable sRNA tags were used to search the release of version 20.0 of miRBase to identify known miRNAs in B. tabaci. Known miRNAs were defined as sequences that were identical to sequences from Drosophila melanogaster or other insects (Aedes aegypti, Apis mellifera, Tribolium castaneum and Bombyx mori) as outlined previously (Guo et al., 2013). Then, rRNAs, tRNAs, snRNAs, and snoRNAs were removed by mapping the remained sRNA tags to Rfam². Repeat sequences were discarded by using a repeat sequence database³, and protein coding genes were filtered by mapping to the exon and intron of mRNAs of B. tabaci. Finally, novel miRNAs were predicted by exploring the secondary structure, the Dicer cleavage site and the minimum free energy of the former unannotated sRNA tags that could be mapped to the reference sequence by integrating two available software, miREvo and mirdeep2 (Wen et al., 2010; Friedlander et al., 2011).

Differential Expression Analysis of miRNAs Between C ⁺ and C^∗− Strains

miRNA levels (transcripts per million; TPM) were estimated using the following normalization formula: Normalized expression = (mapped read count/total reads)^∗1,000,000 (Zhou et al., 2010). Since there were two biological replicates for each sample, differential expression analysis of miRNAs between libraries was performed using the DESeq R package (1.8.3) (Anders and Huber, 2010). P-values were adjusted using the method of Benjamini and Hochberg (1995). The threshold level of collected P-values was set as 0.05 for significant differential expression.

Target Prediction and Enrichment Analysis

We used the miRanda software to identify genes targeted by miRNAs in the genome of B. tabaci⁴ (Enright et al., 2003). To further reveal functions related to the putative target genes, GO and KEGG enrichment analysis of the predicted target genes was performed using the GOSeq/topGO2.12 and KOBAS 2.0 software (Mao et al., 2005; Kanehisa et al., 2008; Young et al., 2010).

The dataset that was studied in this article is available in the NCBI (SRA) public repository under accession number SRP076077.

Validation of Differentially Expressed miRNAs via qRT-PCR

Total RNA was extracted from the C⁺ and C^∗− whiteflies (taken from the same treatment samples as used for library construction). Reverse transcription was performed by Mir-X^TM miRNA First-Strand Synthesis Kit (Takara catalog no. 638313) under the following conditions: 37°C for 1 h, 85°C for 5 min, and then held at 4°C. Levels of miRNA were assessed in four biological and two technical replicates of each sample, using SYBR Premix EX TaqTM II (Takara catalog no. RR820A). The cDNA was quantified in 20μL reactions, containing 7.2 μL ddH₂O, 10 μL SYBR Advantage Premix (2X), 0.8 μL miRNA-specific Primer (10 uM), mRQ 3′ Primer (supplied with the Kit No. RR820A) and 2 μL template cDNA. The primers used in this study are provided in Supplementary Table S1. The reactions were incubated at 95°C for 10 s, followed by 40 cycles of 95°C for 5 s, 60°C for 20 s. A dissociation curve was obtained to ensure that only one product was amplified after the amplification phase. The 2^−ΔΔCt method for relative quantification of gene expression was used to determine the level of miRNA expression.

Results and Discussion

Overview of the Analysis of sRNA Libraries

In this study, we generated two replicate sRNA libraries for each of the C⁺ and C^∗− strains (Table 1), totaling 44 million raw reads. About 43 million clean reads remained after the filtering step (see section “Materials and Methods” for details). Analysis of the length distribution among the 18–35 nt sRNAs examined here indicated that the highest proportion of small RNAs, in both strain-specific libraries, was in the 21–22 nt range (57.71 and 66.64% in the C⁺ and C^∗− libraries, respectively) (Figure 3). The proportion of small RNAs in the 21–22 nt range is consistent with that observed for miRNAs in animals (Li et al., 2011). But in previous studies examining miRNAs in B. tabaci B/Q and non-viruliferous whiteflies, length distribution showed maximum enrichment for 21–23 nt or 28–30 nt sRNAs (Guo et al., 2013; Wang et al., 2016). This difference may be due to sample treatment prior to library construction. Here, only one day-old adult females were utilized for total RNA isolation, whereas the whiteflies used in previous reports were neither isolated by sex nor segregated by time of emergence.

TABLE 1

TABLE 1. Statistics of small RNA sequences in C⁺ and C^∗− libraries.

FIGURE 3

FIGURE 3. Length distribution of small RNAs in the two libraries. C ⁺, Cardinium-infected B. tabaci Q; C^∗−, Cardinium-uninfected B. tabaci Q.

The mapped sRNAs fell in one of 12 categories, including known miRNAs, rRNAs, tRNA, snRNA, snoRNA, repeat, novel miRNAs, exon, intron and unannotated sRNAs (Figure 4). The percentages of known miRNAs in the C^∗− libraries (19.90%) were higher than those in the C⁺ libraries (13.41%), whereas the corresponding values for novel miRNAs were similar (10.97 and 10.41%, respectively). Thus, differences of both length distribution and composition of two sRNA libraries suggest that infection by Cardinium may inhibit the synthesis of miRNAs in B. tabaci. The same results were found in T. urticae infected with the symbiont Wolbachia (Rong et al., 2014).

FIGURE 4

FIGURE 4. Classification of sRNAs in the C ⁺ and C^∗− libraries. snRNA, small nuclear RNA; snoRNA, small nucleolar RNA; exon: + , same orientation with mRNA exon; exon: -, reverse orientation with mRNA exon; intron: + , same orientation with mRNA intron; intron: -, reverse orientation with mRNA intron; other, unannotated sRNA.

Identification of Known and Novel miRNA

Sequences in our libraries identical to miRNA sequences of D. melanogaster, Aedes aegypti, Apis mellifera, T. castaneum and B. mori were considered to be potentially known miRNAs. Following a BLASTn search against the miRBase 20.0 and subsequent sequence analysis, a total of 125 known miRNAs were identified in the four libraries (Table 1 and Supplementary Table S2). Next, we used the miRNA prediction software miRDeep2 to identify putative novel miRNAs by searching against the B. tabaci Q genome sequence⁵. In total, 100 novel miRNAs were identified from both libraries (Table 1). To determine whether these novel miRNAs are known in other species, we used these miRNA sequences to search against all nucleotide sequences in miRBase Version 20.0. The results showed that 79 novel miRNA sequences of B. tabaci were similar to those in other species, but differed in some nucleotide positions, while 21 novel miRNA sequences were not found in other species (Supplementary Table S3). Thus, 21 novel miRNAs are putatively specific to whiteflies and could thus have a species-specific functions.

Among all identified known and novel miRNAs, three miRNAs [novel-122 (efu-miR-9216), novel-156 (hsa-miR-4762-3p) and tca-miR-263b-5p] were detected only in C⁺ and 11 miRNAs were found only in C^∗−, while 157 miRNAs were found in both libraries (Figure 5). miR-263a/b affects expression of genes contributing to cellular and humoral immunity in Manduca sexta and regulates immunity-related signal transduction by affecting the expression of genes related to the Galleria mellonella tumor necrosis factor receptor superfamily, suggesting that the specific tca-miR-263b-5p in Cardinium-infected B. tabaci Q may be a regulator of immune processes (Mukherjee and Vilcinskas, 2014; Zhang et al., 2014). An earlier study also reported that miR-263a/b negatively regulates apoptosis, chaeta development and compound eye morphogenesis in D. melanogaster, suggesting that tca-miR-263b-5p may play a role in cell apoptosis in C⁺ strain (Hilgers et al., 2010).

FIGURE 5

FIGURE 5. Venn chart of miRNA expression in the C ⁺ vs. C^∗− strains. The blue circle represents C ⁺, Cardinium-infected B. tabaci Q; the red circle represents C^∗−, Cardinium-uninfected B. tabaci Q.

Differentially Expressed miRNAs Between C⁺ and C^∗− Strains

The expression levels of known and novel miRNAs in C⁺ and C^∗− libraries were compared after normalization (normalized expression = actual miRNA count/total count of clean reads × 100000). We identified 23 differentially expressed miRNAs, including 13 up-regulated miRNAs and 10 down-regulated miRNAs in the C⁺ strain relative to the C^∗− strain (Table 2 and Figure 6). The strongest up-regulation was observed for novel-140, novel-27 (cte-miR-87b), novel-256 (ppc-miR-2235a-9-3p) and ame-bantam, while the strongest down-regulation was observed for novel-19, tca-miR-993-5p, bmo-miR-993b-5p and novel-193 (pma-miR-133a).

TABLE 2

TABLE 2. List of all significantly differentially expressed miRNAs in C⁺ vs. C^∗− strains.

FIGURE 6

FIGURE 6. Expression of miRNAs in C⁺ vs. C^∗−. The x-axis represents the fold change of miRNA expression in different samples; the y-axis represents the significance of the variances on transcription levels of miRNAs. Each point represents one miRNA. Red points represent the up-regulated miRNAs. Blue points indicate equally expressed miRNAs. Green points represent down-regulated miRNAs. P < 0.05.

For up-regulated miRNAs, previous studies showed that bantam miRNA simultaneously stimulates cell proliferation and prevents apoptosis in Drosophila, indicating that an up-regulated level of ame-bantam in the C⁺ strain may help resist apoptotic cell death (Brennecke et al., 2003). miR-87 has been reported to be associated with anti-pathogen and immune responses (Liu et al., 2015). Thus up-regulation of novel-27 (cte-miR-87b) in the present study suggests it may play a role in immune responses in Cardinium-infected B. tabaci Q. With respect to down-regulated miRNAs, lmi-miR-133 was shown to mediate phenotypic plasticity and behavioral changes between social and solitary phases of henna and pale locusts by targeting key players in the dopamine synthesis pathway, indicating that novel-193 (pma-miR-133a) may modulate phenotypic plasticity in B. tabaci Q (Yang et al., 2014).

Prediction of Target Genes for Differentially Expressed miRNAs

To further clarify the biological functions of the differentially expressed miRNAs, we analyzed the targets of 23 of them. A total of 12760 target genes for these 23 sRNAs listed in Table 2 were predicted. Among them, the number of target genes for each differentially expressed miRNAs ranged from 339 to more than one thousand. These results are consistent with previous findings indicating multiple miRNAs can be used to regulate a single gene and, conversely, that a single miRNA can target multiple genes, confirming the complexity of the miRNA-gene regulation network (Thomson et al., 2011).

To further assess the presumptive functions of the genes predicted to be targeted by differentially expressed miRNAs, GO (Gene Ontology) annotation enrichment was performed. GO ontologies containing four different biological processes and 20 molecular functions were predicted (Figure 7 and Supplementary Table S4). The enrichment of each GO term within biological processes was compared, and most were involved in metabolic processes. As for molecular functions, the molecular function terms and most related binding terms comprised most of the targets, which is consistent with a regulatory role in transcription and translation for these miRNAs (Hobert, 2008). In addition, a KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis was preformed to elucidate the biological interpretation of the genes targeted by differentially expressed miRNAs (Figure 8). A total of 90 highly diversified biochemical pathways were identified that involved the miRNA-targeted genes, after considering a corrected p-value and the number of genes involved (Supplementary Table S5). Furthermore, a total of 1739 genes were enriched in the top 20 pathways. The “Metabolic pathways” category was the most significantly enriched term with 812 genes, followed by “Carbon metabolism” (98 genes) and “Starch and sucrose metabolism” (82 genes). The “rich factor” (x axis of Figure 8) represents the ratio of the number of differentially expressed genes to the total number of unigenes in the pathway. The most significant degree of enrichment was observed for the glycolysis “Glycolysis/Gluconeogenesis,” “Inositol phosphate metabolism” and “Nitrogen metabolism” pathways, all related to energy metabolism.

FIGURE 7

FIGURE 7. GO analysis results of the target genes of C⁺ vs. C^∗− strains. The x-axis is the GO category; the y-axis is the number of genes.

FIGURE 8

FIGURE 8. The 20 most enriched KEGG pathways based on target genes of differentially expressed miRNAs in C⁺ vs. C^∗− strains. The x-axis shows the rich factor. The y-axis shows the pathway names. The size of each point represents the number of genes enriched in a particular pathway. The larger the value of the rich factor and the smaller the value of q-value indicate a more significant degree of enrichment.

Cardinium-Responsive miRNAs May Be Associated With Apoptosis Process in B. tabaci Q

The target genes of several significantly differentially expressed miRNAs in B. tabaci Q have been identified to be associated with apoptosis (Table 2). The gene BTA013163.1, targeted by down-regulated novel_19, encodes glucosylceramidase, whose deficiency can cause neurodegeneration and apoptosis in the brain of mice (Enquist et al., 2007). The observed down-regulation of novel_19 in C⁺ strain may indicate its role in inhibiting apoptosis. Moreover, the gene BTA021458.1 targeted by novel_140 codes for Apoptosis-stimulating of p53 protein 2. Therefore, the observed up-regulation of novel_140 in C⁺ strain may also indicate its role in decreasing apoptosis. It is well known that Wolbachia and Cardinium infection can lead to decreases in apoptosis in their respective hosts, Drosophila mauritiana and Encarsia suzannae (Fast et al., 2011; Mann et al., 2017). Female-biased genes encoding ribosomal proteins indicate an increase in general translational activity of Cardinium in female wasps (Mann et al., 2017). In this study, the gene BTA027909.1 targeted by down-regulated novel_19, encodes 30S ribosomal protein S11, indicating a potential increase in general translational activity of Cardinium. We, therefore, speculated that Cardinium in B. tabaci Q probably could inhibit apoptosis in its host, which in turn, would be beneficial to the maintenance of Cardinium itself.

Cardinium-Responsive miRNAs May Regulate Reproduction and Development in B. tabaci Q

Cardinium can induce a number of reproductive disruptions in arthropods to optimize its transmission in a manner similar to that described for Wolbachia, including CI and feminization (Gotoh et al., 2007; White et al., 2011). CI was expected to an increase in the production of symbiont-infected females in the population (Poinsot et al., 2003; Werren et al., 2008). In the present study, some target genes of differentially expressed miRNAs were found to be related to female reproduction. For example, the target gene BTA017001.1 (targeted by novel_193) codes for Histone H3.3, which is believed to play an essential role in early developmental manifestations of CI (Serbus et al., 2008). In addition, the gene BTA022535.1, which is targeted by down-regulation of novel_36, codes for “meiosis arrest female protein 1 homolog” and is involved in oogenesis (Table 2).The miRNA tca-miR-13a-3p was up-regulated, which is expected to result in the down-regulation of BTA019714.1, a gene involved in female sterility. Moreover, target genes BTA013707.1 (targeted by dme-miR-305-5p), BTA017354.1 (targeted by tca-miR-993-5p), BTA023288.2 (targeted by tca-miR-305-5p) and BTA004944.1 (targeted by bmo-miR-993b-5p) code for cuticle protein in specific developmental stages. The down-regulated miRNAs should result in enhanced cuticular protein production, which may be related to growth and development. Thus, our results appear to support previous observations to the effect that Cardinium may manipulate host female reproduction and development, although this hypothesis will need to be further tested (Rong et al., 2014).

Cardinium-Responsive miRNAs May Be Essential in B. tabaci Q Under Stress Conditions

Accumulating evidence indicates that miRNAs play an important role in drug resistance, as revealed by miRNA expression profiling (Xin et al., 2008; Ma et al., 2010; Zheng et al., 2010). In the present study, target gene BTA004618.1 (targeted by novel_24) codes for multidrug resistance-associated protein 7, and target genes such as BTA017976.1 (targeted by novel_140), BTA002500.1 (targeted by novel_27) and BTA029742.1 (targeted by novel_256) are related to glutathione S-transferase and cytochrome P450. In addition, target gene BTA026695.3 (targeted by novel_256) and BTA028304.1 (targeted by aae-miR-277-3p) are related to drug resistance as well. Those up-regulated miRNAs may reduce their targeted gene expression. Additionally, the target gene BTA011400.1 (targeted by novel_140) and BTA019100.1 (targeted by novel_256), encoding “Heat shock protein (HSP) 70 B2” and “Activator of 90 kDa HSP ATPase homolog 1.” The up-regulated miRNAs may result in reduced HSP protein production, which may depress tolerance. Taken together, our observations support our earlier speculation that the C⁺ strain may have lower thermotolerance, insecticide resistance and detoxification than the C^∗− strain.

qPCR Validation of Differentially Expressed miRNAs

To further validate the expression patterns of miRNAs identified in this work, we randomly selected nine differentially expressed miRNAs to measure the expression levels via quantitative real time polymerase chain reaction (qRT-PCR). The results showed that seven miRNAs had qPCR expression patterns similar to those revealed by our RNA-seq analysis in the C⁺ and C^∗− strains. For example, novel-140 had been found to be highly expressed in the C⁺ strain, displaying a 2.8 log2-fold change compared with the C^∗− strain (Table 2 and Figure 9A), while similar results were obtained through a qPCR assessment, where approximately 4 times greater expression was found in the C⁺ strain (Figure 9B). However, two miRNAs, novel_1 and ame-miR-87, exhibited a decrease in expression in C⁺ as compared to C^∗− strain, which is the opposite of what the Illumina sequencing revealed. The reason for these differences is currently unknown. Overall, the qPCR results tend to indicate that the Illumina sequencing data accurately reflected the expression level of miRNAs.

FIGURE 9

FIGURE 9. Expression patterns of differentially expressed miRNAs identified based on Illumina sequencing (A) and qRT-PCR (B). X-axis, name of the miRNAs that were selected for qRT-PCR assessment; column above the X-axis: miRNAs that were up-regulated in the C⁺ strain; column below the X-axis: miRNAs that were upregulated in the C^∗− strain.

Conclusion

Using high-throughput sRNA sequencing, we screened out 23 miRNAs exhibiting differential expression in response to Cardinium infection in B. tabaci Q. These miRNAs are involved in several biological processes, including cell apoptosis, reproduction and development. The major finding of this study is the identification of several miRNAs overexpressed in infected whiteflies that target genes involved in thermotolerance and insecticide resistance, pointing to a compromised resistance to heat and xenobiotic stresses in C⁺ individuals. These findings may help explain fitness variation among Cardinium-infected B. tabaci Q whiteflies in response to environmental stress. In general, our findings lay a solid foundation for further functional study of the interactions between Cardinium and its host whitefly, B. tabaci Q.

Ethics Statement

The present research complies with all laws of the country (China) in which it was performed and was approved by the Department of Science and Technology of the Qingdao Agricultural University, China (Permit Number: 20110712).

Author Contributions

DC contributed to experimental design and management. HL carried out the data analysis and drafted the manuscript. XW and TD participated in data analysis. All authors read and approved the final manuscript.

Funding

This research was supported by the National Natural Science Foundation of China (31572064) and the Taishan Mountain Scholar Constructive Engineering Foundation of Shandong.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Novogene Company (Beijing) for helping with small RNA sequencing.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2018.01580/full#supplementary-material

Footnotes

1. ^ http://gigadb.org/dataset/100286
2. ^ http://rfam.xfam.org/
3. ^ http://www.repeatmasker.org/cgi-bin/WEBRepeatMasker/
4. ^ http://gigadb.org/dataset/100286
5. ^ http://gigadb.org/dataset/100286

References

Anders, S., and Huber, W. (2010). Differential expression analysis for sequence count data. Genome Biol. 11:R106. doi: 10.1186/gb-2010-11-10-r106

PubMed Abstract | CrossRef Full Text | Google Scholar

Asad, S., Hussain, M., Hugo, L., Osei-Amo, S., Zhang, G., Watterson, D., et al. (2018). Suppression of the pelo protein by Wolbachia and its effect on dengue virus in Aedes aegypti. PLoS Negl. Trop. Dis. 12:e0006405. doi: 10.1371/journal.pntd.0006405

PubMed Abstract | CrossRef Full Text | Google Scholar

Asgari, S. (2013). MicroRNA functions in insects. Insect Biochem. Mol. Biol. 43, 388–397. doi: 10.1016/j.ibmb.2012.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartel, D. P. (2009). MicroRNAs: target recognition and regulatory functions. Cell 136, 215–233. doi: 10.1016/j.cell.2009.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. 57, 289–300.

Google Scholar

Brennecke, J., Hipfner, D. R., Stark, A., Russell, R. B., and Cohen, S. M. (2003). bantam encodes a developmentally regulated microRNA that controls cell proliferation and regulates the proapoptotic gene hid in Drosophila. Cell 113, 25–36. doi: 10.1016/S0092-8674(03)00231-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Bühler, M., Verdel, A., and Moazed, D. (2006). Tethering RITS to a nascent transcript initiates RNAi- and heterochromatin-dependent gene silencing. Cell 125, 873–886. doi: 10.1016/j.cell.2006.04.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Bushati, N., and Cohen, S. M. (2007). microRNA functions. Annu. Rev. Cell Dev. Biol. 23, 175–205. doi: 10.1146/annurev.cellbio.23.090506.123406

CrossRef Full Text | Google Scholar

Carlsbecker, A., Lee, J. Y., Roberts, C. J., Dettmer, J., Lehesranta, S., Zhou, J., et al. (2010). Cell signalling bymicroRNA165/6 directs gene dose-dependent root cell fate. Nature 465, 316–321. doi: 10.1038/nature08977

PubMed Abstract | CrossRef Full Text | Google Scholar

Chu, D., Gao, C. S., De Barro, P., Zhang, Y. J., Wan, F. H., and Khan, I. A. (2011). Further insights into the strange role of bacterial endosymbionts in whitefly, Bemisia tabaci: comparison of secondary symbionts from biotypes B and Q in China. Bull. Entomol. Res. 101, 477–486. doi: 10.1017/S0007485311000083

PubMed Abstract | CrossRef Full Text | Google Scholar

Enquist, I. B., Bianco, C. L., Ooka, A., Nilsson, E., Månsson, J. E., Ehinger, M., et al. (2007). Murine models of acute neuronopathic Gaucher disease. Proc. Natl. Acad. Sci. U.S.A. 104, 17483–17488. doi: 10.1073/pnas.0708086104

PubMed Abstract | CrossRef Full Text | Google Scholar

Enright, A. J., John, B., Gaul, U., Tuschl, T., Sander, C., and Marks, D. (2003). MicroRNA targets in Drosophila. Genome Biol. 5:R1. doi: 10.1186/gb-2003-5-1-r1

PubMed Abstract | CrossRef Full Text | Google Scholar

Esau, C., Davis, S., Murray, S. F., Yu, X. X., Pandey, S. K., Pear, M., et al. (2006). miR-122 regulation of lipid metabolism revealed by in vivo antisense targeting. Cell Metab. 3, 87–98. doi: 10.1016/j.cmet.2006.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, Y. W., Liu, L. Y., Zhang, H. L., Jiang, D. F., and Chu, D. (2014). Competitive ability and fitness differences between two introduced populations of the invasive whitefly Bemisia tabaci Q in China. PLoS One 6:e100423. doi: 10.1371/journal.pone.0100423

PubMed Abstract | CrossRef Full Text | Google Scholar

Fast, E. M., Toomey, M. E., Panaram, K., Desjardins, D., Kolaczyk, E. D., and Frydman, H. M. (2011). Wolbachia enhance Drosophila stem cell proliferation and target the germline stem cell niche. Science 334, 990–992. doi: 10.1126/science.1209609

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedlander, M. R., Mackowiak, S. D., Li, N., Chen, W., and Rajewsky, N. (2011). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40, 37–52. doi: 10.1093/nar/gkr688

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaelen, R. B., Benjamin, B. N., Colin, F., and Nancy, A. M. (2009). Evolution and diversity of facultative symbionts from the aphid subfamily Lachninae. Appl. Environ. Microb. 75, 5328–5335. doi: 10.1128/AEM.00717-09

PubMed Abstract | CrossRef Full Text | Google Scholar

Gotoh, T., Noda, H., and Ito, S. (2007). Cardinium symbionts cause cytoplasmic incompatibility in spider mites. Heredity 98, 13–20. doi: 10.1038/sj.hdy.6800881

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, Q., Tao, Y. L., and Chu, D. (2013). Characterization and comparative profiling of miRNAs in invasive Bemisia tabaci (Gennadius) B and Q. PLoS One 9:e59884. doi: 10.1371/journal.pone.0059884

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, L., Kelly, S., Hunter, M., and Perlman, S. (2010). Population dynamics and rapid spread of Cardinium, a bacterial endosymbiont causing cytoplasmic incompatibility in Encarsia pergandiella (Hymenoptera: Aphelinidae). Heredity 104, 239–246. doi: 10.1038/hdy.2009.130

PubMed Abstract | CrossRef Full Text | Google Scholar

Hilgers, V., Bushati, N., and Cohen, S. M. (2010). Drosophila microRNAs 263a/b confer robustness during development by protecting nascent sense organs from apoptosis. PLoS Biol. 2010:e1000396. doi: 10.1371/journal.pbio.1000396

PubMed Abstract | CrossRef Full Text | Google Scholar

Hobert, O. (2008). Gene regulation by transcription factors and microRNAs. Science 319, 1785–1786. doi: 10.1126/science.1151651

PubMed Abstract | CrossRef Full Text | Google Scholar

Hussain, M., Frentiu, F. D., Moreira, L. A., O’Neill, S. L., and Asgaria, S. (2011). Wolbachia uses host microRNAs to manipulate host gene expression and facilitate colonization of the dengue vector Aedes aegypti. Proc. Natl. Acad. Sci. U.S.A. 108, 9250–9255. doi: 10.1073/pnas.1105469108

PubMed Abstract | CrossRef Full Text | Google Scholar

Iliopoulos, D., Hirsch, H. A., and Struhl, K. (2009). An epigenetic switch involving NF-kappaB, Lin28, Let-7 MicroRNA, and IL6 links inflammation to cell transformation. Cell 139, 693–706. doi: 10.1016/j.cell.2009.10.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa, M., Araki, M., Goto, S., Hattori, M., Hirakawa, M., Itoh, M., et al. (2008). KEGG for linking genomes to life and the environment. Nucleic Acids Res. 36, 480–484. doi: 10.1093/nar/gkm882

PubMed Abstract | CrossRef Full Text | Google Scholar

Langmead, B., Trapnell, C., Pop, M., and Salzberg, S. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10:R25. doi: 10.1186/gb-2009-10-3-r25

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H. R., Harwood, J. D., Liu, T. X., and Chu, D. (2018). Novel proteome and acetylome of Bemisia tabaci Q in response to Cardinium infection. BMC Genomics 19:523. doi: 10.1186/s12864-018-4907-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, T. T., Wu, R. M., Zhang, Y., and Zhu, D. H. (2011). A systematic analysis of the skeletal muscle miRNA transcriptome of chicken varieties with divergent skeletal muscle growth identifies novel miRNAs and differentially expressed miRNAs. BMC Genomics 12:186. doi: 10.1186/1471-2164-12-186

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y. X., Zhou, Y. H., Wu, J. Y., Zheng, P. M., Li, Y. J., Zheng, X. Y., et al. (2015). The expression profile of Aedes albopictus miRNAs is altered by dengue virus serotype-2 infection. Cell Biosci. 5:16. doi: 10.1186/s13578-015-0009-y

PubMed Abstract | CrossRef Full Text | Google Scholar

López, S. N., and Andorno, A. V. (2009). Evaluation of the local population of Eretmocerus mundus (Hymenoptera: Aphelinidae) for biological control of Bemisia tabaci biotype B (Hemiptera: Aleyrodidae) in greenhouse peppers in Argentina. Biol. Control 50, 317–323. doi: 10.1016/j.biocontrol.2009.05.002

CrossRef Full Text | Google Scholar

Luan, J. B., Ruan, Y. M., Zhang, L., and Liu, S. S. (2008). Pre-copulation intervals, copulation frequencies, and initial progeny sex ratios in two biotypes of whitefly, Bemisia tabaci. Entomol. Exp. Appl. 129, 316–324. doi: 10.1111/j.1570-7458.2008.00777.x

CrossRef Full Text | Google Scholar

Ma, J., Dong, C., and Ji, C. (2010). MicroRNA and drug resistance. Cancer Gene Ther. 17, 523–531. doi: 10.1038/cgt.2010.18

PubMed Abstract | CrossRef Full Text | Google Scholar

Mann, E., Stouthamer, C. M., Kelly, S. E., Dzieciol, M., Hunter, M. S., and Schmitz-Esser, S. (2017). Transcriptome sequencing reveals novel candidate genes for Cardinium hertigii-caused cytoplasmic incompatibility and host-cell interaction. mSystems 21:e00141-17. doi: 10.1128/mSystems.00141-17

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, X., Cai, T., Olyarchuk, J. G., and Wei, L. (2005). Automated genome annotation and pathway identification using the KEGG orthology (KO) as a controlled vocabulary. Bioinformatics 21, 3787–3793. doi: 10.1093/bioinformatics/bti430

PubMed Abstract | CrossRef Full Text | Google Scholar

Mukherjee, K., and Vilcinskas, A. (2014). Development and immunity-related microRNAs of the lepidopteran model host Galleria mellonella. BMC Genomics 15:705. doi: 10.1186/1471-2164-15-705

PubMed Abstract | CrossRef Full Text | Google Scholar

Penz, T., Schmitz-Esser, S., Kelly, S. E., Cass, B. N., Müller, A., Woyke, T., et al. (2012). Comparative genomics suggests an independent origin of cytoplasmic incompatibility in Cardinium hertigii. PLoS Genet. 8:e1003012. doi: 10.1371/journal.pgen.1003012

PubMed Abstract | CrossRef Full Text | Google Scholar

Poinsot, D., Charlat, S., and Merçot, H. (2003). On the mechanism of Wolbachia induced cytoplasmic incompatibility: confronting the models with the facts. Bioessays 25, 259–265. doi: 10.1002/bies.10234

PubMed Abstract | CrossRef Full Text | Google Scholar

Rong, X., Zhang, Y. K., Zhang, K. J., and Hong, X. Y. (2014). Identification of Wolbachia-responsive microRNAs in the two-spotted spider mite, Tetranychus urticae. BMC Genomics 15:1122. doi: 10.1186/1471-2164-15-1122

PubMed Abstract | CrossRef Full Text | Google Scholar

Serbus, L. R., Casper-Lindley, C., Landmann, F., and Sullivan, W. (2008). The genetics and cell biology of Wolbachia-host interactions. Annu. Rev. Genet. 42, 683–707. doi: 10.1146/annurev.genet.41.110306.130354

CrossRef Full Text | Google Scholar

Thomson, D. W., Bracken, C. P., and Goodall, G. J. (2011). Experimental strategies for microRNA target identification. Nucleic Acids Res. 39, 6845–6853. doi: 10.1093/nar/gkr330

PubMed Abstract | CrossRef Full Text | Google Scholar

Turelli, M., and Hoffmann, A. A. (1991). Rapid spread of an inherited incompatibility factor in California Drosophila. Nature 353, 440–442. doi: 10.1038/353440a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, B., Wang, L. L., Chen, F. Y., Yang, X. L., Ding, M., Zhang, Z. K., et al. (2016). MicroRNA profiling of the whitefly Bemisia tabaci middle East-Aisa Minor I following the acquisition of Tomato yellow leaf curl China virus. Virol. J. 13, 183–197. doi: 10.1186/s12985-016-0469-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z. Z., Ye, X. Q., Shi, M., Li, F., Wan, Z. H., Zhou, Y. N., et al. (2018). Parasitic insect-derived mi RNAs modulate host development. Nat. Commun. 9:2205. doi: 10.1038/s41467-018-04504-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Weeks, A. R., Velten, R., and Stouthamer, R. (2003). Incidence of a new sex-ratio-distorting endosymbiotic bacterium among arthropods. Proc. R. Soc. Lond. B 270, 1857–1865. doi: 10.1098/rspb.2003.2425

PubMed Abstract | CrossRef Full Text | Google Scholar

Wen, M., Shen, Y., Shi, S., and Tang, T. (2010). miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinformatics 13:140. doi: 10.1186/1471-2105-13-140

PubMed Abstract | CrossRef Full Text | Google Scholar

Werren, J. H., Baldo, L., and Clark, M. E. (2008). Wolbachia: master manipulators of invertebrate biology. Nat. Rev. Microbiol. 6, 741–751. doi: 10.1038/nrmicro1969

PubMed Abstract | CrossRef Full Text | Google Scholar

White, J. A., Kelly, S. E., Cockburn, S. N., Perlman, S. J., and Hunter, M. S. (2011). Endosymbiont costs and benefits in a parasitoid infected with both Wolbachia and Cardinium. Heredity 106, 585–591. doi: 10.1038/hdy.2010.89

PubMed Abstract | CrossRef Full Text | Google Scholar

Xin, F., Li, M., Balch, C., Thomson, M., Fan, M., Liu, Y., et al. (2008). Computational analysis of microRNA profiles and their target genes suggests significant involvement in breast cancer antiestrogen resistance. Bioinformatics 25, 430–434. doi: 10.1093/bioinformatics/btn646

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, J., Xu, X. X., Li, S. Z., Wang, S., Xu, X. J., Zhou, X. Q., et al. (2017). Genome-Wide Profiling of Plutella xylostella immunity-related miRNAs after Isaria fumosorosea infection. Front. Physiol. 8:1054. doi: 10.3389/fphys.2017.01054

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, M., Wei, Y., Jiang, F., Wang, Y., Guo, X., He, J., et al. (2014). MicroRNA-133 inhibits behavioral aggregation by controlling dopamine synthesis in locusts. PLoS Genet. 10:e1004206. doi: 10.1371/journal.pgen.1004206

PubMed Abstract | CrossRef Full Text | Google Scholar

Young, M. D., Wakefield, M. J., Smyth, G. K., and Oshlack, A. (2010). goseq: gene ontology testing for RNA-seq datasets. Genome Biol. 11:R14. doi: 10.1186/gb-2010-11-2-r14

PubMed Abstract | CrossRef Full Text | Google Scholar

Zchori-Fein, E., Perlman, S. J., Kelly, S. E., Katzir, N., and Hunter, M. S. (2004). Characterization of a ‘Bacteroidetes’ symbiont in Encarsia wasps (Hymenoptera: Aphelinidae): proposal of ‘Candidatus Cardinium hertigii’. Int. J. Syst. Evol. Microbiol. 54, 961–968. doi: 10.1099/ijs.0.02957-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Zheng, Y., Jagadeeswaran, G., Ren, R., Sunkar, R., and Jiang, H. (2014). Identification of conserved and novel microRNAs in Manduca sexta and their possible roles in the expression regulation of immunity-related genes. Insect Biochem. Mol. Biol. 47, 12–22. doi: 10.1016/j.ibmb.2014.01.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, T., Wang, J., Chen, X., and Liu, L. (2010). Role of microRNA in anticancer drug resistance. Int. J. Cancer 126, 2–10. doi: 10.1002/ijc.24782

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, L., Chen, J., Li, Z., Li, X., Hu, X. D., Huang, Y., et al. (2010). Integrated profiling of microRNAs and mRNAs: microRNAs located on Xq27.3 associate with clear cell renal cell carcinoma. PLoS One 5:e15224. doi: 10.1371/journal.pone.0015224

PubMed Abstract | CrossRef Full Text | Google Scholar