piRNA Profiling of Dengue Virus Type 2-Infected Asian Tiger Mosquito and Midgut Tissues

The Asian tiger mosquito, Aedes albopictus, is a competent vector for the majority of arboviruses. The mosquito innate immune response is a primary determinant for arthropod-borne virus transmission, and the midgut is the first barrier to pathogen transmission. Mosquito antiviral immunity is primarily mediated by the small interfering RNA pathway. However, the roles that the P-element induced wimpy testis (PIWI)-interacting RNA (piRNA) pathway play in antiviral immunity in Ae. albopictus and its midgut still need further exploration. This study aimed to explore the profiles of both viral-derived and host-originated piRNAs in the whole body and midgut infected with Dengue virus 2 (DENV-2) in Ae. albopictus, and to elucidate gene expression profile differences of the PIWI protein family between adult females and their midguts. A deep sequencing-based method was used to identify and analyze small non-coding RNAs, especially the piRNA profiles in DENV-2-infected Ae. albopictus and its midgut. The top-ranked, differentially-expressed piRNAs were further validated using Stem-loop qRT-PCR. Bioinformatics analyses and reverse-transcription PCR (RT-PCR) methods were used to detect PIWI protein family members, and their expression profiles. DENV-2 derived piRNAs (vpiRNA, 24–30 nts) were observed in both infected Ae. albopictus and its midgut; however, only vpiRNA in the whole-body library had a weak preference for adenine at position 10 (10A) in the sense molecules as a feature of secondary piRNA. These vpiRNAs were not equally distributed, instead they were derived from a few specific regions of the genome, especially several hot spots, and displayed an obvious positive strand bias. We refer to the differentially expressed host piRNAs after DENV infection as virus-induced host endogenous piRNAs (vepiRNAs). However, we found that vepiRNAs were abundant in mosquito whole-body tissue, but deficient in the midgut. A total of eleven PIWI family genes were identified in Ae. albopictus; however, only AalPiwi5–7 and AalAgo3(1–2) were readily detected in the midgut. The characteristics of piRNAs in DENV-2-infected Ae. albopictus adult females were similar to those previously described for flavivirus infections but were not observed in the midgut. The reduced levels of vepiRNAs and incomplete expression of PIWI pathway genes in midgut samples from DENV-2-infected Ae. albopictus suggests that viral regulation of host piRNAs may not be an important factor in the midgut.


Introduction
Small RNA (sRNA) regulatory pathways (SRRPs) control essential aspects of antiviral immune responses to arboviruses in mosquitoes. There are three major SRRPs: canonical small interfering RNA expression of host and/or specific viral mRNAs [31,32]. The origin of vpiRNA is not completely clear; however, in some viruses, this pathway could also be involved in the ping-pong mechanism, such as with the host cell native piRNA [6]. Furthermore, Li et al. were the first to characterize Ultrashort RNA (usRNA, also unusually-small RNA, 13-19 nt in length) in Kaposi's sarcoma-associated herpesvirus (KSHV)-infected cells [33]. These usRNAs were also characterized in dengue virus (DENV)-infected Aedes mosquitoes [20,25], as well as in bunyavirus La Crosse virus (LACV)-and Sindbis virus (SINV)-infected mosquito cells [19]. However, the role of usRNAs in the antiviral immune response still needs to be further elucidated.
Aedes albopictus (Asian tiger mosquito), considered one of world's worst invasive alien species, which also has an aggressive biting behavior, has great potential to transmit multiple virus infections. As a secondary vector of dengue viruses, Ae. albopictus-borne dengue fever (DF) and dengue hemorrhagic fever (DHF) are still major threats to human health in tropical urban areas of Asia [34,35]. Transmission of DENV through Aedes requires a virus to effectively invade and replicate in midgut epithelial cells, disseminate into the hemocoel from the midgut, and move to the salivary gland acini [36]. Thus, the midgut is the first barrier to pathogen transmission. Furthermore, SRRPs have also been found in the midguts of West Nile virus (WNV) orally-exposed Culex pipiens [37]. Therefore, further investigation is required to determine whether the piRNA pathways are broadly involved in the antiviral response in Ae. Albopictus; their role in the midgut also needs to be investigated.
In this study, we used a deep sequencing approach to characterize that Ae. albopictus, infected with dengue virus type-2 (DENV-2), produces both vpiRNAs and vepiRNAs. We then determined whether specific piRNAs are produced in mosquito whole body and midguts in response to infection. This will significantly enhance our understanding of the relationship between the piRNA pathway and the mosquito innate immune response to virus infections.

Mosquitoes
The experiments were performed using the Ae. albopictus Foshan (Guangdong, China) strain, which was kindly provided by the Center for Disease Control of Guangdong (CDC, Guangdong), China, and established in the laboratory in 1981. All mosquitoes were reared in a bioclimatic insectary at 25 ± 1 • C, with 70-80% relative humidity, under a 10:14 h daily light-dark (LD) cycle. Larvae were fed with a mixture of dry yeast powder and fish food (Yee ® , Shandong, China) at a ratio of 1:3, and adults were fed with a sugar solution (10% glucose). Adult female mosquitoes were allowed to take blood by feeding on an anesthetized mouse, which was then returned to its cage.

