Transcriptional Response of Wolbachia to Dengue Virus Infection in Cells of the Mosquito Aedes aegypti.

Aedes aegypti transmits one of the most significant mosquito-borne viruses, dengue virus (DENV). The absence of effective vaccines and clinical treatments and the emergence of insecticide resistance in A. aegypti necessitate novel vector control strategies. A new approach uses the endosymbiotic bacterium Wolbachia pipientis to reduce the spread of arboviruses. However, the Wolbachia-mediated antiviral mechanism is not well understood. To shed light on this mechanism, we investigated an unexplored aspect of Wolbachia-virus-mosquito interaction. We used RNA sequencing to examine the transcriptional response of Wolbachia to DENV infection in A. aegypti Aag2 cells transinfected with the wAlbB strain of Wolbachia. Our results suggest that genes encoding an endoribonuclease (RNase HI), a regulator of sigma 70-dependent gene transcription (6S RNA), essential cellular, transmembrane, and stress response functions and primary type I and IV secretion systems were upregulated, while a number of transport and binding proteins of Wolbachia, ribosome structure, and elongation factor-associated genes were downregulated due to DENV infection. Furthermore, bacterial retrotransposon, transposable, and phage-related elements were found among the up- and downregulated genes. We show that Wolbachia elicits a transcriptional response to virus infection and identify differentially expressed Wolbachia genes mostly at the early stages of virus infection. These findings highlight Wolbachia's ability to alter its gene expression in response to DENV infection of the host cell. IMPORTANCE Aedes aegypti is a vector of several pathogenic viruses, including dengue, Zika, chikungunya, and yellow fever viruses, which are of importance to human health. Wolbachia is an endosymbiotic bacterium currently used in transinfected mosquitoes to suppress replication and transmission of dengue viruses. However, the mechanism of Wolbachia-mediated virus inhibition is not fully understood. While several studies have shown mosquitoes' transcriptional responses to dengue virus infection, none have investigated these responses in Wolbachia, which may provide clues to the inhibition mechanism. Our results suggest changes in the expression of a number of functionally important Wolbachia genes upon dengue virus infection, including those involved in stress responses, providing insights into the endosymbiont's reaction to virus infection.

INTRODUCTION

Arthropod-borne viruses (arboviruses) account for a significant number of viral infections in humans and livestock, representing a considerable global public health burden (1, 2). Dengue virus (DENV) is overwhelmingly the most important mosquito-borne flavivirus, responsible for an estimated 390 million infections annually with a further 2.5 billion people at risk of infection worldwide (1, 3). DENV consists of four distinct serotypes (DENV-1 to -4) and is a member of the Flaviviridae family that encompasses lipid-enveloped, single-stranded positive-sense RNA viruses (4). The main transmission vectors, the mosquitoes Aedes aegypti and Aedes albopictus, are distributed in high abundance worldwide, especially in regions with tropical and subtropical climates such as Asia, mid-Africa, Central America, and most of South America (5, 6). DENV has fully adapted to a human-mosquito-human transmission cycle, with A. aegypti being the primary vector of DENV in high-density urban environments (7, 8). A. albopictus is also capable of transmitting DENV; however, it predominantly colonizes suburban, rural, and forest areas (8). Dramatic increases in population density accompanied by rapid urbanization and cross-linked global transport systems are responsible for the rise of mosquito vectors throughout the tropical regions of the world (1, 7–9).

In light of the absence of an effective vaccine or clinical treatment, sustainable alternatives to control the global spread of DENV are being investigated. Historically, vector control programs against A. aegypti, reliant on chemical and biological agents, have been the methods of choice to reduce DENV transmissions (10). However, with the emergence of insecticide resistance in A. aegypti across many regions of the world, the effectiveness of chemical-based vector control interventions has declined dramatically (10–13). As a result, novel vector control strategies are needed to stop the expansion of mosquito-borne arboviruses.

A new effective approach to reduce the spread of arboviruses such as DENV, Zika virus (ZIKV), and chikungunya virus (CHIKV) uses the endosymbiotic bacterium Wolbachia pipientis as a biological control agent (14–16). It is estimated that various strains of Wolbachia infect 40% to 60% of all arthropod species (17–19). Wolbachia are Gram-negative, obligate intracellular alphaproteobacteria that have evolved to invade, influence, and manipulate their arthropod host’s reproductive system in order to establish within a host population (18, 20, 21). Wolbachia’s ability to proliferate through a population is mainly achieved by inducing cytoplasmic incompatibility (CI) in its host (20, 22). The phenomenon of CI occurs when a Wolbachia-infected male insect mates with a female insect that is uninfected or not infected with a compatible strain of Wolbachia, generating nonviable offspring (20, 22). Consequently, CI promotes the unaided dispersion of the exclusively maternally inherited Wolbachia infection into the host population (20, 22). Additionally, in Wolbachia-transinfected female mosquitoes, Wolbachia conveys a fitness advantage by protecting the mosquitoes against infectious microbes, including viruses (20, 23, 24). The ability of certain strains of Wolbachia to provide protection against RNA virus infections was originally discovered in Drosophila melanogaster by two independent groups of researchers in 2008 (25, 26). The discovery of this insect protection against RNA viruses revealed the association between virus resistance and Wolbachia infection, paving the way for an effective mechanism to control mosquito-borne diseases.

