Article

# Small RNA Response to Infection of the Insect-Specific Lammi Virus and Hanko Virus in an *Aedes albopictus* Cell Line

Pontus Öhlund <fnref n="1" />*, Juliette Hayer <fnref n="2" /> ID, Jenny C. Hesson <fnref n="3" /> ID and Anne-Lie Blomström <fnref n="1" />

<figure><img src="image_1.png" /></figure>

check for
updates

**Citation:** Öhlund, P.; Hayer, J.; Hesson, J.C.; Blomström, A.-L. Small RNA Response to Infection of the Insect-Specific Lammi Virus and Hanko Virus in an *Aedes albopictus* Cell Line. *Viruses* **2021**, *13*, 2181. https://doi.org/10.3390/v13112181

Academic Editor: A. Lorena Passarelli

Received: 22 September 2021

Accepted: 27 October 2021

Published: 29 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_2.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** RNA interference (RNAi)-mediated antiviral immunity is believed to be the primary defense against viral infection in mosquitoes. The production of virus-specific small RNA has been demonstrated in mosquitoes and mosquito-derived cell lines for viruses in all of the major arbovirus families. However, many if not all mosquitoes are infected with a group of viruses known as insect-specific viruses (ISVs), and little is known about the mosquito immune response to this group of viruses. Therefore, in this study, we sequenced small RNA from an *Aedes albopictus*-derived cell line infected with either Lammi virus (LamV) or Hanko virus (HakV). These viruses belong to two distinct phylogenetic groups of insect-specific flaviviruses (ISFVs). The results revealed that both viruses elicited a strong virus-derived small interfering RNA (vsiRNA) response that increased over time and that targeted the whole viral genome, with a few predominant hotspots observed. Furthermore, only the LamV-infected cells produced virus-derived Piwi-like RNAs (vpiRNAs); however, they were mainly derived from the antisense genome and did not show the typical ping-pong signatures. HakV, which is more distantly related to the dual-host flaviviruses than LamV, may lack certain unknown sequence elements or structures required for vpiRNA production. Our findings increase the understanding of mosquito innate immunity and ISFVs' effects on their host.

**Keywords:** insect-specific flaviviruses; Hanko virus; Lammi virus; small interfering RNA; Piwi-interacting RNA

## 1. Introduction

Most of our current knowledge on virus interactions with the mosquito vector has come from studies on human pathogenic arboviruses such as West Nile virus (WNV), Dengue virus (DENV), and Chikungunya virus (CHIKV). However, mosquitoes are often naturally and persistently infected with a group of viruses that are known as insect-specific viruses (ISVs). This group of viruses are unable to infect vertebrates and are maintained in the mosquito population through vertical transmission from mother to offspring [1-4]. Moreover, phylogenetic analyses have shown that many ISVs belong to viral families associated with pathogenic arboviruses such as *Flaviviridae*, *Togaviridae*, and *Peribunyaviridae* [5,6]. Insect-specific flaviviruses (ISFVs) can further be divided into two distinct phylogenetic subgroups. Those of the first subgroup are phylogenetically distinct from all other known flaviviruses and are referred to as classical ISFVs. The ISFVs of the other subgroup are phylogenetically affiliated with the medically important dual-host flaviviruses but have the insect-restriction phenotype and are referred to as dual-host affiliated ISFVs [7]. Two members of the insect-specific flavivirus group are Lammi virus (LamV) and Hanko virus (HakV), which were both isolated from mosquitoes in Finland. LamV was isolated from *Aedes cinereus* mosquitos and belongs to the dual-host affiliated

<footnote n="1">¹ Department of Biomedical Sciences and Veterinary Public Health, Swedish University of Agricultural Sciences, P.O. Box 7028, 750 07 Uppsala, Sweden; anne-lie.blomstrom@slu.se</footnote>

<footnote n="2">² Department of Animal Breeding and Genetics, Swedish University of Agricultural Sciences, SLU-Global Bioinformatics Centre, P.O. Box 7023, 750 07 Uppsala, Sweden; juliette.hayer@slu.se</footnote>

<footnote n="3">³ Department of Medical Biochemistry and Microbiology, Zoonosis Science Center, Uppsala University, P.O. Box 582, 751 23 Uppsala, Sweden; jenny.hesson@imbim.uu.se</footnote>

<footnote>* Correspondence: pontus.ohlund@slu.se; Tel.: +46-18-672-409</footnote>
ISFV category [8], and HakV was isolated from *Aedes caspius* mosquitoes and belongs to the classical ISFV category [9]. Little is known about how these ISFVs interact with the mosquito and whether they trigger an antiviral immune response in the host.

The main antiviral immune response in mosquitoes is believed to be the RNA interference (RNAi) pathway, which includes three major classes of regulatory small RNAs (sRNAs): Small interfering RNA (siRNA), microRNA (miRNA), and P-element-induced wimpy testis (Piwi)-interacting RNA (piRNA) [10]. The siRNA pathway is possibly the most important one with regards to antiviral defense in the mosquito. Upon a viral infection, double-stranded RNA (dsRNA) will form as replicative intermediates in the cytoplasm. These foreign (exogenous) dsRNAs are recognized and processed by the Dcr2-R2D2 complex into siRNAs of 21 nucleotides (nt) in length. Furthermore, the siRNA is incorporated into the RNA-induced silencing complex (RISC), where one of the strands is degraded and the remaining guide strand is used as a template to recognize complementary target RNAs, which are subsequently cleaved and degraded by Ago2 [11–13]. Functional Dcr2, R2D2, and Ago2 are crucial for limiting virus replication and dissemination in the mosquito, as demonstrated by many studies [14–18].