Virus Strains and Infection
The DENV-2 New Guinea C strain (NGC) was used in this study and was obtained from Sun Yat-sen University Zhongshan School of Medicine. DENV-2 stocks were obtained by inoculating a nearly confluent monolayer of C6/36 cells in a 25-cm 2 tissue culture flask, with the virus diluted to 1:5 in 1 mL RPMI-1640 medium (Gibco; Thermo Fisher Scientific, Inc., Waltham, CA, USA) containing 2% heat-inactivated fetal bovine serum (FBS; Gibco; Thermo Fisher Scientific, Inc.). After five days in culture, the cells and supernatant were subsequently harvested via gentle pipetting. Then, the cells were pelletized using centrifugation at 3000 rpm for 15 min, and the viral supernatant was harvested and adjusted to 20% FBS, aliquoted, and then stored at −80 • C. The viral titers were titrated using a viral plaque formation on BHK-21 cells. Results showed that the virus titer was 10 6.85 plaque-forming units per ml (PFU)/mL.
Ae. albopictus adult females (4-6 days old) were starved for 24 h and then moved to 2.5 L disposable paper feeding containers. Adult females were fed with defibrinated sheep's blood (Applied Biological Products Management, Aldinga Beach, Australia) mixed with DENV-2 viral supernatant (5:1). Prior to feeding, blood-virus mixtures were incubated at 37 • C for 1 h using glass membrane feeders covered by a porcine intestinal membrane [38]. After feeding, mosquitoes were briefly anesthetized using CO 2 ; non-engorged and partially-fed mosquitoes were discarded. The control mosquitoes were exposed to uninfected bloodmeal.

Mosquito Tissue Dissection and RNA Isolation
DENV-2-positive mosquitoes were identified by reverse-transcription PCR (RT-PCR) with DENV-2-specific primers at 9 days post-infection (dpi). Mosquito forelegs were removed, and RNA was extracted using TRIzol Reagent (Invitrogen; Thermo Fisher Scientific, Inc.). RT-PCR was employed to identify infected mosquitos with DENV-2-specific primers [39] (Table S1). The RT process was performed using the SuperScript ® III first-strand synthesis system (Thermo Fisher Scientific, Invitrogen, Carlsbad, CA, USA), and only DENV-2-positive mosquitoes were used for subsequent sRNA library construction.
Pools of three groups of infected/uninfected midguts (n = 100 per group) were collected at 1 dpi by dissecting anesthetized DENV-2-positive mosquitoes onto a drop of pre-chilled physiological saline, and the slide was then placed on a dissecting microscope. Midguts were placed into RNase-free 1.5 mL microcentrifuge tubes containing 1 mL of RNAlater (Life Technologies, Thermo Fisher Scientific, Walthman, MA, USA) and were then flash frozen in a dry ice/ethanol bath and stored at −80 • C until subsequent RNA extraction. Total RNA from each group was extracted with Trizol ® Reagent, according to the manufacturer's protocol. For the DENV-2-infected/uninfected adult female samples, pools of three groups (n = 50 per group) of whole mosquitoes were combined, frozen in liquid nitrogen, and then ground to a fine powder for RNA extraction using TRIzol ® reagent.

sRNA Library Construction and Sequences Analysis
The primary objective of this study was to profile the small RNA in DENV-2 infected mosquitos and midguts. Instead of constructing individual libraries for each replicate, 6 ng of total RNA per biological replicate in each sample (3 replicates) was pooled together for subsequent sRNA library construction. The small RNAs were purified using 15% denaturing polyacrylamide gel electrophoresis to enrich the desired size range (13-30 nt), and were sequentially ligated with 5 -(5GUUCAGAGUUCUACAGUCCGACGAUC) and 3 -(5PUCGUAUGCCGUCUUCUGCUUGUidT) end RNA oligonucleotide adaptors; then, they were subsequently reversely transcribed with the small RNA RT-Prime (5 -CAAGCAGAAGACGGCATACGA-3 ), and the adapter sequence specific primers were used to amplify the inserted DNA using PCR. After cluster generation of purified PCR products (forward primer: 5 -AATGATACGGCGACCACCGACAGGTTCAGAGTTCTACAGTCCGA-3 ; reverse primer: 5 -CAAGCAGAAGACGGCATACGA-3 ), library preparations were sequenced using an Illumina Genome Analyzer (Illumina, San Diego, CA, USA). All sequencing was carried out at the Beijing Genomics Institute (BGI) (Shenzhen, China). Raw sequence reads were submitted to the NCBI Short Read Archive (SRA) under the accession numbers SRA129084 and SRP049501.

Sequence Analysis
The adaptor sequences were removed, low-quality tags were filtered out, and the contamination of adaptor-adaptor ligation was cleaned up with Trimmomatic-0.30 using default parameters [40]. Unique reads of 24-30 nt were selected for further analysis and were mapped to the genome assembly of the Ae. albopictus Foshan strain (version AaloF1), as provided in VectorBase [41], using the Short Oligonucleotide Analysis Package 2 (SOAP2) [42]. Furthermore, reads matching noncoding tRNAs, snoRNAs, rRNAs, small nuclear RNAs (snRNAs), smiRNAs, and exons of protein-coding genes were removed; additionally, reads containing poly-A/T/C/G nucleotides (a minimum of eight homopolymer repeat nucleotides) were excluded. The remainder of the genome-matched reads were used as piRNA-like sRNA candidates for subsequent analyses.

Ping-Pong Signature
Ping-pong signals were defined as a typical pair of sense and antisense piRNAs, overlapping by 10 nt at their 5 ends [43]. To confirm the fraction of piRNAs involved in ping-pong pairs, we first counted all the uniquely-matched piRNAs; each piRNA was then matched to the database of piRNAs in the pool. Complimentary piRNAS matching along the 5 ends were considered a pair. We counted the length of paired bases (the last paired nucleotide position relative to the 5 ends of the piRNAs) of every partner. Finally, the paired piRNAs were classified according to the length of overlapped paired bases, as previously described, and a 10 bp overlap was defined as a typical ping-pong signature [44].

