Characterization of Viral miRNAs during Adenovirus 14 Infection and Their Differential Expression in the Emergent Strain Adenovirus 14p1

Human adenoviruses (HAdV) express either one or two virus-associated RNAs (VA RNAI or VA RNAII). The structure of VA RNA resembles human precursor microRNAs (pre-miRNA), and, like human pre-miRNA, VA RNA can be processed by DICER into small RNAs that resemble human miRNA. VA RNA-derived miRNA (mivaRNA) can mimic human miRNA post-transcriptional gene repression by binding to complementary sequences in the 3′ UTR of host mRNA. HAdV14 is a member of the B2 subspecies of species B adenovirus, and the emergent strain HAdV14p1 is associated with severe respiratory illness that can lead to acute respiratory distress syndrome. Utilizing small RNA sequencing, we identified four main mivaRNAs generated from the HAdV14/p1 VA RNA gene, two from each of the 5′ and 3′ regions of the terminal stem. There were temporal expression changes in the abundance of 5′ and 3′ mivaRNAs, with 3′ mivaRNAs more highly expressed early in infection and 5′ mivaRNAs more highly expressed later in infection. In addition, there are differences in expression between the emergent and reference strains, with HAdV14 expressing more mivaRNAs early during infection and HAdV14p1 having higher expression later during infection. HAdV14/p1 mivaRNAs were also shown to repress gene expression in a luciferase gene reporter system. Our results raise the question as to whether differential expression of mivaRNAs during HAdV14p1 infection could play a role in the increased pathogenesis associated with the emergent strain.