The miRNA pathway is similar to the siRNA pathway in the sense that dsRNAs are processed into smaller dsRNAs that are loaded into an RISC complex and facilitate target-specific cleavage. In the miRNA pathway, Dcr1 and Ago1 are analogous to Dcr2 and Ago2, respectively, in the siRNA pathway. However, the primary function of the miRNA pathway is post-transcriptional gene regulation, and the miRNAs (20–25 nt) come from endogenous hairpin transcripts [19,20]. The miRNA pathway has been shown to play a critical role in modulating host genes to control viral infection, e.g., the Aae-miR-2940-5p miRNA was found to be selectively downregulated in C6/36 cells as a response to WNV infection. This resulted in a lower expression of the metalloprotease m41 FtsH gene, which has been shown to be important for WNV replication [21]. Other studies have observed miRNA-driven gene modulation during the active viral infection of Zika virus (ZIKV) [22], DENV [23], and WNV [24].

The primary roles of the piRNA pathway are to silence transposons and maintain the integrity of the germline [25]. The key proteins of this pathway are members of the PIWI and Argonaut superfamilies, and the *Aedes albopictus* genome is known to encode seven PIWI proteins (Piwi 1–6 and Ago3) [26]. The piRNAs are 24–30 nt in length and usually show specific signatures features such as an uracil (U) bias at the first position of the antisense strand and an adenine (A) at the 10th position of the sense strand. Recent studies on mosquitoes have suggested a possible role of the piRNA pathway in antiviral defense. For example, the production of ping-pong-dependent virus-derived piRNAs (vpiRNAs) by the PIWI-5 and Ago3 proteins has been observed during viral infection, suggesting an antiviral function [27]. Furthermore, the knockdown of the PIWI-4 proteins in mosquito cell lines has been found to increase the viral titers in comparison to controls [28–30]. The exact antiviral mechanisms of the piRNA pathway are still unknown.

In this study, we infected the *Aedes albopictus* U4.4 cell line with either of two ISFVs from two distinct phylogenetic groups, LamV or HakV, to investigate the small RNA response. We found that both LamV and HakV elicited a strong virus-derived small interfering RNA (vsiRNA) response that increased over time and targeted the whole viral genome. However, only the LamV-infected cells produced viral pi-like RNAs (vpi-like RNAs), which were mainly derived from the antisense genome and did not show the characteristic ping-pong signatures.

## 2. Materials and Methods

### 2.1. *Cell Culture*

*Aedes albopictus* U4.4 cells (kindly provided by Associate Professor G. Pijlman, Wageningen University, Wageningen, The Netherlands) and C6/36 cells (Sigma-Aldrich, Darmstadt, Germany) were cultured at 28 °C in a Leibovitz L-15 medium (Gibco, Paisley, UK) supplemented with 10% fetal bovine serum (FBS) (Gibco, Paisley, UK), 10% tryptose
phosphate broth (TPB) (Gibco, Paisley, UK), amphotericin (250 µg/mL) (Gibco, Grand Island, NY, USA), and Pen Strep (100 U penicillin/mL and 100 µg streptomycin/mL) (Gibco, Grand Island, NY, USA).

## 2.2. Virus Stocks and Virus Titration

The ISFVs LamV (Strain: 2009/FI/Original) and HakV (Strain: 2005/FI/UNK) were obtained from the European Virus Archive—Global (EVAg). Virus stocks were propagated in C6/36 *Aedes albopictus* cells until a clear cytopathic effect (CPE) was shown (4 days post infection (d.p.i.)), and then supernatant was collected, centrifuged, and frozen at -80 °C. We did not obtain any reliable virus stock titers with the traditional methods. Thus, to quantify the ISFV stocks, plasmid standards containing the PCR target region of the virus were ordered (GeneScript Biotech, Leiden, The Netherlands). The plasmids were used to construct a qPCR standard curve, and a stock concentration was calculated as RNA copies/µL. RNA extractions and qRT-PCR protocols are described below.

## 2.3. In Vitro Infection

On the day of infection, U4.4 cells grown in 24-well plates to a confluence of 85–90% were used (approximately 350,000 cells per well). The cells were infected as triplicates for each timepoint. In brief, per infection, 35,000 RNA copies of either LamV or HakV were mixed with 200 µL of infection medium (Leibovitz L-15 medium containing 2% FBS and 10% TPB) and added to the respective wells. After one hour of incubation at 28 °C, the inoculum was discarded, and 500 µL of Leibovitz L-15 medium supplemented with 5% FBS, 10% TPB, and PEST were added. Mock infected cells were used as controls. The U4.4 cells were sampled every 24 h until 72 h post-infection (p.i.) for small RNA sequencing, and the supernatant was collected from 24 to 96 h p.i. for qPCRs.

## 2.4. RNA Extractions

The RNA used for growth curves and virus titration was extracted from 200 µL of supernatant/virus stocks in 750 µL of TRIzol™ (Thermo Fisher Scientific, Carlsbad, CA, USA) and homogenized. The aqueous phase obtained after the addition of chloroform and the subsequent centrifugation step was collected, diluted to 1:1 with freshly prepared 70% ethanol, and purified with GeneJet spin columns (Thermo Fisher Scientific, Vilnius, Lithuania). The RNA was eluted in 40 µL of nuclease-free water and stored at -80 °C until further processing.

Small RNA used for high-throughput sequencing was isolated from cells that were collected by adding 750 µL of TRIzol™ to respective wells. The aqueous phase was obtained in the same manner as described above, and the RNA was further purified using the mirVana™ PARIS™ Kit (Thermo Fisher Scientific, Vilnius, Lithuania) according to the protocol to isolate small RNA (<200 nt) provided by the manufacturer.

## 2.5. qPCR

First-strand cDNA was generated using the SuperScript™ III Reverse Transcription kit (Thermo Fisher Scientific, Carlsbad, CA, USA) with Random Hexamers (Thermo Fisher Scientific, Carlsbad, CA, USA) following the manufacturer's instructions, with an input of 5 µL of RNA in a total reaction volume of 20 µL. The qPCR was performed using the iTaq Universal SYBR® Green Supermix (Bio-Rad laboratories Inc., Hercules, CA, USA) with 2 µL of template cDNA and 0.5 µM of each corresponding virus primer (Table 1) in a total volume of 20 µL per reaction. The qPCR was carried out using the Bio-Rad® CFX96 real-time PCR system (Bio-Rad laboratories Inc., Hercules, CA, USA), with amplification conditions consisting of an initial denaturation at 95 °C for 30 s, followed by 40 cycles of denaturation at 95 °C for 7 s, and finally annealing/extension and plate reading at 60 °C for 30 s. A melt curve was generated starting at 60 °C with a 0.5 °C increase up to 96 °C.