Stem-Loop qRT-PCR Analyses
Candidate piRNAs were reverse transcribed to cDNA per biological replicate using specific stem-loop RT primers (Table S2). Quantitation was performed on an ABI 7500 Fast Real-Time PCR system (Applied Biosystems Life Technologies, Foster City, CA, USA) using cDNA specific forward primer and a universal reverse primer, as listed in Table S2. U6 small nuclear RNA (U6 snRNA) was used for normalization in all samples. The specificity of amplification was confirmed by a melting curve analysis and also by running PCR products on agarose gels (3%). The quantitative RT-PCR analysis was performed three times using independent purified RNA samples with three replicates for each sample.

Phylogenetic Analyses
Putative amino acid sequences of full-length PIWI subfamily proteins were aligned using Clustal X software [46] (version 1.81) using the default parameters. Phylogenetic relationships among proteins were estimated using the neighbor-joining (NJ) method using the MEGA 4 package [47]. The reliability of the resulting topologies was tested using the bootstrap method [48] with 1000 pseudoreplicates.

Gene Expression Analyses of the PIWI Protein Family
To explore the expression profiles of the PIWI protein family gene in Ae. albopictus, time series transcriptome data from the midguts of DENV-2-infected adult females and carcasses of Ae. albopictus (three replicates per sample, n = 50-70 per replicate)(GenBank SRA database accession number SRP077936) were used to analyze the gene expression levels [49]. Clean reads were aligned to the Ae. albopictus Foshan strain genome reference with Bowtie2 (version 2.3.4.1) [50] and then the sequencing reads for each RefSeq gene were calculated using RSEM (version 1.273) [51] to estimate the expression levels. The trimmed mean of the M-values normalization method (TMM) was used to normalize the gene expression levels in different libraries [52]. Finally, the transcript expression levels were compared, based on transcripts per million transcripts (TPM). The TPM can be seen as counts and therefore, as Poisson distributed values. A test for differences in TPM was found in two experimental conditions based on Poisson distribution statistics [53].
The gene expression profiles were further confirmed using RT-PCR. Total RNA was extracted from Ae. albopictus adults from different tissues (whole body and midgut), and at different developmental stages (0-24 h embryos, 24-48 h embryos, 1st-2nd larvae, 3rd-4th larvae, pupae, 2 days post-emergence sugar meal adult males and females, 24 h and 48 h post blood feed adult females) using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA). Any residual DNA was removed using a TURBO DNA-free™ Kit (Ambion, Life Technologies, TX, USA). First-strand cDNA was synthesized from the total RNA using Oligo (dT) primers and a RevertAid First Strand cDNA Synthesis Kit (Thermo scientific). PIWI protein family genes were detected by RT-PCR using sequence-specific primers. Ae. albopictus ribosomal protein 7 gene (AalrpS7) (GenBank: JN132168) was used as a control. The sequences of the primers used in the RT-PCR are shown in Table S1. The amplified products of RT-PCR were separated on 1.0% or 2.0% agarose gels, depending on the size of the product, and purified using an E.Z.N.A Gel Extraction Kit (Omega Bio-tek, Doraville, GA, USA). These cDNA fragments were cloned into a plasmid pGEM-T Easy Vector (Promega, Madison, WI, USA) and confirmed via sequencing (Beijing Genomics Institute, BGI, Beijing, China).

sRNA Seq Analyses
In order to explore the DENV-2-induced expression profiles of piRNAs in Ae. albopictus and midguts, we sequenced four libraries from uninfected adult females, DENV-2-infected whole adult females, and midguts from uninfected and DENV-2-infected adult females. After trimming the adaptors, cleaning up contaminations, and excluding low-quality reads, clear reads were obtained from four different libraries (Table S3). Both uninfected female adults ( Figure 1A) and uninfected female midguts ( Figure 1B) showed a major peak at 22 nucleotides, accompanied by a sub-peak at 26-28 nucleotides; in DENV-2-infected female adults, the 22 nt class was the secondary peak, and the primary peak shifted to 27 nt, which may represent PIWI-interacting RNAs (piRNAs) ( Figure 1A). sequence-specific primers. Ae. albopictus ribosomal protein 7 gene (AalrpS7) (GenBank: JN132168) was used as a control. The sequences of the primers used in the RT-PCR are shown in Table S1. The amplified products of RT-PCR were separated on 1.0% or 2.0% agarose gels, depending on the size of the product, and purified using an E.Z.N.A Gel Extraction Kit (Omega Bio-tek, Doraville, GA, USA). These cDNA fragments were cloned into a plasmid pGEM-T Easy Vector (Promega, Madison, WI, USA) and confirmed via sequencing (Beijing Genomics Institute, BGI, Beijing, China).

sRNA Seq Analyses
In order to explore the DENV-2-induced expression profiles of piRNAs in Ae. albopictus and midguts, we sequenced four libraries from uninfected adult females, DENV-2-infected whole adult females, and midguts from uninfected and DENV-2-infected adult females. After trimming the adaptors, cleaning up contaminations, and excluding low-quality reads, clear reads were obtained from four different libraries (Table S3). Both uninfected female adults ( Figure 1A) and uninfected female midguts ( Figure 1B) showed a major peak at 22 nucleotides, accompanied by a sub-peak at 26-28 nucleotides; in DENV-2-infected female adults, the 22 nt class was the secondary peak, and the primary peak shifted to 27 nt, which may represent PIWI-interacting RNAs (piRNAs) ( Figure 1A).

