- Research article
- Open Access
Highly focused transcriptional response of Anopheles coluzzii to O’nyong nyong arbovirus during the primary midgut infection
BMC Genomicsvolume 19, Article number: 526 (2018)
Anopheles mosquitoes are efficient vectors of human malaria, but it is unknown why they do not transmit viruses as well as Aedes and Culex mosquitoes. The only arbovirus known to be consistently transmitted by Anopheles mosquitoes is O’nyong nyong virus (ONNV, genus Alphavirus, family Togaviridae). The interaction of Anopheles mosquitoes with RNA viruses has been relatively unexamined.
We transcriptionally profiled the African malaria vector, Anopheles coluzzii, infected with ONNV. Mosquitoes were fed on an infectious bloodmeal and were analyzed by Illumina RNAseq at 3 days post-bloodmeal during the primary virus infection of the midgut epithelium, before systemic dissemination. Virus infection triggers transcriptional regulation of just 30 host candidate genes. Most of the regulated candidate genes are novel, without known function. Of the known genes, a significant cluster includes candidates with predicted involvement in carbohydrate metabolism. Two candidate genes encoding leucine-rich repeat immune (LRIM) factors point to possible involvement of immune protein complexes in the mosquito antiviral response. The primary ONNV infection by bloodmeal shares little transcriptional response in common with ONNV infection by intrathoracic injection, nor with midgut infection by the malaria parasites, Plasmodium falciparum or P. berghei. Profiling of A. coluzzii microRNA (miRNA) identified 118 known miRNAs and 182 potential novel miRNA candidates, with just one miRNA regulated by ONNV infection. This miRNA was not regulated by other previously reported treatments, and may be virus specific. Coexpression analysis of miRNA abundance and messenger RNA expression revealed discrete clusters of genes regulated by Imd and JAK/STAT, immune signaling pathways that are protective against ONNV in the primary infection.
ONNV infection of the A. coluzzii midgut triggers a remarkably limited gene regulation program of mostly novel candidate genes, which likely includes host genes deployed for antiviral defense, as well as genes manipulated by the virus to facilitate infection. Functional dissection of the ONNV-response candidate genes is expected to generate novel insight into the mechanisms of virus-vector interaction.
O’nyong-nyong virus (ONNV) and Chikungunya virus (CHIKV) are closely related alphaviruses of the Semliki Forest virus complex . Both viruses cause febrile illness in humans and the symptoms are hard to distinguish from other alphavirus infections, Dengue fever or malaria . Therefore, the ONNV burden on human populations has likely been underestimated due to misdiagnosis.
Indeed, a recent serodiagnostic study showed that in coastal Kenya, both CHIKV and ONNV circulate at high levels in human populations , with a high possibility of co-infections. Anopheles mosquitoes transmit ONNV, while CHIKV is transmitted by Aedes mosquitoes. More research effort has been focused on Aedes mosquitoes, due to their ability to also transmit multiple arboviruses such as dengue or Zika (Flaviviruses), which has led to a relative neglect about Anopheles interactions with arboviruses.
Virus survival of the bloodmeal digestive environment and initial infection of the midgut epithelium is thought to be the first infection bottleneck, called the midgut infection barrier (MIB) [4, 5]. The MIB is probably a distinct step from virus escape into the systemic compartment (the midgut escape barrier, MEB), followed by dissemination to other tissues including the salivary glands. Consistent with that prediction, we previously identified distinct antiviral immune signaling responses of the Anopheles primary midgut infection by bloodmeal, as compared to the subsequent systemic disseminated infection . Earlier studies in which the virus was introduced by intra-thoracic injection [7,8,9,10] modeled only the disseminated systemic infection, even if it includes infection of midgut cells by dissemination, and do not measure the primary antiviral responses of the midgut.
To deepen the understanding of Anopheles interactions with RNA viruses, we sequenced messenger RNAs and small RNAs of A. coluzzii midguts three days after oral infection with ONNV, a time when the virus has not yet escaped to the systemic compartment . Analysis of the datasets identified differentially regulated transcripts for messenger RNA (mRNA) and microRNA (miRNA) during ONNV infection of A. coluzzii, and predicted novel potential miRNAs. Transcriptional responses of the primary midgut ONNV infection were compared to published transcriptional datasets of A. coluzzii infected with ONNV by needle injection, and to the primary midgut infection with the eukaryotic protozoan malaria parasite, Plasmodium. In addition, joint analysis of mRNA and miRNA expression levels predicted putative regulatory networks of immune signaling.
Differential A. coluzzii mRNA transcript abundance during ONNV infection
The primary infection of ONNV in the mosquito midgut after an infectious bloodmeal is controlled by antiviral mechanisms largely distinct from those active in the subsequent disseminated systemic infection . To identify candidate antiviral factors and immune pathways solicited during the primary midgut infection, we fed mosquitoes on an ONNV-containing bloodmeal or control normal bloodmeal, purified RNA three days post-bloodmeal before virus dissemination to the systemic compartment, and transcriptionally profiled mRNA abundance by Illumina RNA sequencing (RNAseq). Virus bloodmeal titer was optimized to produce midgut infection rates of 100% among bloodfed mosquitoes, as previously determined , which augments experimental and statistical power. Although ONNV infection is limited to the midgut at 3 days (d) post-bloodmeal, we purified RNA from entire mosquitoes in order to also detect any potential systemic signaling prior to virus escape from the midgut.
Only 30 candidate genes were significantly regulated during ONNV primary infection of A. coluzzii as compared to mosquitoes fed a normal bloodmeal, generating a discrete and highly focused transcriptome response (Fig. 1; Additional file 1: Table S1, column F: ONNV log2FoldChange, column G: ONNV adjusted p-value). Of the 30 regulated candidates, just six have assigned gene names. Thus, the majority of ONNV-regulated genes carry little or no functional information, and are difficult to interpret. Functional analysis using Gene Ontology (GO) terms detected significant enrichment for predicted carbohydrate metabolic processes (GO:0005975, carbohydrate metabolic process, 201 annotated, 0.67 expected, 6 observed, p = 3.06e-05; Additional file 1: Table S1, column H: Associated GO-terms), suggesting that ONNV replication and/or the host immune response may impose metabolic costs on the mosquito host. The ONNV-regulated candidate genes include known or predicted immune-related genes LRIM4, LRIM10, CLIPB8, FREP63 and AGAP000376, although the different GO terms associated with immunity were not significantly enriched. The latter gene, AGAP000376, displays orthology with Drosophila melanogaster tsf1, a midgut immune factor in the Imd immune pathway (Additional file 2: Figure S1) [11,12,13].
Comparison with other Anopheles transcriptional datasets
The current Illumina RNAseq dataset was compared to previously reported microarray gene expression data, also from A. coluzzii, infected by intra-thoracic injection with ONNV . In that study, 211 transcripts were differentially regulated at four days post-injection. Only three genes overlap with the transcripts regulated by the normal bloodmeal-infection route used in our study (Fig. 2): LRIM4 (AGAP007039) and AAPP (AGAP009974) are upregulated both in blood-infected midgut and needle-infected mosquitoes at one day post-injection, while LRIM10 is downregulated in blood-infected midgut but was upregulated in needle-infected mosquitoes at four days post-injection.
The small number of overlapping regulated genes between the studies, and the opposite regulation of LRIM10, is consistent with the dichotomy in immune signaling previously described for the two compartments, that is, the primary midgut infection by bloodmeal, and the secondary disseminated systemic infection . Needle infection is a model for the secondary disseminated viral infection, bypassing the primary midgut infection, including the processes of the MIB and MEB. After needle-infection, most cell types of the mosquito are infected rapidly, which could explain the large differential transcript set, while infection of the midgut epithelium by bloodmeal is kinetically different because of the retarding effect of the midgut infection barrier, and may even produce systemic signals that induce qualitatively different responses in the secondary disseminated infection . However, an influence of technical differences between the studies for infection methods (bloodmeal without wounding, or needle infection with wounding), and transcript profiling platforms (RNAseq or microarray) also cannot be excluded.
We compared differential regulation induced by ONNV infection of the midgut with previously reported expression data from A. coluzzii infected with the midgut stage of Plasmodium falciparum . Unfortunately, a comprehensive comparison is not possible, because the Plasmodium infection data are from 2006, and used a now-obsolete reference genome annotation that contained different gene models. Many gene IDs from the annotation used in the Plasmodium study were retired, some gene IDs were split, and many new gene IDs were created since then. Technical differences would also limit the value of completely remapping the 2006 microarray probe sequences to current candidate gene IDs, because the microarray coverage of transcripts was much lower than the current RNAseq coverage. Nevertheless, we remapped the recognizable gene ID matches to the current A. gambiae PEST strain genome (assembly AgamP4) in order to compare the two biological datasets to the extent possible.
Among the 30 current ONNV-regulated transcripts, a total of 12 gene IDs could be linked to old gene IDs from the Plasmodium dataset (Additional File 1: Table S1, column U: ID_Old-ENSANGT). All 12 shared transcripts were induced by ONNV infection (Additional file 3: Figure S2; Additional file 1: Table S1, column V: ONNV, 1 = induced, 0 = no change, − 1 = repressed). Seven of the transcripts were also induced by mosquito feeding on a viable but invasion-defective P. falciparum mutant (Additional file 1: Table S1, column Y: PfCTRP), while only four and three genes were regulated by midgut infection with P. falciparum and P. berghei, respectively (Additional file 1: Table S1, column W, PfGUT, column X, PbGUT). There is no clear pattern of gene function among the shared genes, but this is not surprising given the small numbers as well as the low level of gene functional annotation. A comprehensive comparison would require repeating the Plasmodium study using RNAseq analysis mapped to the current reference genome and gene IDs.
Differential miRNA abundance during ONNV infection
The small RNA fraction was sequenced from the same above A. coluzzii samples fed on an ONNV-containing or normal bloodmeal, and was analyzed using miRDeep2. 118 known miRNAs were detected, and 182 potential novel miRNAs were predicted in A. coluzzii (Additional file 4: Table S2). These figures are similar to the total number of 256 known miRNAs in Drosophila (miRbase BDGP5.0).
Differential analysis of miRNA abundance after ONNV infection detected one miRNA (aga-chr2L_23370) that was regulated during the primary infection of A. coluzzii with ONNV (Additional file 5: Figure S3, Additional file 4: Table S2). We observe a log2 fold change of 1.949 of miRNA aga-chr2L_23370 in ONNV-infected mosquitoes as compared to normal bloodfed controls (adjusted p = 2.52e-04). This miRNA was identified in a recent catalog of Anopheles miRNAs . That analysis identified four miRNAs significantly regulated by bloodmeal, and six by P. berghei-infected blood, but levels of aga-chr2L_23370 were not altered by any condition tested. These two sets of results suggest that neither viral nor Plasmodium infection exerts a large impact on Anopheles miRNA expression. However, one miRNA displayed a functional effect on the efficiency of Plasmodium development in Anopheles [17, 18].
Thus, to date ONNV infection is the only condition known to regulate the expression level of miRNA aga-chr2L_23370. The potential functional effect of miRNA aga-chr2L_23370 on the host response to ONNV infection, and viral replication, will require further work to understand.
Correlation of miRNAs with regulation of mRNA expression
Sequence-based target prediction methods of miRNAs are based on the prediction of complementarity between the seed region of mature miRNAs (nucleotide positions 2–7) and the 3’ UTR region of messenger RNAs [19, 20]. A drawback of these methods is a high rate of false positive predictions due to the small size of the seed region of mature miRNAs . Consequently, it has been recommended to filter these sequence based predictions by using the joint analysis of expression levels of miRNAs and mRNA targets .
In order to generate a large dataset to empirically test the correlation between miRNA expression and regulation of mRNAs, we silenced the key immune signaling nodes, Rel2 and Stat-A, with and without ONNV infection, to perturb large networks of gene expression in the signaling pathways Imd and JAK/STAT, respectively. These two immune pathways were chosen because we previously showed that their activity is required for full A. coluzzii antiviral protection from initial midgut ONNV infection . Perturbing these pathways caused differential expression of a total of 11 miRNAs, with expression levels of five miRNAs altered by silencing of Rel2, four by silencing of Stat-A, and two by control treatment with dsRNA for LacZ (Additional file 4: Table S2). The behavior of these regulated miRNAs clustered according to the different experimental treatments (Additional file 6: Figure S4). Consistent with the results above indicating a small and focused ONNV-dependent regulation of abundance of mRNAs (Fig. 1) and miRNAs (Additional file 5: Figure S3, Additional file 4: Table S2), the hierarchical clustering of miRNA expression response was driven largely by the pathway silencing treatments, with little effect due to the presence or absence of ONNV infection.
Here, we use deep sequencing of mRNAs and miRNAs in the same samples in order to empirically test the accuracy of sequence-based miRNA target predictions, and to ascertain candidate miRNA-mRNA pairs potentially involved in A. coluzzii antiviral immunity. Co-expression levels were determined for miRNA-mRNA target pairs, under the general mechanistic model assuming that miRNAs repress the expression of mRNA from their target genes . When miRNA-mRNA target pairs from sequence-based predictions were filtered for those having a significant negative correlation for expression, the number of potential target genes was 336, with just one gene (AGAP001650) potentially targeted by both Rel2 and Stat-A-induced miRNAs (Fig. 3a, Additional file 7: Table S3). Principal component analysis (PCA) of the Pearson correlation matrix for the 11 miRNAs and the 336 filtered candidate mRNA targets demonstrates a clear clustering of target genes of miRNAs induced under different treatments (Fig. 3b).
We further explored the discovery of candidate target genes for regulation by specific miRNAs by imposing only a significance threshold for anti-correlation between expression levels of miRNA-mRNA pairs, without in silico target site prediction. When filtered strictly by significant anti-correlation, the number of potential target genes of the 11 miRNAs increased more than 10-fold to 4614 mRNA candidates induced among the different treatment conditions. This result is consistent with observations that computational miRNA target site prediction can have limited accuracy , and in particular can under-predict functional miRNA-mRNA interactions . In the current case, filtering only for anti-correlation of abundance, clear clustering of target genes according to treatment is again observed (Fig. 4a, Additional file 7: Table S3), but now with greater density of candidate target genes within the clusters (Fig. 4b).
These results strongly suggest that nine miRNAs are involved in inhibiting specific clusters of target genes of the Imd and JAK/STAT immune pathways previously implicated in A. coluzzii midgut antiviral protection to ONNV . In addition, another cluster of genes responding to two miRNAs induced after treatment with control LacZ dsRNA appear to be involved in wounding and/or exposure to dsRNA. Moreover, the results suggest that computational target site prediction may under-predict target genes in this system, although it cannot be ruled out that the additional fraction of anti-correlated mRNAs not predicted as target sites may also include genes indirectly regulated by direct miRNA targets.
Here we show that ONNV infection of Anopheles, the natural vector of this virus, causes a highly focused and limited impact on the host transcriptome during the primary midgut blood-induced infection. The most striking observation is that just 30 host candidate genes and one miRNA were regulated by ONNV infection when controlled for the effect of the bloodmeal alone. The midgut epithelial cells are the first point of contact between host and pathogen at the portal of initial infection. The host genes regulated during this initial interaction are likely to be crucial for the subsequent infection, probably including both host genes deployed for antiviral defense as well as genes manipulated by the virus to facilitate infection. A majority of the ONNV-regulated candidate genes have unknown molecular function. Little is known about Anopheles interaction with viruses, and functional dissection of these candidates would likely generate novel biological insights into the host-virus interaction, and the structure of Anopheles immunity.
Two main categories of observed regulated candidate genes are associated with carbohydrate metabolism and immunity. A simple interpretation of the carbohydrate metabolism category could be that viral infection of the midgut epithelium imposes a metabolic burden on the infected cells because of the demands of viral replication, simultaneous with the high energetic demands of bloodmeal digestion. Induction of immune processes probably limits viral replication and consequent cytopathic effects, but also imposes the metabolic costs of immunity. The net host response probably balances the fitness costs of viral replication with the costs and benefits of digestion and immunity. Converting the bloodmeal to egg production is probably prioritized by the host because it is directly linked to host reproductive fitness. Conversely, viral fitness is linked not only to replication but ultimately to transmission, which requires a vector healthy enough to fly and bite a vertebrate host. This need should modulate overt arbovirus pathogenicity to the mosquito host.
Mosquitoes are persistently exposed in nature to members of the natural virome, including insect-specific viruses (ISVs) [26,27,28,29], which have probably been an important evolutionary force shaping mosquito antiviral immunity. In general, prevalent members of the natural virome should evolve towards lower pathogenicity in their habitual host lineages. Consistent with that expectation, a symbiotic Anopheles densovirus (studied for its potential use as a paratransgenesis tool) elicited low transcriptional response upon infection, and neither mosquito age at infection nor feeding regime affected viral titers [30, 31]. However, maintenance of bacterial microbiome commensals in the non-pathogenic commensal state requires active policing by basal host immunity . By analogy, the maintenance of adapted ISVs as non-pathogenic may also result from a dialog with host immunity. Presumably, the same antiviral mechanisms used in basal maintenance are also deployed against arboviruses when encountered, which are often in the same families as members of the insect virome .
There was only one previous study of gene expression in Anopheles infected with ONNV, in our knowledge, which displayed little similarity of gene regulation to the current study . The difference between the study results can likely be explained by the experimental designs. The previous study infected the hemocoel by injection of ONNV, thus bypassing the midgut epithelium and physiological influence of virus delivery in a bloodmeal. The Anopheles antiviral response to ONNV is highly compartmentalized, and involves the activity of distinct signaling pathways and molecules in the primary midgut infection as compared to the secondary disseminated infection .
There have been several transcriptomic studies of the Aedes aegypti mRNA response to virus infection. The results of one study  are not comparable to ours because tissues were analyzed 10 d after an infective bloodmeal, which detected mainly the secondary systemic dissemination into all tissues including midgut. Another study used needle infection of viruses for transcriptome discovery . A recent study delivered chikungunya virus in an oral saline meal rather than bloodmeal, and detected 78 differentially regulated transcripts at 2 d post-infection .
Two studies with comparable experimental design to detect the primary midgut infection response were reported in Ae. aegypti. One study after a dengue virus oral bloodmeal detected 30–89 differential transcripts at, respectively, 1 d and 4 d post-infection . More recently, analysis of the transcriptomic response of Ae. aegypti miRNAs  and mRNAs  after Zika virus infection detected 17 differentially regulated miRNAs, and 54 or 101 differentially expressed mRNAs at 2 d and 7d post-infection, respectively. The results of these two studies in Ae. aegypti both appear broadly similar to the current Anopheles and ONNV study in the magnitude of mRNA response to viral infection, while the Ae. aegypti miRNA response may appear more robust than the miRNA regulation in Anopheles.
A ternary immune complex protects A. coluzzii against Plasmodium infection, and is comprised of the leucine-rich repeat (LRR) immune factor LRIM1, either APL1A or APL1C, and a complement cofactor Tep1, Tep3 or Tep4 [39,40,41]. We previously found that activities of APL1A and APL1C are required for antiviral protection against ONNV midgut infection in A. coluzzii, but that the cofactors in the anti-Plasmodium response, LRIM1, Tep1, and Tep3, are not involved . In the current study, we observe that two LRR proteins are regulated during ONNV infection. LRIM4 is upregulated whereas LRIM10 is downregulated, suggesting that LRIM4 could be a partner in a putative antiviral ternary complex in conjunction with APL1A or APL1C, and an unknown Tep factor.
We also identify the clip-domain serine protease CLIPB8 as a downregulated candidate gene during ONNV infection. CLIP proteins modulate diverse immune pathways in insects. CLIPB8 is a component of the Anopheles pro-phenoloxidase activation system , and is required for melanization of P. berghei parasites and Sephadex beads [43, 44]. The repression of CLIPB8 during ONNV infection raises the possibility that either the protective response of Anopheles to ONNV infection requires an inhibition of the melanization cascade, or that the virus suppresses the cascade as host manipulation. That possibility would be consistent with the observed decrease of melanization of P. berghei parasites in mosquitoes co-infected with ONNV .
The extensive small RNA sequence data in this study produced a comprehensive database of known and potential miRNAs in A. coluzzii. We identified one miRNA, aga-chr2L_23370, which is differentially regulated by ONNV expression. Characterizing this miRNA and its target genes could shed light on the RNA virus-Anopheles interaction. No previous miRNA studies have been performed on the Anopheles viral response to our knowledge. Multiple studies in Aedes or Culex, often cell lines, have identified miRNAs that respond to arbovirus infection, and can participate in antiviral protection [37, 45,46,47,48,49,50,51]. A comparable study in Ae. aegypti blood-infected with Zika virus detected ten miRNAs significantly regulated at 2 d post-bloodmeal as compared to normal bloodmeal .
Interestingly, nine A. coluzzii miRNAs were differentially regulated by silencing of the Imd and JAK/STAT immune signaling pathways (with two more miRNAs regulated by the dsLacZ control), suggesting an influence by miRNA expression upon mosquito basal immunity. Anti-correlation of these 11 miRNAs with transcriptome-wide mRNA abundance, either with or without miRNA target site prediction, identified defined clusters of regulated genes that plausibly respond to miRNA control. These results provide novel candidate miRNAs and target gene networks for the study of Anopheles immune and antiviral regulation.
O’nyong nyong virus (ONNV) is the only virus known to be consistently transmitted by Anopheles mosquitoes. Anopheles are efficient vectors of human malaria, but it is unknown why they do not transmit viruses as efficiently as Aedes and Culex mosquitoes. Here, we identified the candidate genes that respond to ONNV infection in the mosquito stomach after virus infection. The cells of the stomach are the first site of antiviral defense, and early events there are likely to be crucial for the course of the subsequent infection. We found that the expression levels of only 30 candidate genes are changed by primary infection with ONNV. The function of most of the genes is unknown, but they probably include a combination of mosquito self-defense functions, and also mosquito genes that are manipulated by the virus in order to promote infection. Nevertheless, the vector must remain fit enough to fly in order to transmit the virus, suggesting that the virus may prioritize viral stealth and immune evasion rather than pathogenicity. Examination of the ONNV-response candidate genes to determine their functions is expected to generate novel insight into the mechanisms of virus-vector interaction.
The A. coluzzii Ngousso strain was originally initiated with mosquitoes collected in Cameroon in January 2006, and belongs to the M molecular and Forest chromosomal forms . Mosquitoes were reared under standard conditions at 26 °C and 80% relative humidity, with a 12 h light/dark cycle as previously described .
Viruses and mosquito infection
The ONN-eGFP infectious clone of the ONNV SG650 strain  was kindly provided by Dr. Brian Foy, Colorado State University. ONNV SG650 was first isolated from human serum in Uganda in 1996. Infectious clone cDNA was linearized by NotI, and viral RNA was transcribed using T7 RNA polymerase (NEB), m7(5′)ppp(5′)A-Cap Analog (NEB) and rNTPs (Promega). Transcribed RNA was electroporated into BHK-21 cells. Virus was recovered after 72 h, and passaged once on A. gambiae 4a3a cells [53, 54] to increase titers.
Mosquitoes were infected by feeding on an infectious bloodmeal containing viral titers at the high end as found at the peak of clinical viremia, as previously described . These conditions ensure that all mosquitoes are positive for virus infection at 3 d post-bloodmeal, which augments experimental and statistical power. Briefly, female mosquitoes deprived of sugar for 12 h were allowed to feed for 15 min through a Hemotek membrane covering a glass feeder containing the blood-virus mixture maintained at 37 °C. The infectious bloodmeal was composed of a virus suspension diluted (1:3) in washed rabbit erythrocytes from arterial blood, resuspended at 50% in dialyzed rabbit serum (Sigma R4505). ATP was added to a final concentration of 5 μM. Fully engorged females were transferred to small cardboard containers and maintained with free access to 10% sucrose at 28 ± 1 °C. Final bloodmeal titer fed to mosquitoes was between 1 to 3x107pfu/ml as assessed by plaque assay.
Double-stranded RNA synthesis and injection
Double-stranded RNAs were synthesized from PCR amplicons using the T7 Megascript Kit (Ambion) as described . The sequences of primers used for synthesis of dsRNA templates are in Additional file 8: Table S4. For each targeted gene, 500 ng of dsRNA (up to 69.9 nl, depending on the concentration) was injected using glass capillary needles into the thorax of cold-anesthetized 1–2 d-old A. coluzzii females using a nano-injector (Nanoject II, Drummond Scientific). dsLacZ was used as control.
The efficiency of the gene silencing effect was monitored 2–3 d after dsRNA injection (at the time of infection) on a pool of 5 whole mosquitoes. One step RT-qPCR (Power SYBR Green RNA-to-Ct 1 Step Kit, Applied Biosystems) was performed following the manufacturer’s instructions. The sequences of primers used for gene knockdown verification are in Additional file 8: Table S4, the annealing step was performed at 60 °C.
RNAseq and small RNA sequencing
12 mosquitoes per pool, in each of two biological replicate pools, were collected 3 d after infection by oral feeding with ONNV and homogenized in TriReagent (Sigma). The long RNA fraction (> 200 nt) and small RNA fractions (< 200 nt) were extracted from the same biological pool using Nucleospin miRNA (Macherey Nagel) using the Trizol protocol and submitted to a Bioanalyser (Agilent) for quality assessment.
For RNAseq, 10μg of RNA (> 200 nt fraction) was used to prepare directional libraries using the TruSeq Stranded mRNA Library sample prep kit with index sets A and B (Illumina), with the following protocol modifications. Chemical fragmentation of polyA RNA was done using Ambion reagent (AM8740), followed by purification on RNeasy columns (Qiagen, #74204). Phosphatase and PNK treatments were performed and followed by purification on RNeasy columns (Qiagen, #74204). The fragmented RNA was then ligated with 3′- and 5′- TruSeq adapters, as described in the original protocol. Synthesis of cDNA was performed by reverse transcription. The cDNA product was then specifically amplified by 11 cycles of PCR and products purified on Agencourt AMPure XP beads (Beckman Coulter Genomics, # A63881).
For small RNA sequencing, the 10- to 30-nt RNA fragments were purified from the small RNA fraction (< 200 nt) on a 15% urea PAGE (BioRad). Small RNA libraries were then prepared using TruSeq Small RNA Library reagents (Illumina). The purified small RNAs were ligated with 3’ RNA adaptor (RA3 5’-TGGAATTCTCGGGTGCCAAGG). The 5’ RNA adaptor (RA5 5’-GUUCAGAGUUCUACAGUCCGACGAUC) was then ligated to the 5′ phosphate end of the small RNAs. Reverse transcription was done to convert RNA to cDNA, which was then selectively enriched by 11 cycles of PCR. The PCR products were purified on 5% PAGE (BioRad) and checked on a Bioanalyser DNA1000 chip (Agilent).
The resulting libraries were quality-controlled on a Bioanalyser DNA1000 chip (Agilent). Libraries were sequenced using the Illumina Hiseq 2000 in a multiplexed 51 + 7 bases single read using a TruSeq SR cluster kit v3 cBot HS, and a TruSeq SBS kit v3 HS 50 cycles (Illumina). Primary analysis of the sequences was performed with Casava software (v1.7 Illumina). Library preparation and sequencing was performed by the Transcriptomics and Epigenomics Core Facility of Institut Pasteur.
Analysis of RNAseq sequencing data
Raw fastq files were clipped of their adaptor and filtered for quality using Cutadapt . Sequence reads were then mapped against the A. gambiae PEST strain genome (assembly AgamP4) using BWA mem version 0.7.7 with default parameters . Read counting in genes/features was done with HTSeq-count version 0.6.1  and the AgamP4.3 gene set. Differential gene expression analysis was done under R version 3.3.1  with the package DESeq2 version 1.14.1  The raw p-values were adjusted for multiple testing according to the Benjamini and Hochberg procedure  and only the genes with an adjusted p-values under 0.05 were considered as differentially expressed. Venn Diagrams were drawn using the R package VennDiagram version 1.6.0, and GO term analysis was performed using the R package topGO .
Analysis of small RNA sequencing data
Raw fastq files were clipped of their adaptor and filtered for quality with Cutadapt . Next, a database of miRNA was build using reference miRNA from miRBase  and from the literature [16, 63]. For this purpose, genomic coordinates of reference miRNAs in gff format were integrated with bedtools v2.26  in order to generate a non-redundant set of 167 A. coluzzii miRNA genes. Genomic coordinates of miRNA genes are based on A. gambiae PEST strain genome (assembly AgamP3).
The miRNA database was complemented using our small RNA datasets to predict potential novel miRNA using miRdeep2 version 188.8.131.52  (miRNA_ONY_pre.fasta + mature), which yielded a prediction of 300 miRNAs, of which 118 were reference miRNAs and 182 were novel miRNAs detected in the present study (Additional file 4: Table S2, miRNA databases).
The quantifier.pl module of miRdeep2 was used to map filtered and clipped small RNA sequence reads over reference and novel miRNAs, quantifying the number of reads falling into an interval 2 nt upstream and 5 nt downstream of the mature/star sequence. Raw counts were analyzed for differential expression using DEseq 2  using same procedure followed with RNAseq data above.
Sequence-based prediction of target genes of reference and new miRNAs were performed with MiRanda . Pairwise comparison between log 10 transformed normalized expression matrices of mRNAs and miRNAs were performed using Pearson correlation with rcorr package in R version 3.3.1.
Midgut escape barrier
Midgut infection barrier
O’nyong nyong virus
Principal component analysis
Illumina RNA sequencing
Powers AM, Brault AC, Shirako Y, Strauss EG, Kang W, Strauss JH, Weaver SC. Evolutionary relationships and systematics of the alphaviruses. J Virol. 2001;75(21):10118–31.
Atkins GJ. The pathogenesis of alphaviruses. ISRN Virology. 2013;2013:1–22.
LaBeaud AD, Banda T, Brichard J, Muchiri EM, Mungai PL, Mutuku FM, Borland E, Gildengorin G, Pfeil S, Teng CY, et al. High rates of o'nyong nyong and chikungunya virus transmission in coastal Kenya. PLoS Negl Trop Dis. 2015;9(2):e0003436.
WCt B, Bennett KE, Gorrochotegui-Escalante N, Barillas-Mury CV, Fernandez-Salas I, de Lourdes Munoz M, Farfan-Ale JA, Olson KE, Beaty BJ. Flavivirus susceptibility in Aedes aegypti. Arch Med Res. 2002;33(4):379–88.
Hardy JL, Houk EJ, Kramer LD, Reeves WC. Intrinsic factors affecting vector competence of mosquitoes for arboviruses. Annu Rev Entomol. 1983;28:229–62.
Carissimo G, Pondeville E, McFarlane M, Dietrich I, Mitri C, Bischoff E, Antoniewski C, Bourgouin C, Failloux AB, Kohl A, et al. Antiviral immunity of Anopheles gambiae is highly compartmentalized, with distinct roles for RNA interference and gut microbiota. Proc Natl Acad Sci U S A. 2015;112(2):E176–85.
Keene KM, Foy BD, Sanchez-Vargas I, Beaty BJ, Blair CD, Olson KE. RNA interference acts as a natural antiviral response to O'nyong-nyong virus (alphavirus; Togaviridae) infection of Anopheles gambiae. Proc Natl Acad Sci U S A. 2004;101(49):17240–5.
Sim C, Hong YS, Vanlandingham DL, Harker BW, Christophides GK, Kafatos FC, Higgs S, Collins FH. Modulation of Anopheles gambiae gene expression in response to o'nyong-nyong virus infection. Insect Mol Biol. 2005;14(5):475–81.
Sim C, Hong YS, Tsetsarkin KA, Vanlandingham DL, Higgs S, Collins FH. Anopheles gambiae heat shock protein cognate 70B impedes o'nyong-nyong virus replication. BMC Genomics. 2007;8:231.
Waldock J, Olson KE, Christophides GK. Anopheles gambiae antiviral immune response to systemic O'nyong-nyong infection. PLoS Negl Trop Dis. 2012;6(3):e1565.
Yoshiga T, Georgieva T, Dunkov BC, Harizanova N, Ralchev K, Law JH. Drosophila melanogaster transferrin. Cloning, deduced protein sequence, expression during the life cycle, gene localization and up-regulation on bacterial infection. Eur J Biochem. 1999;260(2):414–20.
Levy F, Bulet P, Ehret-Sabatier L. Proteomic analysis of the systemic immune response of Drosophila. Mol Cell Proteomics. 2004;3(2):156–66.
Buchon N, Broderick NA, Poidevin M, Pradervand S, Lemaitre B. Drosophila intestinal response to bacterial infection: activation of host defense and stem cell proliferation. Cell Host Microbe. 2009;5(2):200–11.
Karlikow M, Goic B, Saleh MC. RNAi and antiviral defense in Drosophila: setting up a systemic immune response. Dev Comp Immunol. 2014;42(1):85–92.
Dong Y, Aguilar R, Xi Z, Warr E, Mongin E, Dimopoulos G. Anopheles gambiae immune responses to human and rodent Plasmodium parasite species. PLoS Pathog. 2006;2(6):e52.
Biryukova I, Ye T, Levashina E. Transcriptome-wide analysis of microRNA expression in the malaria mosquito Anopheles gambiae. BMC Genomics. 2014;15:557.
Winter F, Edaye S, Huttenhofer A, Brunel C. Anopheles gambiae miRNAs as actors of defence reaction against Plasmodium invasion. Nucleic Acids Res. 2007;35(20):6953–62.
Dennison NJ, BenMarzouk-Hidalgo OJ, Dimopoulos G. MicroRNA-regulation of Anopheles gambiae immunity to Plasmodium falciparum infection and midgut microbiota. Dev Comp Immunol. 2015;49(1):170–8.
Ghildiyal M, Zamore PD. Small silencing RNAs: an expanding universe. Nat Rev Genet. 2009;10(2):94–108.
Filipowicz W, Bhattacharyya SN, Sonenberg N. Mechanisms of post-transcriptional regulation by microRNAs: are the answers in sight? Nat Rev Genet. 2008;9(2):102–14.
Ritchie W, Flamant S, Rasko JE. Predicting microRNA targets and functions: traps for the unwary. Nat Methods. 2009;6(6):397–8.
Muniategui A, Pey J, Planes FJ, Rubio A. Joint analysis of miRNA and mRNA expression data. Brief Bioinform. 2013;14(3):263–78.
Flynt AS, Lai EC. Biological principles of microRNA-mediated regulation: shared themes amid diversity. Nat Rev Genet. 2008;9(11):831–42.
Didiano D, Hobert O. Perfect seed pairing is not a generally reliable predictor for miRNA-target interactions. Nat Struct Mol Biol. 2006;13(9):849–51.
Martin HC, Wani S, Steptoe AL, Krishnan K, Nones K, Nourbakhsh E, Vlassov A, Grimmond SM, Cloonan N. Imperfect centered miRNA binding sites are common and can mediate repression of target mRNAs. Genome Biol. 2014;15(3):R51.
Carissimo G, Eiglmeier K, Reveillaud J, Holm I, Diallo M, Diallo D, Vantaux A, Kim S, Menard D, Siv S, et al. Identification and characterization of two novel RNA viruses from Anopheles gambiae species complex mosquitoes. PLoS One. 2016;11(5):e0153881.
Fauver JR, Grubaugh ND, Krajacich BJ, Weger-Lucarelli J, Lakin SM, Fakoli LS 3rd, Bolay FK, Diclaro JW 2nd, Dabire KR, Foy BD, et al. West African Anopheles gambiae mosquitoes harbor a taxonomically diverse virome including new insect-specific flaviviruses, mononegaviruses, and totiviruses. Virology. 2016;498:288–99.
Colmant AMG, Etebari K, Webb CE, Ritchie SA, Jansen CC, van den Hurk AF, Bielefeldt-Ohmann H, Hobson-Peters J, Asgari S, Hall RA. Discovery of new orbiviruses and totivirus from Anopheles mosquitoes in eastern Australia. Arch Virol. 2017;162(11):3529–34.
Nanfack Minkeu F, Vernick KD. A systematic review of the natural Virome of Anopheles mosquitoes. Viruses. 2018;10(5)
Ren X, Hughes GL, Niu G, Suzuki Y, Rasgon JL. Anopheles gambiae densovirus (AgDNV) has negligible effects on adult survival and transcriptome of its mosquito host. PeerJ. 2014;2:e584.
Barik TK, Suzuki Y, Rasgon JL. Factors influencing infection and transmission of Anopheles gambiae densovirus (AgDNV) in mosquitoes. PeerJ. 2016;4:e2691.
Buchon N, Broderick NA, Lemaitre B. Gut homeostasis in a microbial world: insights from Drosophila melanogaster. Nat Rev Microbiol. 2013;11(9):615–26.
Xi Z, Ramirez JL, Dimopoulos G. The Aedes aegypti toll pathway controls dengue virus infection. PLoS Pathog. 2008;4(7):e1000098.
Colpitts TM, Cox J, Vanlandingham DL, Feitosa FM, Cheng G, Kurscheid S, Wang P, Krishnan MN, Higgs S, Fikrig E. Alterations in the Aedes aegypti transcriptome during infection with West Nile, dengue and yellow fever viruses. PLoS Pathog. 2011;7(9):e1002189.
Dong S, Behura SK, Franz AWE. The midgut transcriptome of Aedes aegypti fed with saline or protein meals containing chikungunya virus reveals genes potentially involved in viral midgut escape. BMC Genomics. 2017;18(1):382.
Bonizzoni M, Dunn WA, Campbell CL, Olson KE, Marinotti O, James AA. Complex modulation of the Aedes aegypti transcriptome in response to dengue virus infection. PLoS One. 2012;7(11):e50512.
Saldana MA, Etebari K, Hart CE, Widen SG, Wood TG, Thangamani S, Asgari S, Hughes GL. Zika virus alters the microRNA expression profile and elicits an RNAi response in Aedes aegypti mosquitoes. PLoS Negl Trop Dis. 2017;11(7):e0005760.
Etebari K, Hegde S, Saldana MA, Widen SG, Wood TG, Asgari S, Hughes GL. Global transcriptome analysis of Aedes aegypti mosquitoes in response to Zika virus infection. mSphere. 2017;2(6)
Baxter RH, Steinert S, Chelliah Y, Volohonsky G, Levashina EA, Deisenhofer J. A heterodimeric complex of the LRR proteins LRIM1 and APL1C regulates complement-like immunity in Anopheles gambiae. Proc Natl Acad Sci U S A. 2010;107(39):16817–22.
Povelones M, Waterhouse RM, Kafatos FC, Christophides GK. Leucine-rich repeat protein complex activates mosquito complement in defense against Plasmodium parasites. Science. 2009;324(5924):258–61.
Mitri C, Bischoff E, Takashima E, Williams M, Eiglmeier K, Pain A, Guelbeogo WM, Gneme A, Brito-Fravallo E, Holm I, et al. An evolution-based screen for genetic differentiation between Anopheles sister taxa enriches for detection of functional immune factors. PLoS Pathog. 2015;11(12):e1005306.
Zhang X, An C, Sprigg K, Michel K. CLIPB8 is part of the prophenoloxidase activation system in Anopheles gambiae mosquitoes. Insect Biochem Mol Biol. 2016;71:106–15.
Paskewitz SM, Andreev O, Shi L. Gene silencing of serine proteases affects melanization of Sephadex beads in Anopheles gambiae. Insect Biochem Mol Biol. 2006;36(9):701–11.
Volz J, Muller HM, Zdanowicz A, Kafatos FC, Osta MA. A genetic module regulates the melanization response of Anopheles to Plasmodium. Cell Microbiol. 2006;8(9):1392–405.
Varjak M, Maringer K, Watson M, Sreenu VB, Fredericks AC, Pondeville E, Donald CL, Sterk J, Kean J, Vazeille M, et al. Aedes aegypti Piwi4 is a noncanonical PIWI protein involved in antiviral responses. mSphere. 2017;2(3)
Shrinet J, Jain S, Jain J, Bhatnagar RK, Sunil S. Next generation sequencing reveals regulation of distinct Aedes microRNAs during chikungunya virus development. PLoS Negl Trop Dis. 2014;8(1):e2616.
Yan H, Zhou Y, Liu Y, Deng Y, Chen X. miR-252 of the Asian tiger mosquito Aedes albopictus regulates dengue virus replication by suppressing the expression of the dengue virus envelope protein. J Med Virol. 2014;86(8):1428–36.
Slonchak A, Hussain M, Torres S, Asgari S, Khromykh AA. Expression of mosquito microRNA Aae-miR-2940-5p is downregulated in response to West Nile virus infection to restrict viral replication. J Virol. 2014;88(15):8457–67.
Skalsky RL, Vanlandingham DL, Scholle F, Higgs S, Cullen BR. Identification of microRNAs expressed in two mosquito vectors, Aedes albopictus and Culex quinquefasciatus. BMC Genomics. 2010;11:119.
Hussain M, Torres S, Schnettler E, Funk A, Grundhoff A, Pijlman GP, Khromykh AA, Asgari S. West Nile virus encodes a microRNA-like small RNA in the 3′ untranslated region which up-regulates GATA4 mRNA and facilitates virus replication in mosquito cells. Nucleic Acids Res. 2012;40(5):2210–23.
Varjak M, Donald CL, Mottram TJ, Sreenu VB, Merits A, Maringer K, Schnettler E, Kohl A. Characterization of the Zika virus induced small RNA response in Aedes aegypti cells. PLoS Negl Trop Dis. 2017;11(10):e0006010.
Harris C, Lambrechts L, Rousset F, Abate L, Nsango SE, Fontenille D, Morlais I, Cohuet A. Polymorphisms in Anopheles gambiae immune genes associated with natural resistance to Plasmodium falciparum. PLoS Pathog. 2010;6(9):e1001112.
Holm I, Lavazec C, Garnier T, Mitri C, Riehle MM, Bischoff E, Brito-Fravallo E, Takashima E, Thiery I, Zettor A, et al. Diverged alleles of the Anopheles gambiae leucine-rich repeat gene APL1A display distinct protective profiles against Plasmodium falciparum. PLoS One. 2012;7(12):e52684.
Muller HM, Dimopoulos G, Blass C, Kafatos FC. A hemocyte-like cell line established from the malaria vector Anopheles gambiae expresses six prophenoloxidase genes. J Biol Chem. 1999;274(17):11727–35.
Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. 2011. 2011;17(1).
Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. ARXIV. 2013;
Anders S, Pyl PT, Huber W. HTSeq--a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31(2):166–9.
R: A language and environment for statistical computing. http://www.R-project.org/.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Benjamini Y, Hochberg Y. Controlling the false discovery rate - a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B-Methodological. 1995;57(1):289–300.
Alexa A, Rahnenfuhrer J. topGO: enrichment analysis for gene ontology. R package version. 2016:2300.
Kozomara A, Griffiths-Jones S. miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 2014;42(Database issue):D68–73.
Castellano L, Rizzi E, Krell J, Di Cristina M, Galizi R, Mori A, Tam J, De Bellis G, Stebbing J, Crisanti A, et al. The germline of the malaria mosquito produces abundant miRNAs, endo-siRNAs, piRNAs and 29-nt small RNAs. BMC Genomics. 2015;16:100.
Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26(6):841–2.
Mackowiak SD. Identification of novel and known miRNAs in deep-sequencing data with miRDeep2. Curr Protoc Bioinformatics. 12:10.
Stark A, Brennecke J, Russell RB, Cohen SM. Identification of Drosophila MicroRNA targets. PLoS Biol. 2003;1(3):E60.
We thank the Institut Pasteur core facility, the Center for the Production and Infection of Anopheles (CEPIA), for mosquito rearing, and Corinne Genève and Emma Brito-Fravallo of the GGIV research unit for rearing mosquitoes and laboratory support.
This work received financial support to KDV from the European Research Council, Support for Frontier Research, Advanced Grant 323173 AnoPath; European Commission, Horizon 2020 Infrastructures 731060 Infravec2; National Institutes of Health NIAID USA, AI121587; and Agence Nationale de la Recherche Laboratoire d’Excellence “Integrative Biology of Emerging Infectious Diseases” ANR- 10-LABX-62-IBEID. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Availability of data and material
All RNAseq datasets are available at the EMBL ArrayExpress Archive (https://www.ebi.ac.uk/arrayexpress/), accession number E-MTAB-6443.
All small RNAseq datasets are available at the EMBL ArrayExpress Archive (https://www.ebi.ac.uk/arrayexpress/), accession number E-MTAB-6444.
Ethics approval and consent to participate
There were no human subjects. The protocol for the ethical treatment of the animals used in this study was approved by the research animal ethics committee of the Institut Pasteur, “C2EA-89 CETEA Institut Pasteur” as protocol number B75–15-31. The Institut Pasteur ethics committee is authorized by the French Ministry of Higher Education and Research (MESR) under French law N° 2001–486, which is aligned with Directive 2010/63/EU of the European Commission on the protection of animals used for scientific purposes.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. A. coluzzii mRNAs differentially expressed 3 d after ONNV infection by bloodmeal as compared to control mosquitoes fed a normal bloodmeal. Table lists differential regulation data for the 30 ONNV-regulated genes shown in Fig. 1, including bioinformatically predicted functional information where available. (XLSX 13 kb)
Figure S1. Vector base ortholog representation of AGAP000376. (PDF 13 kb)
Figure S2. Venn diagram of differentially expressed transcripts in A. coluzzii infected with ONNV or Plasmodium. Transcriptional response of A. coluzzii to ONNV infection 3 d post-bloodmeal as measured by RNAseq in the current study (ONNV-Gut) is compared to a published study of A. coluzzii response to Plasmodium infection as measured by microarray (Dong et al., 2006). Compared conditions were transcripts differentially expressed from P. falciparum–infected midgut (Pb-Gut), P. berghei-infected midgut (Pb-Gut), or midgut after a bloodmeal containing an invasion-incompetent mutant of P. falciparum (Pf-CTRP). (PDF 181 kb)
Table S2. A. coluzzii miRNAs and regulation by treatments. Table includes novel and known miRNAs, sequences, and differential regulation data after dsRNA treatments, in mosquitoes with and without ONNV infection. Genomic coordinates of miRNA genes are based on A. gambiae PEST strain genome (assembly AgamP3). (XLSX 59 kb)
Figure S3. Venn diagram of A. coluzzii differentially expressed miRNAs during ONNV infection. Mosquitoes were either not treated with dsRNA (no dsRNA) or were treated before bloodfeeding, with or without ONNV, with dsRNA for control LacZ (dsLacZ treatment), for Imd pathway factor Rel2 (dsRel2 treatment), or for JAK/STAT pathway factor Stat-A (dsStat−A treatment). Name and number of differentially expressed miRNAs are indicated. Details on miRNAs are in Additional File 4: Table S2. (PDF 147 kb)
Figure S4. Hierarchical clustering of 11 A. coluzzii miRNAs differentially expressed among treatments. Treatments are indicated on the right vertical axis. Mosquitoes were either not treated with dsRNA (No-dsRNA) or were treated before bloodfeeding with dsRNA for control LacZ (dsLacZ), for Imd pathway factor Rel2 (dsRel2), or for JAK/STAT pathway factor Stat-A (dsStat−A). Bloodfeeding was either without ONNV (Non-infected) or with ONNV (Infected). Names of miRNAs are indicated on the x-axis. Details on miRNAs are in Additional File 4: Table S2. (PDF 292 kb)
Table S3. A. coluzzii miRNA-mRNA target gene pairs by sequence-based target prediction and anti-correlation with mRNA expression levels. Table indicates complete transcript table with miRNA data. (XLSX 6821 kb)
Table S4. Primer list. Primers used for synthesis of double-stranded RNAs (prefix T7, T7 RNA polymerase promoter underlined) or qPCR analysis (prefix q) of target genes. Final suffix indicates forward, F, or reverse, R, sense of primers. (DOCX 13 kb)