Characterization and potential role of microRNA in the Chinese dominant malaria mosquito Anopheles sinensis (Diptera: Culicidae) throughout four different life stages

microRNAs (miRNAs) are one kind of small non-coding RNAs widely distributed in insects. Many studies have shown that miRNAs play critical roles in development, differentiation, apoptosis, and innate immunity. However, there are a few reports describing miRNAs in Anopheles sinensis, the most common, and one of the dominant malaria mosquito in China. Here, we investigated the global miRNA expression profile across four different developmental stages including embryo, larval, pupal, and adult stages using Illumina Hiseq 2500 sequencing. In total, 164 miRNAs were obtained out of 107.46 million raw sequencing reads. 99 of them identified as known miRNAs, and the remaining 65 miRNAs were considered as novel. By analyzing the read counts of miRNAs in all developmental stages, 95 miRNAs showed stage-specific expression (q < 0.01 and |log2 (fold change)| > 1) in consecutive stages, indicating that these miRNAs may be involved in critical physiological activity during development. Sixteen miRNAs were identified to be commonly dysregulated throughout four developmental stages. Many miRNAs showed stage-specific expression, such as asi-miR-2943 was exclusively expressed in the embryo stage, and asi-miR-1891 could not be detected in larval stage. The expression of six selected differentially expressed miRNAs identified by qRT-PCR were consistent with our sequencing results. Furthermore, 5296 and 1902 target genes were identified for the dysregulated 68 known and 27 novel miRNAs respectively by combining miRanda and RNAhybrid prediction. GO annotation and KEGG pathway analysis for the predicted genes of dysregulated miRNAs revealed that they might be involved in a broad range of biological processes related with the development, such as membrane, organic substance transport and several key pathways including protein processing in endoplasmic reticulum, propanoate metabolism and folate biosynthesis. Thirty-two key miRNAs were identified by microRNA-gene network analysis. The present study represents the first global characterization of An. sinensis miRNAs in its four developmental stages. The presence and differential expression of An. sinensis miRNAs imply that such miRNAs may play critical roles in An. sinensis life cycle. A better understanding of the functions of these miRNAs will have great implication for the effective control of vector population and therefore interrupting malaria transmission.


Background
Anopheles sinensis is one of the most important disease vectors in China. It is considered as the primary vector of P. vivax malaria due to its wide distribution and high density in central China where several malaria outbreaks occurred in history [1,2]. Like many mosquitoes, An. sinensis goes through four life stages that include egg, larva, pupa and the adult. During these developmental periods, a fine-tuning complex biological process, such as embryogenesis, organ differentiation, and metamorphosis, etc., is completed.
For the past decades, facilitated by high-throughput sequencing technology, a large number of small noncoding RNAs (ncRNAs) have been discovered in many species [3,4]. Amongst them, microRNAs (miRNAs) have become the most popular research topic for its critical role in regulating important biological events at the post-transcriptional level. Mature miRNAs are singlestranded, evolutionarily conserved endogenous small ncRNAs of approximately 22 nucleotides (nt) in length, which could lead to transcriptional decay through translational inhibition, transcript degradation or both [5,6].
The first miRNA repertoire of mosquito species was reported from the African malaria vector Anopheles gambiae [7]. Since then, mosquito miRNAs have been identified and characterized in a number of mosquito species, including Anopheles stephensi [8,9], Aedes aegypti [10][11][12], Culex fatigans [13][14][15], and many others, but not yet in An. sinensis. Several studies suggested that miRNAs may contribute to an important regulatory role in the post-transcriptional level in developmental stages. For instance, temporal and stage specific expression of miR-NAs has been identified in the embryo, larval, pupal and adult stage of Ae. albopictus [16], An. anthropophagus [17], and An. stephensi [9]. Among them, the expression of aal-miR-286b was increased in the embryo, while aal-miR-2942 had the highest level of expression in larvae although it was normally expressed at the embryo, pupal and adult stages [16]. In addition, ast-miR-1891, ast-miR-190-3p, ast-miR-285, ast-miR-988-3p, and ast-miR-989 were absent in the larval stage, while ast-miR-8-3p was the most abundant in both male and female larval stage [9]. These results indicated that dysregulated miRNAs might play pivotal roles in embryogenesis and metamorphosis throughout development. However, the exact functions of these miRNAs remain to be determined.
In this study, we provided a comprehensive repertoire of An. sinensis miRNAs and investigated the differential expression of miRNAs using Solexa sequencing technology, with an aim to identify key regulatory miR-NAs through development stages. Selected candidate miRNAs were then validated by qRT-PCR in samples that have been used in sequencing. To explore the potential gene regulatory networks, 32 dysregulated miRNAs were defined as key miRNAs during development and then subjected to miRNA-gene network analysis in terms of developmental pathway analysis. The data described here served as an important step towards the understanding of gene expression modulation by miRNAs in An. sinensis, which will further advance functional studies for novel mosquito control approaches developing and promote elimination campaign in China.

