Excessive production and extreme editing of human metapneumovirus defective interfering RNA is associated with type I IFN induction

Type I IFN production is one of the hallmarks of host innate immune responses upon virus infection. Whilst most respiratory viruses carry IFN antagonists, reports on human metapneumovirus (HMPV) have been conflicting. Using deep sequencing, we have demonstrated that HMPV particles accumulate excessive amounts of defective interfering RNA (DIs) rapidly upon in vitro passage, and that these are associated with IFN induction. Importantly, the DIs were edited extensively; up to 70 % of the original A and T residues had mutated to G or C, respectively. Such high editing rates of viral RNA have not, to our knowledge, been reported before. Bioinformatics and PCR assays indicated that adenosine deaminase acting on RNA (ADAR) was the most likely editing enzyme. HMPV thus has an unusually high propensity to generate DIs, which are edited at an unprecedented high frequency. The conflicting published data on HMPV IFN induction and antagonism are probably explained by DIs in virus stocks. The interaction of HMPV DIs with the RNA-editing machinery and IFN responses warrants further investigation.

As a first line of defence, the innate immune system protects cells against viral infections through production of type I IFNs. For successful infection and replication, viruses use a variety of strategies to evade the host innate immune system (reviewed by Randall & Goodbourn, 2008). Respiratory syncytial virus (RSV), the mammalian pneumovirus most closely related to HMPV, uses two non-structural proteins (NS1 and NS2) as IFN antagonists (Spann et al., 2005). Homologues of NS1 and NS2 are absent in the genome of HMPV (van den Hoogen et al., 2002), indicating that HMPV must have a different strategy to evade the innate immune system. It has been suggested that HMPV is a strong activator of the retinoic acid-inducible gene 1 (RIG-I) and/or mitochondrial antiviral-signalling protein (MAVS) signalling pathways leading to IFN production (Baños-Lara et al., 2013;Bao et al., 2007;Liao et al., 2008). In contrast, the interaction of the attachment protein (G) of HMPV with RIG-I resulted in inhibition of RIG-I activation . However, using small interfering RNA (siRNA) silencing of G, an antagonistic role of G was not confirmed (Preston et al., 2012). Goutagny et al. (2010) reported that only NL/1/00 induced IFN production, in contrast to NL/1/ 99. Using the same virus strains, we could not reproduce IFN production upon NL/1/00 infection. These discrepancies between different studies led us to hypothesize that defective interfering RNAs (DIs) present in virus stocks may be responsible for the reported activation of the IFN pathway. DIs consist of partially deleted viral genomes that arise spontaneously when virus stocks are generated by passaging at high m.o.i. in mammalian cells, due to errors made by the replicase complex (Huang, 1973). The presence of DIs in virus stocks was reported in the 1960s and their generation by diverse DNA and RNA viruses has been widely documented (Rima et al., 1977;Sekellick & Marcus, 1978;Strahle et al., 2006). DIs with a 'copyback' or 'snapback' structure are strong IFN inducers, which correlates with their ability to self-anneal and form dsRNA (Sekellick & Marcus, 1978;Strahle et al., 2006). 'Copyback' types are generated during generation of the negative-sense strand by the replicase, which leaves its positive-sense template and resumes synthesis at the beginning of the 'daughter' genome negative-sense chain that it is still carrying. The complementary ends of the two chains are responsible for the circular 'panhandle' structures that these genomes form as naked RNA (Kolakofsky, 1976). 'Snapback' DIs are generated in a similar way to the 'copyback' DIs; however, the replicase does not only copy back the 39 end but also, when it crosses to the nascent strand, copies back the complete strand, resulting in a covalently linked positive-and negative-sense strand, forming a duplex (Lazzarini et al., 1981).
Here, we have demonstrated for the first time, to our knowledge, that HMPV accumulates DIs with a snapback structure rapidly upon in vitro passage, causing activation of the IFN pathway upon infection. Strikingly, the genomes of each of the DIs displayed extensive AAG or TAC hypermutation. In addition, our data indicated that the RNA-editing enzyme adenosine deaminase acting on RNA (ADAR) was the most likely enzyme responsible for editing of HMPV DIs.

