miRNAsofAedesaegypti (Linnaeus 1762) conserved in six orders of the class Insecta

Aedes aegypti L. is the most important vector of arboviruses such as dengue, Zika, chikungunya, Mayaro, and yellow fever, which impact millions of people’s health per year. MicroRNA profile has been described in some mosquito species as being important for biological processes such as digestion of blood, oviposition, sexual differentiation, insecticide resistance, and pathogens dissemination. We identified the miRNAs of Ae. aegypti females, males and eggs of a reference insecticide susceptible strain New Orleans and compared them with those other insects to determine miRNA fingerprint by new-generation sequencing. The sequences were analyzed using data mining tools and categorization, followed by differential expression analysis and conservation with other insects. A total of 55 conserved miRNAs were identified, of which 34 were of holometabolous insects and 21 shared with hemimetabolous insects. Of these miRNAs, 32 had differential expression within the stages analyzed. Three predominant functions of miRNA were related to embryonic development regulation, metamorphosis, and basal functions. The findings of this research describe new information on Ae. aegypti physiology which could be useful for the development of new control strategies, particularly in mosquito development and metamorphosis processes.


Results miRNAs annotation.
A total of 21,632,581 reads were obtained; 40,815 showed similarity in the structure and biogenesis of 80 different miRNAs of Aedes aegypti. The number of readings for each biological stage is provided in Table 1. miRNAs taxonomic distribution. Of the 80 miRNAs annotated from Ae. aegypti, 55 showed significant conservation with 713 miRNAs from other insect species, while 25 showed Ae. aegypti-specific conservation, considering as cut-off an E value ≤ 0.005. According to the NCBI taxonomy information of these insect species, we identified that of the 713 miRNAs, 445 miRNAs (62.41%) were present in the order Diptera, 130 in Hymenoptera (18.23%), 84 in Lepidoptera (11.78%), 29 in Coleoptera (4.07%), 19 in Hemiptera (2.66%) and 6 in Orthoptera (0.84%) as described in Fig. 1 and Supplementary Material 1.

