In Yokohama City between January 2013 and June 2016, 1308 clinical specimens (throat swabs, nasal swabs, nasal secretions, and nasal aspirate fluids) were collected from patients suffering upper or lower ARIs in 16 sentinel hospitals and clinics (eight pediatric clinics, four internal medicine clinics, and four hospitals) participating in the National Epidemiological Surveillance of Infectious Diseases (NESID), instituted by the Infectious Diseases Control Law in Japan, and in eight other medical institutions (one pediatric clinic and seven hospitals) (Anon, 2010). There were 113 hospitals and 2962 clinics in Yokohama City, and the population in the city was ∼3.7 million. Thus, the numbers of hospitals and clinics participated in this study were ∼9.7 and ∼0.4%, respectively, in the city. Before collecting the clinical specimens in which to analyze the viruses causing ARIs, the physicians at each medical institution obtained the informed consent of the patients or their guardians. The de-identified clinical specimens were sent to Yokohama City Institute of Public Health and subjected to multiplex RT–PCR with the Seeplex^® RV15 OneStep ACE Detection kit (Seegene, Seoul, South Korea), which identifies 15 major respiratory viruses. The clinical specimens that were positive for HMPV were analyzed further.

The RNAs were purified from clinical specimens (140 μl) using the QIAamp Viral RNA Mini Kit (Qiagen, Hilden, Germany), according to the manufacturer’s instruction, and dissolved in 60 μl of distilled water. The RNA in 5 μl of the purified RNA solution was reverse transcribed, and the G gene was amplified from the cDNA with RT–PCR using the PrimeScript II High Fidelity One Step RT–PCR Kit (TaKaRa Bio, Otsu, Japan) and a pair of G-gene-specific primers, SH7 (5′ primer) (van den Hoogen et al., 2004) and GR (3′ primer) (Banerjee et al., 2011). The reaction mixture was prepared in a 50 μl solution, and the reverse-transcription reaction was performed at 45°C for 10 min. The reverse transcriptase was inactivated at 94°C for 2 min. The G gene cDNA was amplified with PCR with 40 cycles of 10 s at 98°C, 15 s at 56°C, and 10 s at 68°C. The cDNA products were separated electrophoretically in 3% NuSieve^TM 3:1 Agarose gel in 0.5% TBE buffer. When no G-gene-specific band was detected or was only barely detectable, 5 μl of the RT–PCR product was subjected to semi nested PCR. In this assay, the reaction mixture was prepared in a 50 μl solution, and the RT–PCR products were denatured at 94°C for 2 min. The G gene was then PCR amplified with 40 cycles of 10 s at 98°C, 15 s at 50°C, and 30 s at 68°C with Tks Gflex^TM DNA Polymerase (TaKaRa Bio) and a pair of G-gene-specific primers, SH7 (5′ primer) and GNR2 (3′ primer). The sequence of the GNR2 primer was 5′-GGATTCATTAAGAGGATCCATTG-3′. The products of the semi nested PCR were separated electrophoretically to detect the G gene.

The nucleotide sequences of the G gene of the HMPV strains were determined with direct sequencing. The amplified HMPV G genes were purified with the QIAquick PCR Purification Kit (Qiagen). The PCR products were subjected to cycle sequencing with the BigDye Terminator ver. 1.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA) and the SH7, GR, and GNR2 primers. The products of the cycle sequencing reactions were purified through Centri-Sep Spin Columns (Thermo Fisher Scientific), and the nucleotide sequences were determined with a 3500 Genetic Analyzer (Applied Biosystems). The nucleotide sequence data for the HMPV strains were aligned and edited with the BioEdit software (ver. 7.2.5) (Hall, 1999).

A multiple-sequence alignment was constructed with the MAFFT software (ver. 7.304b) using the default settings (Katoh et al., 2002). Phylogenetic analyses were performed with the maximum likelihood method in the MEGA software (ver. 7.0.20) (Kumar et al., 2016), and the statistical significance of the tree topologies was tested with bootstrapping (50 replicates).

Potential acceptor sites for N- and O-linked sugars were predicted with the NetNGlyc 1.0¹ and NetOGlyc 3.1 programs², respectively.