IFN production upon HMPV infection
The ability of HMPV to induce IFN production was tested following inoculation of A549 cells with two different virus stocks of recombinant HMPV strain NL/1/00. One stock was generated with a maximum of two passages in Vero-118 cells at low m.o.i. of 0.01 (hereafter named P2 low ), and a second stock was generated by five passages in Vero-118 cells at a high m.o.i. of 3 (hereafter named P5 high ).
A549 cells were inoculated with these stocks as well as with sendai virus (SeV), measles virus Edmonston strain (MeV-Edm) and parainfluenza type 5 (PIV-5) as controls. PIV-5 does not activate the IFN pathway (Childs et al., 2012), whilst SeV (Cantell strain) and MeV-Edm are known to be strong inducers (Shingai et al., 2007;Strahle et al., 2006;Takaki et al., 2011). At 24 h after inoculation at an m.o.i. 3, the supernatants were tested for IFN content and compared with mock-inoculated cell supernatants. The supernatant of P2 low -inoculated cells displayed similar low amounts of IFN to the supernatants of cells inoculated with PIV-5. In contrast, P5 high -inoculated cells produced nearly as much IFN as SeV-or MeV-Edm-inoculated cells (Fig.  1). The lack of IFN induction by P2 low was not due to a lower infection efficiency, as FACS analysis 24 h after inoculation revealed infection of more than 90 % of the cells for all viruses. HMPV NL/1/00 thus only activated the IFN pathway when virus stocks were prepared at a high m.o.i. The ability of SeV, vesicular stomatitis virus (VSV), MeV-Edm and PIV-5 to strongly activate the IFN pathway has been associated with the presence of DIs in these virus stocks (Baum et al., 2010;Killip et al., 2013;Marcus & Sekellick, 1977;Sekellick & Marcus, 1982). In addition, studies with VSV revealed that the IFN-inducing capacity of DIs was not affected by UV treatment (Marcus & Gaccione, 1989). UV treatment of P5 high , P2 low and SeV virus stocks resulted in reduced infection of A549 cells with P5 high and SeV but did not affect IFN induction by these viruses. This suggested that the IFN induction by P5 high may be attributed to the presence of DIs (Fig. S1, available in the online Supplementary Material).

Deep sequencing demonstrates the presence of DIs
To investigate the possible presence of DIs, RNA from virus stocks of P2 low and P5 high was subjected to arbitrarily primed 454 sequencing. Mapping of the sequence reads against the full-length HMPV genome resulted in 28 121 and 8200 sequence reads for P2 low and P5 high , respectively. The genomes of P2 low and P5 high were covered with a mean of 588 and 168 reads per nucleotide, respectively. No coverage was obtained for the last 50 and 4 nt for P2 low and P5 high , respectively. In addition, the ends of the genomes were covered with fewer than 10 reads from nt 13 164 for P2 low and from nt 13 325 for P5 high . Fig. 2(a), depicting the ratio between the depths of sequence obtained for P2 low and for P5 high per nucleotide position, demonstrates that, at the end of the genome, P5 high displayed a greater depth of sequence than P2 low . From nucleotide 11 726 onwards, P5 high displayed up to 31 times more depth of sequence than P2 low . In fact, for P5 high 16.3 % of the 8200 reads mapped to the region beyond position 11 kb of the genome, compared with only 2.2 % of 28 121 reads obtained for P2 low . Examination of the P5 high reads at this part of the genome revealed the presence of reads that only aligned for approximately 50 % with the reference sequence. Subsequent BLAST searches conducted with all sequence reads mapping for at least 20 % to the reference sequence revealed the presence of reads that contained two copies of (almost) the same genome region: one part of the read mapped to the genomic negative-sense sequence and the other part mapped to the antigenomic positive-sense sequence of the viral genome, with the two parts being complementary. In total, 47 dsRNA structures were identified in P5 high and none in P2 low . All 47 reads aligned to the region upstream of position 11 200 of the viral genome, the same region where P5 high displayed deeper coverage than P2 low . Alignment of these reads revealed 12 unique structures that differed in sequence and position, as shown in Fig. 2(b). As assignment of the orientation of the two strands could not be carried out based on the deepsequencing results, the assignment was based on the results obtained in the BLAST search. The complementarity of the two strands indicated that P5 high contained DIs with a 'snapback' structure. The 12 DIs had variable lengths of the single-stranded parts at the cross-over point (i.e. the position where the polymerase switches from the template to the daughter strand and where the positive-and negative-sense strands of the DI are connected). No specific sequence motif was detected at the position of cross-over points that may have served as a signal for the polymerase to switch strands. Table S1 depicts the exact positions of all DIs compared with the viral genome.