Wolbachia does not naturally occur in A. aegypti; however, stable transinfections can be achieved through microinjections (27). Several studies in A. aegypti investigating artificially introduced Wolbachia strains such as wMelPop, wMel, and wAlbB have demonstrated Wolbachia-induced changes in expression of immune-related genes (15, 28, 29). These findings led the way to two main hypotheses for Wolbachia-mediated virus inhibition, that of host innate immune priming and the competition for limited cellular host resources required for virus replication (14, 30). Mosquito species such as A. aegypti depend on their innate immune system to control arbovirus replication, as they lack the adaptive immune system of vertebrates (31). The interaction between the replicating DENV and its host A. aegypti’s innate immune system affects viral replication and the successive transmission of DENV (32). In this regard, a number of studies have found that the innate immune response and mediated pathways such as the RNA interference (RNAi), Toll, immune deficiency (IMD), and JAK/STAT signaling pathways have antiviral effects in A. aegypti (33, 34). Initially, the innate immune response in mosquitoes is activated through the recognition of pathogens by various receptors (35). The RNAi pathway is activated by the detection of double-stranded RNA (dsRNA) molecules of viral genomes or replication intermediates (34). Within the RNAi pathway, two RNA molecules, microRNA (miRNA) and small interfering RNA (siRNA), are vital for the RNAi pathway’s functionality (36). Interestingly, both miRNAs and siRNAs are modified upon Wolbachia infections (37). The inhibition of certain miRNAs induces a reduction of Wolbachia density, indicating that Wolbachia promotes its own preservation within the host (37). Several studies have confirmed that Wolbachia infections enhance the basal gene expression of the Toll, JAK/STAT, and RNAi pathways (15, 28, 29, 38). A study conducted in an A. aegypti (Aag2)-derived clonal cell line infected with Wolbachia demonstrated upregulation of endonuclease argonaute-2 (Ago2), a core catalytic component within the RNA-induced silencing complex (RISC) of the RNAi pathway (38). The Wolbachia-induced upregulation of Ago2 results in increased activity of the RNAi pathway, subsequently providing better protection against arboviruses (38). Overall, these findings have identified the RNAi pathway as the key insect antiviral pathway capable of limiting arboviruses such as DENV, CHIKV, and Sindbis virus (SINV) (31, 32, 39). However, studies suggest that RNAi may play no or a limited role in Wolbachia-mediated virus inhibition (38, 40).

Competition for limited cellular host resources also has been subject to investigations (14, 41). In particular, blood meal-derived molecules such as cholesterol are essential in mosquito egg development (41, 42). A. aegypti infected with Wolbachia strains wMelPop and wMel showed a 15% to 25% reduction in total cholesterol quantities compared to that in uninfected controls (41). However, subsequent cholesterol supplementation administered via blood meal did not correct the altered embryo morphology phenotype, and it did not improve fecundity or egg viability deficits (41). Understanding the contributions of A. aegypti, DENV, and Wolbachia in the virus inhibition phenotype is challenging and is further complicated by the fact that Wolbachia cannot be genetically modified (30, 43). The majority of knowledge to date has been gained from comparing the durability and level of Wolbachia-mediated virus inhibition in different combinations of vector species, virus genotypes, and Wolbachia strains (30, 44). Dissecting the mechanism of Wolbachia-mediated virus inhibition in A. aegypti proves difficult, and the inhibition is likely multilayered, with no individual mechanism responsible for this phenotype (30, 43). The diverse differences of Wolbachia strains infecting the same A. aegypti host confirm the contribution of the symbiont’s genome to the mosquito-Wolbachia association (30, 45). Remarkably, the ability and magnitude of Wolbachia-mediated antiviral inhibition can differ between Wolbachia strains (15, 46–50).

Thus far, the arthropod host A. aegypti response to viral infections has been extensively investigated; however, to our knowledge, the transcriptional response of Wolbachia to virus infections in the mosquito has not. For this reason, we examined whether Wolbachia’s transcriptional response to DENV infection may shed light on the Wolbachia-mediated viral inhibition mechanism based on the function of differentially expressed genes. To test our hypothesis, we compared the transcriptional response of Wolbachia wAlbB strain to DENV infection in a transinfected A. aegypti (Aag2)-derived cell line. The outcomes of this study provide a basis for follow-up investigations into the mechanism.

RESULTS

RNA-Seq data analysis of Aag2.wAlbB cells after infection with DENV.

To identify differentially expressed genes of Wolbachia wAlbB strain in response to DENV infection in transinfected A. aegypti Aag2 cells, DENV-2-infected and uninfected Aag2.wAlbB cells were collected at 1, 6, and 24 h postinfection (hpi) in three biological replicates for each time point. We chose these early time points, because previous studies have indicated that Wolbachia-mediated virus inhibition is exerted at early hours of infection in insect cells, including mosquito cells (51, 52). Prior to the experiment, wAlbB’s relative density was quantified by quantitative PCR (qPCR), which showed amplification of the wsp gene in Aag2.wAlbB cells, whereas no amplification of the gene was detected in the Aag2.tet cells (Aag2.wAlbB cells previously treated with tetracycline to remove Wolbachia) (Fig. 1A). Wolbachia density is expressed in terms of the number of wsp copies relative to the number of the mosquito’s RPS17 copies. Furthermore, increasing DENV genome copy accumulation was confirmed in the samples at time points 1, 6, and 24 hpi by quantitative reverse transcription-PCR (RT-qPCR) (Fig. 1B). While DENV genome copy numbers did not significantly change between 1 and 6 hpi, a significant increase in genome copy numbers was detected at 24 hpi. Based on a previous study in A. albopictus C6/36 cells, accumulation of viral genomic RNA of the four DENV serotypes started from 13 to 18 hpi (53). Therefore, no detectable viral replication between 1 and 6 hpi is expected.

FIG 1

Relative density of Wolbachia in Aag2.wAlbB cell line and infectivity of the DENV inoculum. (A) The relative Wolbachia density was measured by qPCR analysis using genomic DNA extracted from Aag2.wAlbB cells and tetracycline-cured Wolbachia‐free Aag2.tet cells as a control. Primers targeting the Wolbachia surface protein gene (wsp) and A. aegypti RPS17 to amplify the target and normalizing genes were used, respectively. (B) Validation of virus infection by RT-qPCR analysis confirming DENV-2 genome copy accumulation in infected Aag2.wAlbB samples used for RNA-Seq. Uninfected Aag2.wAlbB cells served as control (mock). Sample 24a-infected was excluded from analysis. The bars and error bars represent means and standard deviations (SDs), respectively, from the three biological replicates. One-way analysis of variance (ANOVA) with Tukey’s post hoc multiple-comparison tests were used for data analysis. ns, not significant; **, P < 0.01.

