The transcriptome of the mosquito Aedes fluviatilis (Diptera: Culicidae), and transcriptional changes associated with its native Wolbachia infection

RNA sequencing and de novo transcriptome assembly

Functional annotation

Annotation was performed through local alignments using BlastX (with an e-value of 1e-10) against the NCBI non-redundant database (NR), which returned annotation for 26,066 contigs. Of these, 14,209 contained unique hits (Table 1). As previous molecular characterization of Ae. fluviatilis has been limited to a few microsatellites [33], our analysis relied entirely on finding similarity to other genomes. Consequently, 31,947 contigs (55% of the entire data set) did not get a significant BlastX hit, which could be explained by misassembly or from lack of representation in the NR database. 920 of these were considered to be protein-coding sequences based on predictions made by Transdecoder, but no further characterization was possible for the remaining 31,027 contigs. For the purposes of our analysis we have assumed that gene function is identical to what was described in the blast hits, but this may not be the case, as genetic homology does not always imply functional homology. However, we are fairly confident given that majority of annotations were previously described in closely related mosquitoes.

The assembled transcriptome of Ae. fluviatilis showed a high degree of similarity to Ae. aegypti, with 18,082 (69.38%) of the annotated contigs most closely matched to that species (Fig. 1). This was not surprising as Ae. aegypti has the most thoroughly annotated genome of the mosquitoes in the Aedes genus, meaning that this observed similarity may have more to do with the prevalence of Ae. aegypti sequences in public databases than the closeness of the relationship between the two species [34]. A further 5,616 (21.55%) of the annotated contigs were homologous to genes from other mosquito species, and overall 94.9% of the matches from BlastX were arthropod in origin. We conducted divergence analyses to better clarify the relationship between Ae. fluviatilis and other mosquito species (Additional file 1), and determined that Ae. fluviatilis diverged from these species approximately 98 million years ago (95% highest posterior density interval: 64.1 to 133.5 million years ago). A further 383 contigs matched to non-arthropod animals.

Across both the Tet and wFlu libraries we identified 751 contigs of bacterial origin (2.86%), and 112 of fungal origin (0.43%), which may potentially represent part of the Ae. fluviatilis intestinal microbiota, or could be the result of environmental contamination during sample collection. The bacterial sequences represented 250 distinct taxa, with the majority associated with a single contig. The diversity of sequence origins may indicate that the microbiota comprises many species, or that the majority of gut flora belonged to species that were not represented in the NR database. It is also possible that some of these contigs actually came from Wolbachia, but were different in composition from previously sequenced Wolbachia genes. Given the relatively low number of contigs of bacterial origin, the low copy numbers of these contigs and the low number of hits for each bacterial taxon in our data set, it was not possible to gain significant insight into the potential influence of wFlu on the composition or role of the microbiota.

Only 42 contigs of Wolbachia origin were identified in our data (Additional file 2). These contigs came from a variety of strains including wAna, wBm, wMelPop, and wPip. 38 were expressed in the wFlu libraries, and 4 in the Tet libraries - potentially via integration into the host genome. Wolbachia genomes typically contain between 900–1400 coding sequences [35]. Of the 42 identified Wolbachia contigs, 22 were hypothetical proteins with unknown function. Interestingly, one contig was a heme biosynthesis protein, which suggests that wFlu promotes heme production, and provides further evidence to support the theory that Wolbachia alters host iron and redox homeostasis, even amongst native associations [36, 37]. We also identified a glycosidase hydrolase, which plays a role in sugar metabolism. The remainder of the Wolbachia contigs fell into two broad functional categories. The first included genes involved in DNA/RNA processing, DNA repair and RNA synthesis, and are likely part of normal Wolbachia replication and transcription processes, as similar genes have been identified in the wMel genome [38]. A further group of contigs included Wolbachia membrane proteins, and ankyrin genes, which are used to attach the bacterial membrane to the host cytoskeleton. These genes likely represent part of the machinery used to mediate Wolbachia-host interactions [39], are typically present in large numbers in Wolbachia genomes [35, 38, 40–42], and may also be involved in CI [43].