Detection and quantification of DIs by Northern blot assays
To replicate, DIs need to contain polymerase initiation sites at the 59 end (genomic/negative-sense strand) and the reverse complement of this trailer at the 39 end (antigenomic/positive-sense strand). The deep-sequencing and bioinformatics approach only identified a read as a DI if the read overlapped with a cross-over point. All DI reads were too short to reach from the cross-over point to the trailer sequences or the reverse complement thereof. To test whether the double-stranded structures detected with deep sequencing were indeed DIs with a 'snapback' structure, we performed electrophoresis of viral RNA and Northern blotting using a probe complementary to the trailer. The Northern blot assay demonstrated the presence of similar amounts of full-length viral genome in P5 high and P2 low (Fig. 3, top band) with an additional band of approximately 7 kb detected in both virus stocks. The identity of the 7 kb fragment is unknown. In addition to these large RNA molecules, smaller RNAs ranging in size from 0.5 to 4 kb were detected but only in P5 high . As the RNA was analysed on denaturing gels, fragments of this size may represent DIs, double-stranded in nature, of 0.25-2 kb. Quantification using a ChemiDoc Imaging system indicated that 80 % of the total RNA, in the size range 0.5-4 kb for P5 high , was probably of DI origin.
In theory, the detection of dsRNA structures could have been a result of the process of deep sequencing, as this method involves ligation of linkers to PCR products. Reverse transcription (RT)-PCR assays directly on the virus stocks conducted with primers designed based on the consensus sequence obtained for DI # 1 (Fig. S2) confirmed the presence of DIs in the original virus stocks. Besides DI # 1, these RT-PCR assays identified additional DIs that were not identified during analysis of the 454 sequencing reads (Fig. S3).

HMPV DI RNAs are hypermutated
Analysis of the sequence variation of the deep-sequencing reads revealed extensive AAG and TAC mutations at the Type I IFN induction by hypermutated HMPV DIs end of the genome of P5 high compared with the reference sequence, and this was not observed in P2 low (Fig. 4a). The location of hypermutations matched with the region where P5 high displayed more sequence depth than P2 low and where DIs were detected. In the region downstream of nt 11 000, 53.6 % of the A residues in the reference sequence displayed mutation to G and 35.7 % of the T residues displayed mutation to C (in at least 2 % of the reads mapping to that position) (Fig. 4b, c). In the region beyond position 12 000, 75 % of the reads displayed hypermutation and the mutation percentages in this region were even higher: 87.3 % AAG and 82.8 % TAC mutation. Individual reads varied in the position and number of mutations, with a range of 5-77 % of the A residues mutated to G and 5.7-57.1 % of the T residues mutated to C. We hypothesized that this kind of hypermutation was probably caused by the host RNA-editing enzyme ADAR. ADAR binds to dsRNA and deaminates A to I, which is recognized as a G during virus replication (Bass, 2002). Three human ADAR genes have been identified, and two of these have editing activity, ADAR1 and ADAR2.