The nucleotide sequence data reported in the present study were deposited in the DDBJ/EMBL/GenBank nucleotide sequence database under accession numbers LC192170–LC192253.

The overall rate of evolutionary change (nucleotide substitutions per site per year) in the HMPV G gene was estimated with the BEAST 2 program³, which uses a Bayesian Markov Chain Monte Carlo (MCMC) approach (BEAST 2: A Software Platform for Bayesian Evolutionary Analysis⁴). An aligned G gene sequence dataset from the 84 HMPV strains detected in Yokohama City from 2013 to 2016 was analyzed with the BEAST 2 program using a strict HKY model. The MCMC chain was run for a sufficient length of time to ensure convergence (all expected sample size values exceeded 200 with a 10% burn-in).

These analyses have been done as a part of the National Epidemiological Surveillance of Infectious Diseases, Japan (NESID) as stipulated under the Infectious Diseases Control Law, and before collecting the clinical specimens, physicians in each medical institution obtained the informed consent of the patients or their guardians. The ethics committee of Yokohama City Institute of Public Health approved this study.

The consensus nucleotide sequence was determined using both the first 180-nucleotide sequence (1st-180nts) and the second 180-nucleotide sequence (2nd-180nts) in the 15 A2b_180nt-dup HMPV strains (in total, 30 sequences of 180 nt) and was assumed to be the corresponding 180 nt of the hypothetical ancestral strain of the 15 A2b_180nt-dup HMPV strains. The number of nucleotide substitutions in the 1st- and 2nd-180nts of each A2b_180nt-dup HMPV strain was determined and compared with the corresponding 180 nt of the hypothetical ancestral strain, and the year in which 180nt-dup occurred in the HMPV genome was estimated based on the mean rate of nucleotide substitutions for the HMPV G gene.

Human metapneumovirus was detected in 91 of 1308 specimens from ARI patients, accounting for 7.0% of the specimens. The basic and clinical characteristics of the patients infected with HMPV are summarized in Table 1. The entire nucleotide sequences of the open reading frames of the G proteins were determined for 84 of the 91 HMPV strains. These 84 HMPV strains had 73 unique G gene sequences, containing at least one nucleotide substitution, whereas the remaining 11 sequences were identical to other G gene sequences. A phylogenetic analysis was performed using the G gene sequences of the 84 HMPV strains and all the HMPV G gene sequences available at the National Center for Biotechnology Information (NCBI) nucleotide sequence database⁵. In total, 574 G gene sequences of HMPV strains were used to construct the phylogenetic tree (Figure 1). Among our 84 strains, 6, 43, 10, and 25 strains were classified in subgroups A2a, A2b, B1, and B2, respectively (Figure 1). No HMPV strain belonged to subgroup A1. These data are similar to previous findings in Japan (Matsuzaki et al., 2008; Mizuta et al., 2010; Toda et al., 2010; Omura et al., 2011; Nidaira et al., 2012; Nakamura et al., 2013). Table 2 shows the numbers of HMPV strains detected between January 2013 and June 2016 in Yokohama City and their subgroup classification. Multiple subgroup strains were detected in each year. Subgroup A2b strains were detected every year, and approximately half the A2b subgroup strains detected in 2014, 2015, and 2016 contained 180nt-dup in the G gene (Figure 2). 180nt-dup is a duplication of the 180 nt at nucleotide positions 371–550 (the first nucleotide of the initiation codon of the G gene is deemed to be nucleotide position 1). No A2b strain detected in 2013 contained 180nt-dup. The 15 A2b_180nt-dup HMPV strains formed a small cluster on the phylogenetic tree (Figure 1), and the cluster also contained four classical A2b strains, which lack 180nt-dup and were detected in 2014 and 2015 in Yokohama City. These data suggest that 180nt-dup occurred in a specific ancestor of the HMPV strains in this cluster (Figure 2).