The genomes of both Wolbachia strains and their insect hosts typically contain large quantities of mobile elements including transposons [38]. These can alter or disrupt the expression of genes, depending on their point of insertion into the host genome [34, 44]. We found 61 contigs related to transposable elements, none of which were matched to a Wolbachia genome. We also found a further 28 of viral origin, which could represent sequences from past or present members of the Ae. fluviatilis viral flora, which have not been well characterised.

As part of the contig annotation process, gene ontology (GO) and KEGG orthology (KO) terms for each contig were identified using Blast2GO and the KEGG database, respectively (Fig. 2). We identified genes associated with a wide range of biological processes, functions and structures. From the whole set of contigs, 23,035 were assigned to at least one GO term. The majority of contigs in the biological process category were associated with cellular or metabolic process. Cell, and cell part were the cellular component terms with the greatest frequency, while binding and catalytic activity were the most common terms associated with the molecular function category. The most common KO terms were related to diseases, and molecular information processing.

Differentially expressed contigs

Each RNA-Seq library was independently mapped to the assembled transcriptome using Bowtie. An average of 95% of the reads from each library were successfully mapped, and these data were used to generate counts for each contig. A list of differentially expressed contigs was generated using the baySeq package from Bioconductor, with 66% of these contigs (data not shown) also determined to be differentially expressed via analysis with the DESeq2 package, again from Bioconductor. Given this high level of concordance between the lists, we chose to proceed with further analysis of the baySeq list, as that method of analysis is known to be more sensitive [45]. Through baySeq we calculated the FPKM (Fragments Per Kilobase of transcript per Million mapped reads) for each contig, which considers the size of the contig in base pairs, and the overall data set size, as a measure of expression.

A total of 257 differentially expressed contigs were identified using baySeq, 159 were associated with the wFlu libraries, and 98 with the Tet libraries. Of these, 50 (19.4%) from the wFlu libraries (upregulated by wFlu) and 32 (32.7%) from the Tet libraries (downregulated by wFlu) had no matches after the BlastX search. This left 109 and 66 possible orthologs to known genes with increased and decreased expression due to Wolbachia infection, respectively (Additional file 3). Excluding multiple matches to genes with the same accession number, there were 95 unique contigs that were upregulated by Wolbachia, and 59 unique contigs that were downregulated, on the whole mosquito level. It is possible that analysis of whole mosquitoes precluded detection of genes that experienced a pattern of tissue-specific differential expression.

To confirm the accuracy of the sequencing data, the expression of six contigs was examined using RT-qPCR (Fig. 3). Two contigs, comp10645_c1_seq2 and comp10453_c0_seq5 were predicted to be upregulated by wFlu, and further analysis confirmed that expression was significantly higher in wFlu mosquitoes (Mann Whitney U tests: 10645 - P < 0.0001; 10453 - P < 0.0001). A further 2 contigs were predicted to be downregulated by wFlu, and the first of these, comp15178_c0_seq1, demonstrated that pattern (Mann Whitney U test: P < 0.0001). However for the second, comp14155_c0_seq1, expression levels were not significantly different due to the presence of Wolbachia (Unpaired t test: t = 0.1196 P = 0.9057). All four of those contigs were predicted to display differential expression in both the bayseq and DESeq2 analyses. Only low levels of comp14155_c0_seq1 were detected during sequencing, which could explain the lack of differential expression observed during RT-qPCR. The final 2 contigs that were examined were predicted to have equivalent expression levels between treatments, and both fit that pattern (Unpaired t tests: comp2025_c0_seq1 - t = 1.500, P = 0.1449; comp2041_c0_seq1 - t = 0.7994, P = 0.4308). Further validation of the differential expression of these genes using different techniques will likely prove valuable for future research.

Immune stimulation and suppression

Transinfection with Wolbachia in mosquito species has typically led to widespread increases in the expression of immune genes, including those involved in the Toll and IMD immune pathways, and a large number of antimicrobial peptides [10, 11, 14, 20]. This immune activation has been linked to interference of dengue viruses in mosquitoes [10, 14, 15], and immune suppression has been hypothesized as a potential cause of the enhancement of Plasmodium infection seen in some native and transient Wolbachia infections [28, 29]. We observed no changes in the transcription of genes in the Toll or IMD pathways, including genes such as rel1, or rel2, which might have been indicative of a change in pathway regulation. We also saw no changes in the expression of any of the antimicrobial peptides commonly affected by Wolbachia. The absence of systemic immune activation is common amongst native Wolbachia infections, and may be symptomatic of increased tolerance on the part of the host, and reduced pathogenicity on the part of the symbiont [46, 47].