Illumina HiSeq sequencing of RNA from DENV-2-infected and uninfected Aag2.wAlbB cells generated a total of 1,629,779,502 raw reads with 150-bp paired-ends from the 18 RNA samples (see Table S2 in the supplemental material). Following the adapter trimming and quality filtering using the trim reads sequence tool in CLC Genomics Workbench (CLC-GWB), a total of 1,629,758,722 high-quality reads remained. Each sample and each biological condition were represented by an average of 92,742,982 and 88,341,320 million reads, respectively. The trimmed reads with a Phred quality score of >Q30 corresponding to an error rate of 0.04% were mapped to the A. aegypti three chromosomes, mitochondrion, wAlbB, DENV-2, and Aag2 cell line-related viruses (cell-fusing agent virus [CFAV] and Aedes albopictus negev-like virus [AaNLV]) reference genomes using the RNA sequencing (RNA-Seq) analysis tool in CLC-GWB as described in Materials and Methods. An average of 90,044,894 and 86,018,515 reads mapped to the A. aegypti, wAlbB, and virus genomes in the DENV-infected and uninfected samples, respectively, equivalent to 97.4% of the total reads (Fig. 2). Sample 24a infected with DENV-2 was excluded from further analysis, as we found that RNA from Aag2.tet cells at 24 hpi was accidently submitted for sequencing instead of Aag2.wAlbB at 24 hpi. However, the data were used for de novo assembly of the DENV-2 ET-300 genome (see below). For each sample and biological condition, averages of 27,005,569 and 29,785,073 reads for DENV-infected and uninfected samples, respectively, equivalent to 29.1% to 33.7%, were mapped to the wAlbB’s reference genome (54). Principal-component analysis (PCA) analysis showed distinct separation at 1, 6, and 24 hpi between DENV-infected and uninfected Aag2.wAlbB samples, indicating a high level of similarity among the biological replicates from different samples (Fig. 3). The 2.6% unmapped reads were used for de novo assembly of contigs using the CLC-GWB de novo assembly tool. BLASTn analysis of the assembled contigs to the nonredundant nucleotide NCBI DNA database produced a variety of long and short contigs. The longest contigs between 1,801 and 13,186 bp in size were identified as additional A. aegypti host-related viruses, Aedes anphevirus, Primus virus, Phasi Charoen-like virus, Aedes albopictus densovirus, Australian Anopheles totivirus, and unplaced scaffolds of A. aegypti genome artifacts (55).

FIG 2

Composition of Aag2.wAlbB normalized reads mapped to the combined reference file and the wAlbB reference genome. Aag2.wAlbB reads were mapped to A. aegypti, wAlbB, and viral reference genomes for DENV-infected (A) and uninfected (B) Aag2.wAlbB cells. Sample 24a-infected was excluded from analysis. The bars represent total numbers of reads mapped to all the genomes, with the red portion representing reads mapped to the wAlbB genome and the blue portion representing reads that did not map to the wAlbB genome but to the rest (A. aegypti and viral reference genomes).

FIG 3

Principal-component analysis (PCA) of the RNA-Seq samples. The PCA figure represents a two-dimensional scatterplot of the first two principal components of the RNA-Seq data. DENV-infected and uninfected Aag2.wAlbB sample groups at different time points are represented by different colors as indicated by the legend provided within the graph. Each dot represents a biological replicate of an RNA-Seq sample. Sample 24a-infected was excluded from analysis.

Information relevant to the RNA-Seq data. Download Table S2, XLSX file, 0.1 MB.

Full-length virus genome read counts expressed in normalized read counts (trimmed mean of M values [TMM] adjusted) confirmed the increasing DENV abundance in each DENV-infected Aag2.wAlbB sample used for RNA-Seq (Fig. 4).

FIG 4

Accumulation of DENV-2 read counts following infection in RNA-Seq samples. Validation of virus infection: full-length virus genome read counts expressed in normalized read counts (TMM adjusted) confirmed DENV abundance in each infected Aag2.wAlbB sample used for RNA-Seq. Sample 24a-infected was excluded from analysis. The bars and error bars represent means and SDs from the three biological replicates.

Transcriptional response by Wolbachia with respect to time and DENV infection.

A total of 96 genes were significantly differentially expressed across the three time points combined; 13, 31, and 52 differentially expressed genes (DEGs) at 1 versus 6, 6 versus 24, and 1 versus 24 hpi, respectively (Fig. 5; see also Table S3). There was a total of 16 DEGs that overlapped between the time points (Fig. 5). The open reading frames (ORFs) from the identified DEGs featured complete coding with no premature stop codon truncations. Approximately one-quarter of those DEGs, 26.0% (25), were upregulated versus 74.0% (56) being significantly downregulated. Specifically, 4, 13, and 8 genes were upregulated at 1 versus 6, 6 versus 24, and 1 versus 24 hpi in contrast to 9, 18, and 44 downregulated genes at the three time points, respectively (Table S3).

FIG 5

Differentially expressed genes in wAlbB across time points after DENV infection. Venn diagram showing DEGs (fold change of ≥2.0 and an adjusted P value of <0.05) in wAlbB in DENV-infected Aag2.wAlbB mosquito cells when different time points were compared. Sixteen overlapping DEGs were found between the time points. Each colored circle represents a sample collection time point.

Differentially expressed Wolbachia wAlbB genes upon DENV infection over time. Download Table S3, XLSX file, 0.1 MB.