Discussion
Metamorphosis diversity in insects (ametabolous, hemimetabolous and holometabolous) is a fundamental part of their evolution, being a determining factor in the conservation of these organisms 15 . So far, different miRNAs, such as miR-let-7, miR-100 and miR-125, have been conserved and participate in metamorphosis in different organisms such as Blattella germanica Linnaeus, 1767 (hemimetabolous) 16 and Anopheles gambiae (holometabolous) [17][18][19] . We showed that 34 miRNAs corresponded to only one or more holometabolous organism, classified in the orders Diptera, Hymenoptera, Lepidoptera and Coleoptera and that 21 miRNAs shared conservation with some hemimetabolous insects.
Another portion contains 13 miRNAs (Fig. 5A, purple section; Supplementary Table S1); six of these showed overexpression in eggs [ Fig. 5C, second bar "Eggs (Up)", (Supplementary Table S1)], but they are not related to specific biological events as yet. We observed miR-10 and miR-71-5p overexpression in holometabolous and hemimetabolous insects. Only five have been described previously in D. melanogaster embryos playing developmental functions: miR-307 regulated cell death processes through the silencing of wntless 6 ; miR-33 in cell proliferation events by silencing the genes CDK6 y CCND1 31 ; miR-9a in muscle myofibril development by silencing troponin-T 32 ; and miR-9c-5p in the cleaning of unstable maternal transcripts 33 . Finally, miR-iab-4-5p www.nature.com/scientificreports/  www.nature.com/scientificreports/ regulating the distribution and development of the wings and halters through the silencing of Ultrabithorax 34 . Thus, we can suggest embryonic development function for these 5 miRNAs in Ae. aegypti. Five miRNAs (miR-1891, miR-286b, miR-2944a-5p, miR-2944b-5p and miR-996) were found conserved only in the members of the family Culicidae (Fig. 5B, bar 1 "Culicidae", Supplementary Table S1); only two www.nature.com/scientificreports/ (miR-1891 and miR-996) have been reported overexpressed in eggs of Ae. aegypti, Anopheles gambiae and An. stephensi Liston 1901 35,36 , and three (miR-286b, miR-2944a-5p, miR-2944b-5p) were highly expressed in the early development stages of mosquitoes 35,37 . The family formed by miR-996, miR-279, and miR-296 have undergone some duplication events along insect evolution, making it diversify in holometabolous organisms 38 . Thus, we suggested that miRNAs conserved specifically in Culicidae could have undergone similar events that would have led to miRNomic regulation in embryonic stages in a specific way. Discrepancies were found between previous studies and our findings regarding the following miRNAs: miR-9a has been reported in embryo development function in D. melanogaster eggs 32 , but in contrast, no significant changes were found in the biological approaches that we used (eggs, females and males); miR-71-5p with overexpression in Manduca sexta (Linnaeus 1763) eggs 39 was in contrast to the expression in males in our findings; miR-10 was overexpressed in D. melanogaster, T. castaneum, Apis mellifera (Linnaeus 1758) and B. germanica eggs 40,41 and downregulated in females according to our results. Thus, as explained above, we can infer new and different roles for these miRNAs in Ae. aegypti.
Twelve miRNAs conserved between Ae. aegypti and the class Insecta were categorized in diverse basal functions such as neuronal signaling, hypoxia functions, nervous system, and circadian cycle (Fig. 5A, blue section; Supplementary Table S1). The functions of these miRNAs are: miR-1, regulating the architecture and maintenance of D. melanogaster cardiac tissue through DELTA silencing 21 ; miR-137, in neuronal signaling of Drosophila Parkinson disease models 42 ; miR-190, orchestrating the functions related to hypoxia in D. melanogaster 43 ; miR-263b-5p and miR-279, in the circadian cycle in D. melanogaster 44 ; miR-2c and miR-31, in the maintenance of the nervous system of D. melanogaster 45,46 and miR-957, in the interaction Cx. quinquefasciatus/Western Nile virus (WNV) 47 ; miR-9b, silencing serine/threonine kinase protein in Drosophila 48 ; miR-932-5p, in the neuroplasticity mechanisms and long-term memory of A. mellifera 49 ; miR-305-5p, in the intestinal homeostasis of D. melanogaster 50 ; and miR-8-5p, in post-embryonic development 51 and growth factor regulation of body fat in larval stages 40 in D. melanogaster.
Of these 12 miRNAs, 9 showed a basal expression in the life stages of Ae. aegypti analyzed in the present study (Fig. 5C, bar 1 "Basal exp. "; Supplementary Table S1) and the remaining 3, were downregulated in eggs (miR-932-5p) or upregulated (miR-31) and downregulated (miR-305-5p) in females. Conservation profiles of the 12 miRNAs showed 4 with conservation only in holometabolous, 6 in holometabolous and hemimetabolous and 2 specifically in Diptera (Fig. 5B). Of the basal function of the miRNAs with holometabolous and hemimetabolous conservation, miR-2c and miR-8-5p have different biological functions depending on the type of metamorphosis that was explored. For the hemimetabolous functions for these miRNAs, miR-2c was related to metamorphosis processes through the silencing of KR-H1 in B. germanica 28 and miR-8-5p in chitin biosynthesis in Locusta migratoria Linnaeus 1758 52 and motor coordination of B. germanica 53 .
Other biological functions were observed (sex, blood digestion and detoxification) that despite having a low abundance are relevant in the biology of Ae. aegypti. miRNAs were found to be related to sexual development, where 2 were found to be unique for Culicidae (miR-2942 and miR-989) (Fig. 5B, bar 1 "Culicidae"; Supplementary Table S1), which regulates ovary development in Ae. aegypti, Cx. quinquefasciatus and An. stephensi 14,[54][55][56] .
On the other hand, miR-193 was shown to be conserved in the order Diptera (Fig. 5B, bar 2 "Diptera"; Supplementary Table S1), being overexpressed in D. melanogaster 57 . Also, the experimental design did not consider sex as a factor for expression analysis, so our results pointed to be a possible function in feminization. The battery of miR-7 and miR-993 was conserved in the holometabolous organisms (Fig. 5B, bar 3 "Ho"; Supplementary Table S1), where it has been observed that they fulfill masculinization roles in D. melanogaster 58,59 . Another battery of miR-184 and miR-100 was conserved in holo-and hemimetabolous insects (Fig. 5B, bar 4 "Ho + He"; Supplementary Table S1), both being highly expressed in males of An. anthropophagus Xu & Feng 1975 60 . miR-184 have been associated with regulation of ovary development in D. melanogaster 21 and with high expression in females of the 6 th nymphal stage in B. germanica 30 and miR-100 with essential metamorphosis regulation conserved among holometabolous and hemimetabolous insects 28 . Our results suggest that the feminization of mosquitoes may require a miRNomic regulation conserved uniquely in Culicidae.
Finally, only miR-285 was found to be associated with detoxification functions. Previously, it has been described in a resistant strain (91-R) of D. melanogaster in exposure to DDT 61 . This agrees with our results of overexpression in females because resistant strains were not analyzed in this study, and we used the susceptible reference strain New Orleans. However, expression findings infer that miR-285 plays an essential role in the development of insecticide resistance in a specific way in Diptera.
In conclusion, the sequencing of miRNAs in eggs and adults (males and females) of Ae. aegypti showed a total of 55 miRNAs with relevant conservation with other members of the class Insecta. Of these miRNAs, 32 showed DE in the stages analyzed, showing to be related to development, metamorphosis, and sex functions, among others. On the other hand, 12 miRNAs with basal expressions were shown to be mainly related to diverse biological functions. While a lower abundance of conservation was expected in specific functions, such as detoxification, blood digestion, and immune response, there was a predominance of highly conserved miRNAs related to embryonic development and metamorphosis processes or a diversity of conservation clusters in miRNAs related to sexuality. These results suggest that the development of wings, with miR-iab-4-5p and masculinization, with the battery of miR-7 and miR-993, could be ancestral miRNomic mechanisms.
We propose to coin the word "phylorthology" to refer to the phylogenetic relationship that species have with relation in a conserved sequence in 100%. www.nature.com/scientificreports/

