We collected 1,064 MuV sequences obtained from 644 clinical samples and 420 MuV isolates. These samples were obtained from 19 local governmental health institutes nationwide and five local base hospitals in 24 prefectures during 1986 to 2017. The clinical samples (including 362 buccal swabs, 433 throat swabs, 119 spinal fluids) were collected from patients with clinically suspected mumps in the hospitals by doctors and sent to the local government laboratories or the National Institute of Infectious Diseases (NIID) and stored at −70°C or lower in each laboratory until they have served for analyses. The case definition for mumps includes clinical symptoms of rapid swelling of uni- or bilateral salivary glands with pain. A confirmed case was defined by the clinical symptoms plus a laboratory confirmation by detecting MuV RNA by RT-PCR from clinical samples. Among the 1,064 specimens collected, 623 were sequenced in the local government laboratories. The remaining 441 samples that were collected in the hospitals were sequenced at NIID. The phylogenetic analyses were conducted at NIID and National Center for Global Health and Medicine.

The MuV isolation and propagation were performed as previously reported (Kidokoro et al., 2011). Briefly, the supernatants of centrifuged clinical samples which were soaked in Eagle’s Minimum Essential Medium (MEM), which were supplemented with 3% fetal bovine serum, 200 μg/m streptomycin, and 200 U/ml penicillin, were inoculated onto Vero cells. The cells were cultured at 37°C and the supernatants were harvested when those cytopathic effects (CPE) became prominent. Two blind passages were performed on all CPE negative tissue cultures.

Viral RNA was extracted from clinical samples or the culture fluids of mumps isolates by using the QIAamp Viral RNA Mini Kit (Qiagen, Tokyo, Japan) following the manufacturer’s instructions. The SH gene (the genome positions including from 6,076 to 6,803) was amplified with a set of primers SH-1 (5′-GCRACYAAAGARATCAGRAGRATC-3′) and SH-2 (5′-AGCCTTGATCATTGATCATCC-3′), using PrimeScript™ One Step RT-PCR Kit Ver.2 (TaKaRa, Kyoto, Japan) according to the manufacturer’s instructions. When the RT-PCR products were not detected, we performed a second PCR using nested primer set designed to amplify the genome positions from 6,130 to 6,689 with a set of SH-3 (5′-TCAAGYAGTGTCGAYGATCTC-3′) and SH-4 N (5′-AGCCTTGATCATTGATCATCC-3′) primers. The PCR products were sequenced using the BigDye Terminator v3.1 Cycle Sequencing Kit and 3130xl Genetic Analyzer (Applied Biosystems Japan, Tokyo, Japan).

We applied the viral RNAs extracted from culture fluids of 82 MuV isolates of the current study (of these, 77 were genotype G, 5 were genotype L) to the NEBNext Ultra RNA Library Prep Kit for Illumina, followed by the NEBNext rRNA Depletion Kit (Human/Mouse/Rat; NEW ENGLAND BioLabs Japan, Inc. Tokyo, Japan) according to the manufacturer’s instructions. Thereafter, sequences were generated using paired-end sequencing on the MiSeq platform (Illumina K.K. Tokyo, Japan). The 77 strains of genotype G isolates were constituted by 21 strains which represented the subclades based on SH sequence, and the phylogenetically related 56 isolates. Of those, 31 isolates share the identical SH gene sequence. Five of the 6 genotype L strains also share the identical SH gene sequence. We used CLC Genomics Workbench version 8.0.1(CLC Bio, Qiagen, Germany) for assembly, reference mapping, and generation of consensus sequences of mumps isolates. The sequences were deposited in the DNA Data Bank of Japan (DDBJ) and the accession numbers are provided in Supplementary Table S1.

