Profiling of Transcriptome-Wide N6-Methyladenosine (m6A) Modifications and Identifying m6A Associated Regulation in Sperm Tail Formation in Anopheles sinensis

Recent discoveries of reversible N6-methyladenosine (m6A) methylation on messenger RNA (mRNA) and mapping of m6A methylomes in many species have revealed potential regulatory functions of this RNA modification by m6A players—writers, readers, and erasers. Here, we first profile transcriptome-wide m6A in female and male Anopheles sinensis and reveal that m6A is also a highly conserved modification of mRNA in mosquitoes. Distinct from mammals and yeast but similar to Arabidopsis thaliana, m6A in An. sinensis is enriched not only around the stop codon and within 3′-untranslated regions but also around the start codon and 5′-UTR. Gene ontology analysis indicates the unique distribution pattern of m6A in An. sinensis is associated with mosquito sex-specific pathways such as tRNA wobble uridine modification and phospholipid-binding in females, and peptidoglycan catabolic process, exosome and signal recognition particle, endoplasmic reticulum targeting, and RNA helicase activity in males. The positive correlation between m6A deposition and mRNA abundance indicates that m6A can play a role in regulating gene expression in mosquitoes. Furthermore, many spermatogenesis-associated genes, especially those related to mature sperm flagellum formation, are positively modulated by m6A methylation. A transcriptional regulatory network of m6A in An. sinensis is first profiled in the present study, especially in spermatogenesis, which may provide a new clue for the control of this disease-transmitting vector.