Among the upregulated genes at 1 versus 6 hpi were those encoding transposases, conjugal transfer protein TraJ from the type IV secretion system, and MSF transporter. Downregulated genes identified in that comparison were those encoding an ATPase, ankyrin repeat domain-containing protein, RNase HI, ATP synthase subunit C, and five hypothetical proteins (Table S3). DEGs, when 6 and 24 hpi samples were compared, represented transposases, ankyrin repeat domain-containing proteins, RuvC, O-methyltransferase, a number of ribosomal proteins, maturase, ferredoxin, MspII family outer membrane protein, elongation factor Tu, a reductase, TrbC from the type IV secretory system, iron-sulfur cluster assembly accessory protein, sodium:proton antiporter, and a number of hypothetical proteins (Table S3). When data from 1 versus 24 hpi were compared, similar types of genes were differentially expressed, including those encoding transposases, O-methyltransferase, maturase, epimerase, ribosomal proteins, ion antiporters, ATPases, elongation factor Tu, and ankyrin repeat domain-containing proteins plus a number of other proteins that are all listed in Table S3.

Transcriptional response by Wolbachia to DENV infection.

When DENV-infected and uninfected samples at each time point were compared, a total of 22 genes were significantly differentially expressed across the three time points combined, considering a fold change of ≥2.0 and an adjusted P value of <0.05: 4, 5, and 13 DEGs at 1, 6, and 24 hpi, respectively (Fig. 6; see also Table S4). These included genes encoding DNA mismatch repair protein MutS, RNase HI, and transposases at 1 hpi, glutathione synthase, group II intron reverse transcriptase, and transposase at 6 hpi and O-methyltransferase, SPFH domain-containing protein, group II intron reverse transcriptases, sodium:proton antiporter, ribosomal proteins, and iron-sulfur cluster assembly accessory protein at 24 hpi (Table S4). There was only one DEG that overlapped between the time points, encoding DNA mismatch repair protein MutS (Fig. 6).

FIG 6

Differentially expressed genes in wAlbB due to DENV infection. Venn diagram showing DEGs (fold change of ≥2.0 and an adjusted P value of <0.05) in wAlbB in DENV-infected versus uninfected Aag2.wAlbB mosquito cells at each time point. One overlap of DEGs was found between the time points. Each colored circle represents a sample collection time point.

Differentially expressed Wolbachia wAlbB genes upon DENV infection at each time point. Download Table S4, XLSX file, 0.1 MB.

However, considering a fold change of ≥1.5, a total of 148 genes were significantly differentially expressed across the three time points combined: 39, 65, and 44 DEGs at 1, 6, and 24 hpi, respectively (Fig. S2; Table S4). There was a total of 11 overlapping DEGs between the time points (Fig. S2). These included genes encoding MutS, ATP-binding protein, DnaJ, two group II intron reverse transcriptases, ankyrin domain-containing protein, TatD hydrolase, cytidylyltransferase, ribosomal protein S9, and a hypothetical protein. A third of those DEGs, 25.0% (37), were upregulated versus 75.0% (57) being significantly downregulated. Specifically, 21, 1, and 15 genes were upregulated at 1, 6, and 24 hpi in contrast to 18, 64, and 29 downregulated genes at the three time points, respectively (Table S4).

Differentially expressed genes in wAlbB due to DENV infection. Venn diagram showing DEGs (fold-change of ≥1.5 and an adjusted P value of <0.05) in wAlbB in DENV-infected versus uninfected Aag2.wAlbB cells at each time point. Eleven overlapping DEGs were found between the time points. Each colored circle represents a sample collection time point. Download FIG S2, TIF file, 0.1 MB.

De novo assembly of DENV-2 ET-300 genome.

Despite the usage of DENV-2 ET-300 strain in several studies, the complete genome of the virus was not available on GenBank at the time this study took place. To assemble the viral genome, we used reads from Aag2.wAlbB sample 24a. In Aag2.wAlbB sample 24a, 678,422 of 100,812,616 (0.67%) reads mapped to a 10,696-nucleotide long contig, achieving an average coverage of 9,476×. The largest ORF from this contig showed complete coding with no premature stop codon truncations present. The ORF is 10,176 nucleotides in length and encodes a 3,392-amino-acid polyprotein. BLASTn analysis identified the ORF as 99.28% identical with the highest pairwise nucleotide identity (10,091/10,164) to dengue virus 2 genomic RNA, strain D2/Hu/OPD030NIID/2005 (GenBank identifier [ID] LC111438.1; E value, 0.0; query cover, 99%). BLASTp analysis of the ORF revealed that it is most similar to the polyprotein of dengue virus 2 (GenBank ID BAU36333.1; E value, 0.0; identity, 99.71%; query cover, 100%) and contains the three classic structural (C, prM, and E) and seven nonstructural (NS1, NS2A, NS2B, NS3, NS4A, NS4B, and NS5) proteins of the DENV polyprotein (58–60). The complete genome sequence obtained in this study is 100% identical to a recently deposited sequence for DENV-2 ET-300 (accession number MT921572).

Validation of DEGs by quantitative reverse transcription-PCR.

The representative selection of differentially expressed Wolbachia genes identified by RNA-Seq analysis was validated by RT-qPCR. For this, the DENV infection experiment was repeated in Aag2.wAlbB cells with identical conditions to the preparation of samples for RNA-Seq. DENV infection was confirmed in cells collected at 1, 6, and 24 hpi (see Fig. S3). The results showed an overall consistency between RNA-Seq and RT-qPCR, when DEGs in DENV-infected Aag2.wAlbB cells with respect to time and DENV infection were considered (Fig. 7A). Similarly, we found consistency between the two approaches when DENV-infected and uninfected samples were compared (Fig. 7B).

FIG 7

Validation of differentially expressed Wolbachia genes after DENV-2 infection using RT-qPCR. (A) The bar graphs represent RNA-Seq normalized gene reads and RT-qPCR relative expression results of the DEGs in DENV-infected Aag2.wAlbB cells with respect to time and DENV infection. (B) The bar graphs represent RNA-Seq normalized gene reads and RT-qPCR relative expression results of the DEGs in DENV-infected versus uninfected Aag2.wAlbB cells at 1 hpi.