Molecular Characterization of DENV-2-Derived sRNAs
DENV-2-derived small RNAs (DENV-2 vsRNAs) were identified by aligning the small RNA library of DENV-2-infected adult females and DENV-2-infected female midguts to the DENV-2 genome through pattern matching using SOAP2 (Table S4). The sRNA segments (13-30 nt), of which sequences showed complete matches to the DENV-2 genome at each base position, were considered to be vsRNAs. vsRNAs are distinguishable from host-derived sRNAs in that they may show imperfect complementarity to sites in the minus-or plus-strand of the viral genome. The analysis results show that the DENV-2-infected adult female libraries had 11,716 reads, representing 6345 unique perfectly-matched DENV-2 vsRNAs, accounting for only 5.60 × 10 −2 percent of the total sRNA reads. The DENV-2-infected female midgut libraries displayed the number of vsRNAs (44,494 reads and 14,885 unique reads), and accounted for 37.41 × 10 −2 percent of the total abundance of sRNA reads, possibly due to the rapidly-increasing levels of viral proliferation in the midgut, or to the

Molecular Characterization of DENV-2-Derived sRNAs
DENV-2-derived small RNAs (DENV-2 vsRNAs) were identified by aligning the small RNA library of DENV-2-infected adult females and DENV-2-infected female midguts to the DENV-2 genome through pattern matching using SOAP2 (Table S4). The sRNA segments (13-30 nt), of which sequences showed complete matches to the DENV-2 genome at each base position, were considered to be vsRNAs. vsRNAs are distinguishable from host-derived sRNAs in that they may show imperfect complementarity to sites in the minus-or plus-strand of the viral genome. The analysis results show that the DENV-2-infected adult female libraries had 11,716 reads, representing 6345 unique perfectly-matched DENV-2 vsRNAs, accounting for only 5.60 × 10 −2 percent of the total sRNA reads. The DENV-2-infected female midgut libraries displayed the number of vsRNAs (44,494 reads and 14,885 unique reads), and accounted for 37.41 × 10 −2 percent of the total abundance of sRNA reads, possibly due to the rapidly-increasing levels of viral proliferation in the midgut, or to the insulation of the dsRNA trigger in cytoplasmic membrane-enclosed vesicles [54][55][56]. A small number of vsRNAs perfectly mapped both the viral and host genomes (Table S4), which matched the Ae. albopictus genome in a random manner, due to the complete nature of sRNA homology for both viral and host genomes; thus, it is difficult to discover the origin of this micro group of sRNAs.
The predominant vsRNA sizes, both in DENV-2-infected adult females ( Figure 2A) and DENV-2-infected female midgets ( Figure 2B), were grouped into 21 nt (vsiRNA), and were considered to be the size characterization for DCL products. However, in contrast to the 21 nt vsRNA percentage (89.3% of vsRNA) in the DENV-2-infected adult female libraries, 21 nt vsRNA only occupied 14.61% of vsRNAs reads in the DENV-2-infected female midgut libraries. Though vsRNAs in DENV-2-infected adult females also displayed a sub-peak at 26 nt, these 24-30 nt length vsRNAs were less than 3% of the abundance of DENV-2-infected adult female vsRNAs. DENV-2-infected female midgut vsRNAs displayed an approximate Gaussian distribution, and the 24-30 nt length sRNAs corresponded to 31.25% of the reads. insulation of the dsRNA trigger in cytoplasmic membrane-enclosed vesicles [54][55][56]. A small number of vsRNAs perfectly mapped both the viral and host genomes (Table S4), which matched the Ae. albopictus genome in a random manner, due to the complete nature of sRNA homology for both viral and host genomes; thus, it is difficult to discover the origin of this micro group of sRNAs. The predominant vsRNA sizes, both in DENV-2-infected adult females ( Figure 2A) and DENV-2-infected female midgets ( Figure 2B), were grouped into 21 nt (vsiRNA), and were considered to be the size characterization for DCL products. However, in contrast to the 21 nt vsRNA percentage (89.3% of vsRNA) in the DENV-2-infected adult female libraries, 21 nt vsRNA only occupied 14.61% of vsRNAs reads in the DENV-2-infected female midgut libraries. Though vsRNAs in DENV-2infected adult females also displayed a sub-peak at 26 nt, these 24-30 nt length vsRNAs were less than 3% of the abundance of DENV-2-infected adult female vsRNAs. DENV-2-infected female midgut vsRNAs displayed an approximate Gaussian distribution, and the 24-30 nt length sRNAs corresponded to 31.25% of the reads. The vsRNAs in the DENV-2-infected adult females ( Figure 3A) and the DENV-2-infected female midgut libraries ( Figure 3B) showed obviously positive (genome) strand biases of 68.14% and 93.56% positive-sense vsRNAs, respectively. In the DENV-2-infected female midgut libraries, the ratio of plus-to minus-strand-derived vsRNAs was relatively stable across the 13-30 nt length. However, in DENV-2-infected adult females, the ratio of positive-sense-derived vsRNAs peaked at 17 nt in length and increased from 22 nt, following the increased length of vsRNAs. In particular, for vsRNAs from 25 to 29 nt in length in the DENV-2-infected adult females libraries, nearly 96.27% of the vsRNAs found were derived from the positive-sense ( Figure 3A). The vsRNAs in the DENV-2-infected adult females ( Figure 3A) and the DENV-2-infected female midgut libraries ( Figure 3B) showed obviously positive (genome) strand biases of 68.14% and 93.56% positive-sense vsRNAs, respectively. In the DENV-2-infected female midgut libraries, the ratio of plus-to minus-strand-derived vsRNAs was relatively stable across the 13-30 nt length. However, in DENV-2-infected adult females, the ratio of positive-sense-derived vsRNAs peaked at 17 nt in length and increased from 22 nt, following the increased length of vsRNAs. In particular, for vsRNAs from 25 to 29 nt in length in the DENV-2-infected adult females libraries, nearly 96.27% of the vsRNAs found were derived from the positive-sense ( Figure 3A).