Role of ADAR in editing of DIs
We analysed whether the editing patterns in HMPV DIs were consistent with those described for ADAR editing by Eggington et al. (2011) who, refining earlier work by Lehmann & Bass (2000), determined that ADARs preferentially target As with a certain 59 and a variable 39 neighbour. They found that ADAR editing of AAI was more likely when the 59 neighbour was U, followed by A, C and G (U.A.C.G). The 39 preference was less specific, with a different preference for each of the four ADARs they considered and equal preference for some nucleotides. We analysed all reads (n5954) mapping downstream of nt 11 000 for the editing probability depending on the 39-and 59-neighbour nucleotide. The mean edit probability for each of the 16 combinations of 59 and 39 neighbours of an edited A or T was calculated. No pattern was identified for the 39-neighbouring nucleotide. In contrast, the 59-neighbour preference was in accordance with the probabilities as described for ADAR (Eggington et al., 2011), and was the same for AAG and TAC mutation (Fig. 4d).
In human cells, two isoforms of the ADAR1 protein are generated through the use of alternative promoters and alternative splicing. The long form (ADAR1-p150) is IFN inducible and localizes mostly in the cytoplasm, whereas a short form (ADAR1-p110) is expressed constitutively and localizes, like ADAR2, predominantly in the nucleus (George et al., 2011). Vero cells are derived from African green monkeys, which might have different ADAR genes compared with human cells. RT-PCR assays performed with primer sets designed based on the human ADAR1 sequences in combination with the available African green monkey ADAR1 sequences detected mRNA of both ADARs in Vero-118 cells (Fig. S4a). Subsequent sequencing of the amplified fragments revealed that the Vero-ADAR1-p150 and Vero-ADAR-p110 mRNA displayed 95 % nucleotide sequence identity with the human ADAR1-p150 and ADAR1-p110 mRNAs and 99 % with those of macaques (Fig. S5). In addition, Western blotting using an antibody against human ADAR1 revealed expression of ADAR1-p110 and ADAR1-p150 in Vero-118 cells, and upon IFN treatment a slight increase in expression of ADAR1-p150 was observed (Fig. S4b).
Given the cytoplasmic replication of HMPV, the editing motifs and the expression of ADAR mRNA, we suggest that ADAR1 is the most likely enzyme responsible for editing of HMPV DIs.