We did observe differential expression of 6 contigs directly involved in mosquito immunity, 4 upregulated (2 cell wall hydrolases, a galactose specific c-type lectin, and a gram negative bacteria binding protein) and 2 downregulated (a mucin like protein, and a galactose specific c-type lectin). These genes are typically associated with bacterial binding and degradation, but many have also been linked to Plasmodium infection and could have contributed to the enhancement of P. gallinaceum infection in Ae. fluviatilis [32]. Gram negative bacteria binding proteins, for instance can have a broader role in immune stimulation, and are involved in the response to Plasmodium infection in An. gambiae [48]. Some galactose-specific c-type lectins are highly differentially expressed by Wolbachia [11], their expression is stimulated by Plasmodium infection, and they have been shown to protect Plasmodium against melanisation by the host immune system [49]. There are also examples of mucin like genes and cell wall hydrolases that are critical to Plasmodium development in mosquitoes [50, 51].

We also identified other differentially expressed genes have also been linked with insect immunity or Plasmodium infection. Plasmodium infection is dependent on glycolysis, and the metabolism of amino acids and lipids [52], all affected by wFlu infection (see below). There was evidence of altered sugar metabolism and transport, which could promote the development of Plasmodium [53]. Additionally, there was increased expression of a β-hexosaminidase, and these genes have a peptidoglycan hydrolase effect, can be bactericidal, and have a broader immune function that could facilitate interaction with Plasmodium [54]. Finally, we also saw upregulation and downregulation of serine proteases, which are involved in both digestion and immunity. Serine proteases have been linked to Plasmodium growth and development, invasion of host cells, and can protect the parasite against the host oxidative stress response [55, 56]. Any of these transcriptional changes could be linked to the enhancement of P. gallinaceum, and therefore merit further investigation.

Redox process

One interesting characteristic of Wolbachia infection is the breakdown of redox homeostasis, which has been seen across multiple infected hosts [57, 58]. This effect often manifests through the induction of reactive oxygen species (ROS), and altered expression of genes involved in oxidative stress response, which occur with both native Wolbachia associations and transinfections [10, 12, 58], but it is unclear if they represent an immune response to Wolbachia, or normalization of the redox processes altered by infection [37]. Increased oxidative stress has been linked to pathogen interference against viruses and Plasmodium [10, 11, 20], and may be part of Wolbachia’s immune evasion strategy [37]. Similarly, higher oxidative stress levels are an important part of the anti-Plasmodium immune response [59].

In Ae. fluviatilis we saw evidence of an altered redox response due to wFlu infection, in the form of higher levels of a short-chain dehydrogenase, an oxidoreductase, which can induce ROS, and 3 cytochromes p450, which act as oxidases [60, 61]. Wolbachia have also been linked with iron metabolism and storage, and this can influence key physiological traits such as fecundity [62–64]. wFlu infection induced higher levels of neuferricin, a protein that binds iron-rich heme, which is interesting given that Wolbachia produce enzymes involved in heme synthesis [38], and may utilise it as an energy source [36]. Likewise, heme is essential to the development of Plasmodium in mosquitoes [65]. Given that wFlu appears to alter redox homeostasis, it would be interesting to see if wFlu infection induces higher ROS levels. It should be noted that tetracycline treatment impacts mitochondrial function, which can lead to changes in insect oxidative stress response [66]. While we did use tetracycline to clear the wFlu infection, our experiments were performed more than 2 years (approximately 30 generations) after antibiotic treatment.

Metabolism

We observed that 33 genes associated with metabolism and digestion were differentially expressed as a result of wFlu infection. Wolbachia demonstrate many nutrition-based physiological changes in different insect hosts, including mutualism through nutritional provision [40, 67], better performance under nutritional stress [62], and competition for nutrients leading to effects on host fecundity and fertility, and on pathogen interference [68, 69]. These processes are linked to a variety of nutrients and micronutrients including amino acids, iron, and flavin adenine dinucleotide, which indicates that Wolbachia interact with a wide range of host metabolic pathways.