Introduction
China was certified as a malaria-free country on 30 June 2021 (https://www.who.int/ news/item/30-06-2021-from-30-million-cases-to-zero-china-is-certified-malaria-free-by-who, accessed on 30 June 2021), but the risk of re-establishment of malaria transmission in China cannot be neglected, mainly due to the fact that nearly 3000 malaria cases imported yearly [1] and the widespread distribution of malaria mosquito vectors, especially the most widely distributed species of Anopheles sinensis [2,3], which is considered to be a competent vector for Plasmodium vivax [4], and a potential vector for P. falciparum according to the positive susceptibility experiments in the laboratories [2,4]. Moreover, introduced vivax malaria cases have been reported in Liaoning province, Hunan province, and Yunnan border areas, respectively, which were speculated to be the secondary transmission caused by imported vivax malaria cases from Southeast Asia [5][6][7]. In addition, An. sinensis is also the vector that transmits pathogens severely impacting global human health, such as lymphatic filariasis [8,9], Japanese encephalitis virus [10], and Rickettsia felis [11]. However, available tools are not sufficient for mosquito control, which is the mainstay approach for 2 of 21 diseases prevention and control. Thus, new strategies instead of traditional measures to control vectors are urgently needed [12].
The function and mechanism of m6A in modifying Drosophila mRNA have been studied in insects [14][15][16][17][18]. For instance, Drosophila retained the conserved writer complex (Ime4, Mettl14, Vir, Nito, Fl(2)d) and YTH-domain family proteins (YT521-B, CG6422), but no corresponding demethylases have been identified so far [30]. Flightless and neuronal function defects of Drosophila resulted from the knock-out of Ime4 were verified that brain functions were regulated by neuronal mRNA with m6A modification [16,17,23]. Furthermore, Ime4 and Mettl14 can also regulate the sex determination of Drosophila by controlling the splicing of the Sxl [15,16,18]. Recently, Drosophila Ime4 was found to be potentially involved in the establishment and maintenance of the somatic cyst-cell permeability barrier in spermatogenesis [20]. In other insects, m6A was found to regulate insecticide resistance in Bemisia tobaci [21] and diapause-related genes in silkworms [22] and impact caste differentiation and larval development in honey bees [20,33]. However, the biological roles of m6A mRNA in other insects, including mosquitoes, remain largely unknown.
To further investigate the functions of m6A on mRNA and to facilitate future studies of m6A in mosquitoes, we first report here transcriptome-wide m6A profiling in female and male An. sinensis. This study confirmed that m6A is a highly conserved RNA modification on mRNA across female and male An. sinensis. Intriguingly, m6A in An. sinensis is enriched not only around the stop codon and within 3 UTRs, such as that in yeast and mammalian systems, but also around the start codon, as that in Drosophila and Arabidopsis thaliana. Furthermore, a positive correlation between m6A deposition and mRNA levels indicates a regulatory role of m6A in the gene expression of female and male mosquitoes. More importantly, m6A is found to play potential roles in spermatogenesis, especially in the sperm tail formation of An. sinensis.

Conserved mRNA m6A Toolkit in Anopheles sinensis and Other Anopheles Species
Putative orthologs for the writer complex and YTH domain-containing proteins (Figure 1a, Supplementary File S1) were found in An. sinensis through a sequence-similarity approach using reciprocal best BLAST. Moreover, the conserved domain and phylogenetic analyses (Supplementary File S1) ascertained these orthology relationships. Anopheles spp. retains a core methyl transferase complex (METTL3, METTL14, and WTAP) and other components, such as HAKAI, RBM15 and VIRMA, and YTHDC and YTHDF as potential m6A readers (Figure 1b). However, no corresponding eraser enzymes (FTO and ALKBH5) were identified in Anopheles spp. genomes (Figure 1a,b, Supplementary File S1).
There was no significant difference in gene expressions of core methyltransferase complex between female and male An. sinensis, while that of VIRMA and YTHDCs were significantly lower in the males compared to the females (Figure 1c). Although no ALKBH5 and FTO were identified in Anopheles spp., the gene expression levels of other ALKBHs demethylases, which can potentially compensate for the absence of ALKBH5 and FTO, were significantly lower in the males than the females (Figure 1c). There was no significant difference in gene expressions of core methyltransferase complex between female and male An. sinensis, while that of VIRMA and YTHDCs were significantly lower in the males compared to the females (Figure 1c). Although no ALKBH5 and FTO were identified in Anopheles spp., the gene expression levels of other ALKBHs demethylases, which can potentially compensate for the absence of ALKBH5 and FTO, were significantly lower in the males than the females (Figure 1c). The existence of m6A modification on mRNA of female and male An. sinensis was preliminarily identified by mass spectrometry (Figure 2a). To further obtain the transcriptomewide m6A map, both female and male An. sinensis were interrogated using m6A-targeted antibody coupled with high-throughput sequencing. In total, 6756 m6A peaks (fold change ≥ 2.0) representing the transcripts of 4981 genes in females, and 6289 m6A peaks representing the transcripts of 4664 genes in males, were identified, respectively (Figure 2b). Among them, 4014 m6A tagged genes were detected within both female and male An. sinensis, most of which were related to the oxidation-reduction process, membrane, and protein binding. Moreover, there were 967 female-specific m6A peak genes related to metal ion binding and ATP binding and 650 male-specific m6A peak genes involved in the regulation of transcription and DNA template (Figure 2b). Additionally, in An. sinensis, there were abundant genes with sex-specific m6A peak without known functions, indicating their species-and sex-specific functions associated with m6A modification. Above all, it was estimated that the An. sinensis transcriptome contained 0.7-0.8 m6A peaks per expressed transcript (active expressed transcript FPKM ≥ 2) (Table S1). After analyzing the distribution of m6A in the whole transcriptome for both female and male m6A-IP and non-IP (input) samples, it was found the m6A-IP reads were

m6A Peaks Positions in Female and Male An. sinensis
After analyzing the distribution of m6A in the whole transcriptome for both female and male m6A-IP and non-IP (input) samples, it was found the m6A-IP reads were highly enriched around the 5 UTR, start codon, stop codon, and within 3 UTR (Figure 2c). It was also found the peak callings of female samples were distributed almost equally in transcription start (TSS) and end (TES) sites, while more peak callings were detected in the TES regions than in TSS regions in male An. sinensis (Figure 2c). Moreover, the global hypermethylation of m6A was induced in male An. sinensis compare to females. The subgroups of genes were classified according to their detailed m6A peaks positions (Figure 2d). Most of the m6A sites were found at 3 -UTR (35.08-38.76%), exon (32.92-36.67%), and 5 -UTR (7.12-8.67%). In addition, there were m6A peaks at both 5 UTR and 3 UTR in the females (8.67%) and the males (7.12%).

m6A Peaks Motifs in Female and Male An. sinensis
Top 1000 significant peaks were used to determine whether the m6A peaks in female and male An. sinensis contained the m6A consensus sequences. A core motif sequence of GGACAAGGAGG was identified in the females, showing the conserved pattern of the RRACH motif ( Figure 2e). Meanwhile, a specific GAA/CGAAGA/CAG motif was observed in the male m6A peak genes ( Figure 2e).

Gene Expression and Important Biological Pathways of m6A-Containing mRNAs Genes
The expression levels of genes with m6A modification were significantly higher than that of genes without m6A modification in both female and male An. sinensis (Kruskal-Wallis test, p < 0.001) (Figure 3a,b). Moreover, the expression of m6A modified genes at UTR regions had higher levels than those at exon regions or without m6A modification, respectively (Kruskal-Wallis test, p < 0.001) (Figure 3a   (c) GO enrichment of sex-specific m6A genes in female and male An. sinensis. Female-m6A, m6A genes in female; Male-m6A, m6A genes in male; female high m6A genes, genes with m6A sites ≥ 3 in females; Male high m6A genes, genes with m6A sites ≥ 3 in males.
Then, the enriched GO (gene ontology) was further predicted using the genes containing m6A in female and male mosquitoes to uncover potentially functional insights about m6A in An. sinensis (Figure 3c). It was found that these genes were enriched in tRNA wobble uridine modification, regulation of GTPase activity, basement membrane, signal recognition particle, and phospholipid-binding in the females, while they were enriched in peptidoglycan catabolic process, regulation of GTPase activity, basement membrane, signal recognition particle, exosome and signal recognition particle, endoplasmic reticulum target-ing and RNA helicase activity in the males ( Figure 3c). Importantly, in female An. sinensis, the high m6A-containing genes were associated with the establishment or maintenance of cell polarity and cohesion complex. In male An. sinensis, the high m6A-containing genes were involved in negative regulating of cyclin-dependent protein serine/threonine kinase activity, actin filament capping, spectrin, cohesion complex, and cyclin-dependent protein serine/threonine kinase inhibitor activity (Figure 3c).

mRNA m6A Peaks with Different Abundance between Male and Female An. sinensis
A different intensity was also found in m6A peaks in female and male An. sinensis ( Figure 4a). Among them, 88 down-regulated m6A peaks in 86 genes and 55 up-regulated m6A peaks in 45 genes were found in the males (p < 0.05. FC ≥ or ≤ 2). The peaks in these genes were mostly identified in the exon, especially in other exons rather than the 1st exon ( Figure 4b). GO analysis indicated that the genes with dynamic m6A peaks were enriched in several categories of fundamental biological functions; for example, cilium, cell projection, and cilium assembly were regulated by genes with up-regulated m6A peaks, and cysteine biosynthesis, cystathionine beta-synthase activity were regulated by genes with down-regulated m6A peaks ( Figure 4c).

Different Gene Expression in Male and Female An. sinensis
A total of 3369 differentially expressed genes were identified between male and fe-

Different Gene Expression in Male and Female An. sinensis
A total of 3369 differentially expressed genes were identified between male and female An. sinensis, and 2091 genes were up-regulated, and 1278 were down-regulated, respectively (male vs. female) (Figure 5a,b). Among them, some sex-specific genes were identified, such as Heme peroxidase, Envelysin, Yippee-like, microtubule-associated protein, RP/EB family, Abhydrolase domain-containing protein, and Niemann-Pick Type C-2 were only actively expressed in female mosquitoes, while adenylate cyclase, leucine-rich repeatcontaining protein, serine kinase, homeobox protein Nkx, invertebrate, dynein light chain 1, were specifically identified in male mosquitoes (Figure 5a). Up-regulated differentially expressed genes were dominated in the pathways of Indole alkaloid biosynthesis, betalain biosynthesis, ribosome and ribosome biogenesis in eukaryotes, while down-regulated differentially expressed genes were in the ribosome, ribosome biogenesis, and histidine metabolism ( Figure 5c). In addition, up-regulated expressed genes were highly enriched in biological processes of rRNA procession, microtubule-based movement, neuropeptide signaling pathway, ribosomal large subunit biogenesis, and regulation of ion transmembrane transport; in cellular components of the ribosome, cilium, cell projection, motile cilium; and in molecular functions of transmembrane transporter activity and ion channel activity, respectively ( Figure 5d). The down-regulated expressed genes were highly enriched in biological processes of DNA replication, methylation, pseudouridine synthesis, and protein import into the nucleus; in cellular components of small ribosomal subunit, small-subunit processome, exosome (RNase complex), mismatch repair complex, and nuclear pore; and in molecular function of nucleotide-binding (Figure 5d).  After analyzing the transcriptome sequencing with or without an m6A marker, 54 differentially expressed genes with different m6A abundances in female and male An. sinensis (male vs. female) were identified (Figure 6a). Among them, 22 genes displayed m6A hypermethylation and up-regulated expression (hyper-up), 9 genes had hyper-methylation and down-regulated expression (hyper-down), 6 genes displayed m6A hypo-methylation and up-regulated expression (hypo-up), and 17 genes had hypo-methylation and downregulated expression (hypo-down) (Figure 6a, Table 1). Notably, hyper-up genes were mainly involved in the oxidation-reduction process, cilium assembly, cell projection, cytoskeleton, membrane and cilium, ATP/protein/nucleotide-binding (Figure 6b), while hypo-down genes were enriched in DNA repair-related process in the nucleus, with binding to DNA and ATP (Figure 6c).  The 54 differentially expressed genes with different m6A peaks were carefully annotated with corresponding information from model organisms (human, D. melanogaster, and C. elegans). Interestingly, it was noticed that 14 up-regulated genes involved in spermatogenesis, especially in the development of sperm flagellum, were potentially regulated by m6A modification on mRNA (Figure 6a,d).
m6A Associated with Spermatogenesis of Male An. sinensis Compared to the mammalian spermatogenesis with three stages, mosquitoes retained the conserved mitosis and meiosis process but had a unique spermiogenesis process undergoing eight stages to form the mature sperm (Figure 7a) [34]. Recently, m6A was reported to be involved in the entire process of spermatogenesis, and many m6A-target transcripts have been identified in male germ cells [35][36][37][38][39]. In the present study, the orthologs of genes coding these transcripts were also identified in An. sinensis genome. Importantly, many of these genes were expressed differentially between female and male mosquitoes, and there were m6A modifications on their mRNAs (Table S2, Figure 7a). However, no significant difference between male and female An. sinensis was found in the m6A peaks abundances in most genes.    Unlike mammalians', mosquitoes' sperm is longer and more slender, with its head (nucleus) as wide as the tail (flagellum) [40,41]. The flagellum consists of two mitochondrial derivatives (MD), which extend the length of the flagellum, and an axoneme, a microtubular structure responsible for motility. Mosquito sperm axonemes have 19 microtubules arranged in two concentric circles of nine around a long central tubule (9 + 9 + 1) [40,41] (Figure 7b). Many mammalian sperm-tail-associated proteins in the main structures of sperm, such as axoneme, center pair (CP), outer dense fibers (ODFs), mitochondrial sheath (MS), HTCA (head tail coupling apparatus), and annulus, have been extensively documented, while no protein in axoneme and mitochondrial derivatives of mosquitoes have been systematically studied (Figure 7b).
Here, the orthologs of these functional proteins were identified in An. sinensis genome. Most of these genes were highly expressed in male mosquitoes and contained m6A modification on their mRNA. Particularly, significantly different m6A peak abundances in CCDC40, CFAP157, TEKTINS, TSSK2, Dynein fα, Tctex2, and DRC11 in the males were correlated with their high expression levels (Figure 7b,c, Tables 1 and S3). Moreover, many male-specific or highly expressed genes identified using the sperm proteomics of Aedes aegypti [42] were also tagged with m6A modification in An. sinensis (Table S4). Among them, arginine kinase responsible for maintaining constant ATP level to support the movement of sperm flagellum, and mitochondrial isocitrate dehydrogenase subunit without known functions in mosquitoes, were potentially regulated by the higher m6A modified mRNA in the males, while CFAP57, DUF4746 domain-containing protein, and BTB/POZ domaincontaining protein KCTD12 showed lower m6A modification levels (Figure 6a, Table S4). In addition, it is noticed that kelch-like protein 10 (Klhl10), as a substrate-specific adapter of a CUL3-based E3 ubiquitin-protein ligase complex [43], was also highly expressed in male mosquitoes and tagged with m6A modification (Figure 6a).

Discussion
Mosquito-borne diseases such as malaria, dengue, and Zika are spread by the bite of an infected female mosquito rather than the harmless nonbiting males; thus, control methods biasing the sex ratio of their offspring by reducing the number of blood-sucking females or converting them into harmless males, have long been sought [44][45][46][47]. Genetic elements such as sex-chromosome drives can distort sex ratios to produce unisex populations that eventually collapse, and advanced technological breakthroughs at driving maleness to achieve efficient sex-separation, the ultimate disease refractory phenotype [48][49][50], become possible and may represent efficient and self-limiting methods instead of the traditional vector control measures to effectively reduce and control the target mosquito populations, but the underlying molecular signals and mechanisms in illustrating the entire sex development pathway in mosquitoes remain unknown [50]. The identification of m6A mechanism components, the existence of m6A on mRNA by mass spectrometry, and the transcriptomewide map of a mosquito species An. sinensis presented here provides a starting roadmap for uncovering m6A functions that may affect/control mosquito development in the future. The m6A players and m6A target genes, especially the ones that play important roles in sperm tail, might be the potential target for further mosquito and disease control.
The m6A mRNA modification originated from the last eukaryotic common ancestor and has been considered a conserved mechanism for regulating the gene expression and translation in eukaryotes [14,[25][26][27]. However, the existence and biological roles of m6A on mRNA in most insects, including mosquitoes, remain largely unknown. Here, this modification was first revealed in mosquitoes. Firstly, the identification of putative orthologs of m6A writer complex and YTH domain-containing proteins in Anopheles spp. suggested that mRNA m6A machinery was conserved and important in Anopheles spp. Noticeably, An. sinensis retains most of the m6A players, including writer complex and readers. Consistent with Drosophila [30], no ALKBH5 and FTO were identified in Anopheles spp., but the presence of other ALKBHs demethylases might potentially compensate for the absence of ALKBH5 and FTO. Importantly, the existence of m6A RNA modification was further confirmed by LC-MS in both male and female An. sinensis. Moreover, it was estimated that the An. sinensis transcriptome contains 0.7-0.8 m6A peaks per expressed transcript, which is similar to that obtained in mammals and plants [25,26]. Above all, the m6A mechanism is highly conserved in mosquitoes and potentially plays a role in the biological functions of mosquitoes.
Also, the profiling of the transcriptome-wide m6A distributions in An. sinensis further uncovered the m6A modification characteristics, including m6A positions on mRNA, the relationships between m6A distribution and gene expressions, and regulation of the sexspecific functions. In contrast to mammals and yeast [25,27], m6A on An. sinensis mRNA was not only highly enriched around the stop codon and at 3 UTR, but also around the start codon and at 5 UTR, as that has been found in Drosophila [15] and A. thaliana [26], and they are bona fide m6A enrichment peaks rather than m6Am at the cap of mRNA [25][26][27]. It was surprising to find that the m6A pattern in An. sinensis was distinct from mammals and yeasts but resembled A. thaliana, which is evolutionarily distinct from Amorphea [51]. This evidence shows diverse m6A patterns in different insects [15,52], as well as those in plants [26,53]. It also indicated that although m6A is a conserved mechanism for regulating the gene expression and some aspects of biological processes, the development and realization of this mechanism are different in different species. Therefore, it is necessary to reveal the fine-tuning mechanism of m6A in distinct organisms in the future.
The present study showed that the expression levels of genes with m6A modification were significantly higher than genes without m6A modification in both female and male An. sinensis, indicating m6A is a way to positively regulate the gene expression. Particularly, similar to Bombyx mori [52] and A. thaliana [26], the m6A peaks at 5 UTR and 3 UTR were positively correlated with the overall up-regulation of mRNA expression levels in both female and male An. sinensis. However, this relationship contradicted the findings of a negative correlation between gene expression and m6A methylation around the stop codon and at 3 UTRs in mammalian, yeast, and honey bee systems [19,25,27]. Hence, it is suggested that m6A reader proteins may play roles in recognizing m6A at 5 UTRs or 3 UTRs and subsequently affect the stability of the target mRNA and then directly impact translation through the methylation itself or through the reader proteins [26].
Furthermore, the correlation of the sex specific functions with the distribution and density of m6A methylation in An. sinensis transcriptome was also analyzed. Firstly, the m6A peak density was more abundant in male than in female An. sinensis. Given that no differences in gene expression levels of m6A core writer complex were found in female and male An. sinensis (Figure 1c), the distinct m6A abundances may be contributed by other differentially expressed m6A components, such as VIRMA and potential ALKBH demethylases, because their lower expression levels were positively correlated with the higher abundance of m6A in male mosquitoes. Secondly, different m6A motif patterns in female and male An. sinensis add the complexity of m6A modification in regulating sexrelated functions. GO and KEGG analysis of the genes with m6A modification indicated its critical roles in regulating mosquito development, especially in the sex-related respect. For instance, the transcripts with abundant m6A were highly enriched in cilium-related genes in the males, while m6A in genes associated with the cysteine biosynthetic process were down-regulated. Intriguingly, after crosslink analysis of m6A peak abundances and gene expression levels, it was found that, out of 54 target genes, there were 14 genes associated with spermatogenesis, strongly indicating the relationship between m6A modification and up and down-regulated gene expressions involved in spermatogenesis.
Spermatogenesis is a highly sophisticated and complex process that can be divided into three stages: mitosis, meiosis, and spermiogenesis [35]. Mosquitoes retained the conserved mitosis and meiosis process but have a unique spermiogenesis process undergoing eight stages to form the mature sperm (Figure 7a) [34]. m6A modification has been reported to be involved in mammalian spermatogenesis [27,31,32,35,37,38,[54][55][56][57][58][59][60]. For instance, m6A methyltransferase complexes have been verified to play important roles in mammalian spermatogonial stem cell differentiation, meiosis, and spermiogenesis [27,35,37,54,55]. Sperm motility in the mice was significantly reduced after double-knockout of METTL3 and METTL14, accompanied by flagella defects and abnormal sperm head [38]. Moreover, the highly expressed FTO in the mammalian testis and its relationship with male fertility indicates FTO plays an important role in spermatogenesis [31,32,56,57]. Additionally, YTH domain proteins can also affect spermatogenesis by controlling germline transition into meiosis, participating in the normal function of spermatogenic tubules, and regulating spermatogonia migration and proliferation [58,59]. Recently, many genes were identified to be under the regulation of m6A modification in spermatogonial differentiation, meiosis, spermiogenesis, and other processes [32,37,38,56,60]. These genes with m6A modification were also identified in the An. sinensis transcripts, suggesting their potential roles in mosquito spermatogenesis. However, the m6A peak abundances were not significantly different between male and female adult mosquitoes, as well as the expression levels. It may be attributed to the physiological process that spermatogenesis is mainly developed in mosquito larvae rather than in adult males (Figure 7a), the development stage of which we focused on in this study. However, in adult male An. sinensis, m6A modification was found to be abundant in sperm tail associated proteins, which are involved in sperm tail formation and specifically expressed in mature sperm. Even though the sperm structures between mammals and mosquitoes are different [40,41], the sperm-related genes in An. sinensis identified by bioinformatics methods showed the conservation and complexity of sperm genes in mosquitoes. Meanwhile, the identification of spermatogenesis-related m6A peaks in genes indicates that their transcript, splicing, translation, and storage are potentially regulated by m6A modification. Additionally, the high expression levels of some functional genes in mature sperm tails are correlated with their abundant m6A modification on mRNA.

Mosquito Rearing and Sample Preparation
Anopheles sinensis (China strain) were reared at 28 ± 2 • C, 70-75% humidity in the insectary laboratory at the National Institute of Parasitic Diseases, Chinese Center for Disease Control and Prevention. Adult mosquitoes were maintained in screened cages and provided constant access to water and glucose-soaked sponges.
Male and female adults were identified morphologically and collected after emerging 3 to 4 days. All samples were flash-frozen in liquid nitrogen immediately following collection, and then stored at −80 • C until RNA isolation.

Phylogenetic Analysis of m6A Associated Proteins in Anopheles spp.
A total of 22 human m6A-associated proteins sequences (Table S5) were used as initial queries to search homologues in 22 Anopheles spp. genomes (https://vectorbase. org/vectorbase, accessed on 24 September 2021) with reciprocal Blast method, with an E-value cut-off < 0.05. The amino acid residues were visualized with Jalview Version 2 [61]. Sequences were aligned using MUSCLE (https://www.ebi.ac.uk/Tools/msa/muscle/, accessed on 26 September 2021) with the default parameters. Human and Drosophila were used as the out-group. These alignments were trimmed with TrimAI using the heuristic automated1 method [62]. Maximum likelihood (ML) analyses were performed using online IQ-TREE platform [63]. The best-fitting model was defined by IQ-TREE with Bayesian Information Criterion. The tree branches were tested with ultrafast bootstrapping (1000) and SH-like approximate likelihood ratio test (SH-aLRT, 1000 replicates). The final trees were visualized with TreeGraph2 [64]. The conserved domains were identified with NCBI Conserved Domain Search [65].

RNA Extraction and qPCR
Total RNA of male and female adult Anopheles sinensis was extracted separately using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions. The RNA amount and purity of each sample were quantified using NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA). The integrity of the RNA was determined using the Agilent BioAnalyzer 2100 (Agilent Technologies, USA). Reverse transcription for qPCR was performed with PrimeScript RT master mix (Takara, Japan). Quantitative PCR was carried out using TB Green Premix Ex Taq II (Takara, Japan). The primer sequences specific to target genes were listed in Table S6. The cycling protocol was as follows: 95 • C for 30 s, followed by 35 cycles of 95 • C for 5 s, 60 • C for 10 s, and 72 • C for 30 s. A melting curve was generated by cooling the products to 65 • C and then heating to 95 • C at a rate of 0.1 • C/s while simultaneously measuring fluorescence. Samples were run in triplicate for technical repeats and triplicate for biological repeats. Relative enrichment levels were determined by comparison with An. sinensis 18 s rRNA [66], using the 2−∆∆Ct method.

Analysis of RNA Modification by Quantitative Mass Spectrometry
The Poly (A) RNA was purified from 200 µg total RNA using one round Dynabeads Oligo (dT) (Thermo Fisher, Carlsbad, CA, USA) purification. The mRNA samples (1 µg) were denatured at 95 • C for 5 min, digested by S1 nuclease (1 U), alkaline phosphatase (0.1 U), and phosphodiesterase (0.01 U) at 37 • C for 2 h. After RNA was digested into nucleosides completely, the mixture was extracted with chloroform to remove the enzymes.  Cleaved RNA fragments were incubated with m6A-specific antibody (Synaptic Systems GmbH, Goettingen, Lower Saxony, Germany) for 2 h at 4 • C in IP buffer (50 mM Tris-HCl, 750 mM NaCl, and 0.5% Igepal CA-630). Both IP RNA and total Poly (A) RNA (input) were reverse-transcribed to create the cDNA by SuperScript™ II Reverse Transcriptase (Invitrogen, Carlsbad, CA, USA), which were next used to synthesize U-labeled secondstranded DNAs with E. coli DNA polymerase I (NEB, Ipswich, MA, USA), RNase H (NEB, Ipswich, MA, USA) and dUTP Solution (Thermo Fisher, Carlsbad, CA, USA). An A-base was then added to the blunt ends of each strand and ligate to the T-base overhang of adapters. The size selection for ligated fragments was performed with AMPureXP beads (NEB, Ipswich, MA, USA). After digesting the U-labeled second-stranded DNAs with UDG enzyme (NEB, Ipswich, MA, USA), the ligated products were amplified with PCR by the following conditions: initial denaturation at 95 • C for 3 min, 8 cycles of denaturation at 98 • C for 15 s, annealing at 60 • C for 15 s, extension at 72 • C for 30 s, final extension at 72 • C for 5 min. The cDNA libraries were sequenced at LC-Bio (Hangzhou, Zhejiang, China) using an illumina Novaseq™ 6000 instrument. Libraries from the samples were sequenced in 2 × 150 bp paired-end sequencing (PE150) mode, with average insertions length of 300 ± 50 bp. Raw data of RNA-seq and m6A-seq have been uploaded to GEO database (accession number GSE193379).
Read processing. Fastq files for male and female mosquitoes were trimmed using Fastp software [67] to remove the reads that contained adaptor contamination, low-quality bases, and undetermined bases with default parameter. RNA-seq data from control mRNA and m6A mRNA were mapped to Anopheles sinensis reference genome sequence (Version: v49) using HISAT2 package [68].
Peak calling and motif analysis. Mapped reads of IP and input cDNA libraries from male and female mosquitoes were used to identify peaks by R package exomePeak [69]. The m6A peaks were annotated using R package ChIPseeker [70] and visualized using IGV software [71]. The de novo and known motif predicting were performed using MEME [72] and HOMER [73] to locate the peaks of the motif.
Expression level of mRNA and analysis. The read counts per gene were then normalized to obtain FPKM values (total exon fragments /mapped reads (millions) × exon length (kB)) using StringTie [74]. Differentially expressed genes and differences in abundance of m6A peaks between male and female mosquitoes were calculated as fold change (FC), with absolute log2 FC ≥ 1 and p < 0.05 by R package edgeR [75]. Gene ontology (GO) and KEGG enrichment were carried out with online tools of LC-BIO company (https://www.omicstudio.cn/index, accessed on 22 October 2021).

Conclusions
An. sinensis is the vector for several parasites and viruses threatening human health, and novel tools targeting mosquito physiological developments and physical behaviors in the life cycle are urgently needed for vector control. In the present study, the profile of the transcriptome-wide m6A in male and female An. sinensis is characterized for the first time, and it is confirmed that m6A is a highly conserved modification in mosquitoes but with distinct characteristics compared to other species, including mammals and yeast. Moreover, unique distribution patterns of m6A in An. sinensis are associated with mosquito sexspecific pathways, especially in regulating spermatogenesis and the sperm tail formation of An. sinensis. This evidence paves the way for revealing the detailed regulation mechanism in mosquito physiological development and is a starting roadmap for uncovering m6A functions that may affect/control mosquito development in future, contributing to the control of this disease-transmitting vector.

Conflicts of Interest:
The authors declare no conflict of interest.