Mosquito rearing and sample preparation
Anopheles sinensis (China strain) from insectary were reared at 28 ± 2 °C, 70-75% humidity in the Vector Control Reference Laboratory at the National Institute of Parasitic Diseases (NIPD), Shanghai, China. Adult mosquitoes were maintained in screened cages, and provided constant access to water and glucose-soaked sponges.
Sample collections from different developmental stages of An. sinensis for high throughput sequencing are briefly described below. Embryos were collected within 24 h after placing a damp collection filter paper within the cage and were pooled to represent the embryonic stage. In the larval sample collection, we used mixed larval sample (I-IV instars) including early and late larval samples. Pupal samples were collected from a pool of varied ages. Female adults were collected after emerging one to 3 days. All samples were flash frozen in liquid nitrogen immediately following collection, and then stored at − 80 °C until for RNA isolation.

RNA extraction and small RNA sequencing
Briefly, total RNAs were extracted separately from the four different life cycle stages of An. sinensis (eggs, larvae, pupae and unfed adult females) using TRIzol reagent (Invitrogen, USA) according to the manufacturer's protocol. RNA purity was checked using the NanoPhotometer ® spectrophotometer (IMPLEN, CA, USA). RNA degradation and contamination were monitored on 1% agarose gels. RNA concentration was measured using Qubit ® RNA Assay Kit in Qubit ® 2.0 Fluorometer (Life Technologies, CA, USA). RNA quality assessment was determined by using Agilent 2100 Bioanalyzer RNA Nano 6000 kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). No rRNA depletion was performed. RNA extracts of the four life stages were sent to Newgenes Inc (Shanghai, China) for small RNA sequencing. A total amount of 3 µg of total RNA was ligated with 3′ and 5′ adaptors followed by reverse transcription using RT primers. The libraries were size-selected for sequencing of RNA fragments of 18-30 nucleotides. These small RNA libraries were then sequenced using Illumina Hiseq 2500 platform (Illumina Inc., USA).

Small RNA sequencing data processing and analysis
Following sequencing, raw data (raw reads) of fastq format were firstly processed through custom perl and python scripts. In this step, clean data (clean reads) were obtained by removing reads containing ploy-N, with 5′ adapter contaminants, without 3′ adapter or the insert tag, containing ploy A or T or G or C and low quality reads from raw data. At the same time, Q20, Q30, and GC-content of the raw data were calculated. Then, a certain range of clean reads from each stage was combined to do all the downstream analyses. All reads were mapped to the An. sinensis strain China genome [Gen-Bank: GCA_000441895.2] by Bowtie [18] without mismatch to analyze their expression and distribution on the reference.
For identification of known and novel miRNA sequences present in the four life stage datasets. Firstly, after removing tags originating from protein-coding genes, repeat sequences, rRNA, tRNA, snRNA, and snoRNA, small RNA tags were mapped to RepeatMasker, Rfam database or those types of data from Anopheles gambiae. Mapped small RNA tags were used to look for known miRNA using miRBase21 [19] as reference. Custom scripts adapted from quantitative script of miRD-eep2 were used to obtain the miRNA counts as well as base bias on the first position of identified miRNAs with certain length and on each position of all identified miR-NAs. The reads that did not yield a match were used to predict novel miRNAs. The available software miREvo [20] and mirdeep2 [21] were integrated to predict novel miRNA from the small RNA tags unannotated in the former steps. A number of stringent criterions were used for evaluating whether novel miRNA candidates were genuine miRNAs, such as exploring the secondary structure by RNAfold, calculating minimum free energy (MFE) and avoiding isolated base pairs. In each library, same sequences were pooled together to generate expression files (fasta) with their read count and unique ID used for further analysis.
All identified known miRNAs were named according to their most similar miRNAs in the miRNA database (miR-Base) release 21 (June 2017). All identified novel miRNAs were named sequentially following the asi-mir-nov1. Raw sequencing data were prepared to submit to the National Center for Biotechnology Information (NCBI).