Wolbachia infections in Ae. aegypti show differential expression of large numbers of digestive enzymes, particularly serine proteases and trypsins, which may indicate cannibalisation of host resources to promote Wolbachia propagation. Infection with wFlu increased the expression of 3 serine proteases, and 2 trypsins, and decreased levels of two other serine proteases and 2 zinc metallopeptidases. While the number of affected genes were fewer than seen for the Ae. aegypti transinfections it is sufficient to suggest that a similar process operates in Ae. fluviatilis. Likewise, wFlu infection elevated levels of brain chitinase and chia, an enzyme involved in breakdown of the peritrophic matrix, which may be further evidence of an altered digestive process, and could compromise the integrity of the peritrophic matrix and potentially facilitate Plasmodium invasion.

We observed that wFlu infection had a broad effect on several aspects of host metabolism, including carbohydrate and lipid metabolism, both areas where Wolbachia is lacking key biosynthesis genes, and where Wolbachia transinfection alters transcription in Ae. aegypti [11, 38, 70]. Previous work indicates that wFlu infection leads to elevated levels of glycogen [71], a major carbohydrate reserve, in developing embryos. While it is unclear if this effect occurs in adults, we did see elevated levels of β-glucosidase, which encodes an enzyme involved in the breakdown of complex carbohydrates including glycogen. Similarly, β-hexosaminidase b levels were also higher, with this enzyme involved in the hydrolysis of amino-sugars. Wolbachia heavily utilise host lipids and sterols, and alter the cellular lipid profile [68, 72, 73], potentially to serve bacterial propagation [74], with similar processes underlying Plasmodium development [75]. Lipid metabolism genes affected by wFlu had a diverse range of functions and include glucosyl glucuronosyl transferases, enzymes required for glucuronidation, which has a role in the metabolism of fatty acids and other compounds, and also pancreatic triacylglycerol lipases, which are used to digest triacylglycerides and process fats. Interestingly, wFlu also increased levels of phosphatidylcholine-sterol acyltransferase, which is involved in sterol and phospholipid metabolism, and may be utilised by Wolbachia to fabricate its membranous vacuoles. wFlu infection also decreased levels of diacylglyceride (DAG) synthesis protein, suggesting that the infection likely interferes with DAG synthesis, and therefore lipid transport as DAGs are prominent transport lipids. Likewise, changes in the expression of these genes could potentially promote P. gallinaceum development.

Infection also induced expression of chorion peroxidase, an enzyme involved in ovarian follicle maturation, and decreased the levels of a vitellogenic carboxypeptidase, which is a yolk protein produced in the fat body [76]. The effect on these genes likely arises due to the strong presence of wFlu in the ovaries, and may contribute to further differences in egg nutritional content, as is seen with glycogen [71].

Membranes

We observed that infection with wFlu also affected 11 transcripts involved in membrane structure, or transmembrane transport. This included an ompa-like protein that is part of the Wolbachia membrane, which was indicative of Wolbachia replication. We also saw the upregulation of membrane proteins involved in the transport of glucose and cholesterol through niemann-pick type c2 protein, and these could potentially be involved in altering Plasmodium infection dynamics, given the varied metabolic requirements of the parasite [77]. This is further evidence of a co-opting of metabolic machinery, as Wolbachia has been strongly linked to carbohydrate metabolism as an energy source [38], and cholesterol, as a critical component of the vacuoles that surround the bacteria [78].

Infection also caused the upregulation of 3 membrane-bound chemosensory receptors potentially linked to Ae. fluviatilis behaviour. These three include sulfakinin, a digestive regulator in insects [79], and 2 taste receptors, including an ionotropic receptor and a gustatory receptor linked to sugar feeding [80]. Wolbachia infection can influence host olfaction at the behavioural [81, 82] and molecular level [11], and in the case of wMelPop, also dramatically alters feeding behaviour [83]. These effects could possibly contribute to a minor change in olfactory response or potentially behaviour in response to wFlu infection. In contrast, only 3 membrane-related genes were differentially expressed in the Tet dataset, including an ATP-dependent transporter molecule of bacterial origin, and two membrane receptors with putative roles in immune signalling, and membrane-protein interactions.

Physiology