Confirmation of DENV infection in Aag2.wAlbB cells. (A) RT-qPCR analysis confirmed DENV-2 genome copy accumulation in infected Aag2.wAlbB samples used for validation of RNA-Seq data by RT-qPCR. Uninfected Aag2.wAlbB cells served as a control (mock). The bars and error bars represent means and SDs from two biological replicates. Download FIG S3, TIF file, 0.1 MB.

DISCUSSION

The endosymbiotic bacterium Wolbachia has proven to be an effective biocontrol strategy to reduce the transmission of arboviruses by mosquitoes. A number of studies have confirmed that Wolbachia confers antiviral protection in primary arbovirus mosquito vectors (14–16, 61). Nonetheless, the mechanism underlying the Wolbachia-mediated viral protection is not yet well understood. The mosquito host A. aegypti response to viral infections has been extensively investigated (56, 62, 63); however, to our knowledge, the transcriptional response of Wolbachia to virus infections in the mosquito has not. For this reason, we examined whether Wolbachia’s transcriptional response to DENV infection could shed light on the Wolbachia-mediated viral inhibition mechanism based on the function of differentially expressed genes. We used RNA-Seq in order to study the transcriptional response of the Wolbachia wAlbB strain to DENV-2 infection in A. aegypti Aag2 cells.

We found an exponential increase of the transcriptional changes of Wolbachia with respect to time and DENV infection occurring between 1 and 6 hpi (13 DEGs), followed by 6 and 24 hpi (31 DEGs), with most changes between 1 and 24 hpi (52 DEGs). Similarly, the transcriptional response by Wolbachia to DENV infection when comparing DENV-infected versus uninfected Aag2.wAlbB cells at 1, 6, and 24 hpi showed increases in transcriptional changes across the three time points for 4, 5, and 13 genes, respectively. At 1 hpi, we identified a gene encoding an endoribonuclease (RNase HI) as significantly upregulated due to DENV infection. RNase HI plays an essential role in bacterial cellular processes such as DNA replication, gene expression, DNA repair, and degradation of RNA (64–66) and is an endoribonuclease with the capability of degrading RNA/DNA hybrids generated during viral replication in a non-sequence-specific manner (64, 66, 67). Furthermore, studies have suggested that RNase HI could play a vital role in antiviral defense in prokaryotes (68–71).

Using a fold change difference of ≥1.5 provided a more complete insight into Wolbachia’s transcriptional response to DENV infection and uncovered a number of significant differentially expressed genes related to stress response, secretion, and fundamental bacterial cellular functions based on annotation. Among the upregulated genes at 1 hpi were genes encoding the DNA mismatch repair protein MutS, molecular chaperone DnaJ, nucleotide exchange factor GrpE, ATP-binding cassette (ABC) domain-containing protein, and type I (HlyD family periplasmic adaptor subunit) and IV (conjugative DNA transfer family protein VirD4 and P-type conjugative transfer protein VirB9) secretory system elements. mutS encodes an essential DNA mismatch repair protein, MutS, that identifies and corrects errors generated during DNA replication and other biological processes (72, 73). The heat shock molecular chaperone protein-encoding genes dnaJ and grpE participate in responses to heat shock and cellular stress, stabilizing proteins by preventing their denaturation and actively refolding denatured proteins (74–78). Both dnaJ and grpE are part of the major DnaK/DnaJ chaperone system used by prokaryotic symbionts in response to stress (78, 79). The nucleotide exchange factor GrpE functions as a molecular cochaperone promoting the exchange of ADP for ATP from DnaK within the bacterial cytoplasm (75, 76). ABC domain-containing proteins are part of the ABC transmembrane transport system involved in the import of essential nutrients and the export of various molecules such as proteins, lipids, and toxins (80–82). In Wolbachia, ABC transmembrane transporters are embedded in the inner membrane and usually connected to a membrane fusion protein that crosses the periplasmic space (83). The type I secretion system (T1SS) is a basic secretion system spanning both the inner and outer membranes of Wolbachia (83–85). This enables the transport of substrates (e.g., DNA and proteins) from the cytoplasm to outside the cell into the host (83–85). The T1SS is composed of three proteins: an ABC transporter/ATPase (ABC) embedded in the inner membrane, a membrane fusion protein, and an outer membrane protein (83, 85). The T4SS protein VirD4 is an important ATPase known as a coupling protein responsible for substrate recruitment within the T4SS (86, 87). The T4SS transports substrates (e.g., DNA and proteins) from the cytoplasm to outside the cell or into various target recipients, including eukaryotic cells, and spans both membranes and the periplasm (83, 87–89). Essentially, all members of the Rickettsiales harbor T1SS and T4SS implicated in their establishment, pathogenicity, and transmission (83). Despite all of these genes being significantly upregulated or inducted at the very early time point of DENV infection, they were all downregulated over time, i.e., at 24 hpi.

Among the downregulated genes at 1 hpi is an RNase (RNAP), a DNA-dependent RNA polymerase known as the key enzyme in gene expression and regulation in bacteria (90–92). Despite the fact that the gene encoding the RNA polymerase was downregulated at 1 hpi, it was upregulated over time, and the mean gene expression across all biological replicates and time points remained very high. Furthermore, two ankyrin repeat-containing (ANK) proteins and a site-specific recombinase were identified as WO prophage-related genes (54, 93, 94). ANK (DEJ70_02765) and the site-specific recombinase (DEJ70_04945) are located within the wAlbB WO-like island 01 and wAlbB WO-like island 02, respectively. ANK proteins are rare in most bacteria; however, 34 wAlbB proteins contain at the minimum one copy of an ANK repeat domain, with a total of 81 copies of various types of ANK repeat domains (54). More than 30 Wolbachia strains from different insect and nematode hosts feature these ANK repeat regions (83, 94–96). These ANK proteins mediate a wide range of protein-protein interactions, such as transcription initiation and signaling, and are believed to be involved in Wolbachia-host interactions (54, 94–97).