Quantification of miRNA and differential expression analysis
In order to analyze the relative abundance and expression level of miRNAs, tags per million of total RNA reads (TPM) for each miRNA in all four libraries were calculated. TPM was compared between the libraries to identify miRNAs differentially expressed between two consecutive life stages in the mosquito. Differential expression analysis of two samples was performed using the DEGseq [22] R package. p value was adjusted using q value [23]. q < 0.01 and |log2 (fold change)| > 1 was set as the threshold for significantly differential expression by default.

miRNA target prediction, GO and KEGG enrichment analysis
To predict the potential targeted genes of miRNA, the targeted mRNAs of differentially expressed known and novel miRNAs were predicted by combining miRanda [24] and RNAhybrid [25] to select and assess overlapping predicted targets between each program. Briefly, 3′UTR of An. gambiae downloaded from vectorbase were used for the prediction. The targets showing complementarity with miRNA seed region and binding energy ≤ 20 kcal/ mol were selected. The predicted targets were used for further GO and KEGG enrichment analysis.
Gene ontology (GO) enrichment analysis was used on the target gene candidates of differentially expressed miRNAs ("target gene candidates" in the following). GOseq based Wallenius non-central hyper-geometric distribution [26], which could adjust for gene length bias, was implemented for GO enrichment analysis. GO terms with corrected p value less than 0.05 were considered significantly enriched by target genes. KEGG pathway analyses were performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.ad.jp/ kegg/) database [27]. In this study, we used KOBAS [28] software to test the statistical enrichment of the target gene candidates in KEGG pathways, with the threshold of corrected p < 0.05.

miRNA-gene network analysis
To identify the key miRNAs in An. sinensis development, microRNA-gene network analysis was conducted for the commonly dysregulated miRNAs identified from the microRNA expression files, the exclusively expressed miRNAs in certain stages (TPM > 10) and the most significantly over/under-expressed miRNAs (>tenfold) and their predicted target developmental genes. Developmental genes of An. gambiae were accessed based on gene annotation in VectorBase [29] and through intensive literature surveys [30][31][32][33]. Interactive networks between key miRNAs and their target developmental genes were constructed by using Cytoscape [34].

miRNA expression validated by quantitative RT-PCR
To validate the presence and expression level of the identified miRNAs, nine miRNAs with differential expression were selected for qRT-PCR as described previously [35]. The total RNA was extracted using TRIzol Reagent (Invitrogen, USA). The reverse transcription reaction was performed with RevertAid First Strand cDNA Synthesis Kit (Fermentas, USA), and the reverse-transcribed products were used as the template for qRT-PCR with miRNA-specific primers (Additional file 1: Table S1). All reactions were assayed in three biological and technical replications, and performed in an ABI-7300 (Applied Biosystems, USA) using SYBR Green qPCR kit (Thermo, USA). Quantitative measurements were normalized to the internal control of U6 small nuclear (U6), and relative expression was calculated by using the 2 −ΔΔCt method [14].