We saw that wFlu upregulated several genes with a potential role in neurological function and development. These included calcium and sodium channels, a major allergen, a fasciculation and elongation protein, which is involved in nerve signal transduction, and neuroblast formation protein, and BMP-induced dendrite morphogenesis factor, which are both involved in neural development. While wFlu does not heavily infect many non-reproductive tissues, it is present in the head at low density [15], and so could affect neurological function. Both wMel and wMelPop infect neural tissues [15, 17], and wMelPop has a pronounced effect, causing neurological degradation in Drosophila melanogaster [84]. Wolbachia also affect hormone levels and contribute to host behavioural changes [85–87]. The overall physiological effect of wFlu infection is unclear, however given that genes involved in neurological development, signalling and neurotransmitter trafficking and release were all upregulated by wFlu, there could be critical effects.

While wFlu does infect Ae. fluviatilis heads, it is not found in the omatidia cells of the eye, unlike wMelPop [15], yet wFlu infection decreased the expression of 4 photoreceptor proteins involved in phototransduction, and 3 visual receptors, indicating that wFlu may influence host visual perception. There have been no categorised effects of Wolbachia on host visual process, and no similar genes were affected by wMel or wMelPop [11]. Potential physiological consequences of these changes could be a decreased sensitivity to light, which would be disadvantageous to the mosquito, or may indicate decreased activity during low level light conditions [88–90].

Cellular processes

There was also evidence that wFlu infection altered common cellular processes by influencing the expression of genes involved in DNA repair and replication, and DNA packaging (Tables 2 and 3), with similar genes affected by other Wolbachia infection in mosquitoes - particularly wMelPop [11]. Other contigs were also linked to the processes of transcription and translation, and more specifically to mRNA processing, snRNA and rRNA processing, and protein folding, which is not unexpected given that Wolbachia produces and affects the production of small RNAs [91, 92].

We observed that wFlu altered the expression of chromobox protein homolog 1, a heterochromatin protein that could potentially be involved in epigenetic silencing of gene expression, or in chromosome integrity. Likewise, we observed altered levels of histone 2b, which is involved in DNA packaging as a key component of the nucleosome, and in gene silencing. Wolbachia infection in Ae. aegypti mosquitoes is known to affect cytosine methylation and gene silencing [93], and our results could indicate that there could potentially be a similar influence on silencing through histone modification. In that study, Wolbachia strongly altered the methylation profile of genes involved in membrane function including transport, and here we observed that a similar group of genes were affected during wFlu infection in Ae. fluviatilis.

Reproductive manipulation in Wolbachia-infected insects has also been associated with histones and chromatin. For instance, CI has been linked to delays in the deposition of histones in the male pronucleus [94], while male-killing has been linked to defective chromatin packaging, and altered chromosome behaviour [95]. Critically, the molecular effects of these processes are associated with adult males (CI), or in early stage embryos, and while our data were generated from adult females, these results may suggest that there are broader effects of Wolbachia infection on cellular processes related to DNA packaging, chromatin formation, and the regulation of gene expression in mosquitoes.

Comparison to other Wolbachia infections in mosquitoes

As wFlu is a native Wolbachia infection, and has likely undergone a lengthy period of co-evolution with Ae. fluviatilis, we were quite surprised to discover that infection affected the expression of genes involved in a wide range of host biological processes. The scope of the transcriptional changes associated with infection (257 genes), appears to be similar to that of wMel in Ae. aegypti (327 genes), although that study utilized a microarray rather than RNA-Seq, and the difference in sensitivity of the techniques may be a factor [11]. This similarity is interesting given that the latter association is both relatively novel, and is a transinfection that produces a broader range of physiological effects, and has greater bacterial density [15, 17, 32]. In contrast, transinfection of Ae. aegypti with the more virulent wMelPop strain has been demonstrated to alter the expression of 244 genes in one study [14], and 2723 in another [11], with the studies utilizing different methods of analysis. A larger set of differentially expressed genes would likely reflect the pathogenicity associated with wMelPop infection.

While these wMel and wMelPop data come from transinfections in a different species, we sought to determine if there was a core set of genes or functional groups of genes that were differentially expressed in all three associations, as this might explain something fundamental about the nature of Wolbachia infections in mosquitoes. 75 of the upregulated contigs, and 54 of the downregulated contigs were homologous to previously described Ae. aegypti transcripts in VectorBase. However, only 42 of these were also differentially expressed during either wMel or wMelPop infection. Interestingly, the majority of these were upregulated by wMel and wMelPop, even if they were downregulated by wFlu (Additional file 5), with this difference likely due to the relative novelty of the former transinfections.