Primer pairs for the qPCR were designed using the software Primer3 [31] to generate a product between 170 and 200 bp long and a Tm of 60 °C. Virus reference genomes were
obtained from the NCBI database and are listed together with the primer sequences in Table 1.

<table><caption>**Table 1.** Primer pairs used for the qPCR analysis.</caption><thead><tr><td>Primers</td><td>Binding Site</td><td>Sequence (5' → 3')</td><td>Ref</td></tr></thead><tbody><tr><td>LamV-F</td><td>4659–4678</td><td>TGGGTGTTACCGGGTTATGT</td><td>FJ606789</td></tr><tr><td>LamV-R</td><td>4845–4864</td><td>ACGTTCCATTCAGTTTCCAT</td><td></td></tr><tr><td>HakV-F</td><td>4655–4674</td><td>TGTGTTACGGTGGAAACTGG</td><td>NC030401</td></tr><tr><td>HakV-R</td><td>4842–4861</td><td>CAACTGGTTCTCCGTTGACA</td><td></td></tr></tbody></table>

## 2.6. *High-Throughput Sequencing of Small RNA*

Small RNA isolated from the infected U4.4 cells was quantified and quality-controlled with the 4150 TapeStation System using the RNA ScreenTape Analysis kit (Agilent Technologies, Santa Clara, CA, USA). Triplicates for each timepoint (24–72 h) were pooled and submitted to SciLifeLab (Stockholm, Sweden) for library preparation using the QIAseq miRNA low input kit (QIAGEN, Hilden, Germany), and libraries were sequenced on one flow cell of the NextSeq 2000 with a 101nt(Read1)-8nt(Jndex1) setup using the 'P2' flow cell, which generated 10–20 million reads per sample with a read length of 1 × 100 base-pairs (bps).

## 2.7. *Small RNA Sequence Analysis*

The generated small RNA sequence data (FASTQ files) were analyzed with a pipeline including the trimming of adaptors and the removal of bad quality reads using Trim Galore! (v0.6.6), with a setting to discard all reads below 17 nt and above 40 nt. The trimmed reads were mapped to the LamV (FJ606789) genome and the HakV (NC030401) genome using Bowtie (v1.2.3) [32], allowing for one mismatch. We used the -a, -best, -strata, and -all flags, which instructed Bowtie to only report those alignments in the best alignment stratum and to generate a FASTQ file of the mapped reads in addition to the SAM file. The FASTQ files were used to analyze the size distribution of the mapped reads, which was visualized with the R package ggplot2. The SAM output files were further analyzed with the MISIS-2 software (v2.6) to visualize the alignment and polarity distribution of small RNA to the viral genomes [33]. The ping-pong signatures of 27–30 nt piRNA were analyzed with WebLogo [34].

## 2.8. *Data Availability*

All sequencing data have been made publicly available at the NCBI Sequence Read Archive (SRA) under accession number PRJNA761074.

# 3. Results

To investigate the RNAi response to ISFV infection in mosquitoes, we infected the *Aedes albopictus* U4.4 cell line with LamV or HakV at a concentration of 175,000 RNA copies/mL. Mock infections were used as controls. Cells were collected at 24, 48, and 72 h p.i. and processed for small RNA sequencing. The retrieved sequence data were trimmed from adapters and bad quality reads, resulting in 10.4–16.9 million reads per sample, with read lengths between 17 and 40 nt. The size distributions of small RNAs were similar between all samples, with the highest peak observed at 21 nt (22.3–26.3% of total reads) and the second highest peak observed at 22 nt (13.7–15.3% of total reads) (Figure 1). These two peaks most likely corresponded to the populations of siRNAs and miRNAs in the cell. We could also observe an elevation in the proportion of reads between 26 and 30 nt, which is in the size range of piRNAs (Figure 1).
<figure><img src="image_3.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>C</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Small RNA profiles of U4.4 cells at 24, 48, and 72 h. (A) Small RNA profile of total reads in mock infected cells; (B) small RNA profile of total read in LamV-infected cells; (C) small RNA profile of total reads in HakV-infected cells.</figcaption></figure>

To confirm active viral infection and replication, the supernatant was collected, further processed by RNA extraction and cDNA synthesis, and analyzed with qPCR. The data were used to calculate a growth curve over time (Figure 2). Both LamV and HakV showed steep virus growth during the first 48 h before plateauing. Furthermore, LamV showed a higher replication compared to HakV, with approximately one log more RNA copies at 48 h p.i. and the later timepoints recorded (Figure 2).

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Growth curves of LamV and HakV in U4.4 cells, shown as RNA copies over time. Dots show the mean RNA copy numbers with a standard deviation between the biological triplicates.</figcaption></figure>
### 3.1. Viral Small RNA Profiles in U4.4 Infected with LamV

To further investigate the production of virus-derived small RNAs (vsRNAs), the trimmed sequence data were mapped to the viral genomes. Results from U4.4 cells infected with LamV showed that few reads were virus-specific at 24 h p.i.: only 0.001% (n = 270) of the reads mapped to the LamV genome. However, the proportion of virus-specific reads increased over time, with 0.06% (n = 6096) observed at 48 h p.i. and 0.17% (n = 24,179) observed at 72 h p.i. The majority of the mapped reads had a read length of 21 nt and most likely corresponded to vsiRNA (32.2% at 24 h, 36.7% at 48 h, and 32.9% at 72 h) (Figure 3). At 48 and 72 h p.i. a shoulder of reads in the size range of 27–30 nt that resembled vpiRNA was observed (Figure 3c,d). These 27–30 nt long reads were examined for ping-pong signatures, with antisense piRNA having a 1U bias (shown as T in this dataset) and the sense piRNA having a 10A bias. The 1U characteristics for Piwi-5-bound piRNAs were identified; however, the 10A bias was not observed (Figure 4a–d).