At 6 hpi, a regulator of sigma 70-dependent gene transcription ssrS (6S RNA), a highly abundant noncoding RNA (ncRNA), was the only DEG upregulated, whereas all 64 of the other DEGs were downregulated (98–102). Concurrently, DENV genome copy numbers decreased notably at 6 hpi. Studies in Escherichia coli and other bacteria have shown that 6S RNA specifically binds σ70 RNA polymerase capable of inhibiting and downregulating transcription at σ70-dependent promoters (98–102). Furthermore, ncRNAs such as 6S RNA have been identified to perform several essential functions such as RNA processing, gene regulation, and interaction with proteins in addition to mRNAs (57, 98–102). In Wolbachia strains, wMel, wMelPop-CLA, wMelCS, and wAu studies discovered a range of ncRNAs and confirmed their biological functionality (57, 103). However, both studies failed to detect 6S RNA in the Wolbachia strains used. Mayoral et al. used Illumina TruSeq small RNA deep sequencing and detecting read fragments of 40 bp in size (103), whereas Woolfit et al. produced sequence fragments of 300 bp, acknowledging that this size is longer than most known small RNAs (sRNAs) (57). Eminently, 6S RNA in wAlbB is 162 bp in size, containing conserved prokaryotic promoter regions and boxes −35 and −10 indicative of its functionality to regulate bacterial and even host gene transcription, as previously described (103).

At the later time point, 24 hpi, the majority of DEGs were related to ribosomal structure, translation, DNA mismatch repair, transmembrane, and stress response. Furthermore, a number of secretion system subunit-encoding genes were downregulated due to DENV infection over time. Bacterial transposable elements and a number of hypothetical proteins were found among the upregulated and downregulated genes across all time points. IS3, and IS66 family transposase are insertion sequences (IS) which comprise the simplest form of transposable elements (TEs) of a transposase gene, flanked by inverted repeats (54, 104, 105). IS elements comprise 13% of the wAlbB genome, corresponding to a total number of 218 different IS element families (54).

Furthermore, 16 overlapping DEGs were identified between 6- and 24-hpi and 1- and 24-hpi time points. Most of these genes are related to basic bacterial cellular functions and TEs, such as amino acids, isoprenoid biosynthesis, methyltransferase, ribosomes, group II intron bacterial retrotransposon, and transposase.

The functional importance of the differentially expressed genes in response to DENV infection across the three different time points potentially points toward a transcriptional response by Wolbachia to DENV infection and the associated cellular stress experienced. Wolbachia elicits the upregulation of RNase HI, an endoribonuclease capable of degrading RNA-DNA hybrids generated during viral replication in a non-sequence-specific manner, due to DENV infection at 1 hpi (64, 66, 67). This is followed by a stress response using the major DnaK/DnaJ chaperone system (dnaJ and grpE upregulated due to DENV infection at 1 hpi), which is used by prokaryotic symbionts in response to thermal stress (54, 78, 79), while concurrently upregulating genes associated with T1SS and T4SS transmembrane transport functions, known to facilitate the translocation of effector proteins such as ANKs involved in Wolbachia-host interactions (54, 94–97). Despite the upregulation of the major DnaK/DnaJ chaperone system and secretion pathway subunits at the very early stage of DENV infection, they were downregulated over time. Remarkably, the majority of DEGs at 6 hpi were downregulated, while 6S RNA, a ncRNA and known negative regulator of sigma 70-dependent gene transcription (98–102), was upregulated.

Unlike what we found in mosquito cells, a couple of studies on Drosophila did not find major transcriptional changes in Wolbachia following virus infection. In one study, the investigators combined a Wolbachia-infected Drosophila melanogaster cell line JW18 with the mosquito-borne Semliki Forest virus (SFV; alphavirus) to examine the mechanism of antiviral protection (51). RNA-Seq was performed on samples collected at 7 and 24 hpi; however, Rainey et al. reported a lack of differentially expressed Wolbachia genes between cells infected with and without SFV (51). Similarly, when D. melanogaster flies harboring the wMel2 strain of Wolbachia were infected with Sindbis virus (an alphavirus), no transcriptional changes in the wMel2 transcriptome were observed at 6, 24, and 48 h postinjection (106). The discrepancy between our findings and the two studies in D. melanogaster cell lines and flies could be due to differences in the Wolbachia strains, host cells, and the viruses used in the studies.

In conclusion, our findings shed new light on Wolbachia’s transcriptional response to DENV infection of the host cell and identify differentially expressed genes starting as early as 1 hpi, which to our knowledge, have not been examined so far. Further research, especially into RNase HI, 6S RNA, and other proteins and signaling molecules which are involved in Wolbachia’s response to stress and virus infection is necessary to uncover their potential contribution to antiviral protection. Furthermore, considering Wolbachia and cytoplasmic RNA viruses compete for resources, changes in the expression of genes that affect the availability of energy and resources to the virus, imposed by Wolbachia, could contribute to reduction in virus propagation. This line of research is, however, currently hampered due to the unavailability of tools to genetically modify or silence Wolbachia genes.

MATERIALS AND METHODS

Mosquito cell line maintenance.

An Aag2.wAlbB cell line was developed previously (107) by transinfecting Aag2 cells with wAlbB from Aa23 cells. Aa23 cells infected with wAlbB were established from eggs of Aedes albopictus (108). Cells were maintained in 1:1 Mitsuhashi-Maramorosch (Himedia) and Schneider’s Drosophila medium (Invitrogen) supplemented with 10% fetal bovine serum (FBS; Bovogen Biologicals, France) at 27°C and passaged every 3 to 4 days.

Virus infection experiments.

