The expression profile of Aedes albopictus miRNAs is altered by dengue virus serotype-2 infection

Aedes albopictus is an important vector of Dengue virus (DENV) and it has quickly invaded the tropical and temperate environments worldwide. A few studies have shown that, microRNAs (miRNAs) regulate mosquito defense against pathogens. However, there is no systematic analysis of the impact of DENV infection on miRNA expression in Ae. albopictus. We conducted this study to investigate the miRNA expression of Ae. albopictus upon DENV-2 infection using Illumina RNA sequencing. A total of 103 known and 5 novel candidate miRNAs were identified in DENV-2 infected and non-infected adult female Ae. albopictus. Comparative analysis indicated that 52 miRNAs were significantly down-regulated and 18 were up-regulated significantly after infection. Furthermore, RT-qPCR validated the expression patterns of eleven of these differentially expressed miRNAs. Targets prediction and functional analysis of these regulated miRNAs suggested that miR-34-5p and miR-87 might be involved in the anti-pathogen and immune responses. This is the first systematic study on the impact of DENV infection on miRNA expression in Ae. albopictus. Complex changes in miRNA expression suggest a potential role of miRNAs in antiviral responses by regulating immune-related genes. This investigation provides information concerning DENV-induced miRNAs and offers clues for identifying potential candidates for vector based antiviral strategies.


Background
Dengue fever is among the most important global arboviral diseases of humans. An estimated 50 million human cases of dengue virus (DENV) infections occur annually, and over 40% of the world's population is at risk of infection [1]. Currently, there is no effective chemotherapy, specific treatment or effective vaccine for DENV infection owing to safety concerns and potential antibody-dependent enhancement. Therefore, the only effective way to prevent or control dengue transmission is to control the vector mosquito. Aedes albopictus (Ae. albopictus) is an important vector of DENV and is an important invasive species that quickly and aggressively spreads worldwide in tropical and temperate environments [1][2][3].
MicroRNAs (miRNAs) have been shown to regulate apoptosis, viral infection and other critical biological events in animals and plants [4][5][6][7]. In mosquitoes, a large number of miRNAs have been reported to exhibit altered profiles during infection and regulate the host immune responses [8][9][10][11][12]. Knocking down the Dicer1 and Ago1 mRNAs in Anopheles leads to increased sensitivity to Plasmodium infection, suggesting that miRNAs might be involved in defense reactions [11]. Moreover, Wolbachia uses aae-miR-2940 to regulate a methyltransferase gene for blocking DENV replication [13,14]. In addition, aga-miR-2304, aga-miR-2390 and aae-miR-375 have been identified to be involved in the modulation of host immune response in mosquitoes, and aae-miR-375 enhances DENV-2 infection in an Aedes aegypti (Ae. aegypti) cell line [15,16]. These studies indicate that mosquito miRNAs can positively or negatively modulate the host response to pathogen infection.
However, although Ae. albopictus is an important vector of DENV, the relationship between Ae. albopictus miRNAs and DENV infection remains unknown. In this study, to gain a better understanding of DENV infection in Ae. albopictus at the miRNA level, we performed deep sequencing of small RNAs in Ae. albopictus upon DENV-2 infection. We also comparatively analyzed the differentially expressed miRNAs between DENV-infected and control mosquitoes. Furthermore, we predicted the targets of these regulated miRNAs to investigate the roles of these miRNAs during DENV infection in Ae. albopictus.

Results
High-throughput sequencing and annotation of small RNAs in Ae. albopictus To evaluate the impact of DENV infection on miRNA expression in Ae. albopictus, we performed Illumina sequencing on small RNA libraries obtained from DENV-2 infected mosquitoes and compared it with uninfected mosquitoes. Approximately 10 and 20 million reads were obtained from control and infected mosquitoes, respectively. After removal of adaptors, contaminants and low quality reads, the length distribution of clean reads was summarized in Figure 1. We found that the peak in the control was at 21-22 nt with a percentage of 35.79; however, the peak of the infected sample was at 27-28 nt with a percentage of 44.82. As shown in Table 1, 4,805,249 clean reads from control and 3,309,844 clean reads from infected mosquitoes were mapped to Ae. aegypti genome using SOAP. Intron, exon, tRNA, rRNA, miRNA, snRNA and snoRNA reads were annotated, respectively.
The clean reads were also mapped against DENV genome to discover DENV encoded small RNA population. We found only 76 reads mapping against DENV genome (Additional file 1).