Although the 15 A2b_180nt-dup HMPV strains seemed to be derived from a common ancestor, the 180nt-dup sequences of most A2b_180nt-dup HMPV strains differed from one another (Figure 2). Even within the same viral genome, the 1st-180nt differed from the 2nd-180nt by one to six nucleotide substitutions. These observations were not unexpected, because substitutions would have occurred in each viral genome after the ancestral A2b strain had acquired 180nt-dup. RNA viruses are vulnerable to changes in their nucleotide sequences because the viral RNA-dependent RNA polymerase (RdRp) lacks proofreading activity. Therefore, the original 180-nt sequence in the hypothetical ancestral strain was predicted by constructing a consensus sequence of the 1st- and 2nd-180nts of the 15 A2b_180nt-dup HMPV strains (Figure 2). The 360-nt sequence, consisting of the 1st- and 2nd-180nts of each of the 15 A2b_180nt-dup HMPV strains, was then compared with that of the hypothetical ancestral strain (Figure 2). In the 1st-180nt, strain HMPV/Yokohama.JPN/P7916/2015 had the same nucleotide sequence as the ancestral strain, whereas the other 14 strains had one to six nucleotide substitutions relative to the sequence of the ancestral strain (Figure 2). In the 2nd-180nt, two strains (HMPV/Yokohama.JPN/P7450/2014 and HMPV/Yokohama.JPN/P7875/2015) had the same nucleotide sequence as the ancestral strain, whereas the other 13 strains had one to six nucleotide substitutions relative to the sequence of the ancestral strain (Figure 2). We then estimated the evolutionary rate of the HMPV G gene with the BEAST 2 program to predict the year in which 180nt-dup occurred. The evolutionary rate was estimated to be 4.3 × 10^-3/site/year (95% highest probability density: 3.2–5.4 × 10^-3/site/year), consistent with the rate previously estimated by de Graaf et al. (2008). Based on these data, 180nt-dup was predicted to have occurred in the ancestral strain between 2011 and 2013. Figure 3 shows a phylogenetic tree constructed with the 15 1st-180nts and the 15 2nd-180nts. The tree was rooted with the HMPV A2a strain HMPV/Yokohama.JPN/P7011/2013. Theoretically, the 1st- and 2nd-180nts in each strain should be identical at the time when the duplication occurred. However, among the 15 A2b_180nt-dup HMPV strains, the 1st- and 2nd-180nts of 13 strains have diverged randomly and the 1st- and 2nd-180nts of each of these 13 strains were not located in the same cluster. These data suggest that these 13 strains had the same ancestral strain but acquired nucleotide substitutions in both the 1st- and 2nd-180nts independently of each other. In contrast, the 1st- and 2nd-180nts of two A2b_180nt-dup HMPV strains, HMPV/Yokohama.JPN/P7886/2015 and HMPV/Yokohama.JPN/P7929/2015, were located in a distinct small cluster on the phylogenetic tree, suggesting that these two A2b_180nt-dup HMPV strains had a different ancestral strain.

Figure 4 shows a nucleotide sequence alignment of a representative A2b_180nt-dup HMPV strain (HMPV/Yokohama.JPN/P7450/2014) and a classical A2b HMPV strain (HMPV/Yokohama.JPN/P7015/2013), and their deduced amino acid sequences. Figure 5 shows amino acid sequence alignment of the G protein of 15 A2b_180nt-dup HMPV strains, the hypothetical ancestral strain, and the representative classical A2b strain (HMPV/Yokohama.JPN/P7015/2013). Compared with the classical A2b strain, the A2b_180nt-dup strains have a 60-amino-acid duplication (60aa-dup) in the G protein. Because the 180-nt region at positions 371–550 is duplicated, a codon (AGG) for arginine encoded at nucleotide positions 550–552 is replaced with a codon (ACA) for threonine (nucleotides in the duplicated sequence are underlined; the first nucleotide of the initiation codon of the G gene is deemed to be nucleotide position 1), which is followed by the duplicated nucleotide sequence encoding 60aa-dup (Figure 4). The 60aa-dup is located in the C-terminal half of the extracellular ectodomain (Figure 5) and contains many potential acceptor sites for O-linked sugars. The G protein of the classical A2b strains, which lack 60aa-dup, have 58–63 potential acceptor sites for O-linked sugars. With 60aa-dup, the G proteins of the A2b_180nt-dup HMPV strains have acquired 23–25 additional potential acceptor sites for O-linked sugars. However, 60aa-dup does not affect the number of potential acceptor sites for N-linked sugars.