Mosquito material

Two mosquito lines were used in these experiments. The original, Wolbachia-infected Aedes fluviatilis (wFlu) line was derived from mosquitoes originally captured in Belo Horizonte in 1975 [32]. The mosquito colony was maintained in the laboratory until 2013 when a subset (Tet) was treated with tetracycline hydrochloride to remove the native Wolbachia infection, and then had their gut microbiota recolonized, as previously described [32]. Colony larvae were reared at low density in dechlorinated water, and were fed with fish food (Goldfish Colour, Alcon, Camboriú, Santa Catarina, Cat. No. 0504-2). Adults were maintained in low-density cages in a climate-controlled insectary (temperature: 27 ± 1 °C, RH: 70 ± 10%, photoperiod: 12 h light: 12 h dark), and provided 10% sucrose solution ad libitum. Mosquitoes used in experiments were maintained in small cylindrical cages (diameter – 16 cm, height – 18 cm) of approximately 80–90 individuals. Experiments were conducted more than 2 years after tetracycline treatment. In all experiments, the Tet line served as a Wolbachia-uninfected control line for the wFlu line, to study the effects of Wolbachia infection.

RNA extractions & library preparation

Six groups of samples were prepared for sequencing, with each group going on to form an independent library. Three groups each of 16 6-day old whole adult females from the Flu and Tet lines were collected and total RNA extracted using the Trizol® protocol according to manufacturer’s instructions (Invitrogen), for a total of 6 independent samples, with 3 biological replicates per treatment. Mosquitoes were fed only 10% sucrose prior to collection. RNA levels in each sample were quantified using a NanoDrop ND1000 (ThermoFisher Scientific). Sample degradation levels were checked by running a portion of the samples on a standard non-denaturing agarose gel containing bleach [104].