Length distribution of small non-coding RNAs in the four libraries
Among the libraries, size distribution of reads varied between 18 and 35 nt. We found that one peak with 21 nt in all libraries, followed by another peak with 20-23 nt which represents miRNAs (Fig. 1). This was consistent with the typical size of miRNAs in other reports. There was also a peak at 24-27 nt in almost all developmental libraries, which corresponds to piRNA-like small RNAs suggesting a possible role in embryogenesis and maintaining genome stability during development.

Annotation of small non-coding RNAs in the four libraries
In order to annotate other small non-coding RNAs in the four libraries, all clean reads were aligned against Rfam version 11 [36]. As showed in Fig. 2, the most abundant class of small non-coding RNAs in the eggs' library was miRNAs (50.10%). However, rRNAs were most abundant in the larvae (6.61%) and pupae library (5.8%), and tRNAs in the adult females' library (6.87%).

Identification of miRNAs in the four life stages libraries of An. sinensis
Based on previously known mosquito mature miRNAs in miRBase release 21 (June 2017), we discovered 99 known miRNAs in An. sinensis (Table 2) throughout four development stages. Reads that do not match any known miRNA sequences were further employed by miREvo and mirdeep2 to investigate the folding properties of precursors or hairpins for novel miRNAs. Finally, 65 mature miRNAs were identified as novel (default parameters of miREvo and mirdeep2) (Additional file 2: Table S2). Thus, a total of 164 An. sinensis miRNAs were identified in four life stages libraries of An. sinensis.

miRNA expression pattern during An. sinensis development
The expression patterns of miRNAs in different life stages were examined based on the number of reads revised. The result showed that the expression profiles of miRNAs varied from highly specific to ubiquitous during the four stages. One-hundred and seventeen miRNAs were ubiquitously expressed in all developmental stages, while some miRNAs showed stage-specific patterns (Fig. 3). Among the known miRNAs, asi-miR-276-3p, asi-miR-1, asi-miR-100, asi-miR-8-3p, asi-miR-999, and asi-miR-9a showed constant high expression throughout four development stages. asi-miR-276-3p was the most frequently expressed miRNA in the larvae, the pupae, and the adult females' library (TPM > 3.17 × 10 5 ), whereas asi-miR-1 was the dominant miRNA in the eggs' library (TPM > 1.87 × 10 5 ). Nevertheless, asi-miR-219 had the lowest expression level in all developmental stages. asi-miR-2943 was specifically expressed in the embryo stage, and asi-miR-1891 could not be detected in the larval stage. The highly expressed novel miRNA in the eggs and pupae libraries was asi-miR-nov54 (TPM = 0.45 × 10 3 and TPM = 1.50 × 10 3 ), while in the larvae library and adult females' library, it was asi-miR-nov55 (TPM = 0.91 × 10 3 ) and asi-miR-nov1 (TPM = 0.42 × 10 3 ), respectively. However, the expression level of novel miRNAs was generally low throughout all developmental stages when compared with known miRNAs.

MicroRNA regulation from egg to adult during An. sinensis development
Hierarchical heat map revealed different expression patterns of the miRNAs between the two consecutive 15 significantly up-regulated miRNAs (12 for known and 3 for novel) were expressed form the pupal to adult stage (Fig. 4). It is noteworthy that the total number of dysregulated miRNAs decreased from the embryo to the adult stage. We also identified two known miRNAs (asi-miR-2943 and asi-miR-286b) to be significantly over-expressed (>tenfold) form the embryo to the larval stage and one known miRNA (asi-miR-1891) to be significantly under-expressed form the pupal to the adult stage (>tenfold) (Fig. 4). Sixteen miRNAs were differentially expressed in all developmental stages of the mosquito. Based on the analysis of these regulated miRNAs between the two consecutive developmental stages, five miRNAs (asi-miR-1890, asi-miR-193, asi-miR-263a-5p, asi-miR-263b-5p and asi-miR-nov54) were down-regulated in the larval stage whereas get up-regulated in the pupal stage but down-regulated in the adult stage again. Seven miRNAs (asi-miR-1175-3p, asi-miR-1175-5p, asi-miR-281-3p, asi-miR-281-5p, asi-miR-317, asi-miR-34-5p and asi-miR-nov55) were up-regulated in the larval stage whereas get down-regulated in the pupal but up-regulated in the adult stage again. asi-miR-nov23   and asi-miR-2765 were up-regulated in the larval stage whereas get down-regulated in both pupal and adult stage. asi-miR-2944a-5p was down-regulated in the larval and the pupal stages but up-regulated in the adult stage. asi-miR-989 was down-regulated in larval and pupal stages but up-regulated in the adult stages.