Identification of known miRNAs in Ae. albopictus
To identify known miRNAs in Ae. albopictus, we aligned the small RNAs to the known mature miRNAs and their precursors in miRBase 20.0 to obtain the miRNA count as well as the base bias at the first position. Approximately 1,370 unique sequences (3,944,255 reads) in the control library and 1,284 unique sequences (2,261,930 reads) in the infected library were annotated as miRNA candidates (Table 1). A total of 103 known miRNA genes were identified in the Ae. albopictus library (Additional file 2). Most miRNAs ranged in size from 21 to 22 nt, and the base bias at the first position of the identified miRNAs showed a strong preference for ′U′ at the 5′-end, which is coordinated with previous studies (Additional file 3A and B) [17].
We also analyzed miRNA conservations using the available genome assemblies of other insect species, and we confirmed previous study that miR-1889, miR-1890 and miR-1891 were present only in mosquitoes [10]. The conservation of 103 known miRNAs across three other species-Ae. aegypti, Culex quinquefasciatus and Anopheles gambiae-is shown in Additional file 4. Of these 103 miRNAs, all were found in Ae. aegypti; 76 were found in all three species; 8 were found in Ae. aegypti, but not in A. gambiae or Cx. quinquefasciatus; 6 were found in Ae. aegypti and A. gambiae, but not in Cx. quinquefasciatus; and 13 were common to Ae. aegypti and Cx. quinquefasciatus, but were not found in A. gambiae (Additional file 3C).

Identification of novel miRNAs in Ae. albopictus
After removing the snRNAs, snoRNAs, rRNAs, tRNAs and known miRNAs, we aligned the remaining unannotated reads against the Ae. aegypti genome to predict novel miRNAs in Ae. albopictus. A total of 7107 and 23697 potential novel miRNA reads were mapped in control and infected samples. Based on the criteria for miR-NAs, 5 novel candidate miRNAs were identified (4 in control and 4 in infected mosquitoes). The miRNAs are shown in Additional file 2, along with their Dicer cleavage site, minimum free energy, frequency of reads and typical secondary structures of the characteristic stem- loop hairpins for the predicted precursors. The secondary structures of the novel candidate miRNAs in Ae. albopictus are also shown in Additional file 5.
After that, the remaining unannotated reads were mapped to miRBase 20.0 to predict potential miRNA, a total of 3383 and 19743 predict tags mapped in control and infected library, respectively. As shown in Additional file 6, 3289 and 16353 reads matched to Bombyx mori bmo-miR-2796-3p, showed differential expression in control and infected libraries.

Differentially expressed miRNAs between control and infected Ae. albopictus
To compare miRNA abundance in both libraries, we used TPM (tag per million of total RNA reads) to normalize the miRNA expression. And to directly test our sequencing methods' ability to reproduce miRNA expression levels, a linear regression analysis was performed using TPM of all miRNAs in control and infected mosquitoes. This analysis showed a Pearson correlation coefficient (r) of 0.86, indicating that the two libraries were well correlated.

Validation of differentially expressed miRNAs by RT-qPCR
To further validate the expression of differentially expressed miRNAs from high-throughput sequencing, RT-qPCR was performed on five novel candidate miR-NAs and seven randomly selected known miRNAs. All of the miRNAs showed consistent expression profiles with the small RNA sequencing data. Results confirmed the down-regulation of 7 miRNAs (miR-10, miR-1890, miR-263a-5p, miR-263b-5p, miR-281-5p, let-7 and novel-1) and the up-regulation of 4 miRNAs (miR-2945, novel-3, novel-4 and novel-5) in infected mosquitoes compared with the uninfected ( Figure 4A).