Genome Distribution of DENV-2-Derived sRNAs
To explore the origins of vsRNAs, we focused on the vsRNAs' genome locations by matching vsRNAs with the DENV-2 genome. The results indicated that vsRNAs in both DENV-2-infected adult females and DENV-2-infected female midgut libraries were not evenly distributed, but rather, were clustered in a few defined regions of the DENV-2 genome, especially in several hot spots. In DENV-2-infected adult females, several vsRNA hot spots originated from either negative-sense or positivesense vsRNA, whereas, in the DENV-2-infected female midgut libraries, hot spots were mainly produced from the positive-strand virus genome, and only a few lower abundance vsRNAs accumulated in regions in the negative strand. It is notable that there was a high abundance of reads from one concentrated at a point around 8200 nt (Figure 4), located at nonstructural protein 5 (NS5) of the DENV-2 genome, in both libraries. To further detect whether these hot spot-derived vsRNAs had some relation to secondary structures of DENV-2 genomic RNA structures and the complementary strand, and whether these RNA secondary structures could act as dsRNA substrates for Dcr2, we analyzed potential secondary structures in the DENV-2 genomic RNA and its reverse complementary strand using RNAfold [57]; as a result, most of these predicted local double-stranded stem regions lacked precise base pairing over >21 bp. Hence, it was inappropriate to correlate these secondary structures with vsRNA hot spots.

Genome Distribution of DENV-2-Derived sRNAs
To explore the origins of vsRNAs, we focused on the vsRNAs' genome locations by matching vsRNAs with the DENV-2 genome. The results indicated that vsRNAs in both DENV-2-infected adult females and DENV-2-infected female midgut libraries were not evenly distributed, but rather, were clustered in a few defined regions of the DENV-2 genome, especially in several hot spots. In DENV-2-infected adult females, several vsRNA hot spots originated from either negative-sense or positive-sense vsRNA, whereas, in the DENV-2-infected female midgut libraries, hot spots were mainly produced from the positive-strand virus genome, and only a few lower abundance vsRNAs accumulated in regions in the negative strand. It is notable that there was a high abundance of reads from one concentrated at a point around 8200 nt (Figure 4), located at nonstructural protein 5 (NS5) of the DENV-2 genome, in both libraries. To further detect whether these hot spot-derived vsRNAs had some relation to secondary structures of DENV-2 genomic RNA structures and the complementary strand, and whether these RNA secondary structures could act as dsRNA substrates for Dcr2, we analyzed potential secondary structures in the DENV-2 genomic RNA and its reverse complementary strand using RNAfold [57]; as a result, most of these predicted local double-stranded stem regions lacked precise base pairing over >21 bp. Hence, it was inappropriate to correlate these secondary structures with vsRNA hot spots.

Genome Distribution of DENV-2-Derived sRNAs
To explore the origins of vsRNAs, we focused on the vsRNAs' genome locations by matching vsRNAs with the DENV-2 genome. The results indicated that vsRNAs in both DENV-2-infected adult females and DENV-2-infected female midgut libraries were not evenly distributed, but rather, were clustered in a few defined regions of the DENV-2 genome, especially in several hot spots. In DENV-2-infected adult females, several vsRNA hot spots originated from either negative-sense or positivesense vsRNA, whereas, in the DENV-2-infected female midgut libraries, hot spots were mainly produced from the positive-strand virus genome, and only a few lower abundance vsRNAs accumulated in regions in the negative strand. It is notable that there was a high abundance of reads from one concentrated at a point around 8200 nt (Figure 4), located at nonstructural protein 5 (NS5) of the DENV-2 genome, in both libraries. To further detect whether these hot spot-derived vsRNAs had some relation to secondary structures of DENV-2 genomic RNA structures and the complementary strand, and whether these RNA secondary structures could act as dsRNA substrates for Dcr2, we analyzed potential secondary structures in the DENV-2 genomic RNA and its reverse complementary strand using RNAfold [57]; as a result, most of these predicted local double-stranded stem regions lacked precise base pairing over >21 bp. Hence, it was inappropriate to correlate these secondary structures with vsRNA hot spots.

Characterization of DENV-2-Derived vpiRNA
We next grouped these vsRNAs into three classes, based on their lengths. The first group was usRNAs (13-19 nt); the second group had a higher abundance of 20-23 nt vsiRNAs in both libraries, which depended on dicer for their biogenesis; and the third group was approximately 24-30 nt in length, as candidates of vpiRNAs ( Figure 3B). We focused on the 24-30 nt vpiRNAs. Similar to vsRNAs, vpiRNAs also displayed uneven genome distributions and were enriched at several hot spots. In DENV-2-infected adult females, vpiRNAs peaked at the 3 end of NS5 of the viral genomic RNA (Figure 5), whereas the peak was in the middle region of NS5 of the viral genomic RNA in DENV-2-infected female midguts. These predicted vpiRNAs exhibited a stronger strand bias, both in DENV-2-infected females and DENV-2-infected female midgut libraries, with more than 94.99% and 93.43% of the reads mapping to the viral plus-strand genomic RNA, respectively. The vpiRNAs detected in both libraries are listed in Table S5. We also analyzed the genomic distribution of usRNAs and vsiRNAs. Compared to vpiRNAs, vsiRNAs displayed a similar genomic distribution, and the major hot spots were located at the NS5 region in DENV-2-infected adult females, as well as the NS3 and NS5 regions in DENV-2-infected female midguts ( Figure S1). However, usRNAs showed more irregular distribution patterns throughout the genome, especially in DENV-2-infected female midguts ( Figure S1).