Validation of differentially expressed miRNA
To validate the small RNA sequencing results, we performed quantitative real-time PCR analysis by selecting nine miRNAs that showed dramatic changes in expression level at different developmental stages. The result (Fig. 5) showed that six miRNAs (asi-miR-276-3p, asi-miR-2943, asi-miR-1891, asi-miR-286b, asi-miR-1175-3p, and asi-miR-nov54) showed similar expression patterns as those revealed by our sequencing result. However, three miRNAs (asi-miR-263a-5p, asi-miR-nov23, and asi-miR-nov55) expression levels detected by qRT-PCR were inconsistent with that of our sequencing results due to some unknown reasons.

miRNA target gene prediction, GO enrichment, and KEGG pathway analysis
To understand the role of dysregulated miRNAs in physiological functions and biologic processes during development in An. sinensis, miRNA target gene prediction was performed to identify the putative target genes based on miRNA/mRNA complementarity. Multiple/different miRNA target prediction algorithms were applied to minimize the number of putative or false positive target genes. The target gene candidates of the differentially expressed miRNAs between two consecutive developmental stages of An. sinensis were predicted. A total of 7198 annotated mRNA transcripts were predicted as target gene candidates for 95 differentially expressed miRNAs (Additional file 3: Table S3). Among them, 5296 and 1902 target gene candidates were identified for 68 known and 27 novel miRNAs, respectively. GO term annotation and KEGG Pathway analysis were used to further understand biological functions and signaling pathways that could potentially be regulated by these putative target genes. We found that the putative target genes of the differentially expressed miRNAs between two consecutive developmental stages of An. sinensis were involved in a broad range of biological functions and signaling pathways (Fig. 6). Detailed classifications about the biological process, cellular component, and molecular function of the miRNA targets were shown in Fig. 6a. By GO annotation, targets of miRNAs regulated from the embryo to the larval stage in the mosquito were enriched in intracellular, cytoplasm and cellular process. Targets of miRNAs regulated from larva to pupal stage in the mosquito were enriched in membrane part, membrane, cytoplasm and organic substance transport. Targets of miRNAs regulated from the pupa to the adult stage in the mosquito were enriched in structural molecule activity, membrane part and transport. Similarly, by KEGG pathway analysis, we found that 119 enriched pathways were identified in this study. The most enriched pathway included protein processing in endoplasmic reticulum, propanoate metabolism, folate biosynthesis, Citrate cycle (TCA cycle) and valine, leucine and isoleucine degradation (Fig. 6b).

miRNA-gene network analyses identified key miRNAs and developmental genes
Through analysis, 16 commonly dysregulated miRNAs identified from the miRNA expression files, 15 exclusively expressed miRNAs in certain stages (TPM > 10), and three the most significantly over/under-expressed miR-NAs (>tenfold) were defined as the key miRNAs in An. sinensis development. A total of 546 developmental genes were obtained based on gene annotation in VectorBase and intensive literature surveys. After removing duplicates and invalid target genes, miRNA-gene network was built based on 32 key miRNAs and 104 developmental genes (Additional file 4: Table S4). We developed a network to demonstrate the overlapping miRNA targets for miRNAs since miRNAs usually have many mRNA targets. The number of targets for asi-miR-317, asi-miR-2944a-5p, asi-miR-263b-5p, asi-miR-263a-5p, and asi-miR-2943 was 10, 7, 6, 6 and 5, respectively  Fig. 7). This network also indicated that several mRNAs were predicted to be targets of one miRNAs. For example, AGAP001999, AGAP007696, AGAP005245, and AGAP004050 were targeted by two to five miRNAs, respectively (Fig. 7). As expected, the overall pattern suggests a positive correlation between the miRNA expression and developmental genes, such as AGAP000048 (APC), AGAP005449 (Cbl), AGAP001867 (hep), which encoding important components of major developmental signaling pathways.

