Aedes aegypti Vero strain F20 was used in this study. The strain was originally collected in Vero Beach, Florida in 2015. Mosquitoes were reared at 28°C and 60-80% relative humidity in a climate-controlled room with a light: dark cycle of 14:10 hours. Upon hatching, larvae were separated into pans (35cm x 20cm x 6cm) with approximately 200 larvae per pan and fed a mixture of 0.5% brewer’s yeast and 0.5% liver powder (MP Biomedicals, Solon, OH, USA). Adults were fed with 20% sucrose solution-soaked cotton rolls (Carolina Absorbent Cotton, Charlotte, NC, USA). Mosquito eggs were collected to maintain the colony as described in (38), and female mosquitoes were blood-fed from chickens following standard approved protocols (IACUC protocol 201807682).

Four- to five- day old adult female mosquitoes were transferred to 16 oz cardboard cartons (WebstaurantStore, Lancaster, PA, USA) one day before infection and placed in an incubator at 28°C and 60-80% relative humidity. The sucrose solution-soaked cotton rolls were removed and only water-soaked cotton rolls provided. In the virus-infected group, DENV-2 strain New Guinea C (GenBank Accession # KM204118) was incubated at a multiplicity of infection of 0.01 viruses per cell in Vero cells (African green monkey kidney cell line) for 5 days. Supernatant from the virus culture was mixed with defibrinated bovine blood (Hemostat, Dixon, CA, USA) and fed to female mosquitoes with an artificial feeding apparatus (Hemotek, Lancashire, United Kingdom) with membranes made from sausage casings. In the control group, bovine blood was mixed with M199 culture media (Corning, Manassas, VA, USA) from uninfected Vero cells 5 days after splitting following the same procedure as the infected group. Three fully blood-engorged mosquitoes were collected at 1-hour post-feeding to quantify the ingested DENV titer. After the removal of non-blood fed mosquitoes across both feeding group, the remaining mosquitoes were placed into cartons and maintained in the incubator at 28°C and 60-80% relative humidity. 20% sucrose solution-soaked cotton rolls were provided and changed every day. At 60 hours post-infection, 180 mosquitoes from both treatment and control groups were immobilized on ice. Midguts were dissected in sterile saline solution and immediately transferred to an eppendorf tube on dry ice. All the tissue samples were stored at -80°C.

A total of 90 midguts from treatment and control groups, separately, were divided in triplicate sample each contained with 30 midguts. RNA from a total of 6 samples was extracted using the standard Trizol protocol (38, 39). The three biological replications of control and three infected group samples were sent to the Novogene Bioinformatics Technology Co. (Durham, NC, USA) for library generation and Illumina HiSeq6000 paired-end 150bp sequencing.

All RNA samples passed quality control for RNA purity, integrity and potential contamination before library construction, and mRNA was enriched using oligo dT beads. The enriched mRNA was fragmented and cDNA synthesized using random hexamers. After cDNA was synthesized, a custom second-strand synthesis Illumina buffer was added with dNTPs, RNase H, and Escherichia coli polymerase I to generate the second strand through nick-translation. The cDNA library was finalized after the application of terminal repair, poly-A tailing, sequencing adapter ligation, size selection, and PCR enrichment. The cDNA libraries were quantified and placed into Illumina sequencer to generate the sequencing reads (Novogene Bioinformatics Technology Co., Durham, NC, USA).

The data from Illumina sequencing were transformed and recorded in a FASTQ file. Quality control indicated the clean reads were more than 98% and the error rate was below 0.03%. The HISAT2 version 2.1.0-beta was used to map the filtered sequenced reads to Aedes aegypti LVP_AGWG AaegL5 chromosome and transcripts reference files available from VectorBase (https://vectorbase.org/) (40). Gene expression level was measured as the Fragments Per Kilobase of transcript sequence per Millions of base pairs sequenced (FPKM) values. HTSeq software version v0.6.1 was used to calculate FPKM.

To analyze the difference of FPKM between control and infected groups, T-tests were used to determine if average expression levels between the control and infected groups were equal. We then calculated the false discovery rate (FDR) with the Benjamini-Hochberg procedure (41) implemented in Python SciPy version 1.5.4 (https://scipy.org/), with a FDR under 0.2 considered as a significant difference. This method has been routinely used over the classic correction methods such as Bonferroni (42) or Sidak (43) in high-throughput science allowing discovery of novel genes or genic regions with relatively small sample sizes (44, 45).

To identify the function of gene, several approaches were applied. Gene IDs were annotated using VectorBase (https://vectorbase.org/) and Uniprot (https://www.uniprot.org/) to obtain gene names and full-length sequences. Gene ontology (GO) terms were also used to estimate gene function, as well as homology with putative orthologs from other mosquito species. Sequences were then examined using BLASTx (https://blast.ncbi.nlm.nih.gov/Blast.cgi) to match identities and characterize conserved domains. Information on key differentially expressed genes was also collected and predicted through literature review.

To assess the impact of DENV infection on Ae. aegypti immune pathways, we compiled a list of 352 genes that fall into 23 gene families or 4 functional groups, all implicated in the innate immune response of dipteran insects to bacterial, viral, or parasitic infections. This list was derived from the ImmunoDB database (http://cegg.unige.ch/Insecta/immunodb), which categorizes loci into families based on homology among Ae. aegypti, Anopheles gambiae, and Drosophila melanogaster (46). Gene accession numbers and annotations were updated using the Ae. aegypti genome (AaegL5.3) on VectorBase (www.vectorbase.org). Accession numbers of previously identified genes were searched, and their functional roles were confirmed by at least one of three sources: protein domain information, orthology/paralogy with a gene(s) of similar function, or gene ontology (GO) terms. In addition, VectorBase was searched using the name of each gene family with results filtered by the organism (Ae. aegypti, genome release AaegL 5.3). The names of immunity genes and gene families in our list are in accordance with annotations from ImmunoDB. The final list was intended to err on the side of inclusivity and may include genes without a confirmed role in innate immunity but are closely related to loci (i.e., belong to the same gene family) with empirical support for a role in some aspect of the immune response.

RNA from the same samples used to generate the transcriptome was used to validate expression profiles of a subset of genes determined to have differential expression due to DENV infection via RNA-Seq. Gene expression was quantified using a Bio-Rad CFX96™ Real-Time PCR system and the iTaq Universal SYBR Green One-step kit (Bio-Rad, Hercules, CA, USA), with specific primer sets designed for each gene ( Table S3 ). Four candidate genes were selected from both upregulated and downregulated expression groups and from the two major functional categories we identified ( Tables 1 , 2 ). The PCR products from specific primer sets were sequenced to confirm the targets were correct. Ae. aegypti ribosomal protein S7 gene (GenBank Accession # AY380336) was picked as a control for standardizing gene expression level (39). The qRT-PCR conditions were one cycle of 50°C for 10 min, one cycle of 95°C for 1 min, then 40 cycles of 95°C for 10 s, and 60°C for 30 s, followed by melt curve analysis. The 2[-Delta Delta C(T)] method was applied to analyze the relative changes in gene expression from real-time quantitative PCR results (47). Gene expression data were analyzed by the Wilcoxon Method from JMP Pro (www.jmp.com), to comparing group means, calculate the p-value and determine the presence of significant differences between treatment groups.

To determine virus genome equivalents in blood and freshly fed mosquitoes from the infection study, DENV-2 specific primers (38) were used with the iTaq Universal SYBR Green One-step kit and the Bio-Rad CFX96™ Real-Time PCR (Bio-Rad, Hercules, CA, USA). The qRT-PCR conditions were as described above. Viral genome equivalents were calculated based on a standard curve and recorded as log (base 10) plaque-forming unit equivalents per milliliter (log PFUe/ml). The standard curves for DENV genome equivalents were described previously (38). The blood mixture contained 9.3 log PFUe/ml of DENV-2. The freshly fed mosquito contained 5.78 ± 0.07 log PFUe/ml of DENV-2. The virus titer at 60 hours after infection in the midgut was 6.25 ± 0.74 log PFUe/ml of DENV-2.

An average of 48.4 million reads was generated per sample, of which 90% were clean reads. The average sequencing error rate was 0.03 [Qphred=-10log10(e)] and the average GC content was 46.86%. The mapping percentage (total number of reads that could be mapped to the reference genome) was 87.7%, and more than 90% of mapped reads corresponded to known exon regions. An average of 9.25% reads were mapped to multiple locations and 78.42% reads were uniquely mapped. The full dataset of FKPM values and the FDR rate for each of 15,173 genes are provided as supplemental Table S1 .

We identified 31 genes that were significantly differentially expressed (DEGs) between the control and infected groups based on a false discovery rate of 0.2 ( Tables 1 – 3 ). These included 8 genes located on Chromosome 1, 15 genes on Chromosome 2, and 7 genes on Chromosome 3 ( Figure 1 ). Eleven genes were upregulated, and 20 genes were downregulated after 60 hours of virus infection in the midgut.

Based on data from Vectorbase, only 7 of the DEGs had an annotated function, with 24 described as having an ‘unspecified product’ based on AaegL5 annotation. Through BLASTx analysis we identified conserved domains in a further 13 DEGs, making a total of 20 annotated genes. Based on their predicted functions from Vectorbase and BLASTx analysis, 9 DEGs were grouped as having a role in immunity or oxidative stress response ( Table 1 ). A further 10 DEGs had a putative role relating to nucleic acid or protein process ( Table 2 ). One gene, AAEL023615-RA, was associated with both immunity and nucleic acid processes. The 13 DEGs not linked to these functions are recorded in Table 3 .

The 352 genes known to be involved in the dipteran immune response covered 23 families including Toll, IMD, JAK/STAT, and RNAi pathways. The list was split into 4 functional categories, defined as (i) pathogen/parasite recognition, (ii) signal modulation, (iii) signal transduction, and (iv) effector, oxidation-reduction, and reactive oxygen species. However, when the expression level differences between control and DENV-2 infected mosquitoes at 2.5 days post-infection in the data from this study were examined, 315 genes had FPKM data but none of the canonical immune genes showed expression level differences with a FDR below 0.2 ( Table S2 ).

Four genes were validated with qRT-PCR and supported the RNA-Seq results with the same pattern of expression ( Figure 2 ). The relative transcript expression fold change of AAEL014890 was 2.15 ± 0.96 (p = 0.02), AAEL028039 was 3.1 ± 1.68 (p = 0.02), AAEL008841 was 0.87 ± 0.15 (p = 0.38) and AAEL001417 was 1.04 ± 0.41 (p = 0.77). Although two transcripts did not have statistically significant differences, the relative expression patterns were the same as the RNA-Seq FPKM results.

Canonical immune pathways and genes are generally thought to play a crucial role in a mosquito’s ability to restrict viral infection, with many canonical immune transcripts displaying altered activity after infection (16, 30). Several immune pathways such as the Toll (16), JAK/STAT (19) and RNAi pathways (12) have a demonstrated ability to limit DENV infection in the mosquito. However, in our dataset there were no significant differences in the expression of 352 canonical immune genes ( Table S2 ), which suggests there might be a distinct, uncharacterized aspect of the immune response and host-virus interactions at this time point and in the midgut tissue we investigated.

Based on previous midgut transcriptome studies, canonical immune genes were generally not identified as significantly associated with the DENV replication cycle, although there were a few exceptions. The Inhibitor of apoptosis protein-1 (AAEL009074) and SCF ubiquitin ligase Rbx1 component (AAEL004691) were increased in the susceptible population after 48 hours of infection (48). The cecropin gene (AAEL017211), an antimicrobial peptide, was downregulated after 24 hours of infection (27). Several immunity-related genes including serine protease inhibitors (AAEL003182 and AAEL002704), Clip domain serine protease (AAEL002124 and AAEL001098), C-type lectin (AAEL005641), and fibrinogen-related protein genes (AAEL006704) were differentially expressed after 4 days of infection in the midgut (27). However, the canonical immune gene activity described in these studies was noted to occur either in the early infection stage or late replication phase (36) but not during the early replication phase. These observations in conjunction with our results suggest that distinct patterns of gene expression might underlie host-virus interactions during different phases of the DENV replication cycle. Moreover, gene expression likely varies among different Ae. aegypti populations if they have different genetic backgrounds. A previous study showed genome-wide divergence of Ae. aegypti populations (49).

Based on conserved domains, we identified 9 other DEGs with a putative role in immunity or redox-associated processes, which can be immunomodulatory in Ae. aegypti (50–52). Two DEGs, AAEL014372-RA and AAEL014890-RA, both linked to oxidative stress and immunity, were upregulated after infection ( Table 1 ). AAEL014372-RA is a juvenile hormone-inducible protein, that was thought to regulate the expression of many other genes in D. melanogaster (53). The conserved protein domain cl21453 from AAEL014372-RA belongs to the PKc_like superfamily which was mainly composed of the catalytic domains of serine/threonine-specific and tyrosine-specific protein kinases. The Ras/ERK signaling pathway is included in the PKc_like superfamily and was shown to play a role in resistance of Aedes mosquitoes to DENV infection (52). AAEL014890-RA is a known Cytochrome p450 gene associated with metabolic detoxification and insecticide resistance (54). Cytochrome p450 genes as monooxygenases were upregulated in the refractory Ae. aegypti in response to DENV infection (30). Although cytochrome p450s are associated with viral infection-induced inflammation and oxidative stress in mammals (51), the role in virus infection in mosquitoes is still unclear.

Most of the downregulated DEGs related to oxidative stress and immunity were associated with cell death and survival, which are key processes associated with virus infection (55, 56). Seven DEGs, AAEL023151-RA, AAEL008841-RA, AAEL012507-RA, AAEL017315-RA, AAEL001417-RA, AAEL023615-RA and AAEL027736-RA were downregulated in response to DENV infection ( Table 1 ). Based on BLASTx results, AAEL023151-RA is a CCAAT/enhancer-binding protein zeta, which has been implicated in growth arrest and the control of apoptosis in mammals (57). However, its role in mosquitoes still needs further investigation. AAEL008841-RA is acyl-coenzyme A oxidase which produces H₂O₂ and is involved in lipid metabolism (58). AAEL012507-RA has the conserved protein domain pfam01920, indicating it is likely a prefoldin. Prefoldin is a co-chaperone protein that transfers unfolded polypeptides to the chaperonin containing tailless complex polypeptide (59). On the other hand, AAEL017315-RA has the conserved protein domain cd10228 that suggests it is a 70 -kDa heat shock protein 4 (Hsp70) which is expressed under stress and is an important component of the chaperone system (60). Together, both AAEL012507-RA and AAEL017315-RA might be involved in chaperone activity where the chaperone function of Hsp70 protects the cell from stress-induced apoptosis (61). AAEL001417-RA is a leucine-rich immune protein (LRIM) which is activated as a complement-like defense response (62). Most of the studies of LRIMs have focused on their ability to act against malaria parasites in Anopheles mosquitoes (63). However, increased levels of LRIM16 transcripts were observed in Mayaro Virus-infected Aedes mosquitoes (64), although their exact role in arboviral infection remains to be determined. The presence of a conserved protein domain cl37801 in AAEL023615-RA could indicate a function similar to that of a topoisomerase II-associated protein (PAT1), which has conserved roles in the 5→3 mRNA decay pathway in both Drosophila and humans (65). Furthermore, PAT1 might affect cell survival and regulation of viral gene expression (66). AAEL027736-RA is predicted as cytochrome oxidase subunit I which is one of the proteins encoded from mitochondria DNA (67). The cl00275 domain AAEL027736-RA contained in the Heme-copper oxidase superfamily (68) which is involved in the reduction of oxygen.

Based on the potential functions of oxidative stress and immunity, we propose that DENV induces stress in midgut cells after entry. Altered mitochondrial function leads to changes in oxidative stress and oxidative homeostasis. Mass production of viral protein induces the unfolded protein response (69). Although the genes did not have a direct link to mosquito immunity, the higher levels of reactive oxygen species and unfolded protein could stimulate apoptosis (70, 71). Further experiments must be done to establish the connections between the genes we identified in the midgut and apoptosis and, moreover, their potential to moderate virus infection.

We identified 9 DEGs associated with nucleic acid/protein processes ( Table 2 ). Three DEGs, AAEL023671-RA, AAEL028039-RA, and AAEL009892-RA, were upregulated after 2.5 days post-infection and were associated with DNA binding and DNA replication. Cd09272 is the conserved domain of AAEL023671-RA and belongs to the family of Ribonuclease H (RNase H). RNase H cleaves the RNA of RNA/DNA hybrids and is critical during DNA replication (72). AAEL028039-RA maintains the conserved protein domain pfam00751 which is a DNA binding domain (73) and has been shown to dimerize and bind palindromic DNA (74). The AAEL009892-RA conserved protein domain pfam09789 identifies as an uncharacterized coiled-coil protein with homology to the basic leucine zipper domain, a large group of transcription factors (75).

There were a further 7 DEGs related to the nucleic acid/protein process that were downregulated ( Table 2 ). Most of these genes were associated with mRNA and rRNA processing. The mRNA splicing machinery was shown to be hijacked by the DENV NS5 protein to improve viral replication (76). On the other hand, DENV NS1 protein interacted with ribosomal proteins associated with DENV translation and replication (77). The possible involvement of DENV proteins might be a reason we identified the genes related to mRNA and rRNA processing. AAEL001191-RA has the conserved protein domain similar to that of transmembrane protein 62, which is a member of the metallophosphatase (MPP) superfamily. MPPs have diverse functions including Mre11/SbcD-like exonucleases which are critical to genome stability (78) and Dbr1-like RNA lariat debranching enzymes which are key factors of pre-mRNA splicing (79). AAEL017421-RA is orthologous to nucleolar protein 56. Although a specific function has not been identified (80), other nucleolar proteins are important for the maintenance of heterochromatin and ribosomal RNA (81). Moreover, AAEL017421-RA has a conserved domain, COG1498, that indicates similarity to RNA processing factor Prp31. The Nop domain in RNA processing factor Prp31 is required for ribosome biogenesis (82) and pre-mRNA splicing (83). The conserved protein domain COG1163 in AAEL012250-RA suggests ribosome interacting GTPase 1 which is important to ribosomal structure and biogenesis (84). AAEL007573 is the U3 small nucleolar RNA-associated protein 18 homolog which functions in the maturation of the small subunit rRNA (85). AAEL025585-RA contains a conserved protein domain pfam00096 and is the structure of the zinc finger that recognizes and binds DNA (86). AAEL007382-RA has the conserved protein domain cd00200 indicates WD40 which involved in many biological processes such as signal transduction, pre-mRNA processing, and cell cycle control (87, 88).

Since viruses utilize host transcription and translation machinery to replicate (76, 77), the genes we identified in Table 2 could be linked to viral RNA translation and replication. Furthermore, the virus might subvert the host cell cycle or nucleic acid processes to keep the cell in a state where viral replication is promoted (89). However, more studies are needed to understand the DENV replication process and to investigate interactions between the DEGs and virus.

We reported 13 DEGs were not listed in either immunity/oxidative stress response or to nucleic acid/protein processes ( Table 3 ). Twelve DEGs were identified as non-coding RNAs (ncRNA), specifically long non-coding RNAs (lncRNA) (90). Expression of cellular lncRNAs may be altered in response to viral replication and induced antiviral response (91). A previous study has shown that DENV-2 infection in Ae. aegypti increased the abundance of a number of long ncRNAs, which might be involved in the response to viral infection (92).

To determine if any of the 31 DEGs had a previously discovered change in expression due to arboviral infection, we compared our data with previous reports describing transcriptomic response to arboviral infection in Ae. aegypti (22, 93). Four DEGs were previously described as having transcript level differences under DENV or Zika virus infection (AAEL001417-RA, AAEL023122-RA, AAEL028030-RA and AAEL023254-RA). Moreover, if we expand the database to a virus-injected transcriptomic response study (34), a total of 15 of our DEGs were reported to have expression levels altered in response to infection with at least one of DENV, Zika virus, Yellow fever virus, or West Nile virus. However, 16 DEGs from our studies were not reported in previous studies. For the previously reported 15 DEGs, 9 were affected by DENV at 1 day post infection (dpi), 2dpi or 7dpi (AAEL014372-RA, AAEL008841-RA, AAEL001417-RA, AAEL023615-RA, AAEL009892-RA, AAEL023122-RA, AAEL023254-RA, AAEL025431-RA and AAEL007382-RA) and 6 DEGs (AAEL014890-RA, AAEL012507-RA, AAEL001191-RA, AAEL012250-RA, AAEL017315-RA and AAEL028030-RA) were only expressed differentially following infection with other viruses. Elucidation of potential genes involved in the replication cycle revealed altered expression of 5 DEGs at 1dpi or 2dpi after DENV infection (AAEL023615-RA, AAEL023122-RA, AAEL023254-RA, AAEL025431-RA, and AAEL007382-RA).

The AAEL023615-RA (PAT1) and AAEL007382-RA was described above and is crucial for cell survival (66, 88). AAEL025431-RA has the conserved domain Josephin, which might regulate cell motility, and endocytosis (94). The Josephin domain-containing protein is also involved in de-ubiquitination (95). AAEL023122-RA and AAEL023254-RA are potential lncRNAs that could participate in antiviral responses (92). The fact that the five DEGs exhibited differential expression across multiple studies of DENV or ZIKV infection in mosquitoes could indicate that they play important roles during arboviral infection, making them important candidates for future study.

A quantitative trait locus (QTL) is a chromosomal region responsible for variation of a quantitative trait and several QTLs that control Ae. aegypti susceptibility to virus infection and dissemination have been identified (96–98). Across the 31 DEGs we identified, there were 8 genes located on Chromosome 1, 15 genes on Chromosome 2, 7 genes on Chromosome 3, and 1 that was not assigned ( Figure 1 ). To date, six DENV-2 midgut infection barrier (MIB) QTLs and two DENV-2 midguts escape barrier (MEB) QTLs have been identified across all three chromosomes (36, 99). However, none of the 31 genes we discovered overlapped with the predicted QTL locations. On the other hand, Aubry et al. detected five highly significant QTLs and four QTLs with lower statistical support Chromosome 2 that were associated with Zika virus infection status in Ae. aegypti (100). Three of our DEGs, AAEL022994-RA, AAEL012507-RA, and AAEL008841-RA, fell in the highly significant QTLs regions. Additionally, three further DEGs, AAEL023696-RA, AAEL027382-RA, and AAEL007573-RA, were located within the less significant QTLs regions. These findings suggest that Zika virus-associated QTLs might also be involved in the response to DENV infection.

A greater number of differentially expressed genes may have been identified as significant if sample sizes had been larger. For example, if the same expression pattern held but we had 6 control and treatment groups, the number of genes below P-value of 0.05 after Sidak correction would have been 121. The number of genes greater than 95% power using false discovery rate 0.2 would have been 1,480. Nevertheless, the 31 DEGs we have identified here are those that demonstrated the most extensive response to DENV infection.

Our data revealed key time-specific differences in gene expression when compared to previous transcriptomic studies of DENV-infected midguts (27, 48). Given differences in gene expression between carcass and midgut samples after DENV infection (27), our findings highlight the value of time- and tissue-specific transcriptional studies to investigate viral infection in mosquitoes. Further research in this area will likely be crucial to improving understanding of the virus infection process. For example, DENV-2 dissemination to the abdominal fat body was detected in 35% of mosquitoes between 2 and 3 days post-infection (9). The fat body-mediated antiviral response of mosquitoes (101) and transcriptomes also changed after blood feeding (102). On the other hand, expanding the sampling time points for transcriptional response during the virus replication process could be applied with dual RNA-seq and provide a different view of pathogen- host interactions, moreover, the impact from the endogenous viral elements (103). All in all, our results provide a missing piece of mosquito’s early response against DENV and could offer potential targets for developing novel arbovirus control strategies.