Characterization of DENV-2-Derived vpiRNA
We next grouped these vsRNAs into three classes, based on their lengths. The first group was usRNAs (13-19 nt); the second group had a higher abundance of 20-23 nt vsiRNAs in both libraries, which depended on dicer for their biogenesis; and the third group was approximately 24-30 nt in length, as candidates of vpiRNAs ( Figure 3B). We focused on the 24-30 nt vpiRNAs. Similar to vsRNAs, vpiRNAs also displayed uneven genome distributions and were enriched at several hot spots. In DENV-2-infected adult females, vpiRNAs peaked at the 3′ end of NS5 of the viral genomic RNA (Figure 5), whereas the peak was in the middle region of NS5 of the viral genomic RNA in DENV-2-infected female midguts. These predicted vpiRNAs exhibited a stronger strand bias, both in DENV-2-infected females and DENV-2-infected female midgut libraries, with more than 94.99% and 93.43% of the reads mapping to the viral plus-strand genomic RNA, respectively. The vpiRNAs detected in both libraries are listed in Table S5. We also analyzed the genomic distribution of usRNAs and vsiRNAs. Compared to vpiRNAs, vsiRNAs displayed a similar genomic distribution, and the major hot spots were located at the NS5 region in DENV-2-infected adult females, as well as the NS3 and NS5 regions in DENV-2-infected female midguts ( Figure S1). However, usRNAs showed more irregular distribution patterns throughout the genome, especially in DENV-2-infected female midguts ( Figure S1). Except for the vpiRNAs in the DENV-2-infected adult females libraries having a weak preference for adenine at position 10 (10A) (43.31%) as the feature of secondary sense piRNA ( Figure 6A), these 10A sense piRNAs were thought to have originated from the ping-pong amplification cycle in an Aub slicer-dependent manner. However, a uridine at the first position (5U), as a key characteristic of primary antisense piRNAs was not observed in DENV-2-infected adult females ( Figure 6A) and DENV-2-infected female midguts ( Figure 6B). Furthermore, a significant 10 nt overlap (ping pong signal) between sense and antisense vpiRNA were also absent in DENV-2-infected adult females ( Figure 6C) and DENV-2-infected female midguts ( Figure 6D). These results indicated that vpiRNA biogenesis are involved in non-canonical pathways. Except for the vpiRNAs in the DENV-2-infected adult females libraries having a weak preference for adenine at position 10 (10A) (43.31%) as the feature of secondary sense piRNA ( Figure 6A), these 10A sense piRNAs were thought to have originated from the ping-pong amplification cycle in an Aub slicer-dependent manner. However, a uridine at the first position (5U), as a key characteristic of primary antisense piRNAs was not observed in DENV-2-infected adult females ( Figure 6A) and DENV-2-infected female midguts ( Figure 6B). Furthermore, a significant 10 nt overlap (ping pong signal) between sense and antisense vpiRNA were also absent in DENV-2-infected adult females ( Figure 6C) and DENV-2-infected female midguts ( Figure 6D). These results indicated that vpiRNA biogenesis are involved in non-canonical pathways. Indicated above each axis is the number of possible overlapping pairs of small RNAs with the specified overlap size. Indicated below each axis is the relative frequency of the 5′ base identity for overlapping sequences. The color code for the bases is indicated in the centre box.

Characterization of vepiRNAs
To explore whether host piRNAs are generated in the process of host-virus interaction, we compared the uniquely-matched piRNA profiles of DENV-2-infected samples to uninfected controls, based on normalized expression levels (TPM).
A total of 61.44% of gene-derived piRNAs were matched to exon regions, and 38.56% were mapped to introns, which strongly suggests that gene-derived piRNAs are generated from precursor mRNA (pre-mRNA). The sense/antisense ratios of exon-and intron-derived piRNAs were 0.10 and 0.52, respectively. Exon-and intron-derived piRNAs also showed a strong 1U bias (82.42% and 94.37%, respectively), but no obvious elevations for 10A (28.57% and 21.12%, respectively).
Moreover, we introduced a threshold of 70 read counts, and these relatively highly-expressed piRNAs were analyzed further. Based on Poisson distribution statistics, a total of 1851 piRNAs were regarded as differential expressed vepiRNA candidates (p < 0.01, FDR < 0.01, and log2 fold change >1 or <−1) between DENV-2-infected adult females and uninfected female adults. A total of 22 piRNAs were regarded as differential expressed vepiRNA candidates between the DENV-2-infected female midguts and uninfected female midguts (Table S6). From the Poisson analysis of single sRNA library, selected vepiRNA were used for further qRT-PCR analysis using three biological replicates.
3.6. Validation of piRNA Expression Using Stem-Loop qRT-PCR Analysis