The vsRNAs from the LamV-infected cells mapped across the entire LamV genome, with a few predominant hotspots with a higher coverage. These hotspots were most distinct at 72 h p.i. (Figure 3d) and included four positions in the genome. Starting from the 5' end of the LamV genome, the first hotspot was at position 423, which was in the gene encoding for the capsid protein. The second hotspot was at position 3854, i.e., in the NS2A gene. The third hotspot at position 6808 was located in the NS4A protein, and the fourth hotspot at position 10,437 was positioned in the 3'-UTR of the genome. Moreover, the fourth hotspot had a coverage of over 900 reads, mainly at the size of 39 nt, which could been seen as a small hump in the read length distribution (Figure 3d). These reads were extracted and analyzed with WebLogo, and it was found that the majority consisted of the same sequence (Figure 4e). This suggested that these 39 nt long reads were not degraded products; however, their function and biogenesis are uncertain. Interestingly, a large majority of the vsiRNA mapped to the sense strand of the LamV genome: 79.2% at 48 h and 82.4% at 72 h p.i. (Figure 3, left panel). This suggests that Ago-2 has a preferred bias for the antisense strand as a template and that the RISC complex mainly targets and degrades the sense genomes during LamV infection.

### 3.2. Viral Small RNA Profiles in U4.4 Infected with HakV

The results from HakV-infected U4.4 cells showed a slightly higher production of vsRNAs at 24 h p.i. compared to LamV-infected cells, where 0.01% (n = 1296) of the reads mapped to the HakV genome. The observed proportion of vsRNAs was similar at 48 h p.i., 0.06% (n = 7324), and lower at 72 h p.i., 0.08% (n = 11,179) compared to LamV-infected cells. Moreover, when analyzing the read length distribution, a high proportion of the reads was in the size range of siRNA (21 nt), with 48.4% observed at 24 h p.i, 60.9% observed at 48 h p.i., and 66.5% observed at 72 h p.i. The shoulder of piRNAs (27–30 nt) observed in the LamV-infected cells was absent in those cells infected with HakV (Figure 5). The lack of vpiRNA suggests that HakV does not trigger the proteins responsible for vpiRNA amplification such as Ago3 and PIWI-5. An analysis of the read alignment showed that the vsRNA mapped across the entire HakV genome, with slightly more coverage observed at the 3' end (7840–10,158 nt), which corresponded to the RNA-dependent RNA polymerase (NS5) gene region. The strand polarity showed a more even distribution between sense (68.6%) and antisense (31.4%) compared to LamV-infected U4.4 cells.
A

Lammi virus

<figure><img src="image_7.png" /></figure>

B

LamV 24h

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

C

LamV 48h

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

D

LamV 72h

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

**Figure 3.** Viral small RNA profiles of U4.4 cells after LamV infection. (A) Schematic representation of the LamV genome and vsRNA profiles at 24 h (B), 48 h (C), and 72 h (D) p.i. Figures to the left show vsRNA reads along the LamV genome as a histogram. The positive values (blue) are counts of sRNAs mapped to the sense strand, and the negative values (red) are those mapped to the antisense strand. Figures to the right show the size distribution of reads that mapped to the LamV genome.
<figure><img src="image_14.png" /><figcaption>**Figure 4.** Nucleotide biases of the 27–30 and 39 nt long reads at 72 h p.i. of LamV. The sequences are in DNA format, so U is written as T. Reads derived from the antisense are shown in the left panel. Reads derived from sense strand are shown in the right panel: (A) 27 nt reads, (B) 28 nt reads, (C) 29 nt reads, (D) 30 nt reads, and (E) 39 nt long reads from the fourth hotspot in the LamV genome.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 5.** Viral small RNA profiles of U4.4 cells after HakV infection. (A) Schematic representation of the HakV genome. vsRNA profiles at 24 h (B), 48 h (C), and 72 h (D) p.i. Figures to the left show vsRNA reads along the HakV genome as a histogram. The positive values are counts of vsRNAs mapped to the sense strand, and the negative values are those mapped to the antisense strand. Figures to the right show the size distribution of reads that mapped to the HakV genome.</figcaption></figure>
4. Discussion

RNAi-mediated antiviral immunity is thought to be key in the defense against viral in-
fection in mosquitoes. The production of virus-specific small RNAs has been demonstrated
in mosquitoes and mosquito-derived cell lines for viruses in all of the major arbovirus
families: *Peribunyaviridae*, *Togaviridae*, and *Flaviviridae* [18,35–40]. However, although
mosquitoes and mosquito-derived cell lines are often persistently infected with this group
of viruses, the RNAi response to ISVs is not as well-studied. To our knowledge, only a
handful of studies have investigated the effect of ISVs on the RNAi response in mosquitoes.
One study developed an integrated mosquito small RNA genomics resource and included
data from mosquito cell lines persistently infected with ISVs, including the ISFV cell-fusing
agent virus [41]. Another study investigated the RNAi response of the mosquito U4.4
(*Aedes albopictus*), *Aag2* (*Aedes aegypti*), and CT (*Culex tarsalis*) cell lines, which were shown
to be persistently infected with different ISVs. For example, the U4.4 cells were persistently
infected with Culex Y virus (CYV) and displayed a potent siRNA response against it [42].
In the present study, we investigated the RNAi response to acute ISFV infection over time
in the aforementioned U4.4 cells. At the early time point (24 h p.i.), both the LamV- and
HakV-infected U4.4 cells showed low amounts of vsRNAs (0.01–0.001% of total reads),
which could have been because of the relatively low amount of virus or because the siRNA
and piRNA pathways had not yet responded. Results from the qPCR analysis of the
supernatant showed the replication of both LamV and HakV at this time point (Figure 2).
Moreover, the proportion of vsRNA steadily increased over time, with 0.06% observed at
48 h p.i. and 0.17–0.08% observed at 72 h p.i. Interestingly, the increasing proportion of
vsRNA correlated with the plateau of the virus growth curves (Figure 2) and might have
been a sign of RNAi-induced interference with the virus replication.