Methods
General strategy. The general strategy consisted of obtaining miRNA sequences from 3 life stages in Ae.
aegypti (males, females, and eggs). Total RNA was extracted, and the corresponding miRNA fractions were sequenced. These sequences were annotated and subjected to conservation analysis against the miRNAs available in the class Insecta. Those with significant conservation (E value ≤ 0.005) were subjected to interaction graphics and cladograms, to observe their abundance and diverse patterns of conservation between the organisms of the class Insecta. Finally, a DE analysis of these miRNAs was evaluated using a likelihood ratio statistical analysis.
Sample collection. Ae. aegypti specimens belonging to the New Orleans strain were bred at the insectarium of the Medical Entomology Laboratory of the Facultad de Ciencias Biologicas in the Universidad Autonoma de Nuevo Leon. Breeding of the mosquitoes was performed under laboratory conditions with a temperature of 28 ± 2 °C, a light/dark photoperiod of 12/12 h and relative humidity of 70 ± 2%. A total of 3000 eggs, 100 males and 100 females were separated, and miRNA extraction from the specimens was carried out with the commercial kit miRNeasy Mini, following the provider instructions (Qiagen, Germantown, MD, USA). RNA extracted was treated with RQ1 RNase-Free Dnase (Promega, San Luis Obispo, CA, USA) to remove genomic DNA traces. Standard spectrophotometry (Thermo Fisher, Waltham, MA, USA) and agarose gel electrophoresis were used to evaluate RNA purity and integrity, respectively.
Next generation sequencing. RNA  Differential expression analysis. For the expression analyses, only miRNAs that showed conservation were considered; the total reads were obtained by the mean of each of the approaches explored. A Venn diagram was made to show the dispersion of these miRNAs between the 3 different biological approaches (eggs, females and males) and to discard the miRNAs with specific expression in only one explored biological approach. The mean of the readings was put through depuration and normalization protocols according to Law et al. 73 , where data expression of the readings was transformed into counts per million (CPM). This value was expressed in log2 (log2-CPM) and was normalized using Trimmed Mean of M values (TMM) method. For the determination of the miRNAs with DE in the three contrasted approaches (eggs vs females, eggs vs males and, females vs males), a generalized linear model and a likelihood ratio likelihood (LRT) were used 74  www.nature.com/scientificreports/