Characterization of vepiRNAs
To explore whether host piRNAs are generated in the process of host-virus interaction, we compared the uniquely-matched piRNA profiles of DENV-2-infected samples to uninfected controls, based on normalized expression levels (TPM).
A total of 61.44% of gene-derived piRNAs were matched to exon regions, and 38.56% were mapped to introns, which strongly suggests that gene-derived piRNAs are generated from precursor mRNA (pre-mRNA). The sense/antisense ratios of exon-and intron-derived piRNAs were 0.10 and 0.52, respectively. Exon-and intron-derived piRNAs also showed a strong 1U bias (82.42% and 94.37%, respectively), but no obvious elevations for 10A (28.57% and 21.12%, respectively).
Moreover, we introduced a threshold of 70 read counts, and these relatively highly-expressed piRNAs were analyzed further. Based on Poisson distribution statistics, a total of 1851 piRNAs were regarded as differential expressed vepiRNA candidates (p < 0.01, FDR < 0.01, and log2 fold change >1 or <−1) between DENV-2-infected adult females and uninfected female adults. A total of 22 piRNAs were regarded as differential expressed vepiRNA candidates between the DENV-2-infected female midguts and uninfected female midguts (Table S6). From the Poisson analysis of single sRNA library, selected vepiRNA were used for further qRT-PCR analysis using three biological replicates.
The topology obtained for the PIWI protein family protein phylogeny is shown in Figure 9. Phylogenetic analyses showed high confidence levels in the defined groups. As expected, based on phylogenetic relationships, three well-supported clades were obviously displayed, as described previously [45,58]. Ago3s of three genera of mosquitoes formed a single clade with D. melanogaster. Ae. albopictus PIWI1-6 were grouped into a clade with An. gambiae Ago4, whereas the PIWI7-9 were grouped into a clade with An. gambiae Ago5 (Figure 9).
Results based on RNA-seq ( Figure 10A) and RT-PCR ( Figure 10B) analyses indicated that transcripts of PIWI1-7 and Ago3(1-2) are readily detected in Ae. albopitus adult females, whereas Piwi8/9 is considerably lower. Different from the carcass, in the midgut, only PIWI5-7 and Ago3(1-2) were visibly detected ( Figure 9). Interestingly, PIWI8 and 9 transcripts were not detected in both the midgut and mosquito whole-body using RT-PCR and RNA-Seq analyses, so we continued to explore the temporal expression patterns of PIWI genes further using the RT-PCR method. As a result, Piwi8 and 9 were highly expressed in the embryo stage, whereas Piwi1-4 were shown to have a higher expression level in embryos, adult males, and in post blood meal (PBM) females ( Figure S2). Moreover, DENV-2 infection showed no obvious effects on the transcription levels of PIWI5-7 and Ago3. The PCR products were run on 3% agarose gel in 1X TBE stained with ethidium bromide. Note: vpiRNA2191731, vpiRNA2191732, vpiRNA0705792, vpiRNA0705793, vpiRNA3130226 and vpiRNA0705787 differed by only a few nucleotide base lengths at the 5 and 3 ends. A similar pattern was also found in vpiRNA0676933 and vpiRNA1747127. **: p < 0.01; ***: p < 0.001.
The topology obtained for the PIWI protein family protein phylogeny is shown in Figure 9. Phylogenetic analyses showed high confidence levels in the defined groups. As expected, based on phylogenetic relationships, three well-supported clades were obviously displayed, as described previously [45,58]. Ago3s of three genera of mosquitoes formed a single clade with D. melanogaster. Ae. albopictus PIWI1-6 were grouped into a clade with An. gambiae Ago4, whereas the PIWI7-9 were grouped into a clade with An. gambiae Ago5 (Figure 9).