DISCUSSION
As HMPV causes an acute respiratory infection in mammals it must -like other respiratory viruses -employ a mechanism to counteract the IFN response of the innate immune system. Several studies have reported interactions between HMPV and the innate immune system, although with conflicting conclusions. We did not observe induction of IFN production for HMPV NL/1/00 and NL/1/99 upon passage at low m.o.i. Only after a few passages at high m.o.i. was HMPV-induced IFN production observed. We showed that this induction was caused primarily by the presence of DIs with a 'snapback' structure in the virus stocks, providing a plausible explanation for the discrepancies in the data reported so far. The current study was conducted in Vero-118 cells. HMPV virus stocks have been generated in LLC-MK2 cells by other groups (Ren et al., 2014;Velayutham et al., 2013). However, in our hands, these cells do not produce sufficiently high viral titres to allow passage at high m.o.i. Future research needs to elucidate whether HMPV does generate DIs in LLC-MK2 cells or any other cells.
In theory, the detection of dsRNA structures could have been a result of the process of deep sequencing, as this method involves ligation of linkers to PCR products. However, no double-stranded structures were observed in P2 low , which was analysed using the same protocol. In addition, both RT-PCR and Northern blots were performed using RNA isolated directly from the virus stock.
Of note, the deep-sequencing analyses, Northern blots and RT-PCR assays performed using several independent P5 high stocks consistently revealed the presence of DIs and similar patterns of hypermutation. Our results indicated that different P5 high stocks contained various DIs.
Northern blot analyses revealed that the DIs ranged in size from 0.5 to 2 kb and represented more than 80 % of the RNA. Based on the Northern blot assays, the virus stock contained a larger number of DIs than detected with deep sequencing. As smaller RNAs transfer more efficiently compared with larger RNAs in Northern blotting, 80 % might be an overestimation. However, the deep sequencing may be less efficient in obtaining reads mapping to the 59 end of the genome, which might result in an underestimation of the number of DI reads. In addition, many more reads in the deep-sequencing analysis might have been of DI origin but not identified as such as they lacked the cross-over point. and nt 13 075-13 125 (bottom) for P5 high . Green bars, mutation to A; black bars, mutation to T; blue bars, mutation to C; red bars, mutation to G. (c) Mutation matrices (%) for deep-sequencing reads obtained for the region downstream of nt 11 000 of the genome. Mutation was calculated when at least 5 % of reads displayed a mutation for that nucleotide. (d) Analysis of edited reads for the edit probability of an A or T, given a particular 59-neighbour nucleotide. The highest probability for AAG or TAC editing was detected where the 59-neighbour nucleotide was a U. Red bars indicate the AAG edit probability for the 59 neighbour, and blue bars the TAC edit probability for the 39 neighbour after complementing the neighbour with its base pair.
Analysis of mutations in the deep-sequencing reads obtained for P5 high and P2 low revealed the presence of AAG and TAC mutation only in P5 high and only in the region where DIs were identified. In addition, all reads that were identified as DI reads and the DI-PCR products displayed hypermutation. Although we cannot exclude hypermutation of the viral genome, the location of the observed hypermutation (in the polymerase gene) would, in that case, probably result in decreased virus replication. We observed similar replication of P2 low and P5 high , arguing against hypermutation of the viral genome itself (Fig. S6).
So far, two mammalian RNA-editing enzymes are known: apolipoprotein B mRNA editing enzyme, catalytic polypeptide-like (APOBEC) and ADAR (Hamilton et al., 2010). APOBEC deaminates C to U on ssDNA or RNA. APOBECinduced specific CAU editing has been described for retroviruses and DNA viruses but not for RNA viruses (Bransteitter et al., 2009). Commercial gene expression assays did not detect expression of APOBEC mRNA in our Vero-118 cells (data not shown), whilst we did detect ADAR mRNAs. ADAR binds to dsRNA and deaminates A to I, which, during replication, is recognized as a G (Bass, 2002). Three ADAR genes have been identified, of which two have editing activity, ADAR1 and ADAR2. The IFNinducible ADAR1-p150 resides in the cytoplasm, whilst the constitutively expressed ADAR1-p110 and ADAR2 reside in the nucleus (George et al., 2011). For viruses that replicate in the cytoplasm of the infected host cell, such as HMPV, ADAR1-p150 is the most likely protein responsible for AAI editing (Samuel, 2011). We did not observe a clear upregulation of ADAR1-p150 mRNA in IFN-treated Vero-118 cells, and only a low level of upregulation of protein in IFN-treated A549 and Vero cells. In accordance, in studies detecting expression of endogenous ADAR1-p150, only low expression levels were observed (Cavarec et al., 2013). Vero-118 cells are IFN deficient (Desmyter et al., 1968), which would argue against a role for ADAR1-p150. However, despite the absence of IFN genes in Vero cells, it is well known that expression of IFN-stimulated genes (e.g. ISG56) can be upregulated in an IFN-independent manner (Gélinas et al., 2011). We did observe upregulation of ISG56 mRNA upon SeV or P5 high inoculation of Vero-118 cells (Fig. S7). Tripathi et al. (2014) found indications for IFN-independent induction of ADAR expression through the concerted action of IRF3 and IRF7, although the induction levels were low. However, the promoter that initiates the ADAR1-p150 transcript has a basal constitutive activity, which may be sufficient for editing activity (Patterson & Samuel, 1995).
ADAR-induced hypermutation of MeV has been studied in detail, but it is unclear why this hypermutation was not observed upon inoculation of Vero cells with MeV (Suspène et al., 2011;Wong et al., 1989). Although we have strong indications for the role of ADAR in HMPVinfected Vero cells, studies in ADAR-deficient cells and RNA-binding studies should elucidate the actual role of ADAR in editing HMPV (DI) RNA. The restricted host range of HMPV does not allow studies in ADAR-deficient MEFs or HeLa cells, and the low expression of ADAR will be a challenge for RNA-binding studies.
DIs have been demonstrated to play a role in the pathogenesis of some virus infections, such as in acute dengue virus and influenza virus infections (de Chassey et al., 2013;Li et al., 2011;Saira et al., 2013). Given the relative ease with which HMPV accumulates DIs and the role demonstrated for DIs, the role that DIs may play in the pathogenesis of HMPV is clearly of interest.
In addition, DIs present in virus stocks have been shown to enable persistent infections both in vivo and in vitro (Cave et al., 1985;Roux et al., 1991;Weiss et al., 1983). This is thought to be a result of activation of the innate immune system by DIs (Sekellick & Marcus, 1978). For HMPV, prolonged shedding of virus after experimental infection of mice has been reported (Alvarez et al., 2004;Hamelin et al., 2006;Liu et al., 2009). It would be interesting to investigate whether the virus stocks used in these persistently infected animals contained DIs, as described for RSV Mejías et al., 2008).
Extensive AAI editing has been described in RNAs from MeV isolated from patients with persistent central nervous system infection (Cattaneo et al., 1988;Wong et al., 1989). Subsequent studies demonstrated that ADAR1 acts as a proviral host factor in the context of MeV infection (Toth et al., 2009). In contrast, ADAR has an antiviral function in other virus infections (Samuel, 2011;Taylor et al., 2005). Future research needs to elucidate the roll of ADAR in HMPV pathogenesis.
In conclusion, upon high m.o.i. passaging, HMPV accumulated DIs rapidly, and these DIs were responsible for a robust induction of IFN production. Our results should redirect research aiming to elucidate how HMPV subverts the innate immune system. Based on the high rate and the patterns of hypermutation, the preference of the 59-neighbouring nucleotide and the expression of ADAR in Vero-118 cells, we suggest that HMPV DIs are edited by ADAR. Based on the fact that HMPV (DIs) replicates in the cytoplasm of infected cells, we speculate that HMPV (DI) genomes are hypermutated by ADAR1-p150.