The siRNA pathway is regarded as the most important antiviral defense mecha-
nism, and in line with this, our experiments showed that it was the most abundant
population of sRNAs between all groups and time points. Hence, this was observed
in the general sRNA distribution—including the mock infected cells, most likely due to
the persistent CYV infection (Figure 1)—and in the distribution of ISFV-specific sRNAs
(Figures 3 and 5). The LamV-infected cells had proportions of between 32.2 and 36.7% of
vsiRNAs, and the HakV-infected cells had slightly higher proportions of between 48.4 and
66.5%. The vsiRNAs from LamV-infected cells mapped along the whole genome, with no
particular coldspots but one predominant hotspot that correspond to the NS4AB proteins
(Figure S1). We also observed that most of the vsiRNAs mapped to the sense strand of
the genome, which could be a strategy for the efficient restriction of virus replication.
The vsiRNAs of HakV-infected cells also mapped to the whole genome, with a slightly
higher coverage at the 3' end (7840–10,158 nt) corresponding to the gene coding for the
NS5 protein (Figure S1). The strand polarity distribution was close to even, with slightly
more vsiRNAs mapping to the sense strand. Although the siRNA response is regarded as
the main antiviral immune response, it has also been postulated that siRNA is necessary
for a persistent viral infection in mosquitoes [43,44]. The hypothesis behind this is that the
siRNA response keeps the viral load in the host cell at a tolerable level, thereby sustaining
a persistent infection. Many of the ISVs have been shown to not only persistently infect
mosquitoes but also vertically transmit from mother to offspring [2–4]. Our data showed
that the mosquito U4.4 cell line is able to elicit a strong virus-specific siRNA response to
the studied ISFVs, which could support lifelong infection in mosquitoes.

In the LamV-infected cells, we observed a hotspot positioned in the 3'-UTR with a
coverage of 900 reads, mainly at the size of 39 nt. These reads mainly consisted of the same
sequence (5'GGGAGTCAGGCCTAAATGCCACCGGATGATAGTAGACGG), suggesting
that the reads were not degraded products (Figure 4e). The BLASTn search of the consensus
sequence showed that this sequence could also be found in the 3'-UTR of another ISFV
named Chaoyang virus (NC 017086), indicating that it may be a conserved sequence that
can exist in the UTR of other ISFVs as well. However, the origin, function, and biogenesis
of these 39 nt reads are uncertain.
Apart from the detection of vsiRNAs, putative vpiRNAs were identified in the LamV-infected cells. Earlier studies in an *Aedes aegypti* cell line have shown that the production of ping-pong-dependent vpiRNA relies on the PIWI-5 and Ago-3 proteins. The knockdown of either of these proteins in the Aag2 cell line infected with Sindbis virus or DENV-2 showed a massive decrease in vpiRNAs [27,45]. The ping-pong amplification of piRNAs is a two-step amplification mechanism where the PIWI-5 protein, loaded with a primary piRNA, cleaves a complementary target RNA. The 3' cleavage product is transferred to the protein Ago-3 and used as a template to cleave a complementary target RNA. This generates a new piRNA precursor for the Piwi-5 protein, which gives rise to the same primary piRNA sequence that initiated the amplification. Therefore, this type of amplification predicts an even distribution of piRNAs derived from both strands [27,46].

In the data from our U4.4 cells infected with LamV, we did not observe equal amounts of vpiRNA from both strands. The majority of the putative vpiRNAs (26–30 nt) were derived from the antisense strand, with the ping-pong characteristic nucleotide bias of an U at the first position. However, the few putative vpiRNAs derived from the sense strand did not show the nucleotide bias of an A in the tenth position characteristic for Ago-3 bound piRNA. Moreover, because of the unequal strand distribution, we could not analyze the overlap probability and look for the 10 nt overlap, which is significant for ping-pong-amplified piRNA (Figures 3 and 4 and Figure S2). Similar observations have been made in studies where U4.4 cells were infected with WNV [42] and where Aag2 cells were infected with ZIKV [47]. In these studies, reads in the size range of 25–30 nt mainly mapped to the sense strand but with no 1U bias. These observations suggest that the identified putative LamV vpiRNAs are not ping-pong-dependent, and we hypothesize that they could be primary vpiRNAs cleaved by an ortholog to the zucchini proteins in *Drosophila* or some unknown endonuclease using either exogenous viral RNA or pre-primary RNA transcribed from viral-derived cDNA (vDNA) [30,48]. Further evaluation is needed to understand their function and biogenesis, as well as whether they interact with any of the PIWI proteins. Surprisingly, U4.4 cells infected with HakV did not elicit any production of these putative vpiRNAs (Figure 5), which suggests that the production of vpiRNAs is virus-dependent. HakV belongs to the group of classical ISFVs that are very distinct from all known flaviviruses [7]. It would be interesting to further investigate whether other viruses in the classical ISFV group show similar results, which could reveal factors important for putative vpiRNA production. We speculate that virus-specific sequence elements or structures trigger the production of vpiRNAs by an unknown Piwi-dependent amplification mechanism.

**5. Conclusions**

In conclusion, we showed that two ISFVs from two distinct phylogenetic groups can
trigger the RNAi response during an acute infection of the *Aedes albopictus* U4.4 cell line.
Both viruses elicited a strong vsiRNA response that increased over time and targeted the
whole viral genome. Furthermore, infection with LamV (which is closely related to the
pathogenic dual-host flaviviruses) triggered the production of putative primary piRNAs,
while infection with HakV, which is more distantly related, did not. This suggests that the
mosquito piRNA pathway is virus-specific and might need specific sequence elements or
structures. These results contribute to our understanding of mosquito antiviral immunity,
small RNA machineries, and how ISVs affect mosquitoes.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/
10.3390/v13112181/s1: Figure S1: Mapping of viral siRNA 72 h after ISFVs infection; Figure S2:
Distribution of virus-derived piRNAs 72 h after LamV infection.

**Author Contributions:** P.Ö., J.H., J.C.H. and A.-L.B. conceived and designed the experiments; P.Ö. performed the experiments; P.Ö., J.H. and A.-L.B. analyzed the data; P.Ö. wrote and prepared the original draft; J.H., J.C.H. and A.-L.B. reviewed and edited the manuscript. All authors have read and agreed to the published version of the manuscript.
**Funding:** This study was supported by the Swedish research Council VR (2016-01251), the Swedish University of Agricultural Sciences (SLU), the Vice-chancellor junior career grant awarded to A.L.B., and the Kungl. Skogs- och Lantbruksakademinen (GFS2016-0033).

**Data Availability Statement:** All sequencing data have been made publicly available at the NCBI Sequence Read Archive (SRA) under accession number PRJNA761074.

**Acknowledgments:** The authors would like to acknowledge the support of the National Genomics Infrastructure in Stockholm funded by Science for Life Laboratory, the Knut and Alice Wallenberg Foundation, the Swedish Research Council, and the SNIC/Uppsala Multidisciplinary Center for Advanced Computational Science for assistance with massively parallel sequencing and access to the UPPMAX computational infrastructure. The SNIC through Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) under project SNIC 2017/7-290 is acknowledged for providing computational resources. This publication was supported by the European Virus Archive goes Global (EVAg) project, which has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No 653316.

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