Discussion
Our current understanding of the primary role of the piRNA pathway fulfils a vital function in the preservation of genomic integrity in animal germ cells against invasive transposable elements (TEs) [59]. Furthermore, some non-repetitive and non-transposon-related piRNAs may perform crucial functions in post-transcriptional regulation. For example, piRNA is considered to regulate the transcription of non-repetitive, protein-coding genes in D. melanogaster, and is involved in the sex determination hierarchy in Bombyx mori.
Ae. albopictus is ranked at number four on the list of 100 worst invasive species on Earth [57]. Ae. albopictus exhibits excellent adaptability to a broad spectrum of environmental conditions, is phenotypically polymorphic, and displays variations in its vectorial capacity for mosquito-borne viruses (MBV) [34][35][36]; it has a very large genome (1967 megabase pairs, Mbp) and is 1.43 times larger than the Ae. aegypti genome (1376 Mbp) [37]. Only virus-derived vsRNA profiles of infected Ae. albopictus C6/36 cells have been previously reported [21][22][23], but both adult female whole bodies and their midguts, which act as the first barrier for pathogen transmission, still lack supporting data.
The extrinsic incubation period (EIP) corresponds to the time necessary for a virus to reach mosquito saliva after an infectious blood meal [63]. For DENV, viral particles started to be detected in the salivary glands of Aedes mosquitos at approximately 9 dpi [68]. On the other hand, DENV had not completely passed through the midgut barrier at 1 dpi. Furthermore, vpiRNAs were visible at 1 dpi in SINV-infected Aag2 cells, when infection was fully established. Thus, we analyzed the sRNA transcriptome at these time points.
Three pathways were thought to be involved in the antivirus immune responses in mosquito midguts. First, the cytokine-activated Janus kinase (JAK)-signal transducer and activator of the transcription (STAT) signalling pathway in mosquitoes-was primarily observed in the Ae. aegypti midgut [69]. Transient impairment of the JAK/STAT pathway via RNAi-mediated gene silencing led to increased DENV titres in Ae. aegypti midguts. Second, the exogenous siRNA pathway showed no evidence of involvement in antiviral responses in the midgut of ONNV-infected An. Gambiae [70], and downregulation of SRRP key factors, Ago2, or Dicer2, promoted viral duplication in the midguts of flavivirus-or alphavirus-infected Ae. Aegypti [71][72][73]. Third, certain miRNAs were also confirmed to definitely influence DENV replication in Aedes mosquito midguts [38,39]. The abovementioned results indicate that the exogenous (exo-) siRNA pathway, JAK-STAT pathway, and miRNA pathway play important roles in the regulation of arbovirus replication in the midgut.
In our work, three types of vsRNA, usRNAs, vsiRNAs, and vpiRNA were detected in both DENV-2-infected female Ae. albopictus and midgut libraries. As expected, the exogenous vsiRNAs were more abundant in vsRNA, suggesting that vsiRNA-mediated RNA still played an important role in anti-DENV-2 defense function, both for the whole body and midgut of Ae. albopictus. In particular, vpiRNA was second in abundance, next to vsiRNA, in virus-infected midguts. However, the role of these DENV-2-specific vpiRNAs in the anti-DENV-2 response still need further investigation.
In most +ssRNA mosquito-borne viruses (alphaviruses and flaviviruses), the plus viral genomic RNA strand that corresponds to the viral transcript is over-represented compared to the minus strand [18][19][20][21][22][23]. Similar results were obtained in the present study; both vsRNA and vpiRNA showed obvious strand bias, and were mainly produced from the positive strand of the DENV-2 genome, which was in accordance with the relative amounts of plus and minus strands of DENV-2. vpiRNA production was mostly confined to specific hot spot regions in the DENV-2 genome, especially in the NS5 gene region, which showed similar results to piRNA hot spots in DENV-2 and cell-fusing agent virus (CFAV)-infected Ae. aeagyti and Aag2 cells [24]. In alphaviruses, such as SINV, Chikungunya virus (CHIKV), and Semliki Forest virus (SFV), vpiRNAs also displayed extreme strand bias towards sequences from the viral plus strands, but these alphavirus-derived vpiRNAs were predominantly produced from a subgenomic RNA that was transcribed from an internal promoter sequence [20]. The origins of these piRNA spikes remain unclear, but this difference strongly indicates that the mechanisms underlying alphavirus and flavivirus vpiRNA biogenesis are essentially variant.
The propensity for a "5U", as a distinct characteristic of primary piRNAs, was not observed in vpiRNA for both the DENV-2-infected adult females and the DENV-2-infected female midgut libraries. However, in the DENV-2-infected adult female libraries, vpiRNA showed a weak preference for 10A as a feature of secondary piRNA. Moreover, these vpiRNAs did not exhibit a typical ping-pong signature. These results were similar to those previously described in DENV and cell fusing agent virus (CFAV)-infected mosquito cell lines and Ae. Aegyti [6,[24][25][26]. However, the major characteristics of 5U and 10A bias as well as ping-pong signature, were observed in SINV-, CHIKV-, and SFV-infected mosquitoes and cell lines [18][19][20][21][22][23]. Thus, it is noted that these flaviviruses vpiRNAs may use a unique biogenesis pathway that is distinct from that of alphaviruses. Although flavivirus vpiRNAs exhibit the secondary piRNA features, whether a different, amplification-independent mechanism is responsible for flavivirus vpiRNA production remains to be determined.
As previously described [22], the vpiRNA-associated proteins, Piwi5 and Ago3, were predominantly expressed in the cytoplasm of Aag2 cells. Thus, as +ssRNA virus RNA replication occurs in the cytoplasm, and generally does not enter the nucleus, the key biogenesis factors of vpiRNA and viral RNAs are co-expressed in the cytoplasm. In fact, the vast majority of vpiRNAs were also present in the cytoplasmic fraction, indicating that vpiRNA biogenesis occurs in the cytoplasm of mosquito cells. Thus, it is supposed that vpiRNA production is a purely cytoplasmic event, and other unknown mechanisms must be involved in vpiRNA biogenesis. It has also been previously shown that vpiRNA-like small RNAs are capable of modulating the pathogenesis of alphavirus infections in dicer-2 null mutant mosquito cell lines, and are defective in viral siRNA production [20]. In addition to vpiRNA, host endogenous piRNAs were thought to also be involved in the anti-viral response. For example, the knockdown of Piwi4 (a member of the PIWI gene family) does not cause a reduction in vpiRNA levels [21], but, rather, enhances the replication of SFV in Aag2 cells [21].
We identified 1479 vepiRNAs displaying an up-regulated expression in the DENV-2-infected adult female library; however, no up-regulated expression piRNAs were found in DENV-2-infected female midguts. Thus, we analyzed the PIWI protein family members of adult females and their expression profiles. In the midgut, only Piwi5-7 and Ago3(1-2) were readily detected in every midgut sample. As previously reports note, the biogenesis of TE-derived host piRNA, either direct or indirect, requires nearly all PIWI protein activity, which suggests the lack of necessary PIWI proteins as a reason for vepiRNA deficiency in the midgut [22]. Furthermore, the 10A bias, as a feature of secondary piRNA, does not appear in midgut vpiRNAs, which indicates that the lack of certain PIWI proteins could lead to the loss of the vpiRNA biogenesis pathway in the midgut, compared to other tissues. Whether these 24-31 nt DENV-2-derived small RNA originated from the PIWI process or not still needs to be further explored. As previously described, Piwi1 and Piwi3 were highly expressed specifically in the ovaries of mosquitoes, a tissue that was generally not infected by SINV [74,75]. The ovaries were also less susceptible to DENV-2 infection than other tissues [76], Piwi1-3 may be associated with viral dissemination in mosquitoes. The role of Piwi1-3 in antivirus immunity requires further investigation.
Except for intergenic regions, major vepiRNAs are derived from Class I TEs (retrotransposons); moreover, the majority of them are derived from the antisense strand. In addition to retroviruses, RNA-only viruses have also given rise to transposons in many mosquito genomes, although the specific mechanisms involved in this horizontal RNA-to-DNA information flow are less clear. The exogenous non-self nucleic acids from infecting viruses can be targeted by piRNA-like RNAs, but this requires genetic information flow in an unexpected retrotransposon-dependent manner [77]. Thus, we deduced that vepiRNAs derived from TEs are antisense relative to the proposed ancient viral RNA, originally a response to interactions between infectious viruses and transposons.