Profiling Transcriptional Response of Dengue-2 Virus Infection in Midgut Tissue of Aedes aegypti

Understanding the mosquito antiviral response could reveal target pathways or genes of interest that could form the basis of new disease control applications. However, there is a paucity of data in the current literature in understanding antiviral response during the replication period. To illuminate the gene expression patterns in the replication stage, we collected gene expression data at 2.5 days after Dengue-2 virus (DENV-2) infection. We sequenced the whole transcriptome of the midgut tissue and compared gene expression levels between the control and virus-infected group. We identified 31 differentially expressed genes. Based on their function, we identified that those genes fell into two major functional categories - (1) nucleic acid/protein process and (2) immunity/oxidative stress response. Our study has identified candidate genes that can be followed up for gene overexpression/inhibition experiments to examine if the perturbed gene interaction may impact the mosquito’s immune response against DENV. This is an important step to understanding how mosquitoes eliminate the virus and provides an important foundation for further research in developing novel dengue control strategies.


INTRODUCTION
Dengue virus (DENV) is a mosquito-borne virus that has greatly impacted human health and is a global issue in tropical and sub-tropical areas (1). The dengue disease burden worldwide is estimated to range from 100 to 400 million infections each year (2), and its occurrence closely mirrors that of its primary vector, Aedes aegypti. According to the World Health Organization, the largest number of dengue cases reported in 2019 were 3.1 million cases in the American region and more than 1 million cases around southeast Asia. Dengue virus is described as having four serotypes, DENV-1, DENV-2, DENV-3, and DENV-4, all of which cause dengue fever. In some individuals, dengue fever can progress to severe dengue with life threatening complications (3,4).
The time required for a virus to replicate and become transmissible in the mosquito is called the extrinsic incubation period (EIP) (5). To become transmissible, the virus must pass through several major barriers in the mosquito: the midgut infection barrier (MIB), the midgut escape barrier (MEB), the salivary gland infection barrier (SGIB), and the salivary gland escape barrier (SGEB) (6). After a female Ae. aegypti ingests virus particles contained in a blood meal, the virus binds to receptors on the midgut-cell surface and enters the cell through endocytosis (7). The initiation of infection and replication in the midgut during the first two days is called the eclipse period, in which the virus titer decreases on day 1 and then continues to decline until day 2 post-ingestion (8). At 2 days post-infection, around 30% of midgut epithelial cells are infected (9). After 3 days post-infection, the virus can spread laterally to infect neighboring midgut epithelial cells (9). Infected midgut cells often become degenerated and cell-cell adhesion is often disrupted which may be the path for viral escape from the midgut epithelium (6). The duration of the EIP varies depending on the genotypes of both the mosquito and the virus, as well as on environmental conditions (10).
These processes rely on complex chains of virus-host interactions, many of which are still poorly characterized. The virus must evade the host immune system and co-opt host machinery to complete its replication cycle. There are a wide range of different host immune processes that can prevent the virus from using resources or disturbing critical processes essential for cell function. This includes passive processes and processes triggered by recognition of pathogen-associated molecules, which have been shown to limit pathogen infection and transmission (11). The RNA interference (RNAi) pathway, which limits viral replication by targeting viral RNAs for degradation has been found to play a major role in the antiviral response (12)(13)(14). Three other innate immunity pathways involving pattern recognition and signal transduction have an identified antiviral role: the Toll (15,16), IMD (17,18), and JAK/STAT (19,20) pathways. Activation of these pathways leads to the induction of antimicrobial peptides and other immune effectors. A further pathway, apoptosis, the process of programmed cell death, is also crucial to limiting viral titer (21). Infection with DENV, Zika, and chikungunya viruses often induces the expression of other immune genes, including cysteine-rich venom proteins, clip-domain serine proteases, C-type lectins, and serine proteases, with the effect varying between mosquito populations (22)(23)(24)(25).
The early phase of the viral replication cycle is critical, as it shapes the course of infection and ultimately determines whether the mosquito becomes competent to transmit the virus (36). Most transcriptomic assays in mosquitoes that examined arbovirus infection were conducted within 48 hours of infection (27,30,32,34) or 7-10 days post-DENV infection (16,17,19,37). However, the time between 2 and 3 days after virus infection may be a crucial time point corresponding to the end of an infection bottleneck before viral titers spike and the virus particles escape the midgut epithelium (8,9,36). Mosquito host-virus interactions at this stage are still not well defined, so examining transcript differences associated with infection at this key time point could improve our overall understanding of DENV infection in mosquitoes, and produce candidate genes that might be used to develop for new strategies aimed at preventing arbovirus transmission.
Here we investigate transcriptomic changes in mosquito midguts after 60 hours of DENV infection, which corresponds to the end of the eclipse phase when virus replication rates and DENV titers begin to increase (36). Through RNA-Seq analysis, we evaluate the impact of DENV infection in the Ae. aegypti Vero Beach strain to identify differentially expressed genes associated with the virus replication cycle. A total of 31 genes had transcripts that showed differential expression following DENV infection. Our results provide important insight into DENV replication and indicate differential regulation of genes in two major functional categories -(1) nucleic acid/protein process and (2) immunity/oxidative stress response, during this stage of infection.

Mosquitoes
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).

Virus Infection Experiments
Four-to fiveday 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.

Transcriptome Sequencing
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 pairedend 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).

Differential Gene Expression Analysis
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.

Impact on Canonical Immune Genes
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.

Validation of Gene Expression Differences and Viral Titer Measurement
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) plaqueforming 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.

Transcriptome
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.

Differentially Expressed Genes
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.

Effect on Canonical Immune Genes
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).

Validation
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.

DISCUSSION
We investigated the transcriptional response in midgut tissue of Ae. aegypti at 2.5 days post-infection with DENV-2. Although DENV is detectable in mosquito midgut epithelial cells at 2 days post-infection, the subsequent 24 hours are crucial for the virus to replicate and infect surrounding cells. This makes 2.5 days post-infection a critical timepoint for studying viral replication. We identified 31 transcripts with significantly altered expression levels following DENV infection, and based on their putative functions, grouped the majority as being linked either to immunity/oxidative stress response or to nucleic acid/protein processes. This finding provides new information on the nature of mosquito-virus interactions during the replication and assembly stages of the DENV replication cycle, and offers 7 new target genes with their potential function that could interfere with the virus infection process in a mosquito host.

Canonical Immune Genes Expression
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 Cell maintenance and division, potential role in apoptosis Chr, chromosome; Ctrl, control group; DVI, DENV infected group; Ctrl/DVI columns represented the average of FPKM value ± standard deviation; D, change in expression level in DENV infected group vs control group, ↑ indicated the expression level was significantly higher in the DENV infected group and ↓ showed the expression level was significantly higher in the control group; FDR, false discovery rate; cd_id, conserved domain id on NCBI conserved domain database (CDD); *indicated the transcript was identified in previous arbovirus studies.
(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)

The Immunity/Oxidative Stress Response-Related Genes
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)(51)(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/threoninespecific 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 2 O 2 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.

Nucleic Acid/Protein Processes-Related Genes
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.

Other Genes and the Potential Genes Had Conserved Response Associated With DENV Replication Stage
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) 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.

Gene Location and Potential Connection to QTLs
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)(97)(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.

Study Caveats, Future Directions and Conclusions
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 postinfection (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.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here in NCBI under accession number PRJNA729510.

ETHICS STATEMENT
The animal study was reviewed and approved by Institutional Animal Care and Use Committee University of Florida.