Prediction of targets for differentially expressed miRNAs
To investigate the regulated miRNAs functions, we predicted the potential targets of the 66 differentially expressed known miRNAs in Ae. albopictus upon DENV-2 infection. A total of 1421 genes were predicted, including 35 immune transcripts, such as Toll-like receptor  (TOLLs), clip domain serine protease (CLIPs), serine protease inhibitors (SRPNs), scavenger receptors (SCRs), and so on. The network of interaction between miRNAs and immune targets was shown in Figure 5. Nine immune target genes were randomly selected to perform RT-qPCR in control and DENV2-injected mosquitoes. As a result, 6 targets (AAEL000028, AAEL004540, AAEL005064, AAEL009192, AAEL003697, AAEL003857) were upregulated, AAEL009692 was down-regulated, and 2 target genes (AAEL007613, AAEL010427) were nonregulated ( Figure 4B). To further study the potential roles of the regulated miRNAs in DENV-Ae. albopictus interactions, we performed GO enrichment and pathway analysis on the top three targets of the differentially expressed miRNAs (Additional file 7) using DAVID Bioinformatics software [18]. Ten GO terms were found to be enriched, including three immune terms ("innate immune response", "immune response", and "defense response"). Based on KEGG analysis, many target genes of miRNAs were involved in Jak-STAT signaling pathway and Toll receptor signaling pathway. Taken the KEGG and GO analysis together, miR-34-5p and miR-87 were found to be related to immune pathway (Toll-like receptor signaling pathway) and immune GO terms (defense response, immune response and innate immune response) ( Figure 6).