## References

1. Nasar, F.; Gorchakov, R.V.; Tesh, R.B.; Weaver, S.C. Eilat virus host range restriction is present at multiple levels of the virus life cycle. *J. Virol.* **2015**, *89*, 1404–1418. [CrossRef]

2. Saiyasombat, R.; Bolling, B.G.; Brault, A.C.; Bartholomay, L.C.; Blitvich, B.J. Evidence of efficient transovarial transmission of Culex flavivirus by Culex pipiens (Diptera: Culicidae). *J. Med. Entomol.* **2011**, *48*, 1031–1038. [CrossRef]

3. Bolling, B.G.; Eisen, L.; Moore, C.G.; Blair, C.D. Insect-specific flaviviruses from Culex mosquitoes in Colorado, with evidence of vertical transmission. *Am. J. Trop. Med. Hyg.* **2011**, *85*, 169–177. [CrossRef]

4. Haddow, A.D.; Guzman, H.; Popov, V.L.; Wood, T.G.; Widen, S.G.; Haddow, A.D.; Tesh, R.B.; Weaver, S.C. First isolation of Aedes flavivirus in the Western Hemisphere and evidence of vertical transmission in the mosquito Aedes (Stegomyia) albopictus (Diptera: Culicidae). *Virology* **2013**, *440*, 134–139. [CrossRef] [PubMed]

5. Marklewitz, M.; Zirkel, F.; Kurth, A.; Drosten, C.; Junglen, S. Evolutionary and phenotypic analysis of live virus isolates suggests arthropod origin of a pathogenic RNA virus family. *Proc. Natl. Acad. Sci. USA* **2015**, *112*, 7536–7541. [CrossRef]

6. Li, C.X.; Shi, M.; Tian, J.H.; Lin, X.D.; Kang, Y.J.; Chen, L.J.; Qin, X.C.; Xu, J.; Holmes, E.C.; Zhang, Y.Z. Unprecedented genomic diversity of RNA viruses in arthropods reveals the ancestry of negative-sense RNA viruses. *eLife* **2015**, *4*. [CrossRef]

7. Blitvich, B.J.; Firth, A.E. Insect-specific flaviviruses: A systematic review of their discovery, host range, mode of transmission, superinfection exclusion potential and genomic organization. *Viruses* **2015**, *7*, 1927–1959. [CrossRef] [PubMed]

8. Huhtamo, E.; Putkuri, N.; Kurkela, S.; Manni, T.; Vaheri, A.; Vapalahti, O.; Uzcategui, N.Y. Characterization of a novel flavivirus from mosquitoes in northern europe that is related to mosquito-borne flaviviruses of the tropics. *J. Virol.* **2009**, *83*, 9532–9540. [CrossRef]

9. Calzolari, M.; Ze-Ze, L.; Ruzek, D.; Vazquez, A.; Jeffries, C.; Defilippo, F.; Osorio, H.C.; Kilian, P.; Ruiz, S.; Fooks, A.R.; et al. Detection of mosquito-only flaviviruses in Europe. *J. Gen. Virol.* **2012**, *93*, 1215–1225. [CrossRef]

10. Blair, C.D. Mosquito RNAi is the major innate immune pathway controlling arbovirus infection and transmission. *Future Microbiol.* **2011**, *6*, 265–277. [CrossRef] [PubMed]

11. Galiana-Arnoux, D.; Dostert, C.; Schneemann, A.; Hoffmann, J.A.; Imler, J.L. Essential function in vivo for Dicer-2 in host defense against RNA viruses in drosophila. *Nat. Immunol.* **2006**, *7*, 590–597. [CrossRef]

12. Matranga, C.; Tomari, Y.; Shin, C.; Bartel, D.P.; Zamore, P.D. Passenger-strand cleavage facilitates assembly of siRNA into Ago2-containing RNAi enzyme complexes. *Cell* **2005**, *123*, 607–620. [CrossRef] [PubMed]

13. Rand, T.A.; Petersen, S.; Du, F.; Wang, X. Argonaute2 cleaves the anti-guide strand of siRNA during RISC activation. *Cell* **2005**, *123*, 621–629. [CrossRef] [PubMed]