cDNA libraries were constructed using Illumina Truseq RNA Sample Preparation Kit (Illumina Inc.) according to the manufacturer instructions starting with 4 μg of total RNA. The library products were then sequenced using an Illumina MiSeq platform on a paired-end 300 bp run. After cleaning reads from adaptor sequences, the quality of the reads was assessed using the FastQC program (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Cleaned reads are available for download at the National Center for Biotechnology Information - Sequence Read Archive under the BioProject ID PRJNA320882.

De novo transcriptome assembly and contig annotation

Since no reference genome was available for Aedes fluviatilis, a de novo transcriptome assembly was built with Trinity (https://github.com/trinityrnaseq/trinityrnaseq/wiki) using default parameters. All six Illumina RNAseq datasets were combined in order to assemble a more reliable transcriptome, with a total of 19,919,299 paired-end, high quality reads. Contig sequences were searched for candidate proteins with TransDecoder (https://transdecoder.github.io/), again using standard parameters. The assembled contigs were annotated through local alignments with BlastX (http://blast.ncbi.nlm.nih.gov/Blast.cgi) to the NCBI non-redundant (NR) and KEGG databases. BlastX parameters were set with an -e value of 1e-10. Blast2GO (https://www.blast2go.com/) was used to retrieve Gene Ontologies to annotated transcripts. Phylogenetic and divergence analyses were conducted using sequence data obtained during this study, or from UniProt (www.uniprot.org) or VectorBase (https://www.vectorbase.org). Methods and references for these analyses are described in Additional file 1.

Read mapping and differential gene expression

All Illumina paired-end reads libraries were mapped separately against the Ae. aegypti predicted transcriptome, available at VectorBase, and the Trinity assembled contigs, both with Bowtie2 (http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) using the default parameters while configuring fragment length. The Integrative Genomics Viewer (http://www.broadinstitute.org/igv/) was used to visualize the reads that were mapped back to the assembled transcriptome. Read counts mapped to each transcript were acquired with a custom Pearl script (available upon request).

Two R packages from Bioconductor (www.bioconductor.org), baySeq and DESeq2, were selected in order to identify the contigs that were significantly differentially expressed due to the presence of Wolbachia. The baySeq method is more sensitive, but also carries a greater false positive call rate, hence differential expression was confirmed using DESeq2 [45, 105]. The baySeq FDR (false discovery rate) P value for multiple tests was set to 0.05. The DESeq2 adjusted P value was set to 0.01, and transcripts were filtered according to their log2FoldChange (higher than +1 or lower than -1). We then obtained an estimate of the average fold change for each of the differentially expressed contigs using FPKM (fragments per kilobase of exon per million fragments mapped). Contigs that were differentially expressed but had no BlastX hit were excluded from analysis.

GO terms for the remaining contigs were generated through Blast2GO (https://www.blast2go.com) and VectorBase, or, when data could not be discovered from these sources, through FlyBase (http://flybase.org) or InterPro (http://www.ebi.ac.uk/interpro/), based on the analysis of homologous genes. GO information from contigs that matched to the same BlastX hit was considered only once, to avoid double counting. These GO lists were compiled and used to create two lists, one for the Flu libraries with the native wFlu Wolbachia infection, and the other for the Tet libraries, where the Wolbachia had been removed. Contigs were grouped according to putative function, and these lists were then compared to develop a profile of the transcriptomic effects of wFlu on its mosquito host.

Differentially expressed contig functions and GOs were then compared against the transcriptomic data from wMel and wMelPop-infected Aedes aegypti arrays [11, 14], in order to develop a profile of types of transcriptional changes held in common between the three Wolbachia strains and their mosquito hosts. Contigs that had BlastX matches against Ae. aegypti were compared directly with Ae. aegypti expression data. Those contigs that matched to another species were compared directly against the Ae. aegypti genome in NCBI using BlastN (http://blast.ncbi.nlm.nih.gov/Blast.cgi). Any hits with a substantial match percentage (>80%) and a significant e-value were used for further analysis, using the first hit in a comparison in VectorBase while those that did not were not considered. This information was used to determine if any of the specific genes affected by wFlu infection were also affected by wMel or wMelPop infection.

Confirmation of differential expression

To assess the accuracy of the transcriptomic data set, six contigs were selected at random for expression analysis with RT-qPCR. Two of these contigs were indicated to have higher expression for Flu mosquitoes (AF10645; comp10645_c1_seq2 and AF10453; comp10453_c0_seq5), two had higher expression in Tet mosquitoes (AF15178; comp15178_c0_seq1 and AF14155; comp14155_c0_seq1), while two had equivalent expression between Flu and Tet mosquitoes (AF2025; comp2025_c0_seq1 and AF2041; comp2041_c0_seq1). The first four of these contigs were predicted to be differentially expressed through both the bayseq and DESeq2 analyses. Primers for these contigs were designed from the sequences generated during sequencing, which meant that they were suitable for cDNA. Primer sequences were designed using Primer 3 V0.4.0 (http://bioinfo.ut.ee/primer3-0.4.0/) to have a Tm of 55–60 °C and a product size range of 80–120 bp (Additional file 6).

Mosquitoes from both the Flu and Tet lines were reared as described above, and then females were collected individually at 6 days post eclosion. Total RNA was extracted and quantified as described above. First strand cDNA synthesis was conducted with 1 μg of total RNA from each sample using the M-MLV reverse transcriptase assay according to manufacturer’s instructions (Promega cat: C118A). cDNAs were then diluted 1:10 in nuclease-free water and stored at -30 °C. SYBR-based PCR was used to confirm the expression of each of the test genes, with 15 samples tested per mosquito line. Each gene was quantified in duplicate for all samples using the following mix: SYBR -5 μL, forward and reverse primers (10 μM) - 0.5 μL each, sterile RNase free water - 2 μL, sample 2 μL). RT-qPCR for samples was run on a LightCycler® 96 System (Roche) using the following profile: 10 min pre-incubation at 95 °C, 40 cycles of 15 s at 95 °C, 60 s at 60 °C, melt curve - 95 °C for 15 s, ramp from 60 °C to 95 °C at 1.6 °C/s. Expression values for each gene were normalised against actin1 expression, which had previously been demonstrated to be a good control gene for Ae. fluviatilis [32]. Mean normalised expression values for each gene were calculated using Q-Gene [106], and were compared statistically between Flu and Tet mosquitoes using Mann Whiney U tests (Prism v 6.0 g, Graphpad).