A total of 316 nts of the sequences containing the entire SH gene (the genome positions from 6,218 to 6,533) were used for genotyping. To verify the results of SH gene-based analysis for genotype G, we performed phylogenetic analyses by using whole-genome sequencing (WGS) through whole-genome sequence datasets (n = 101), excluding both terminal sequences (WGS-dt, 15,292 nts between genomic positions 56 and 15,347). In order to compare our genotype G strains with as many SH gene and WGS data as possible, we searched for the National Center for Biotechnology Information (NCBI) Genbank and retrieved 7,824 and 177 sequences of SH gene and WSG, respectively, and selected nonredundant sequences (884 strains for SH gene and 30 strains for WGS-dt) based on their sequences, isolated areas and isolated years. Regarding genotype L, we adopted the dataset of HN gene sequences for the analyses instead of the WGS sequences, because WGS data of genotype L were not available in NCBI except only one strain. Then we aligned them with the sequence data of reference strains proposed by the World Health Organization (WHO; WHO, 2012) using MUSCLE (Edgar, 2004). Phylogenetic trees of SH gene, WGS-dt and HN gene sequences were inferred by the maximum-likelihood (ML) method using the Tamura-Nei model with the gamma-distributed site (TN93 + G) for SH gene and general time reversible model with gamma-distributed site (GTR + G) for WGS-dt and Tamura 3-parameter model (T92) for HN gene as the nucleotide substitution model. The topology of the phylogenetic tree was validated by Felsenstein’s bootstrap test with 1,000 resamplings, and > =75% of the probability was considered as a significant monophyletic group or clade. The multiple alignment, ML tree inference and the substitution model selection were conducted using Molecular Evolutionary Genetics Analysis (MEGA) software version 10.1.1 (Kumar et al., 2016). To confirm the evolutionary process of MuV transmissions in Japan along with time, WGS-dt of genotype G and HN gene of genotype L were also analyzed by Bayesian Markov Chain Monte Carlo (MCMC) coalescence method using BEAST 1.9 or BEAST 2.4.7. We used the GTR + G substitution model for the WGS-dt dataset and the HKY substitution model for the HN gene dataset. Models for the clock and population were evaluated in the actual data by log marginal likelihood estimation using path sampling and stepping-stone sampling (Baele et al., 2012), and adopted a combination of relaxed clock gamma distribution model and Bayesian Skyline Plot (BSP) for the WGS-dt, and a combination of relaxed clock log normal model and Coalescent Exponential Population for the HN gene. To ensure convergence at a minimum effective sample size (ESS) of >200, MCMC chain set in 300 million and were sampled at every 30,000 generations, and removing 300 million steps as burn-in. The MCMC output was analyzed using Tracer v.1.6.0 software¹ and the evolutionary parameters, mean times of the most recent common ancestor (tMRCAs) of the clades identified by ML tree, and the BSP-based population size changes in time were estimated for a most major clade of genotype G, JPC-1. Bayesian maximum clade credibility (BMCC) chronological tree inference was constructed using TreeAnnotator v.2.4.7 software and visualized in DensiTree v.2.0 and FigTree v.1.4.4. The list of the representative MuV isolates in this study and their accession numbers are provided in Supplementary Tables S1, S2.