14. Campbell, C.L.; Keene, K.M.; Brackney, D.E.; Olson, K.E.; Blair, C.D.; Wilusz, J.; Foy, B.D. Aedes aegypti uses RNA interference in defense against Sindbis virus infection. *BMC Microbiol.* **2008**, *8*, 47. [CrossRef]

15. Sanchez-Vargas, I.; Scott, J.C.; Poole-Smith, B.K.; Franz, A.W.; Barbosa-Solomieu, V.; Wilusz, J.; Olson, K.E.; Blair, C.D. Dengue virus type 2 infections of Aedes aegypti are modulated by the mosquito's RNA interference pathway. *PLoS Pathog* **2009**, *5*, e1000299. [CrossRef] [PubMed]

16. Keene, K.M.; Foy, B.D.; Sanchez-Vargas, I.; Beaty, B.J.; Blair, C.D.; Olson, K.E. RNA interference acts as a natural antiviral response to O'nyong-nyong virus (Alphavirus; Togaviridae) infection of Anopheles gambiae. *Proc. Natl. Acad. Sci. USA* **2004**, *101*, 17240–17245. [CrossRef]

17. Franz, A.W.; Sanchez-Vargas, I.; Adelman, Z.N.; Blair, C.D.; Beaty, B.J.; James, A.A.; Olson, K.E. Engineering RNA interference-based resistance to dengue virus type 2 in genetically modified Aedes aegypti. *Proc. Natl. Acad. Sci. USA* **2006**, *103*, 4198–4203. [CrossRef]
18. Brackney, D.E.; Beane, J.E.; Ebel, G.D. RNAi targeting of West Nile virus in mosquito midguts promotes virus diversification. *PLoS Pathog.* **2009**, *5*, e1000502. [CrossRef]

19. Lee, Y.S.; Nakahara, K.; Pham, J.W.; Kim, K.; He, Z.; Sontheimer, E.J.; Carthew, R.W. Distinct roles for Drosophila Dicer-1 and Dicer-2 in the siRNA/miRNA silencing pathways. *Cell* **2004**, *117*, 69–81. [CrossRef]

20. Puthiyakunnon, S.; Yao, Y.; Li, Y.; Gu, J.; Peng, H.; Chen, X. Functional characterization of three MicroRNAs of the Asian tiger mosquito, Aedes albopictus. *Parasit Vectors* **2013**, *6*, 230. [CrossRef]

21. Slonchak, A.; Hussain, M.; Torres, S.; Asgari, S.; Khromykh, A.A. Expression of mosquito microRNA Aae-miR-2940-5p is downregulated in response to West Nile virus infection to restrict viral replication. *J. Virol.* **2014**, *88*, 8457–8467. [CrossRef] [PubMed]

22. Saldana, M.A.; Etebari, K.; Hart, C.E.; Widen, S.G.; Wood, T.G.; Thangamani, S.; Asgari, S.; Hughes, G.L. Zika virus alters the microRNA expression profile and elicits an RNAi response in Aedes aegypti mosquitoes. *PLoS Negl. Trop. Dis.* **2017**, *11*, e0005760. [CrossRef] [PubMed]

23. Campbell, C.L.; Harrison, T.; Hess, A.M.; Ebel, G.D. MicroRNA levels are modulated in Aedes aegypti after exposure to Dengue-2. *Insect Mol. Biol.* **2014**, *23*, 132–139. [CrossRef]

24. Skalsky, R.L.; Vanlandingham, D.L.; Scholle, F.; Higgs, S.; Cullen, B.R. Identification of microRNAs expressed in two mosquito vectors, Aedes albopictus and Culex quinquefasciatus. *BMC Genom.* **2010**, *11*, 119. [CrossRef] [PubMed]

25. Rozhkov, N.V.; Hammell, M.; Hannon, G.J. Multiple roles for Piwi in silencing Drosophila transposons. *Genes Dev.* **2013**, *27*, 400–412. [CrossRef] [PubMed]

26. Marconcini, M.; Hernandez, L.; Iovino, G.; Houe, V.; Valerio, F.; Palatini, U.; Pischedda, E.; Crawford, J.E.; White, B.J.; Lin, T.; et al. Polymorphism analyses and protein modelling inform on functional specialization of Piwi clade genes in the arboviral vector Aedes albopictus. *PLoS Negl. Trop. Dis.* **2019**, *13*, e0007919. [CrossRef]

27. Miesen, P.; Girardi, E.; van Rij, R.P. Distinct sets of PIWI proteins produce arbovirus and transposon-derived piRNAs in Aedes aegypti mosquito cells. *Nucleic Acids Res.* **2015**, *43*, 6545–6556. [CrossRef] [PubMed]

28. Schnettler, E.; Donald, C.L.; Human, S.; Watson, M.; Siu, R.W.C.; McFarlane, M.; Fazakerley, J.K.; Kohl, A.; Fragkoudis, R. Knockdown of piRNA pathway proteins results in enhanced Semliki Forest virus production in mosquito cells. *J. Gen. Virol.* **2013**, *94*, 1680–1689. [CrossRef] [PubMed]

29. Varjak, M.; Maringer, K.; Watson, M.; Sreenu, V.B.; Fredericks, A.C.; Pondeville, E.; Donald, C.L.; Sterk, J.; Kean, J.; Vazeille, M.; et al. Aedes aegypti Piwi4 Is a Noncanonical PIWI Protein Involved in Antiviral Responses. *mSphere* **2017**, *2*, e00144-17. [CrossRef]

30. Tassetto, M.; Kunitomi, M.; Whitfield, Z.J.; Dolan, P.T.; Sanchez-Vargas, I.; Garcia-Knight, M.; Ribiero, I.; Chen, T.; Olson, K.E.; Andino, R. Control of RNA viruses in mosquito cells through the acquisition of vDNA and endogenous viral elements. *eLife* **2019**, *8*, e41244. [CrossRef]

31. Rozen, S.; Skaletsky, H. Primer3 on the WWW for general users and for biologist programmers. *Methods Mol. Biol.* **2000**, *132*, 365–386.