Discussion
In the present study, the miRNA expression profile of An. sinensis in four different life stages was analyzed by NGS technology, and approximately 100 million high-quality reads were obtained. The mapping rates of high-quality reads to the An. sinensis genome were between 77.61 and 85.75%. This might be subject to incomplete reference genome. The size distribution of small RNA sequence reads showed peaks between 20 AND 30 nucleotides comparing to the average miRNA length. Annotation of the sRNAs in the four libraries indicated that the most abundant class of sRNAs in the eggs, larvae, pupae and adult females' libraries was miRNAs, rRNAs (larvae and pupae) and tRNAs respectively. This result is consistent with previous research in Anopheles funestus [37], indicating that protein synthesis during development is of great importance. In total, 99 known miRNAs and 65 novel miRNAs expressed in the four life stages of the An. sinensis were identified. As expected, most of the known miRNAs identified in An. sinensis were highly conserved across diverse mosquito species, suggesting that the conserved miRNAs may involve in many critically fundamental Fig. 4 Fold-change of miRNAs during the development of An. sinensis. a Log2 fold change of miRNAs between embryo and larval stage, b Log2 fold change of miRNAs between larva and pupal stage, c Log2 fold change of miRNAs between pupal and adult stage. X axis represents miRNAs while Y axis represents fold-change of miRNAs regulatory pathways. Novel miRNAs found in this study are not conserved in the Anopheles genus, and we cannot find homologs in other insect species by searching the miRBase, suggesting that these might be species-specific miRNAs. However, the predicted novel miRNAs exhibited much lower expression levels when compared with conserved miRNAs, which is also consistent with previous reports [9,37]. Previously, Dritsou and his colleagues [6] computationally predicted 115 miRNAs (after removing 33 duplicates and 11 putatively new miRNA, 71 left) in An. sinensis from the genome level. We compared our miRNA dataset with theirs and found that 33 miRNAs overlap with our predictions. Among these miRNAs, 19 miRNAs could be found in Ae. aegypti, An. gambiae, and Cu. quinquefasciatus. Six miRNAs (miR-193, miR-2765, miR-277-5p miR-2b, miR-315-5p and miR-9b) were only found in Ae. aegypti. In contrast, we validate more miR-NAs existence and detect their expression pattern at different developmental stages.
To obtain insight into possible stage-specific miR-NAs and their roles in An. sinensis, we investigated the expression profile of the miRNAs identified in the four life stage libraries. One-hundred and seventeen miRNAs were ubiquitously expressed in all developmental stages. The top five most highly expressed miRNA across four developmental stages were asi-miR-276-3p, asi-miR-1, asi-miR-100, asi-miR-9a and miR-8, which is basically in line with previous report in other mosquito species [15,38]. The most highly expressed miRNA in female An. sinensis was asi-miR-276-3p. In contrast, Skalsky et al. [15] reported that miR-184 was the most highly expressed miRNA in C7/10 Ae. albopictus cells and blood-fed female Cx. quinquefasciatus mosquitoes [15]. afu-miR-8 was the dominant miRNA in the unfed adult An. funestus females [37]. In addition, miR-1, miR-184, and miR-263 ranked among the most highly expressed miRNAs in both Ae. aegypti and An. stephensi [27]. aga-bantam, aga-miR-263a, aga-miR-8, aga-miR-10, aga-miR-184 and aga-miR-281 were the most abundantly expressed miR-NAs in An. gambiae [38]. Thus, despite the differences in expression level, most of these reported abundant miR-NAs in different mosquito lineages are still among the top 10 most highly expressed miRNA in An. sinensis.
Comparisons between egg and larva, larva and pupa, and pupa and unfed adult female are shown in Figs. 3 and 4. A majority of miRNAs were differentially expressed in the four stages. For example, 41 miRNAs were up-regulated and another 33 down-regulated between the embryo and larva stage. Especially, asi-miR-2943 and asi-miR-286b were up-regulated more than tenfold change. While between the larva and the pupa, 22 miRNAs were up-regulated and 20 other miR-NAs were down-regulated. We also identified 15 significantly up-regulated miRNAs and 22 significantly down-regulated miRNAs from the pupal to the adult stage. And asi-miR-1891 was significantly down-regulated in this process. Dramatic alternation of miRNA expression during these three stages signifies that these miRNAs may be involved in the regulation of metamorphosis during development. What's more, asi-miR-2943 was exclusively expressed in the embryo stage, which is also highly expressed in the embryos in both Ae. aegypti, An. stephensi [39] and An. anthropophagus [17], indicating its important roles in embryogenesis across many mosquito species. However, miR-2943 was present at low levels in both Cx. quinquefasciatus and Ae. albopictus C7/10 cells [15]. Similarly, miR-1891 was most abundantly expressed in Ae. aegypti and An. stephensi pupal stage [39], however, asi-miR-276-3p was the highest expressed miRNAs in An. sinensis pupae, although asi-miR-1891 also exhibited a high expression. These results indicate that, on the one hand, miRNA profiles of mosquito lineage may differ from each other, one the other hand, specifically expressed miRNAs or regulatory miR-NAs during developmental stages may also be different. In addition, sample preparation, detection methods and sensitivity may also lead to inconsistency between different studies.
Moreover, we observed that the expression level of miRNAs did not correlate with the fold change. Some miRNAs showed higher expression but with low fold changes. For instance, asi-miR-276-3p, asi-miR-1, asi-miR-100 and asi-miR-7 were among the top 10 most highly expressed miRNAs in all developmental stages but showed low log2-fold changes. However, asi-miR-193 and asi-miR-nov54 were low expressed miRNAs in An. sinensis developmental stages but showed high log2-fold changes. Additionally, asi-miR-309a, asi-miR-2943, asi-miR-315-3p and asi-miR-iab-4-3p could not be detected in the pupae and larvae libraries, asi-miR-275-5p and asi-miR-2491-3p was not identified in the pupae and adult females' libraries, and asi-miR-980-3p was not found from the pupae library. These results indicate that the expressions of some miRNAs are probably stage specific. We assumed that stage-specific miRNAs might be involved in regulation of embryogenesis, organ differentiation, growth, and reproduction during a specific developmental stage. Investigation on stage-specific miRNAs may detail understanding of mosquito biology and could provide mosquito-specific targets for disease control. However, the function of majority of the miR-NAs remains to be determined.
To confirm our sequencing results, we performed quantitative real-time PCR analysis using the RNA samples used for small RNA NGS. Thus, nine miRNAs that showed dramatic changes in expression level at different developmental stages were selected for validation. All nine miRNAs (asi-miR-276-3p, asi-miR-2943, asi-miR-1891, asi-miR-286b, asi-miR-1175-3p, asi-miR-263a-5p, asi-miR-nov23, asi-miR-nov54, and asi-miR-nov55) could be detected in the sample. Six miRNAs showed similar expression patterns as those revealed by our sequencing analysis, which indicated a low false discovery rate of our sequencing data and supported the validity of the profiling. However, the expression levels of asi-miR-263a-5p, asi-miR-nov23, and asi-miR-nov55 detected by qRT-PCR were inconsistent with that of our sequencing results. We inferred that this might be caused by the experimental conditions, such as base bias, amplification efficiency, and analytical error in both small RNA NGS and the qRT-PCR, but some other unknown reasons cannot be excluded.
Mosquitoes grow through four phases of metamorphosis: egg, larva, pupa and adult. Knowing the different stages of the mosquito's life will help us to prevent mosquitoes and disease transmission. To further understand the role of miRNAs in An. sinensis, miRNA target gene prediction was performed based on miRNA/mRNA interactions to provide molecular insight into their function by using two computational methods. A total of 7198 annotated mRNA transcripts were predicted as putative target genes for 95 differentially expressed miRNAs. The enriched biological functions and signaling pathways with significant role at a specific stage of development were predicted by GO term annotation and KEGG Pathway analysis. Metamorphosis from embryo to larval and from pupal to adult stage involves biological functions such as cellular component, membrane, cytoplasm and organic substance transport. In addition, enriched pathways involved in Protein processing in endoplasmic reticulum, Propanoate metabolism, Folate biosynthesis, Citrate cycle (TCA cycle) and Valine, leucine and isoleucine degradation were found. These biological functions and pathways might be important for degeneration and formation of tissue/organ, which happens during immature stages to complete its transition from egg to adult.
Network of miRNA-mRNA interactions was predicted by using results from miRNA analysis and mRNA retrieved from developmental genes of An. gambiae and then visualized by Cytoscape. Further analyses indicated that a number of miRNAs (asi-miR-317, asi-miR-2944a-5p, asi-miR-263b-5p, asi-miR-263a-5p, and asi-miR-2943) target a cluster of mRNA and vice versa. The candidate targets mentioned here have been reported to be key components in developmental pathways [31]. For example, AGAP002352, one of the targets of asi-miR-263b-5p, is an annotated apoptosis-related gene (p53) in Ae. aegypti and many other mosquitoes. p53 plays an important role in programmed cell death in relation to development [40]. Also, AGAP003108 (htl), AGAP001867 (hep) and AGAP005245 (lola) were the common target of asi-miR-281-5p. htl is one of two FGF (Fibroblast Growth Factor) receptors, and FGF signaling regulates a variety of biological processes during development, such as cell differentiation and cell proliferation [41]. hep is one component of Ras signaling pathways which also could regulate many developmental processes [42]. lola, which is called axon guidance gene that regulate wiring of the olfactory system in D. melanogaster [43], have been reported to have many orthologs in mosquito [31]. In addition, knock-in and knock-down of specifically and temporally expressed miRNAs in Ae. albopictus by microinjection have improved our understanding of regulated miRNAs in mosquito development [16]. Intervention on aal-miR-286b led to a profound decrease in the hatching rate of embryos. In addition, eclosion rate of larvae in aal-miR-2942 knocked-down group showed a significant decrease compared to control group. A marked reduction in longevity and fecundity was also observed in the miR-1891 inhibitor group. These research findings have established the exact functional roles these developmental genes and miRNAs in regulating the stage-specific development. However, additional extensive experimental analyses are required to confirm the interactions of miRNA-gene regulatory role in many phenotypic changes during developmental stages in further ongoing research work.
Although miRNAs are reported to play a role in development in C. elegans, fruit flies and other insects, and significant stage-specific expression was observed for a number of miRNAs in various mosquito species. However, the role of miRNAs in metamorphosis and development is far less understood and remains an enigma. For some specific or novel miRNAs, their functional roles still need to be confirmed by biological experiments. Furthermore, only very limited specific target genes have been identified for small part of miRNAs. It is important to note that the result in this work is a key step towards improving our understanding of the complexity and regulation mode of miRNAs in An. sinensis embryogenesis and metamorphosis. Further analysis of stage-specific miRNA expression and functions would be helpful to decipher the complex genetic network that controls mosquito at a crucial stage.
Anopheles sinensis is the most common and wide-distributed malaria vector exhibiting a variety of ecological behaviors and feeding habits throughout the China. It is also of great medical and veterinary importance. Thus, identification of key miRNAs required for gene regulation throughout the life cycle will not only help to a better understanding of its developmental biology but also may help to conceive novel approaches to control this vector. However, the function and molecular mechanism of these miRNAs are still need to be further investigated.

Conclusions
In conclusion, this study provides a comprehensive account of miRNA profile and its dysregulated expression patterns across four different stages of An. sinensis development using the high throughput sequencing. Several miRNAs show stage-specific expression pattern in particular developmental stages. We also investigate targets of these dysregulated miRNAs and their related functions in important life activities in order to provide a better understanding of mosquito development. In turn, such information would provide the basis to conceive novel approaches to control this vector in future.