Multiple alignments of the deduced amino acid sequences for the mumps virus proteins, N, V, P, M, F, SH, HN and L, of the 77 isolates of Japanese genotype G (JPC-1, 64 strains; JPC-2, 2 strains; JPC-3, 8 strains; JPC-4, 1 strain; JPC-5, 2 strains) were constructed by using GENETYX-MAC software v16.0.7 (GENETYX CORPORATION, Japan). The alignments were constructed from nonredundant amino acid sequences from each JPC. Excepting SH protein, two current Japanese mumps vaccine strains, Hoshino (AB470486) and Torii (MD217222) were aligned for reference sequences. In regard to SH protein, five additional isolates in our study (MuVi/Mizushima#06.JPN/43.06[G], MuVi/Mie#B11.JPN/34.11[G], MuVs/Kumamoto#175.JPN/53.12[G], MuVi/Mie#D10.JPN/12.14[G], MuVi/Kanagawa259.JPN/27.15[G]), that harbored truncated forms of SH protein, were aligned with 77 isolates.

The age distribution of the 827 cases, whose epidemiological backgrounds are evident, range from 0 to 49 years old (median 6.0-year-old [IQR, 4.0–8.0 year-old]; Figure 2A). Most of the patients (89.8%) were children aged 10 years or younger. Unexpectedly, although the distribution also showed a minor peak during 11 and 15 years old, the proportion of the adolescents and adults aged 16 to 49 years was relatively low (2.5%). One of the reasons for this is that most specimens might have been collected from pediatric clinics. Among the 672 patients whose symptoms were clear, the most common symptoms were parotitis (499 cases, 74.3%) and aseptic meningitis (166 cases, 24.7%). There were five cases (0.7%) of encephalitis or encephalopathy. In this study, orchitis was reported in only two cases (0.3%) as the vast majority of patients were children under 10 years of age. Although the mumps vaccine histories of over 80% of the patients were not described (887 cases, 84.2%), the previously vaccinated or unvaccinated cases were 96 (9.1%) and 81 (7.7%), respectively. The male-to-female ratio of patients was 58.2% versus 41.8%. These results show similar patterns with the data on age- and gender-specific incidences of mumps in Japan, which were from the Infectious Agents Surveillance Report, whose data were collected from approximately 3,000 sentinel pediatric clinics nationwide.² The time intervals between the onset of mumps and sampling of the specimens ranged from 0 to 16 days (median 1.0 d [IQR, 0.0–2.0 d]; Figure 2B). Notably, >75% MuV-positive specimens (454 [76.3%]) were sampled within 2 d of onset. These results highlight that mumps clinical samples should be collected within 2 d of mumps symptom onset, and support previous findings (Rota et al., 2013; Golwalkar et al., 2018; Shah et al., 2018).

Through phylogenetic analyses based on 1,064 SH gene sequences, six genotypes, namely B, F, G, H, J, L, (110 [10.3%], 1 [0.1%], 900 [84.6%], 3 [0.3%], 41 [3.9%], and 9 strains [0.8%], respectively) were identified (Figure 3). The 353 selected representative MuV strains are listed in Supplementary Table S1. Chronological changes and geographical distribution of genotypes identified in the study are shown in Figures 1, 4. In the 1980s, genotype B was in exclusive circulation in Japan (Figure 1A,B) because it is a classic Japanese indigenous genotype, as with all Japanese vaccine strains, whose parental strains were isolated in the late 1960s or early 1970s (Figure 5A). Therefore, the genotype B isolates in the 1980s and 90s were classified into distinct clades against Japanese vaccine strains, except for Urabe strains (Figure 5A). Genotype B viruses continued to circulate until 2001 (Figure 1A,B).

Genotype J first emerged in 1993 and co-circulated in Japan with genotype B until the early 2000s (Figure 1A,B) in Japan. Genotype J strains have also been reported in several other countries, including Singapore, United Kingdom, Taiwan, Malaysia, Spain, Ireland and Thailand (Jin et al., 2015). However, the Japanese strains are genetically distinct from the foreign strains, whose inter-clade divergences ranged from 3.5 to 7.6% (Figure 5B). Additionally, the Japanese strains are highly homologous to each other, with a maximum intra-clade divergence of 2.8%, compared to that of the foreign strains (4.7%).

Genotype L circulated in Japan from 1999 to 2002 and disappeared after a short period (Figure 1A,B). However, the genotype L viruses caused nationwide epidemics, because they were isolated from various parts of Japan, such as Tokyo, Mie, Okayama, Fukuoka prefectures and, Yokohama city in Kanagawa prefecture (Figure 5C), even though only two prefectures, namely Mie and Okayama, were registered in the study (Figure 4). The sequences formed a single clade with high homology, whose divergences ranged from 0.0 to 1.3%, and were evidently divergent from Netherlands strains with the high bootstrap values (Figure 5C). To confirm the result of ML phylogenetic analysis based on the SH genes, we performed time-scaled Bayesian evolutionary analyses along with ML phylogenetic analysis based on the datasets of 11 HN gene sequences (6 Japanese strain and 5 reference strains) of genotype L strains (Supplementary Figures S1A,B). The resulting trees described that genotype L strains were divided into two phylogenetically distant clades, Japanese clade and Netherlands clade with high bootstrap values (100%), which showed similar tree topologies with that of the SH based ML tree (Supplementary Figures S1A,B). We also estimated mean times to the most common recent ancestor (tMCRA) of these two clades, which were 1996.70 (95% highest posterior density (HPD): 1999.10–1993.38) for Japanese clade, and 1955.92 (95% HPD: 1957.50–1951.91) for Netherlands clade, respectively (Supplementary Table S3). The tMCRA of all genotype L strains was 1954.83 (95%HPD: 1957.49–1944.41). These data suggest that the Japanese genotype L viruses could not be directly related with European genotype L strains.

In our study, Japanese genotype H viruses were identified only three strains, one of which was isolated from Kanagawa Prefecture in 2000 (MuVi/Kanagawa.JPN/45.00). The others were isolated from Nara Prefecture from 2016 to 2017 (MuVs/Nara.JPN/22.16, MuVs/Nara.JPN/21.17; Figures 1A,B, 4, 5D). These isolates were classified into two distinct clades (Figure 5D). Two isolates from Nara Prefecture that shared an identical sequence were classified into a clade, including isolates from Asian countries, such as Korea, Mongolia, and Japan, despite the fact that the most closely related strain was the isolate from Sweden (KF840222; Figure 5D). Almost all Japanese genotype H isolates were classified in this clade. Moreover, MuVi/Kanagawa.JPN/45.00 was classified into distinct clades, which primarily included strains isolated from the Western and the Middle Eastern countries, including United Kingdom, Spain, United States, Turkey, and Iran. This is the only isolate of the clade in Japan thus far.

The singular isolate of genotype F in this study was isolated from Aichi Prefecture in 2016 (MuVi/Aichi.JPN/24.16), whose sequence was identical to that of the Chinese isolate in 2014 (KY680444, Figure 5E). The results suggest that this could be an imported case from China, even though epidemiologic links were unclear.

Genotype G was the most frequently isolated genotype in our study (900 isolates out of 1,064 isolates, 84.6%) as it was predominant in the last two decades in Japan, and our surveillance network began to work at full scale during this period. Although previous studies have reported that the Japanese genotype G consists of two major clades (Akiyoshi and Suga, 2014; Aoki et al., 2016; Kuba et al., 2017), we identified at least seven clades in this study using SH-based classification on the basis of the topology of the phylogenetic tree, designated Csh-1, -2, -3, -4, -5, -6, and -7 (Figure 6). However, as most bootstrap values of each clade were <75%, the topology of the phylogenetic tree could be inaccurate. In effect, topologies of phylogenetic trees varied depending on strain combinations of genotype G (Figures 3, 6).

To establish the phylogenetic relationships among the seven clades based on SH sequences, we performed WGS using NGS on 77 strains of genotype G isolates, which contained 21 strains representing each clade by SH (Figure 6) and 56 related strains isolated in the study. They were analyzed by using both ML phylogenetic method and Bayesian MCMC method based on WGS data excluding leader and trailer sequences (WGS-dt, 15,292 nts between genomic positions 56 and 15,292; Figure 7; Supplementary Figure S2; Supplementary Table S2). Through the MCMC analyses, we also estimated the tMCRA of each clades of Japanese genotype G (Supplementary Table S3).

The bootstrap values in ML tree of the WGS-dt were much higher (> = 97%) than those obtained in SH gene detaset. The posterior probabilities of the Bayesian tree were also sufficiently high (> = 0.99). The both resulting trees inferred from the WGS-dt dataset exhibited similar topologies, which showed that Japanese genotype G strains were divided into two major clades and each clade was further divided into two and three minor clades, respectively. Consequently, the seven clades based on SH sequences were consolidated into five clades, which were designated Japanese clade (JPC)-1, -2, -3, -4, and -5, numbered in order of emerging times of the first isolate of each clade for convenience. As a result, the three SH-based clades, Csh-1 to -3, were classified as JPC-1 (Figure 7; Supplementary Figure S2).

In order to clarify the phylogenetic relationships between the JPCs and globally distributed genotype G strains, we constructed phylogenetic tree from the dataset of SH sequences from 867 genotype G sequences retrieved from NCBI and 30 representative strains this study, 28 of which were from WGS-dt data (Supplementary Figure S3).

The JPC-1, tentatively called “Gw,” was the first identified and the most abundant clade in Japan (726 isolates, 80.7% of genotype G). The first identified JPC-1 strain was isolated from Ehime Prefecture in 1997 (MuVi/Ehime#1342.JPN/0.97). The clade subsequently became the great majority of the MuV isolates during the next 5 years. The tMCRA of JPC-1 was 1992.84 (95% HPD: 1996.26–1988.61; Supplementary Figure S2). Some strains, which isolated in foreign regions, Mainland China, Taiwan, North America, were also classified into the JPC-1, all of which were isolated after the tMCRA of JPC-1 (Supplementary Figure S3; Supplementary Table S3). Of these, three Taiwan isolates (MF281016, MF281017, MF281023) were isolated from a Japanese, a Thailander, and a Korean, respectively (Cheng and Liu, 2018). Two isolates of United States (FJ959106, FJ959107) were isolated from outbreaks occurred at the University of Virginia (UVA) in late 2006, and in New Hampshire in June of 2006, one of which, FJ959107, was a imported case from Japan (Rota et al., 2009). The foreign JPC-1-related strains were phylogenetically distant from isolates of the 2006 large outbreak in Iowa, USA (JN012242), and related outbreaks in USA, Canada and Europe (Figure 7; Supplementary Figures S2, S3; Watson-Creed et al., 2006; Centers for Disease and Prevention, 2006a; Rota et al., 2009; Stapleton et al., 2019; Bodewes et al., 2020).

We also estimated the BSP-based effective population sizes for the JPC-1 viruses (Supplementary Figure S4; Supplementary Table S3). The results of the BSP analysis showed that the JPC-1 viruses have rapidly increased their population in the end of 1990s when right after the first JPC-1 isolate was detected, then have been gradually increasing until 2017.

JPC-2 was the second clade that emerged in Japan in 2001. In this study, it was identified in only two cases from Aichi and Mie prefectures (MuVi/Aichi#10879.JPN/6.01[G] and MuVi/Mie#W12.JPN/20.01[G]), the central region of Japan. However, as closely related strains were isolated in Hokkaido, the northernmost prefecture in 2000 or 2001 (AB105477, KF876710), and one of these, AB105477, shares a 100% identical sequence to MuVi/Aichi#10879.JPN/6.01[G], JPC-2 may have caused a nationwide epidemic in Japan within 2 years. The most closely related strain was isolated from a Burmese with recent travel to Myanmar, in Taiwan in 2007(MF281013; Supplementary Figure S3; Cheng and Liu, 2018). The tMCRA of JPC-2 was 1999.35 (95% HPD:2000.91–1996.41; Supplementary Figure S2; Supplementary Table S3).

JPC-3, tentatively called “Ge,” was the third-emerged and second major clade of Japanese genotype G (144 isolates, 16.0% of genotype G), which has been co-circulating with JPC-1 for a decade. JPC-3 was first isolated from Niigata Prefecture at the end of 2009 (MuVi/Niigata#1604.JPN/49.09[G]). In some regions and periods, JPC-3 viruses predominantly circulated, such as Shizuoka City in 2010–2014, Aichi Prefecture in 2013–2016, and Okinawa main island in 2015–2016 (Kuba et al., 2017; Figure 8). In contrast, in some areas such as Osaka, Ibaraki, and Gifu prefectures, the two clades co-circulated in the same area during the same period. Closely related strains of the JPC-3 were also reported from Canada, United States, Taiwan, and China (MN910022, MN910184, JX455939, KF143753, MF281022, KF031051; Supplementary Figure S3). The tMCRA of JPC-3 was 2007.97(95% HPD: 2009.58–2006.07; Supplementary Figure S2; Supplementary Table S3).

The JPC-4 was singly isolated in Okinawa Prefecture in 2015 (LC167469; Kuba et al., 2017), and it was tentatively called “Ghk” because the isolate was phylogenetically most closely related to the Hong Kong isolate (KM234040) in 2014. In this study, two strains of JPC-4 were isolated successively in Kitakyushu City in 2016 (MuVi/Kitakyushu#062.JPN/31.16[G] and MuVi/Kitakyushu#082.JPN/34.16[G]), whose SH gene sequences shared 100% sequence similarity with that of the Okinawa isolate. Closely related foreign strains were isolated in the US and Canada along with China in 2014–2016 (KY013252, MN911765, MN910077, KM234040; Supplementary Figure S3). Of these, KY013252, the isolate from California, USA in 2016, was most homologous with JPC-4 strain (99.7% similarity). The remaining 3 strains, which shared 100% identical sequence, have 99.4% identity to JPC-4 strain. The tMCRA of JPC-4 could not be estimated because of the single WGS-dt sequence.

The JPC-5 was the latest clade of the Japanese genotype G (26 isolates, 2.9% of genotype G), which was first identified in 2016 in Mie, Gifu, and Aichi prefectures (MuVi/Mie#02.JPN/7.16[G], MuVs/Gifu#147.JPN/17.16[G], and LC424348). In 2017, in addition to Gifu Prefecture, the closely related isolates were identified in Ibaraki Prefecture (MuVs/Ibaraki.JPN/15.17[G] and MuVs/Ibaraki.JPN/25.17[G]) and Tokyo (LC325485; Figure 8). The Tokyo isolate shared a 100% identical SH gene sequence with Mie, Gifu, and Aichi isolates. Therefore, JPC-5 might have circulated over an extensive area of Japan in 2 years. Closely related strains were isolated in Canada and United States in 2015–2017 (MN926205, MF138973, KU508690, MH618522, MH618544; L’Huillier et al., 2018). Of these, a Canadian isolate (MN926205) shared the completely identical sequence with three Japanese isolates, MuVi/Mie 02.JPN/7.16[G], MuVs/Aichi.JPN/16.16[G], and LC325485 (Supplementary Figure S3). The tMCRA of JPC-5 was 2014.02(95% HPD: 2015.65–2011.83; Supplementary Figure S2; Supplementary Table S3).

In order to clarify the diversities of amino acid sequences within Japanese genotype G strains, the multiple alignments of all MuV proteins were constructed from the WGS sequences of 77 Japanese genotype G isolates (JPC-1, 64 strains; JPC-2, 2 strains; JPC-3, 8 strains; JPC-4, 1 strain; JPC-5, 2 strains; Supplementary Figure S5; Supplementary Table S4).

N protein is a major component of helical viral nucleocapsid and plays an important role in translation, transcription, and budding, where the characteristic “546DWD548” motif near the C-terminal of N protein plays an essential role for interaction with M protein and viral budding (Ray et al., 2016). N protein is also known as a target of the human CD8⁺ T cells (de Wit et al., 2019, 2020). The alignment of N protein shows that the DWD motif and human T cell epitope, N_115–122, IPNARANL, are conserved among all Japanese genotype G strains along with three vaccine strains (Supplementary Figure S5; Supplementary Table S4).

The V protein blocks interferon (IFN) production and signaling pathway via its 69aa C-terminal cysteine-rich domain (CTD), which interacts with the host STAT (signal transducer and activator of transcription) molecules (Nishio et al., 2002), or the MDA5 (melanoma differentiation-associated gene 5; Andrejeva et al., 2004; Ramachandran and Horvath, 2010), the key molecules in IFN signal transduction pathway. This characteristic cysteine-rich motif, 170H-C-C-C-C-C-C217, was conserved in all strains along with three vaccine strains (Supplementary Figure S5).

The V/P gene gives rise to additional mRNAs by co-transcriptional addition of two non-templated G residues in the specific positions of mRNAs (referred to as RNA editing) resulted to generate the P proteins (Paterson and Lamb, 1990). Therefore, although V and P proteins share a common N-terminal 155-aa sequences, their C-terminal regions are different. The N-terminal common regions were known to harbor 5 experimentally characterized human T cell epitopes, V/P_17-25, V/P _111–122, V/P _142–150, V/P _142–152, V/P _144–152 (Supplementary Figure S5; Supplementary Table S4). Of these, V/P_17-25, which was one of the most broadly recognized by some human haplotypes, was conserved in all Japanese genotype G strains along with two Japanese vaccine strains, except for Jeryl Lynn vaccine strain (Supplementary Table S4). The remaining 4 epitopes were 100% conserved in JPC-2, -3, -4 and two Japanese vaccine strains, however, partially conserved in JPC-1 (84–97%) and JPC-5 (0 or 100%). The T cell epitope in V protein specific C-terminal region, V_194–202, was conserved in the Japanese vaccine strains and JPCs, excepting JPC-1 (75%).

The P specific C-terminal region contains two T cell epitopes, P_235–243, P_382–391. P_235–243 was 100% conserved in all Japanese genotype G strains along with all vaccine strains. P _382–391 was conserved in the Japanese vaccine strains and JPCs, excepting JPC-1 (97%) (Supplementary Figure S5; Supplementary Table S4).

The M protein plays an essential role in the virion assembly and budding, by means of associating with the inner surface of the MuV envelope and the cytoplasmic tails of the envelope proteins, as well as the N proteins in the RNP complex. The sequence 24FPVI27 near the amino-terminal end of M protein plays an essential role for proper budding of MuV (Li et al., 2009). In the alignment of M protein, this motif was conserved in all 3 vaccine strains and JPC-3, JPC-4, JPC-5 strains (Supplementary Figure S5; Supplementary Table S4). However, interestingly, in three strains of JPC-1 and all 2 strains of JPC-2, this motif sequence was disrupted by V26I or I27V mutations. It is not clear how these mutations influence on budding process and pathogenicity of these mutant viruses. However, since physico- chemical properties of valine and isoleucine are similar, the influence of the mutations on viral replication could be restricted. M protein contains 10 T cell epitopes, M_108-116, M_120-128, M_161-169, M_187-195, M_303-312, M_309-316, M_337-345, M_350-358, M_364-372, M_365-372, all of which were conserved in the Japanese vaccine strains and isolates of all JPCs excepting for JPC-1. In JPC-1 strains, 8 of 10 epitopes, excepting M_108-116, M_350-358 were conserved. Meanwhile, in Jeryl Lynn stain, 5 of 10 epitopes, M_120-128, M_161-169, M_337-345, M_350-358, M_365-372, were not conserved.

F protein, one of two envelope glycoproteins of MuV, is a target antigen of neutralization antibodies. The seven potential glycosylation sites and proposed three B-cell epitopes, 221I, 323 N, 373D (Santak et al., 2015) were conserved in all strains (Supplementary Figure S5; Supplementary Table S4). The polybasic cleavage site (99RxKR102, where x is any amino acid) recognized by cellular protease were conserved in JPC-1, -2, -5, and vaccine strains. Interestingly, in all strains of JPC-3 and -4, the basic amino acid, lysine (K) at position 101 were substituted to an analogous basic amino acid arginine (R) like as other paramyxoviruses, SV5 and RSV (Klenk and Garten, 1994). F protein contains 4 experimentally characterized T-cell epitopes, F_112-121, F_152-161, F_253-262, F_445-453, of which, all epitopes were conserved in two Japanese vaccine strains and JPCs excepting for JPC-1, even though its conserved ratios were high (84–100%; Supplementary Table S4). In Jeryl Lynn strain, a half of the epitopes were not identical.

HN protein, another envelope glycoprotein of MuV, plays a key role in the viral entry to the host cells via binding to the cellular receptor molecules, sialic acids (Kubota et al., 2016). After that, HN protein promotes a conformational change of F protein, resulted in a membrane fusion (Bose et al., 2014). The16 amino acid residues associated with receptor binding and 3 residues associated with fusion promotion by HN protein were conserved in all strains (Supplementary Figure S5). The eight potential glycosylation sites on HN ectodomain were all conserved. HN protein is the most important target antigen for neutralization antibodies, and at least three B-cell epitopes, HN_265–288, HN_329–340, HN_352–360, were experimentally characterized (Orvell, 1984; Cusi et al., 2001; Kulkarni-Kale et al., 2007). The averages of the sequence homologies of the epitopes between Japanese vaccine strain Torii and JPC isolates were relatively higher (96–100%) than that between Jeryl Lynn strain and JPC isolates (88–92%). Notably, 4 experimentally identified T-cell epitopes in genotype G HN protein, HN_88-96, HN_157-165, HN_352–360, HN_505-513, were conserved in two Japanese vaccine strains and also highly conserved in JPCs (98–100%), while 2 of 3 epitopes were not conserved in Jeryl Lynn strain (Supplementary Figure S5; Supplementary Table S4).

L protein is the largest protein in the MuV viral components which has multifunctional enzymatic activities, in which 9 T-cell epitopes, L_249-256, L_336-344, L_441-449, L_591-600, L_740-748, L_1299-1307, L_1336-1343, L_1806-1813, L_1983-1992, L_1986-1994, were experimentally identified (de Wit et al., 2019; Kaaijk et al., 2021). All of them were conserved in two Japanese vaccine strains and almost all of JPC isolates (98–100%; Supplementary Figure S5; Supplementary Table S4). Meanwhile, 5 of 9 epitopes were not conserved in Jeryl Lynn strain.

In the study, truncated forms of mutant SH protein were identified in 7 strains, which coded three different length of SH proteins, 8-, 44-, and 50-aa (Supplementary Figure S5). These mutant SH proteins were generated through three types of nucleotide substitutions, C to A, C to T, or G to A, in different sites that generated new internal stop codons. These mutant viruses were all classified into JPC1, even though their phylogenetic distance and their isolation years (from 2006 to 2015) were not necessarily close (Figure 6; Supplementary Figure S5).