METHODS
Cells and viruses. 293-T cells and Vero-118 cells were grown as described previously . A549 cells were grown in Ham's F-12 medium (Invitrogen) supplemented with 10 % FCS, High Clone (HC-FCS; Greiner Bio-One), with 100 IU penicillin ml 21 , 100 mg streptomycin ml 21 and 2 mM glutamine (PSG). The generation of recombinant HMPV NL/1/00 has been described previously  Bioassays for IFN. The IFN content in supernatants was measured using the ISRE-firefly luciferase reporter construct as described previously (Patel et al., 2013).
Arbitrarily primed PCR and virus genome sequencing were performed as described previously (van Boheemen et al., 2012). GS Junior sequence reads were trimmed by 25 nt at both the 59 and 39 ends of the reads to remove primer sequences and aligned to the HMPV genome (GenBank accession no. AF371337) using CLC Genomics software v.4.6.1 (CLC bio). Using this software, singlenucleotide polymorphism analysis was conducted with a minimum of 2 % variation frequency and a minimal coverage of five reads per position.
Northern blot analysis. RNA was isolated from purified virus stock (10 7 TCID 50 per sample) using TRIzol (Invitrogen) according to the manufacturer's instructions. After electrophoresis of the RNA on a 1 % agarose gel, RNA was transferred to a positively charged nylon membrane (BrightStar-Plus; Ambion) and incubated with the probe. The probe, complementary to the 59 end of genome, was generated with a Strip-EZ RNA system (Ambion) using a PCR fragment based on the full-length cDNA construct of NL/1/00  with PCR primers 59-TAATACGACTCACTATAGGGGAAAATGA-TAAAATG-39 (forward primer, T7 promoter underlined) and 59-GGTTTTTTTGCCGT-39 (reverse complement to end of trailer). The blot was analysed with a NorthernMax-Gly system (Ambion), exposed to a phosphorimager screen and analysed with a Storm 850 PhosphorImager and ChemiDoc MP Imaging System and Image Lab v.4.0.1 software (Molecular Dynamics).