Introduction
Human adenoviruses (HAdVs) are non-enveloped viruses with a linear doublestranded DNA genome of~36,000 bp. HAdVs are grouped into seven species (A-G) and consist of over 100 different types (consecutively numbered) based on genomic sequencing. HAdVs can infect the respiratory, gastrointestinal, and urinary tracts as well as the conjunctiva. Overall, HAdV usually results in mild, self-limited infections in immunocompetent individuals. HAdV14 is a member of the B2 subgroup of adenovirus and, as such, predominately causes kidney and urinary tract infections but can also be associated with respiratory infections. Outbreaks of emergent strains of Ad have resulted in severe and sometimes fatal infections in otherwise healthy people. In the past fifteen years, an emergent strain of HAdV14, HAdV14p1, emerged first in the U.S. and subsequently throughout the world [1][2][3][4][5][6]. HAdV14p1 was first identified as the causative agent of an outbreak of acute respiratory distress syndrome (ARDS) in U.S. military populations [7]. HAdV14p1 is 99.7% identical to HAdV14 yet displays increased lung pathogenesis in the Syrian hamster that models HAdV14p1 pathogenesis in humans [8,9]. Why HAdV14p1 has increased pathogenesis is still unclear.
The HAdV genome encodes roughly 40 proteins. In addition to protein encoding RNAs, non-coding RNAs (ncRNA) are also present in HAdV-infected cells [10][11][12]. Viralassociated RNA (VA RNA) are encoded by the HAdV genome and are the most abundant viral ncRNA found in infected cells [13,14]. All HAdV encode at least one VA RNA transcript~160 nt in length, and approximately 80% of HAdV serotypes encode at least 2 VA RNAs (VA RNA-I and VA RNA-II), with VA RNA-I encoded by all HAdV serotypes [15,16]. VA RNAs are required for efficient viral replication, as the deletion of VA RNA-I reduces viral titer by 20-fold and the deletion of both VA RNA-I and VA RNA-II results in a 60-fold decrease in viral titer [17,18]. Unlike protein encoding HAdV RNAs, VA RNAs are transcribed by cellular RNA polymerase III and are expressed throughout the viral replication cycle [13].
VA RNA forms a complex secondary structure that is critical for its biological function [19][20][21]. Overall, despite their low sequence conservation among HAdV, VA RNAs share similar secondary structure and organization of the apical stem, central domain, and terminal stem [16,21]. The crystal structure of HAdV-2 VA RNA-I reveals a sharply bent, coaxially stacked, wobble-enriched, and pseudoknot-anchored molecule to allow VA RNA to interfere with nearly all host systems that interact with double-stranded RNA (dsRNA) [19]. The best characterized activity of VA RNA is its inhibition of protein kinase R (PKR) through the ability of VA RNA-I to bind and inhibit PKR dimerization [22][23][24]. PKR is part of the cellular innate immune response to viral infection, and, by sensing dsRNA produced during viral infection, PKR activation results in the global inhibition of cap-dependent protein synthesis. VA RNA also binds to 2 -5 Oligoadenylate Synthetase 1 (OAS1), which results in the activation of OAS1 [25]. While the interaction of full-length VA RNA activates OAS1, the binding of cleaved VA RNAI (produced following DICER cleavage) to OAS1 inhibits OAS1 activity [11,26,27]. Interestingly, the cleaved VA RNAI has a higher affinity for OAS1 than full-length VA RNAI, suggesting that the inhibitory activity of cleaved VA RNAI may be the predominant activity of VA RNAI on OAS1 [28]. VA RNAI binds to the retinoic acid-inducible gene I (RIG-I). The triphosphorylyated 5 -end nucleotide of VA RNA triggers RIG-I signaling, which results in an increased type I interferon response [29][30][31][32]. These roles in modulating the cellular innate immune response to HAdV have been long understood.
Recently, new roles for VA RNAI in regulating cellular miRNA processes have been elucidated. Following transcription, VA RNAI is exported from the nucleus to the cytoplasm by Exportin 5 (Exp5). Exp5 exports RNA molecules (e.g., pre-miRNA and tRNAs) that contain a 'mini-helix' from the nucleus [33][34][35][36]. VA RNAI contains a 'mini-helix' at the end of the terminal stem, which results in VA RNA being a direct competitor with cellular small RNAs for interaction with Exp5 [37]. In addition to competing with pre-miRNA for nuclear export through Exp5, once in the cytoplasm, VA RNAs are cleaved by DICER, like pre-miRNA, and as VA RNA expression increases, more VA RNAs are associated with DICER than host pre-miRNA, resulting in decreased host miRNA expression [27,38,39]. DICER processes VA RNA into small RNAs that resemble host miRNA (called mivaRNA) that can be incorporated into the RNA-induced silencing complex (RISC) [37,38,[40][41][42]. RISC association allows for mivaRNA to bind and repress target mRNA expression [38,[40][41][42][43].
Group B2 HAdV members encode only the VA RNAI gene. Bioinformatic analysis of the HAdV14 genome showed that the VA RNAI gene of Ad14 is~98% identical to the other group B2 members HAdV11, HAdV34, and HAdV35 [44]. In this report, we further characterize the VA RNA gene, identifying the intragenic promoter elements, transcriptional stop sites, potential OAS1 consensus sites, and GGGU/ACCC sites. Small RNA-seq was used to identify the mivaRNA produced during infection. Despite 100% sequence identity of the VA RNAs, differential expression of the mivaRNA exists between HAdV14-and HAdV14p1-infected A549 cells. We discuss the possibility that this differential mivaRNA expression could play a role in the pathogenesis of HAdV14p1.

Infection of A549 Cells and Isolation of Total Cellular RNA
A549 cells were infected with either HAdV14 or HAdV14p1 at a multiplicity of infection (MOI) of 10 pfu/cell in suspension for 1 h at 37 • C, after which cells were plated and allowed to adhere until collected. Adherent and non-adherent cells were collected at 6, 12, 24, 36, and 48 h post-infection. Total RNA was isolated using miRNeasy kit (Qiagen, Germantown, MD, USA) with on column DNase treatment. The total RNA in each sample was quantified using the Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA), and quality was measured using the RNA6000 Nano chip on the Agilent 2100 Bioanalzyer (Agilent Technologies, Santa Clara, CA, USA). All samples had an RNA integrity number greater than 7.

Small RNA Library Preparation, Sequencing and Data Analysis
TruSeq Small RNA library prep kit (Illumina, San Diego, CA, USA) was used to create sequencing libraries. Specifically, adapters were ligated using the 5 phosphate and 3 hydroxyl groups common to most mature miRNAs. After adapter ligation, samples were reverse transcribed and amplified. Finally, the libraries were size selected using a 6% polyacrylamide gel and concentrated using ethanol precipitation. Purified libraries were normalized and pooled to create a double-stranded cDNA library and were sequenced on the Illumina MiSeq to render 50 base pair single end reads at the Loyola Stritch School of Medicine Genomics Facility. Adapter sequences were removed, and low-quality reads were trimmed from raw sequencing reads using Cutadapt (v. 1.11). The resulting reads were mapped to the HAdV14 deWit genome (GenBank accession number AY803294) in CLC Genomics Workbench (Qiagen).

Differential Expression Analysis
Differential expression was conducted with CLC Genomics Workbench. The VA RNA sequence was added to miRbase 22 and annotated to contain either total 5 and 3 reads or mivaRNA 5 A, 5 G, 3 A and 3 C seeds. Quantify miRNA 1.2 (CLC) was used with the following settings: allow for length-based isomiRs, no additional upstream bases, 5 additional downstream bases, no missing or mismatch bases, minimum sequence length = 18, maximum sequence length 25. Differential Expression for RNA-Seq 2.6 (CLC) was used for differential expression analysis with normalization set to Trimmed Mean of Means (TMM).

HAdV-14 mivaRNA 3 UTR Luciferase Reporter Assays
The complement of the either the 5 or 3 seed sequences of Ad14 VA RNA were cloned into the multiple cloning site of pmirGLO Dual Luciferase Vector (Promega, Madison, WI, USA). Constructs were sequenced to confirm the correct orientation of the complimentary HAdV14 VA RNA seed sequences. Plasmids were purified using PureLink HiPure Plasmid Maxi Prep (Invitrogen) and transfected into A549 cells using Lipofectamine LTX (Invitrogen). Twenty-four hours after transfection, cells were infected with Ad14 at an MOI of 10. After an additional 24 h, luciferase activity was determined using a dual-luciferase assay (Promega). Firefly luciferase activity was normalized to Renilla luciferase activity. One-way ANOVA was followed by the Holm-Sidák test with p < 0.05 considered significant (Prism 9; GraphPad Software, San Diego, CA, USA).

Prediction of mivaRNA Host Gene Targets and Functional Enrichment Analysis
TargetScan and miRDB were used to predict cellular targets for mivaRNAs by seed sequences and filtered to remove duplicate targets. KEGG (https://www.genome.jp/kegg/, accessed on 2 May 2021) was used to predict the cellular pathways and processes potentially regulated by mivaRNA-target mRNA interaction.
2.7. RNA-Seq Library Preparation, Sequencing, Differential Expression Analysis, and Bioinformatic Analysis Illumina stranded mRNA libraries were made from the total RNA preps from A549infected cells. Sequencing was performed at University of Oregon Genomics & Cell Characterization Core (Eugene, OR, USA) on an Illumina HiSeq4000 set for pair ends 100 bp reads and demultiplexed. Reads were trimmed with Trim Reads 2.4 (CLC) before RNA-seq Analysis 2.21 (CLC) was performed, mapping to the human genome (Hg38) with a mismatch cost of 2, insertion cost of 3, deletion cost of 3, and length/similarity fractions set to 0.8. Paired reads were counted as a single read and expression value was set to reads per kilobase million (RPKM). Differential expression for RNA-Seq 2.4 (CLC) using gene expression values from RNA-seq analysis was performed and the results were filtered to the predicted mivaRNA targeted genes. Ingenuity Pathway Analysis (IPA) (Qiagen) was used for functional enrichment with the differential RNA-seq analysis from HAdV14-infected cells at 36 hpi vs. uninfected cells. FDR was set to 0.05 with absolute expression fold changes being >1.5 and expression intensity being >1.

Analysis of HAdV14/HAdV14p1 VA RNA Gene Sequence and Secondary Structure
Like other HAdV group B2 members, HAdV14 and HAdV14p1 genomes encode only 1 VA RNA sequence that resembles VA RNA-I at nucleotides 10,452-10,613 [44]. Houng and colleagues have shown that the sequences of the VA RNA genes of HAdV14 and HAdV14p1 are 100% identical [9]. Analysis of the VA RNA-I sequence reveals two intragenic promoter elements (Box A (+13-23) and Box B (+55-65)) that are required for RNA polymerase III transcription ( Figure 1A). VA RNA-I transcription in many serotypes has been shown to be initiated at one of two potential transcription start sites, either a G(+1) or three nucleotides upstream at A(−3) [16,27]. The HAdV14 VA RNA-I gene contains both possible start sites. Termination sequences (T1A and T1B) are found at +152-161. Additionally, there is a backup termination sequence that is found at +192-198. The VA RNA gene contains the conversed GGGU (+33-36) and ACCC (+120-123) that form the central domain ( Figure 1A,B). Figure 1B shows a predicted secondary structure for the HAdV14 VA RNA molecule consisting of an apical stem, central domain, and terminal stem. HAdV VA RNA I contains two OAS1 consensus sites (WW(N 9 )WG) found in the central domain [45,46]. HAdV14 contains two OAS1 consensus sites at +36-48 and +93-105; both are found in the central domain.

Identification and Temporal Expression of HAdV14 Small RNAs
Previous studies have shown that the majority of small RNAs produced during HAdV infection map to the VA RNA gene(s). To examine the viral small RNAs produced during HAdV14 or HAdV14p1 infection, A549 cells were infected with either HAdV14 or HAdV14p1 ( Figure 2). Four independent small RNA libraries from infection with either HAdV14 or HAdV14p1 per time point were sequenced. At 6 hpi, less than 1% of the reads mapped to the HAdV14 genome, then increased to ~20% at 24 hpi, then remained at ~15% at 36 and 48 hpi (Table 1). At all time points, >92% of all small RNAs that mapped to the HAdV14 genome mapped to VA RNA gene ( Table 1), showing that the VA RNA gene is responsible for the vast majority of all HAdV14-encoded small RNAs. The diversity of the reads that aligned to the VA RNA gene increased during infection mainly by extensions or deletions at either the 5′ or 3′ ends of the small VA RNAs, with the total number of unique small VA RNAs similar between HAdV14 and HAdV14p1 (Supplementary Materials Table S1). The majority of the reads that aligned to the VA RNA gene mapped to the 5′ and 3′ ends of the terminal stem, with each small VA RNA comprising more than 1% of the total small RNA pool sequenced ( Table 2). The number of unique small RNAs that mapped to the VA RNA gene increased from 6 to 24 hpi and then stabilized through 48 hpi ( Figure 3A). Examination of the read counts showed temporal expression of smal RNAs produced from the 5′ or 3′ ends of the VA RNA gene. At 6 and 12 hpi, the 3′ smal RNAs are predominant, while the 5′ small RNAs show increased expression 36 and 48 hp ( Figure 3B-F). The same temporal expression patterns are seen in HAdV14p1-infected cells ( Figure 3B-F). The differences in read counts between HAdV14-and HAdV14p1 infected cells suggest that there may be differential expression of 5′ and 3′ small VA RNAs produced during infection. Differential expression analysis showed that, for the most part 5′ small VA RNAs during HAdV14p1 infection were not differentially expressed compared to HAdV14 infection, except at 12 hpi. In contrast, 3′ small VA RNAs were differentially expressed during HAdV14p1 infection compared to HAdV14 infection except at 24 hpi.

Identification and Temporal Expression of HAdV14 Small RNAs
Previous studies have shown that the majority of small RNAs produced during HAdV infection map to the VA RNA gene(s). To examine the viral small RNAs produced during HAdV14 or HAdV14p1 infection, A549 cells were infected with either HAdV14 or HAdV14p1 ( Figure 2). Four independent small RNA libraries from infection with either HAdV14 or HAdV14p1 per time point were sequenced. At 6 hpi, less than 1% of the reads mapped to the HAdV14 genome, then increased to~20% at 24 hpi, then remained at~15% at 36 and 48 hpi (Table 1). At all time points, >92% of all small RNAs that mapped to the HAdV14 genome mapped to VA RNA gene ( Table 1), showing that the VA RNA gene is responsible for the vast majority of all HAdV14-encoded small RNAs. The diversity of the reads that aligned to the VA RNA gene increased during infection mainly by extensions or deletions at either the 5 or 3 ends of the small VA RNAs, with the total number of unique small VA RNAs similar between HAdV14 and HAdV14p1 (Supplementary Materials Table S1). The majority of the reads that aligned to the VA RNA gene mapped to the 5 and 3 ends of the terminal stem, with each small VA RNA comprising more than 1% of the total small RNA pool sequenced ( Table 2). The number of unique small RNAs that mapped to the VA RNA gene increased from 6 to 24 hpi and then stabilized through 48 hpi ( Figure 3A). Examination of the read counts showed temporal expression of small RNAs produced from the 5 or 3 ends of the VA RNA gene. At 6 and 12 hpi, the 3 small RNAs are predominant, while the 5 small RNAs show increased expression 36 and 48 hpi ( Figure 3B-F). The same temporal expression patterns are seen in HAdV14p1-infected cells ( Figure 3B-F). The differences in read counts between HAdV14-and HAdV14p1-infected cells suggest that there may be differential expression of 5 and 3 small VA RNAs produced during infection. Differential expression analysis showed that, for the most part, 5 small VA RNAs during HAdV14p1 infection were not differentially expressed compared to HAdV14 infection, except at 12 hpi. In contrast, 3 small VA RNAs were differentially expressed during HAdV14p1 infection compared to HAdV14 infection, except at 24 hpi. MiSeq, and reads were mapped to the HAdV14 genome to identify HAdV14 mivaRNAs. Quantitation of reads was performed with CLC Genomics Workbench. TargetScan and miRDB were used to predict human mRNA targets based on seed sequences of the identified mivaRNAs. To examine gene expression in A549 cells during infection, Illumina stranded mRNA-seq libraries were constructed, sequenced on a NextSeq 2000, reads mapped to the human genome, and differential gene expression determined. Ingenuity Pathway Analysis was used to explore differential expression of mivaRNA predicted target genes. Created with BioRender.com.

Identification of Potential HAdV14/14p1 mivaRNAs
While there is diversity in the small RNAs that map to the VA RNA gene, repression of target gene expression involves the binding of miRNA to the target mRNA through the miRNA seed sequence, nucleotides 2-8 of the miRNA, and its complement on the target gene mRNA. Therefore, the seed sequences provide a better metric for both identifying and quantitating HAdV mivaRNAs. Unique miRNA seed sequences were identified for each small VA RNA molecule ( Table 3). The total number of unique seed sequences increased during HAdV14 or HAdV14p1 infection and plateaued at 24 hpi ( Figure 4A). The number of unique seed sequences was similar in both HAdV14-and HAdV14p1infected cells ( Figure 4A). To identify potential HAdV14/14p1 mivaRNAs, we identified the seed sequences that had total read counts >1000 and were found in small RNAs 20-33 nucleotides long. Eleven seed sequences were identified that met those criteria, and all of them mapped to either the 5 or 3 terminal ends of the VA RNA ( Figure 4B). Based on the read counts, there are four main mivaRNAs produced during HAdV14/14p1 infection, two from both the 5 and 3 ends. The mivaRNAs from the 5 end start at either the A (−3) or G (+1) and have been named mivaRNA 5 A or mivaRNA 5 G. At the 3 end, the predominant mivaRNAs start at either C (138) or A (139) and have been designated mivaRNA 3 C and 3 A, respectively ( Figure 4B). All four of the predominant mivaRNAs average 20 to 24 nucleotides in length, with the 5 mivaRNAs showing the greatest length diversity, out to 50 nucleotides long (Supplementary Materials Table S1). Overall, >95% of the reads are between 20 and 24 nt. Using our read counts ( Table 2) and unique seed sequence data (Table 3), we predict that there are five dominant DICER cleavage sites in HAdV-14/14p1 VA RNA. These predicted cleavage sites generate greater than 98% of the unique seed sequences that map to the VA RNA gene ( Figure 4C). We predict at least four secondary cleavage sites that generate the unique seed sequences in Table 3 that have more than 1000 total counts ( Figure 4C).

Differential Expression of mivaRNA during HAdV14 and HAdV14p1 Infection
Infection with HAdV14 or HAdV14p1 resulted in differential expression of 3′ small VA RNAs (Figure 3). With the dominant mivaRNAs identified, we performed differential expression analysis of each mivaRNA over time in HAdV14p1-infected vs. HAdV14infected cells. There was no differential expression of mivaRNA 5′A during infection

Differential Expression of mivaRNA during HAdV14 and HAdV14p1 Infection
Infection with HAdV14 or HAdV14p1 resulted in differential expression of 3 small VA RNAs (Figure 3). With the dominant mivaRNAs identified, we performed differential expression analysis of each mivaRNA over time in HAdV14p1-infected vs. HAdV14infected cells. There was no differential expression of mivaRNA 5 A during infection ( Figure 5A). However, mivaRNA 5 G ( Figure 5B) showed differential expression that was lower at 6 and 12 hpi in HAdV14p1-infected cells and then was higher in HAdV14p1infected cells at 48 hpi. In HAdV14p1-infected cells, both of the 3 mivaRNAs (3 C and 3 A) showed a lower expression at 6 and 12 hpi compared to HAdV-14-infected cells ( Figure 5C,D). At later times (36 and 48 hpi), mivaRNA 3 C and 3 A were both expressed at higher levels during HAdV14p1 infection than during HAdV14 infection. Overall, the data show that there are temporal mivaRNA expression differences between HAdV14-and HAdV14p1-infected cells.
infected cells at 48 hpi. In HAdV14p1-infected cells, both of the 3′ mivaRNAs (3′C and 3′A) showed a lower expression at 6 and 12 hpi compared to HAdV-14-infected cells ( Figure  5C,D). At later times (36 and 48 hpi), mivaRNA 3′C and 3′A were both expressed at higher levels during HAdV14p1 infection than during HAdV14 infection. Overall, the data show that there are temporal mivaRNA expression differences between HAdV14-and HAdV14p1-infected cells.

Expression of HAdV14 mivaRNA and A549 miRNA Expression during Infection
DICER processes pre-miRNAs and VA RNA into mature miRNA and mivaRNA respectively [27,39]. During HAdV infection, it has been shown that VA RNAs can suppress miRNA biogenesis by competing with host pre-miRNAs for cleavage by DICER [38]. Infection with HAdV14p1 results in nearly a 50% reduction in cellular miRNA at 24 hpi that persists through full CPE (Table 4). To explore the relative abundance of HAdV14 mivaRNAs in relation to A549 miRNAs, the identified HAdV-14 mivaRNAs were added to the miRBase v22 list of human miRNAs to allow for counting of mivRNAs and host miRNAs in the CLC Genomics Workbench. At 6 hpi, only mivaRNA 3′C was one of the top 20 total miRNAs expressed ( Table 5). Beginning at 12 hpi, mivaRNA 3′C, 5′A, and 5′G were found in the top six total miRNAs expressed and remained in the top six through full CPE at 48 hpi. After 6 hpi, all four predominant mivaRNAs were expressed at similar levels to the top cellular miRNAs. In general, increased expression of mivaRNAs resulted in decreased expression of all of the top 20 cellular miRNAs.

Expression of HAdV14 mivaRNA and A549 miRNA Expression during Infection
DICER processes pre-miRNAs and VA RNA into mature miRNA and mivaRNA, respectively [27,39]. During HAdV infection, it has been shown that VA RNAs can suppress miRNA biogenesis by competing with host pre-miRNAs for cleavage by DICER [38]. Infection with HAdV14p1 results in nearly a 50% reduction in cellular miRNA at 24 hpi that persists through full CPE (Table 4). To explore the relative abundance of HAdV14 mivaRNAs in relation to A549 miRNAs, the identified HAdV-14 mivaRNAs were added to the miRBase v22 list of human miRNAs to allow for counting of mivRNAs and host miRNAs in the CLC Genomics Workbench. At 6 hpi, only mivaRNA 3 C was one of the top 20 total miRNAs expressed ( Table 5). Beginning at 12 hpi, mivaRNA 3 C, 5 A, and 5 G were found in the top six total miRNAs expressed and remained in the top six through full CPE at 48 hpi. After 6 hpi, all four predominant mivaRNAs were expressed at similar levels to the top cellular miRNAs. In general, increased expression of mivaRNAs resulted in decreased expression of all of the top 20 cellular miRNAs.

HAdV14 mivaRNA Are Functional miRNAs
Since HAdV VA RNAs are processed by DICER into mivaRNAs, we hypothesized that the mivaRNAs can repress gene expression like cellular miRNAs. To test the ability of HAdV14 mivaRNA to function as miRNA, the complement of each of the four main mivaRNA seed sequences were cloned into the 3 UTR of luciferase in the pmirGLO Dual Luciferase Vector. Each reporter vector was transfected into A549 cells 24 h prior to infection with HAdV-14 at an MOI of 10, and luciferase activity was determined 24 hpi. As shown in Figure 6, HAdV-14 infection resulted in the repression of luciferase activity from all four mivaRNAs at 24 hpi. The strongest repression was seen from the 3 mivaRNAs that are expressed at higher levels at 6 and 12 hpi ( Figure 3B). uses 2022, 14, x FOR PEER REVIEW Figure 6. Regulation of luciferase activity by mivaRNAs during viral infection sequence of the mivaRNA seeds were cloned into pmirGLO Dual Luciferase V into A549 cells. Twenty-four hours after transfection, A549 cells were infected of 10 and luciferase activity was determined 24 h after infection. Control see with no mivaRNA complementary sequence. Luciferase activity is expr Firefly:Renillia luciferase activity. Mean ± SD, n = 8, one-way ANOVA followe ** p < 0.01, *** p < 0.001.

Prediction of mivaRNA Host Targets and Functional Enrichment Analy
Having established that mivaRNAs are expressed at levels miRNAs (Table 5) and are capable of repressing gene expression (Figu that have the complementary seed sequence in the 3′ UTR of their m determine if there was any potential for HAdV14 mivaRNAs to rep TargetScan and mirDB were used to predict potential host target mR sequences identified ( Figure 4B). All but one of the mivaRNAs are uni that are not found in human miRNAs. Th CAAAAATCCAGGATACGGAATCGAGTCGTT encodes the same human miR-584c-3p and was excluded from further analysis. The rem seeds were used to predict potential gene targets. Overall, there were 2 that could be targeted by the Ad14 mivaRNAs (Supplementary Materi 2184, 1919 were recognized by KEGG, and KEGG pathway analysis wa predicted gene targets. The top pathways included many signal tran Figure 6. Regulation of luciferase activity by mivaRNAs during viral infection. The complementary sequence of the mivaRNA seeds were cloned into pmirGLO Dual Luciferase Vector and transfected into A549 cells. Twenty-four hours after transfection, A549 cells were infected with Ad14 at an MOI of 10 and luciferase activity was determined 24 h after infection. Control seed is pmirGLO Vector with no mivaRNA complementary sequence. Luciferase activity is expressed as normalized Firefly:Renillia luciferase activity. Mean ± SD, n = 8, one-way ANOVA followed by Holm-Sidák test ** p < 0.01, *** p < 0.001.

Prediction of mivaRNA Host Targets and Functional Enrichment Analysis
Having established that mivaRNAs are expressed at levels similar to cellular miRNAs (Table 5) and are capable of repressing gene expression ( Figure 6) of target genes that have the complementary seed sequence in the 3 UTR of their mRNA, we sought to determine if there was any potential for HAdV14 mivaRNAs to repress cellular genes. TargetScan and mirDB were used to predict potential host target mR-NAs using the seed sequences identified ( Figure 4B). All but one of the mivaRNAs are unique seed sequences that are not found in human miRNAs. The 3 -mivaRNA-CAAAAATCCAGGATACGGAATCGAGTCGTT encodes the same seed sequence as human miR-584c-3p and was excluded from further analysis. The remaining 10 mi-vaRNA seeds were used to predict potential gene targets. Overall, there were 2184 predicted genes that could be targeted by the Ad14 mivaRNAs (Supplementary Materials  Table S2). Of the 2184, 1919 were recognized by KEGG, and KEGG pathway analysis was performed on the predicted gene targets. The top pathways included many signal transduction pathways and cellular processes involved in viral infection, including metabolism, cytoskeleton/tight junctions/focal adhesion/cell adhesion molecules, mRNA transport/splicing/regulation, and cell death processes (Tables 6 and S3). Table 6. Enrichment of mivaRNA targets by KEGG analysis.

Analysis of Predicted Target Gene Expression by RNA-Seq
Since our in vitro studies showed that the mivaRNAs can functionally repress luciferase gene expression ( Figure 6), we sought to determine if mivaRNA might regulate expression of predicted target genes. RNA-seq analysis was performed on the total RNAs from cells at all time points. Differential expression analysis was performed against uninfected control RNA libraries. The genes were filtered to the 1919 genes recognized by KEGG analysis, and genes that were down regulated by at least 1.5-fold with an FDR of <0.05 were considered significantly down regulated. As seen in Figure 7A, from 6 hpi to 24 hpi, down regulated target genes increased from~200 to~370 for both HAdV14-and HAdV14p1-infected cells. Between 24 and 36 hpi, the number of down regulated genes plateaued (HAdV14) or slightly decreased (HAdV14p1) from the 24 hpi number. At 48 hpi (full CPE), the number of down regulated genes doubled for both HAdV14 and HAdV14p1 from the 36 hpi time point. At both 36 and 48 hpi, there were more down regulated target genes in HAdV14-infected cells than in HAdV14p1-infected cells. From 6 hpi to 24 hpi, ≥78% of down regulated targets were shared between HAdV14-and HAdV14p1-infected cells ( Figure 7B). This similarity to target repression decreased at 36 and 48 hpi. To gain a better understanding of how the down regulated gene expression might alter cellular pathways, we performed IPA with the 36 hpi differential gene expression data sets. This time point was selected, as there should be sufficient repression observed at 36 hpi, and the time course of viral infection has not resulted in 4+ CPE. As seen in Table 7, the top ten pathways enriched from our RNA-seq data are very similar between HAdV14-and HAdV14p1-infected cells. The two pathways that are not the same in the top ten are not far out of the top ten (Supplementary Materials Table S4). Combined with the ability of mivaRNAs to functionally repress luciferase gene expression, these data suggest that mivaRNAs can regulate host gene expression during infection. However, further studies are needed to determine whether other viral genes or indirect effects of viral replication affect host gene expression. seen in Table 7, the top ten pathways enriched from our RNA-seq data are very similar between HAdV14-and HAdV14p1-infected cells. The two pathways that are not the same in the top ten are not far out of the top ten (Supplementary Materials Table S4). Combined with the ability of mivaRNAs to functionally repress luciferase gene expression, these data suggest that mivaRNAs can regulate host gene expression during infection. However, further studies are needed to determine whether other viral genes or indirect effects of viral replication affect host gene expression.

Discussion
The adenoviral VA RNA gene is a multifunctional non-coding RNA molecule. VA RNAs are processed by DICER into mivaRNA that resemble and function like cellular miRNAs. While all human adenoviruses encode at least one VA RNA gene, the amount of sequence similarity between VA RNA genes from all human adenoviruses is limited between all human Ad, but is somewhat conserved within Ad groups. In this study, we sought to identify the mivaRNA produced during prototype HAdV14 deWit infection, assess predicted target gene expression, and compare mivaRNA expression with that of the emergent pathogenic strain HAdV14p1, using A549 cells that are permissive for HAdV14 and HAdV14p1 infection.
HAdV14, like other group B2 HAdV, only encode one VA RNA gene [44]. In HAdVs that only encode one VA RNA, that VA RNA resembles group C HAdV VA RNAI. Here, further bioinformatic analysis of the HAdV14 VA RNA gene has shown that the VA RNA gene contains both RNA Pol III Box A and B promoter sites as well as two internal terminator sequences and a backup terminator sequence. The HAdV14 VA RNA gene also encodes two potential OAS1 consensus activation sequences. Based on the predicted secondary structure of the HAdV14 VA RNA molecule, it includes a terminal stem, central domain, and apical stem. Overall, the bioinformatic results suggest that, like other HAdVs that only encode one single VA RNA gene, HAdV14 VA RNA resembles that of VA RNAI, rather than VA RNAII.
To characterize the mivaRNA produced from the HAdV14 VA RNA gene, we infected A549 cells and sampled mivaRNA expression through all phases of the Ad14 infectious cycle from very early (6 hpi) all the way through to full CPE (48 hpi). At all times, >90% of all small RNA reads that mapped to the HAdV14 genome mapped specifically to the VA RNA gene. These results are consistent with those seen during Ad5 infection of IMR-90 cells [48]. Zhao and colleagues first observed that there was temporal expression of mivaRNA molecules during the infection cycle [48]. We also observed temporal expression of mivaRNAs during HAdV14 infection. Our data showed that early during infection, the 3 mivaRNAs are the most abundant HAdV14 mivaRNAs, while between 24 and 36 hpi, the 5 mivaRNAs begin to be more abundant. This is opposite of what was observed following HAdV2-infected IMR-90 cells. The reason for this difference is unclear. The difference could be virus-specific (HAdV2 vs. HAdV14), cell type-specific (lung fibroblast, IMR-90, vs. lung epithelial, A549), or a combination. Supporting this theory are studies by Kamel et al., which showed that infection of 293 cells with HAdV5, HAdV11 (a group B2 HAdV), and HAdV37 resulted in more 3 mivaRNA expression early during infection [27]. The same temporal expression patterns were seen during HAdV14p1 infection, as was observed by Zhao [48]. When we normalized mivaRNA expression, we observed that there was a significant difference in the expression of total mivaRNAs between HAdV14 and HAdV14p1 over time. Specifically, early in the infection cycle, HAdV14 mivaRNAs showed higher expression than HAdV14p1. After 24 hpi, this pattern is flipped, with HAdV14p1 mivaRNAs expressed at a higher level than HAdV14. Previous studies in our lab and others have shown that there is no difference in the replication rates of HAdV14 or HAdV14p1 in either tissue culture cells or in vivo in the lungs of Syrian hamsters, thus eliminating differences in viral replication as the reason for the differential expression of mivaRNAs [8,49]. The studies here were not designed to test whether this mivaRNA expression pattern has any role in the increased pathogenesis seen during HAdV14p1 infection. Studies are underway to test this possibility using a Syrian hamster model of differences in HAdV14p1 vs. HAdV14 lung pathogenesis [8,50].
Our results identified the core seed sequences produced during Ad14/Ad14p1 infection as coming from both the 5 and 3 terminal stems. Minor variants are also produced. How these arise is unclear. One possibility is that RNA pol III synthesis of VA RNAs might have some minor wobble in its exact initiation site, as these minor mivaRNA species show very slight variations in either the 5 or 3 ends. It is possible that this could have an effect on expression, as small differences in the 5 start can result in a completely different mivaRNA seed sequence. We identified a group of long mivaRNAs produced from the 3 terminal stem that start between 8 and 4 nucleotides upstream of the predominant 3 mivaRNA. Interestingly, the longest mivaRNA produced from the 3 terminal stem has the exact same seed sequencing as hsa-miR-584c-3p. It has been well documented that HAdV mivaRNAs can regulate gene expression, like gene-encoded miRNAs [26,40,42]. Using a luciferase reporter vector, we have shown that both the 5 and 3 HAdV-14 mivaR-NAs can suppress gene expression during viral infection. Our results confirm that the 3 mivaRNAs are more abundantly expressed than 5 mivaRNAs early during infection, as the 3 mivaRNAs showed stronger repression of luciferase activity than the 5 mivaRNAs. mivaRNAs have been shown to associate with RISC complexes, but not all mivaRNAs associate with RISC complexes at the same rate [27,38,40,42]. This could account for the difference in repressive activity seen between the mivaRNAs in our reporter assay. The association of HAdV14 mivaRNAs with the RISC complex was not a focus of these studies, but is under investigation.
The ability of mivaRNAs to target and regulate host gene expression is an attractive theory for potential roles in HAdV replication. As a result, numerous studies have been conducted to determine potential target genes for mivaRNAs. Bioinformatic, microarray studies, and RNA-seq analysis have shown that many of the potential cellular gene targets for mivaRNAs are involved in DNA repair, cell growth, cell death, RNA metabolism, cell signaling pathways, and metabolism [26,42,48]. Our bioinformatic analysis (Table 6) on the predicted HAdV14 mivaRNA targets show similar cellular pathways to be potential targets, despite the lack of sequencing identity between HAdV14 seed sequences and HAdV2/5 seed sequences. Our IPA analysis (Table 7) also shows functional enrichment of the same cellular pathways from our RNA-seq data on HAdV14 mivaRNA target genes. This suggests that there is some conserved biological reason for mivaRNA expression during viral infection. However, studies have raised questions as to the role of mivaRNAs in viral replication. The mutation of the seed sequences of HAdV5 mivaRNAs had no impact on viral replication in HEK293 cells [51]. One underlying problem with studies (including these studies) on HAdV mivaRNAs is that they have been conducted in tissue culture cell lines rather than relevant primary cells or in a permissive animal model, which might clarify the role of mivaRNAs during infection in vivo. HAdV mivaRNAs have also been found in persistently infected lymphoid cells, suggesting a potential role in regulating lytic and persistent infection [52,53]. In addition to mivaRNAs from the VA RNA gene, another viral small RNA that maps to the Major Late Promoter region has been identified that regulates viral gene expression; we are in the process of exploring if there are other HAdV14 small viral RNAs [54].
In this study, we characterized the HAdV14 VA RNA gene and the mivaRNAs produced from it. We identified the predominant seed sequences and the cellular genes that the mivaRNAs can target. Our bioinformatic analysis using RNA-seq data of the expression of mivaRNA target genes showed that HAdV14 mivaRNA can regulate similar pathways seen from HAdV2/5 studies despite the lack of identity of the mivaRNA seed sequences. Despite 100% sequence identity later in infection, HAdV14p1-infected cells express significantly more mivaRNA than HAdV14-infected cells. Whether this difference has an impact on viral replication or pathogenesis is unknown. Due to a lack of studies in immunocompetent animal models, the role of HAdV mivaRNA in viral pathogenesis is largely unexplored. To answer questions about viral factors, such as mivaRNAs, we have developed an immunocompetent Syrian hamster model of HAdV14p1 lung pathogenesis that displays pathology similar to what is seen in the lungs of severely infected human patients [8,50]. We propose that the Syrian hamster model of HAdV14p1 lung pathogenesis could provide the model system to explore the roles of mivaRNA in viral replication and pathogenesis in vivo.