32. Langmead, B.; Trapnell, C.; Pop, M.; Salzberg, S.L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. *Genome Biol.* **2009**, *10*, R25. [CrossRef]

33. Seguin, J.; Otten, P.; Baerlocher, L.; Farinelli, L.; Pooggin, M.M. MISIS-2: A bioinformatics tool for in-depth analysis of small RNAs and representation of consensus master genome in viral quasispecies. *J. Virol. Methods* **2016**, *233*, 37–40. [CrossRef]

34. Crooks, G.E.; Hon, G.; Chandonia, J.M.; Brenner, S.E. WebLogo: A sequence logo generator. *Genome Res.* **2004**, *14*, 1188–1190. [CrossRef] [PubMed]

35. Leger, P.; Lara, E.; Jagla, B.; Sismeiro, O.; Mansuroglu, Z.; Coppee, J.Y.; Bonnefoy, E.; Bouloy, M. Dicer-2- and Piwi-mediated RNA interference in Rift Valley fever virus-infected mosquito cells. *J. Virol.* **2013**, *87*, 1631–1648. [CrossRef] [PubMed]

36. Morazzani, E.M.; Wiley, M.R.; Murreddu, M.G.; Adelman, Z.N.; Myles, K.M. Production of virus-derived ping-pong-dependent piRNA-like small RNAs in the mosquito soma. *PLoS Pathog.* **2012**, *8*, e1002470. [CrossRef]

37. Siu, R.W.; Fragkoudis, R.; Simmonds, P.; Donald, C.L.; Chase-Topping, M.E.; Barry, G.; Attarzadeh-Yazdi, G.; Rodriguez-Andres, J.; Nash, A.A.; Merits, A.; et al. Antiviral RNA interference responses induced by Semliki Forest virus infection of mosquito cells: Characterization, origin, and frequency-dependent functions of virus-derived small interfering RNAs. *J. Virol.* **2011**, *85*, 2907–2917. [CrossRef]

38. Hess, A.M.; Prasad, A.N.; Ptitsyn, A.; Ebel, G.D.; Olson, K.E.; Barbacioru, C.; Monighetti, C.; Campbell, C.L. Small RNA profiling of Dengue virus-mosquito interactions implicates the PIWI RNA pathway in anti-viral defense. *BMC Microbiol.* **2011**, *11*, 45. [CrossRef]

39. Myles, K.M.; Morazzani, E.M.; Adelman, Z.N. Origins of alphavirus-derived small RNAs in mosquitoes. *RNA Biol.* **2009**, *6*, 387–391. [CrossRef] [PubMed]

40. Dietrich, I.; Shi, X.; McFarlane, M.; Watson, M.; Blomstrom, A.L.; Skelton, J.K.; Kohl, A.; Elliott, R.M.; Schnettler, E. The Antiviral RNAi Response in Vector and Non-vector Cells against Orthobunyaviruses. *PLoS Negl. Trop. Dis.* **2017**, *11*, e0005272. [CrossRef] [PubMed]

41. Ma, Q.; Srivastav, S.P.; Gamez, S.; Dayama, G.; Feitosa-Suntheimer, F.; Patterson, E.I.; Johnson, R.M.; Matson, E.M.; Gold, A.S.; Brackney, D.E.; et al. A mosquito small RNA genomics resource reveals dynamic evolution and host responses to viruses and transposons. *Genome Res.* **2021**, *31*, 512–528. [CrossRef]
42. Goertz, G.P.; Miesen, P.; Overheul, G.J.; van Rij, R.P.; van Oers, M.M.; Pijlman, G.P. Mosquito Small RNA Responses to West Nile and Insect-Specific Virus Infections in Aedes and Culex Mosquito Cells. *Viruses* **2019**, *11*, 271. [CrossRef] [PubMed]

43. Goic, B.; Stapleford, K.A.; Frangeul, L.; Doucet, A.J.; Gausson, V.; Blanc, H.; Schemmel-Jofre, N.; Cristofari, G.; Lambrechts, L.; Vignuzzi, M.; et al. Virus-derived DNA drives mosquito vector tolerance to arboviral infection. *Nat. Commun.* **2016**, *7*, 12410. [CrossRef] [PubMed]

44. Poirier, E.Z.; Goic, B.; Tome-Poderti, L.; Frangeul, L.; Boussier, J.; Gausson, V.; Blanc, H.; Vallet, T.; Loyd, H.; Levi, L.I.; et al. Dicer-2-Dependent Generation of Viral DNA from Defective Genomes of RNA Viruses Modulates Antiviral Immunity in Insects. *Cell Host Microbe* **2018**, *23*, 353–365 e8. [CrossRef]

45. Miesen, P.; Ivens, A.; Buck, A.H.; van Rij, R.P. Small RNA Profiling in Dengue Virus 2-Infected Aedes Mosquito Cells Reveals Viral piRNAs and Novel Host miRNAs. *PLoS Negl. Trop. Dis.* **2016**, *10*, e0004452. [CrossRef] [PubMed]

46. Brennecke, J.; Aravin, A.A.; Stark, A.; Dus, M.; Kellis, M.; Sachidanandam, R.; Hannon, G.J. Discrete small RNA-generating loci as master regulators of transposon activity in Drosophila. *Cell* **2007**, *128*, 1089–1103. [CrossRef]

47. Varjak, M.; Donald, C.L.; Mottram, T.J.; Sreenu, V.B.; Merits, A.; Maringer, K.; Schnettler, E.; Kohl, A. Characterization of the Zika virus induced small RNA response in Aedes aegypti cells. *PLoS Negl. Trop. Dis.* **2017**, *11*, e0006010. [CrossRef]

48. Mohn, F.; Handler, D.; Brennecke, J. Noncoding RNA. piRNA-guided slicing specifies transcripts for Zucchini-dependent, phased piRNA biogenesis. *Science* **2015**, *348*, 812–817. [CrossRef]