Prior to the infection of Aag2.wAlbB cells, Wolbachia density was quantified by quantitative PCR (qPCR). Genomic DNA was extracted from cells using the DIY spin column protocol as previously described (109). qPCR was performed using, QuantiFast SYBR green PCR mix (Qiagen), in accordance with the manufacturer’s instructions. Forward and reverse primers (wAlbB-wsp-qF and wAlbB-wsp-qR) (see Table S1 in the supplemental material) targeting the Wolbachia surface protein gene (wsp) and A. aegypti ribosomal protein subunit 17 gene (RPS17) (AeRPS17-qF and AeRPS17-qR) (Table S1) were used to amplify the target and the normalizing gene, respectively (14, 37). All qPCRs were performed in duplicates using a Rotor-Gene Q thermocycler (Qiagen) under the following conditions: 95°C for 5 min and then 40 cycles of 95°C for 10 s and combined annealing/extension at 60°C for 35 s. This was followed by a standard post-qPCR high-resolution melt analysis to assess the specificity of the amplified wsp and RPS17 qPCR products. The relative abundance of wsp and RPS17 was determined using the relative quantification method as described previously (110).

Primers used in this study. Download Table S1, PDF file, 0.1 MB.

For virus infection, Aag2.wAlbB cells were seeded into 6-well culture plates (Greiner Bio-One) at 1.5 × 106 cells per well and allowed to adhere for 1 h at 27°C. Medium was then removed, and cells were inoculated with DENV serotype 2 (DENV-2) East Timor strain (ET-300) at a multiplicity of infection (MOI) of 1. Uninfected cells were treated the same as Aag2.wAlbB virus-infected cells except without the addition of virus. All 6-well plates were incubated on a rocker for 1 h at room temperature followed by the removal of supernatant and subsequent replenishment with fresh medium. The plates were then incubated at 27°C. Aag2.wAlbB virus-infected and uninfected cells were collected at 1, 6, and 24 h postinfection (hpi) in three biological replicates for each time point. In our experimental design, we used matched time points for infected and uninfected samples to take into consideration possible changes of gene expression over time in both Wolbachia and host cells in cell culture. Cells were pelleted by centrifugation at 2,500 × g for 3.5 min. Supernatant was discarded and pellets were frozen at −80°C for future RNA extraction.

Leftover RNA samples were used to confirm DENV-2 infection (see below for RNA extraction). For quantitation of DENV-2 genome copies in Aag2.wAlbB cells, first-strand cDNA synthesis was performed using M-MuLV reverse transcriptase (New England BioLabs) with a gene-specific (NS5) reverse primer (DENV2-qR) (Table S1) to amplify the DENV genomic RNA and A. aegypti RPS17 reverse primer (AeRPS17-qR) (Table S1) for the synthesis of the reference gene cDNA. In each RT reaction, 500 ng of DNase I-treated RNA was used as the template in a total volume of 20 μl. All the RT reactions were performed in a Veriti 96-well thermocycler (Applied Biosystems) under the following conditions: 65°C for 5 min, 50°C for 50 min, and 75°C for 15 min. For qPCR, the produced cDNA was diluted in a 1:5 ratio with Ultrapure DNase/RNase-free water (Invitrogen). Two microliters of the diluted cDNA was used for qPCRs with both forward and reverse gene-specific primers (DENV2-qF and DENV2-qR) (Table S1) used to amplify the DENV‐2 NS5 region of the viral genome (111). A. aegypti RPS17 forward and reverse primers were used to amplify the normalizing gene (14). All qPCRs were performed in duplicates using a Rotor-Gene Q thermocycler (Qiagen) under the conditions specified above.

RNA extraction and sequencing.

Total RNA from DENV-2-infected and uninfected Aag2.wAlbB cells was extracted using the Qiagen RNeasy minikit (Qiagen) according to the manufacturer’s instructions. The obtained RNA pellets were resuspended in 30 μl of Ultrapure DNase/RNase-free water (Invitrogen) and treated with DNase I enzyme (TURBO DNase; Invitrogen) to remove genomic DNA from RNA samples as per the manufacturer’s instructions. RNA sequencing (RNA-Seq) was conducted by GENEWIZ (GENEWIZ, China) under the following conditions. The RNA quality and quantity of each sample were assessed using a Bioanalyzer 2100 (Agilent Technologies), Agilent RNA 6000 Nano kit (Agilent Technologies), and 1% agarose gel. In total, 5 μg of total RNA per sample with an RNA integrity number of >7.0 was used for library preparation according to the manufacturer’s protocol using NEBNext Ultra RNA library prep kit for Illumina (New England BioLabs). After removal of ribosomal RNAs using a Ribo-Zero-rRNA removal kit (Illumina), the mRNA fragmentation and priming were performed using NEBNext first-strand synthesis reaction buffer and NEBNext random primers (New England BioLabs). First-strand cDNA was synthesized using ProtoScript II reverse transcriptase, and the second-strand cDNA was synthesized using a second-strand synthesis enzyme mix (New England BioLabs). The double-stranded DNA was purified by AxyPrep Mag PCR clean-up (Axygen) and treated with End Prep enzyme mix to repair both ends and add deoxyribosyladenine (dA) tailing in one reaction, followed by a T-A ligation to add adaptors to both ends. Size selection of adaptor-ligated DNA was then performed using AxyPrep Mag PCR clean-up (Axygen), and fragments of ∼360 bp (with the approximate insert size of 300 bp) were recovered. Each sample was then amplified by PCR using P5 and P7 primers (Table S1). The PCR products were cleaned up using AxyPrep Mag PCR clean-up (Axygen), validated using an Agilent 2100 Bioanalyzer (Agilent Technologies), and quantified by a Qubit 2.0 fluorometer (Invitrogen). Then, libraries with different indices were multiplexed and loaded on an Illumina HiSeq instrument according to the manufacturer’s instructions (Illumina). Sequencing was carried out using a 2 by 150-bp paired-end (PE) configuration, and image analysis and base calling were conducted by the HiSeq Control Software (HCS)+OLB+GAPipeline-1.6 on the Illumina HiSeq instrument (Illumina). The sequences were processed and quality controlled by GENEWIZ.