Discussions
MiRNAs have been shown to be involved in the mosquito immune responses [19][20][21][22]. In this study, we aim to identify miRNAs involved in the response to DENV-2 in Ae. albopictus using high-throughput sequencing technology. We obtained about 10 and 20 million reads from control and infected mosquitoes, but more than 60% of reads in each of the sequencing libraries were not annotated. It may be attributed to the differences between Ae. albopictus and Ae. aegypti genome.
We identified 103 known and 5 novel candidate miR-NAs in Ae. albopictus. Previous studies showed that, pathogen infection can influence the expression of mosquito miRNAs. In Ae. aegypti, distinct miRNA expression profiles were observed after CHIKV infection, WNV infection and Wolbachia infection [23][24][25]. Likewise, Plasmodium also changed the miRNAs expression in A. gambiae [11]. In this study, we found that 66 known miRNAs displayed different expressions during DENV infection. Previous reports showed that, miR-1889-5p and miR-9c-5p were down-regulated in Wolbachia and DENV infected Ae. aegypti, respectively [24,26]. miR-989, which was down-regulated in WNV infected Cx. quinquefasciatus and Wolbachia infected Ae. aegypti, has been reported to play roles in mediating flavivirus infection in the mosquito host [12,14]. In our study, miR-989, miR-1889-5p and miR-9c-5p were also down-regulated following DENV infection. These results indicate the possibility that, as observed in other mosquitoes, miRNAs likely participate in host-virus interaction in Ae. albopictus.
In mosquitoes, a large number of miRNAs have been reported to regulate the host immune responses. miR-2940-5p, which is highly induced in Wolbachia-infected Ae. aegypti, was previously reported to enhance Wolbachia efficient maintenance and limit replication of DENV in Ae. aegypti [13]. Another study demonstrated that miR-2940-5p is down-regulated in Ae. albopictus C6/36 cells in response to WNV infection and restrict WNV replication [25]. In addition, miR-2940-5p and miR-2940-3p were reported decreased in CHIKV-infected Ae. albopictus [23]. In this study, we found that, miR-2940-5p and miR-2940-3p were significantly downregulated upon DENV-2 infection. These studies indicate that miR-2940-5p and miR-2940-3p may also be involved in interactions between Ae. albopictus mosquito hosts and DENV.
We have further predicted targets of the 66 modulated miRNAs, and note that some differentially expressed miRNAs target a number of immune genes. Then we chose the top three targets of each miRNA to perform GO enrichment and pathways analysis. Taken the GO and KEGG analysis together, we found miR-34-5p and miR-87 to be related to immune GO terms and immune pathway. The down-regulated miR-87, which was reported previously, decreased upon CHIKV infection in Ae. albopictus mosquitoes [23], targeted the TOLL pathway signalling  The relative expression of miRNAs were calculated against control mosquito as a calibrator using 2 -ΔΔCt method and normalized to 5S rRNA. (B) Expression levels of 9 target mRNAs were detected in control and infected mosquitoes. mRNA levels relative to GAPDH were calculated. T-tests were used for the comparison of RT-qPCR data. Error bars show the standard deviation from three independent experiments with three triplicates each. Asterisks indicate the statistical significance: ***P < 0.001, **P < 0.01, *P < 0.05. immune transcripts. Another target of miR-34-5p, the TOLL pathway signalling NF-kappaB Relish-like transcription factor (REL1, AAEL007696), has been reported in another study as the activation of Rel1 leads to a decreased DENV titer, while repression of Rel1 activation leads to a significant increase in DENV replication [27]. Mosquitoes rely on their effective innate immune system to limit pathogen infections. Studies that were mainly conducted in the Ae. aegypti have shown that, immune responses are largely regulated by three major immune signaling pathways: the Toll, Imd, and JAK-STAT pathways [28]. The Toll pathway has been shown to direct antidengue defenses in the A. aegypti [27,[29][30][31]. In this study, based on KEGG analysis, we found miR-34-5p and miR-87 may be related to the Toll-like receptor signaling pathway. And the GO analysis showed miR-34-5p and miR-87 may be involved in the defense response, immune response and innate immune response. Taken together, these studies indicate that miR-34-5p and miR-87 may be associated with the anti-pathogen and immune responses.

Conclusion
In conclusion, we identified miRNAs and evaluated their expression patterns in Ae. albopictus upon DENV-2 infection using Illumina sequencing. Sixty-six known differentially expressed miRNAs and 5 novel candidate miRNAs were identified. Target prediction and functional analysis of these miRNAs suggested that miR-34-5p and miR-78 might be involved in anti-pathogen and immune responses. This study provides important information for studying the interaction between DENV and Ae. albopictus miRNAs, and it also provides clues for further studies on the mechanisms of immune responses in DENV-infected Ae. albopictus.

Ethics statement
All vertebrate animals were housed and handled in strict accordance with the guidelines of the institutional and National Committees of Animal Use and Protection. All
The Aedes albopictus strain was obtained from the CDC of Guangdong Province, PRC, which was isolated from Foshan, Guangdong, China and established in the laboratory since 1981. Mosquitoes were maintained at 27 ± 1°C, with 70-80% relative humidity and a 12:12 hour light-dark cycle. Larvae were fed with yeast powder, and adults were fed with a 10% glucose solution. Adult female mosquitoes were fed with the blood of an anesthetized mouse, which was then returned to the animal room.
For mosquito infection through intra-thoracic inoculation, 0.2 μL of the above virus suspension (10 6.85 PFU/ ml) was used for inoculation into each female; PBS was used as the control. Then, the mosquitoes were maintained at (27 ± 1)°C with 70-80% humidity and supplied with 10% glucose solution. Pools of 14 whole mosquitoes were harvested at 7 dpi from both DENV2-injected and control mosquitoes. Infection rate, as determined by plaque assay [26,37,38], was 100% at 7 dpi.

RNA isolation and deep sequencing
Pools of 14 whole mosquitoes were ground together to a fine powder in liquid nitrogen at 7 dpi. Total RNA was extracted from 14 whole mosquitoes using TRIzol Reagent (Invitrogen) according to manufacturer's instructions. RNAs ranging from 18 to 30 nt were then purified from a 15% denaturing polyacrylamide gel and sequentially ligated to proprietary adapters according to the manufacturer's instructions (Illumina Inc) before sequencing. The gel-purified ligation products were converted to cDNA and amplified using RT-PCR with 18 PCR cycles to produce libraries that were then sequenced using the Illumina Genome Analyzer (Illumina, San Diego, CA, USA) at Huada Genomics Institute Co. Ltd, Shenzhen, China. Raw sequencing data have been submitted to the National Center for Biotechnology Information (NCBI) (accession number SRA129084).

Bioinformatics
After removing the adapter sequences and low-quality sequences, tags with lengths ranging from 18 to 30 nt were screened against the GenBank database (http:// www.ncbi.nlm.nih.gov/genbank/) and Rfam database (http://www.sanger.ac.uk/software/Rfam) to remove rRNA, tRNA, snRNA, snoRNA, repeats, exon sequences and intron sequences. Considering that no miRNA information for Ae. albopictus was in miRBase 20.0, the remaining sequencing reads were aligned to search all known precursor/mature miRNAs of Ae. aegypti in miRBase 20.0 [39]. Only perfectly or near-perfectly (1-2 mismatches) matching sequences were considered to represent conserved miRNAs, and each miRNA at least 15 reads to be considered expressed. Finally, to discover potential miRNAs, all unannotated reads were aligned against the Ae. aegypti genome using SOAP [40]. RNA secondary structures were predicted using RNAfold (http:// rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi) and analyzed using MIPRED and MIREAP (http://sourceforge.net/projects/ mireap/). The identified stem-loop hairpins fulfilled three criteria: mature miRNAs were present in one arm of the hairpin precursors, which also lacked large internal loops or mismatches; the secondary structures of the hairpins were stable, with free energies of hybridization lower than -20 kcal/mol; and hairpins were located in intergenic regions or introns. Subsequently, miRAlign was adopted to filter mature miRNA genes and predict new miRNA genes, and each novel miRNA in the samples was required to have at least 15 reads to be considered as expressed. The whole process is shown as Additional file 8.

Expression profiling of miRNAs in response to DENV
All miRNA expression data were normalized using tag per million of total RNA reads (TPM) to compare miRNA abundance in both libraries. TPM = (The count of the miRNA in a particular sample)/(Total miRNA counts from this sample) × 1 million. Because extremely low abundances might not accurately reflect their true abundance, the known miRNAs with TPM less than 10 in the two libraries were removed from the differential expression analysis. Changes in known miRNA expression in infected versus control samples were considered significant when their p values were below 0.05. miRNAs with log2-fold changes higher than 1 were designated as significantly up-regulated, and those with log2-fold changes less than -1 were designated as significantly down-regulated.

Target prediction of microRNAs and functional analysis
RNAhybrid and SMVLight were used to predict the targets of the miRNAs (p value < 0.05, Prediction-value > 1, a minimum free energy (mfe) < -20.0 kcal/mol, one G:U mismatch was permitted in the seed sequence, and a maximum of three G:U mismatches were allowed in the flanking sequence).
The vectorbase immune-related targets (https://www. vectorbase.org/) were then used to generate the miR-NA:mRNA network. Cytoscape was used for visualizing the networks [41].
Then, the top three targets of the differentially expressed miRNAs were performed GO enrichment and pathway analysis using the DAVID Bioinformatics Resources 6.7 (http://david.abcc.ncifcrf.gov/).

RT-qPCR
Real time PCR analysis was performed on known miR-NAs and target genes using the miScript PCR system, including the miScript Reverse Transcription Kit (Qiagen, Duesseldorf, Germany) and the miScript SYBR Green PCR kit (Qiagen). Amplification was carried out using MX3005P TM Real Time PCR System (Stratagene, CA, USA). Real time PCR analysis was also performed on novel candidate miRNAs using a stem-loop-mediated reverse transcription real-time PCR method. Expression levels of novel miRNAs were analyzed using a 7500 Real-Time PCR system (Applied Biosystems). All reactions were run in triplicate. 5S rRNA and GAPDH were used as endogenous control for miRNA and DENV-2, respectively. All the qPCR primers used in this study are shown in Additional file 9. Expression levels were then calculated against control mosquito as a calibrator using 2 -ΔΔCt method.

Statistical analysis
All quantitative RT-qPCR results were representative of at least three independent experiments, each with three technical replicates. Statistical tests were performed with SPSS 22.0 (SPSS, Chicago, IL). T-tests were used for the comparison of RT-qPCR data. The P values were performed on the data with the significance threshold selected as 0.05. Asterisks indicate the statistical significance: ***P < 0.001, **P < 0.01, *P < 0.05. Error bars show the standard deviation.