Data quality control using demultiplexing was performed by bcl2fastq 2.17, and RNA-Seq raw data were filtered as follows. Low-quality pair-end reads with ambiguous nucleotide content of N bases of more than 10% and a ratio of bases of low quality (Phred quality score Q < 20) of more than 0.5 in either read were filtered out and discarded before data analysis. The paired-end raw reads were obtained as fastq files and then imported into CLC Genomics Workbench v20.0.2 (CLC-GWB; Qiagen) for further processing and analysis. The raw reads were trimmed from the 3′ ends using the trim reads sequence tool in CLC-GWB in order to remove Illumina sequencing adapter P5 and P7. Default trim settings with automatic read-through adapter trimming and removal of low-quality sequences (Phred quality score = 0.05) were applied.

Differential gene expression analysis.

Aag2.wAlbB high-quality trimmed paired-end reads were mapped to a combined reference genome assembly containing the A. aegypti strain LVP_AGWG AaegL5.0 chromosomes 1 to 3 (NC_035107, NC_035108, and NC_035109), mitochondrion complete genome (MF194022), Wolbachia wAlbB genome (CP031221), cell-fusing agent virus strain Galveston genome (CFAV; NC_001564), Aedes albopictus negev-like virus isolate (AaNLV; MK879802), and the de novo assembled DENV-2 ET-300 genome (see below for assembly). It was previously shown that Aag2.wAlbB cells are persistently infected with CFAV and AalNLV (112). The mapping setting “genome annotated with genes and transcripts suitable for Eukaryotes” was used, which accounts for splicing. Mapping parameters with a minimum length fraction of 0.8, minimum similarity fraction of 0.9, and maximum number of hits for a read of 30 to capture multimapping repetitive elements were applied. The CLC-GWB principal-component analysis (PCA) for RNA-Seq tool was used to create a two-dimensional principal component plot to identify outlying samples and variations in the sample data set.

The empirical analysis of differentially expressed genes (DEGs) was performed using the CLC-GWB differential expression for RNA-Seq tool with the default parameters to identify DEGs between DENV-infected and uninfected Aag2.wAlbB samples. TMM (trimmed mean of M values) normalization was used to calculate effective library sizes as developed and described previously (113). The CLC-GWB differential expression for RNA-Seq tool uses multifactorial statistics based on a negative binomial generalized linear model (GLM) for this purpose (114). Differential gene expression in response to DENV infection was evaluated at 1, 6, and 24 hpi between virus-infected and uninfected samples, whereas the differential expression with respect to time and DENV infection was assessed between 1 versus 6, 6 versus 24, and 1 versus 24 hpi (Fig. S1). In the latter comparisons, up- or downregulation of a gene was considered over time. For both sets of analyses, we used CLC-GWB’s “all group pairs” analysis, which tests for differences between all pairs of groups in a factor, including a Benjamini-Hochberg correction for multiple comparisons, and accordingly adjusts P values. Significance of DEGs was assessed using the Wald test, and results were then filtered based on a fold change of ≥2.0, ≥1.5, and adjusted P value of <0.05. The CLC-GWB create Venn diagram for RNA-Seq tool was used to compare and visualize the overlap of DEGs between the three sampling time points. Open reading frame (ORF) integrity of identified DEGs was assessed using the Swiss Institute of Bioinformatics (SIB) ExPasy protein analysis tool. A summary of the procedures from read mapping to the analysis of wAlbB differentially expressed genes in response to DENV infection is shown in Fig. S1.

The flow chart shows the steps taken for the analysis of the RNA-Seq data followed by the analysis of wAlbB’s differentially expressed genes (DEGs) with respect to time and infection with dengue virus. X, Y, and Z, and their variations represent the number of differentially expressed genes between two conditions. Download FIG S1, TIF file, 0.2 MB.

For validation of DENV infection, normalized read counts mapping to the virus genome were used to estimate DENV abundance in each sample used for RNA-Seq.

De novo assembly of DENV-2 ET-300.

After adapter and quality trimming of Aag2.wAlbB, clean reads were de novo assembled using CLC-GWB de novo assembly tool with the default parameters for word length, bubble, and minimum contig length of 200 nucleotides. Assembled contig reads were mapped back to contigs with stringent mapping settings (mismatch cost, 2; insertion cost, 3; deletion cost, 3; length fraction, 0.8; similarity fraction, 0.8) to eliminate false-positive mapping. CLC-GWB find ORFs tool with default settings was used to identify ORFs within the assembled contig. The assembled contig was then queried using NCBI’s BLASTn (BLASTn v.2.10.1+) nonredundant nucleotide and BLASTp (BLASTp v.2.10.1+) nonredundant protein sequence databases.

Validation of DEGs using quantitative reverse transcription-PCR.

For validation of the differentially expressed Wolbachia genes determined through RNA-Seq analysis, a total of 8 up- and downregulated genes were selected. The NCBI Primer-BLAST primer design tool was used to design primers for the selected DEGs and the 16S rRNA normalizing gene (115). The primers are listed in Table S1. New RNA samples were generated by repeating the infection experiment with DENV-2 ET-300 under identical conditions as before by collecting virus-infected and uninfected Aag2.wAlbB cells at 1, 6, and 24 hpi. RNA extraction, reverse transcription, and qPCR were performed as described above, except that first-strand cDNA synthesis was performed using SuperScript III reverse transcription (Invitrogen) and random hexamer primers (New England BioLabs) to synthesize a cDNA strand from extracted RNA as per the manufacturer’s instructions. No reverse transcriptase control reactions were also included in the RT-qPCRs. Wolbachia’s 16S rRNA gene was used as the normalizing gene. The relative abundance of 16S rRNA and selected DEGs was determined using the relative quantification method as described previously (110).

Data availability.

The RNA-Seq data generated in this study were submitted to SRA under the project ID PRJNA669319.