Genotype replacement of the human parainfluenza virus type 2 in Croatia between 2011 and 2017 – the role of neutralising antibodies

Abstract

Previously we reported on the HPIV2 genotype distribution in Croatia 2011–2014. Here we expand this period up to 2017 and confirm that G1a genotype has replaced G3 genotype from the period 2011–2014. Our hypothesis was that the G1a-to-G3 genotype replacement is an antibody-driven event. A cross-neutralisation with anti-HPIV2 sera specific for either G1a or G3 genotype revealed the presence of genotype-specific antigenic determinants. By the profound, in silico analyses three potential B cell epitopic regions were identified in the hemagglutinin neuraminidase (regions 314–361 and 474–490) and fusion protein (region 440–484). The region identified in the fusion protein does not show any unique site between the G1a and G3 isolates, five differentially glycosylated sites in the G1a and G3 genotype isolates were identified in epitopic regions of hemagglutinin neuraminidase. All positively selected codons were found to be located either in the region 314–316 or in the region 474–490 what indicates a strong positive selection in this region and reveals that these regions are susceptible to evolutionary pressure possibly caused by antibodies what gives a strong verification to our hypothesis that neutralising antibodies are a key determinant in the inherently complex adaptive evolution of HPIV2 in the region.

Introduction

Human parainfluenza virus type 2 (HPIV2) is a member of the genus Rubulavirus in the Paramyxoviridae family. The HPIV2 genome is a non-segmented single-stranded negative-strand RNA that contains 15 654 nucleotides. It encodes six genes for seven proteins: the nucleoprotein (NP gene), V protein/phosphoprotein (V/P gene), matrix protein (M gene), fusion protein (F gene), hemagglutinin neuraminidase (HN gene) and large protein (L gene). The major targets of the anti-HPIV2 antibodies are surface glycoproteins hemagglutinin neuraminidase and fusion protein [1, 2]. Hemagglutinin neuraminidase is anchored in the viral envelope as a tetramer (dimer of dimers) together with the fusion protein. The hemagglutinin neuraminidase has a triple role in viral infection: initial attachment to a sialic acid-containing receptor, cleaving the receptor and activating the F protein. The F protein is synthesised as a biologically inactive F₀ precursor from which a prefusion metastable form of disulphide-linked F₂ + F₁ is formed by the proteolityc action of the cellular endoprotease (reviewed in [3]). Binding of the HN protein to the receptor initiates the F protein to undergo conformational changes which ultimately lead to the fusion of the viral and cellular membranes [4, 5] what finally leads to the intrusion of the ribonucleocapsid into the target cell.

The dynamics and evolution of HPIV2 are largely unknown. Most HPIV2 infections are treated symptomatically outside of the hospital setting. Thus the diagnosis is missing making it difficult to follow the genetic variation and evolution of the HPIV2. There are very few studies on the evolution and phylogenetic analysis of this virus [6–8]. These phylogenetic analyses identified four genotypes (G1-4). Our previous study [8] explored hospitalised HPIV2 cases in Zagreb, Croatia during 4-year period (2011–2014) and showed that the dominant genotype in this area was G3 genotype while only one strain isolated in 2014 was positioned in the distant G1a genotype. Therefore now we extended our monitoring of the HPIV2 up to 2017 in order to determine whether G3 genotype remained dominant or it was overruled by G1a genotype.

Respiratory diseases are very often the cause of hospitalisations worldwide. The major cause is infection with a respiratory syncytial virus (RSV), while human parainfluenza viruses (HPIVs) infections are the second major cause of hospitalisation for respiratory tract illness in young children [9]. HPIV are spread by respiratory secretions from infected persons or contact with contaminated surfaces. Infections with these viruses are associated with significant morbidity and at least as much acute respiratory illness (ARI) in infants and young children as [10–13]. HPIV2 causes ARI (croup, pneumonia, bronchitis and bronchiolitis) in infants, young children and immunocompromised elderlies [14] although with lower prevalence rate than HPIV3 and HPIV1. Infections with HPIV2 are also recognised less often because of its mild nature, especially, in otherwise healthy adults. For that reason, in Croatia, the diagnosis for HPIV2 is conducted only for hospitalised children or adults with immunocompromised conditions. Most HPIV2 infection have been detected in autumn or early winter at 1 or 2-year interval. The overall relative contribution to the burden of disease in infants and young children is estimated at approx. 7:3.5:2:1 for HPIV3, HPIV1, HPIV2 and HPIV4, respectively [10, 15–17]. Also, these viruses are proven to be associated with nosocomial ARI in immunocompromised patients [16, 18, 19]. In spite of the impact of these viruses on the public health no licensed vaccine or any other specific anti-viral treatment against any type of HPIV exists because of the incomplete knowledge of the immune response against these viruses. Almost all children encounter HPIVs within the first few years after birth, but immunity seems to be incomplete and re-infections occur throughout life. However, after primary infection, a certain level of protection develops which is able to prevent severe disease later in life but it is unable to completely suppress the infection [9]. The reason for this incomplete protection may be a genetic drift of the virus as we imply in this study.

In this study, we examine the role of neutralising antibodies as the promoter of the inter-seasonal dynamics of the HPIV2 and potentially as s primary cause of the incomplete immunity to re-infections with viruses of distant genotypes.

Materials and methods

Isolation of virus from clinical specimens, amplification of the F and HN genes and sequencing

Clinical specimens were obtained from 38 patients hospitalised for respiratory infections. Samples were collected during 2011–2017 as nasopharyngeal secretions (details in Table 1) and detected as HPIV2 positive by direct fluorescence assays (DFA Light Diagnostics, Chemicon International, Temecula, CA, USA). None of the samples was found to be a co-infection with another respiratory pathogen. This study was approved by the ethics committee of the University Hospital for Infectious Diseases ‘Dr Fran Mihaljević’.

Table 1. HPIV2 isolates collected from patients during 2011–2017 in Zagreb, Croatia

^a Isolate number is given in parentheses.

^b m – month; y – year.

^c P-pneumonia, B-bronchiolitis, C-croup, URI-upper respiratory tract infection.

^d Accession number of the sequence of the F and HN genes submitted to the GenBank.

The procedure to obtain the sequences of the F and HN genes consisted of the virus isolation from each specimen on Vero cells in the presence of trypsin in order to decrease the number of mutations during cultivation in the cell culture as shown in [8]. A number of isolated viruses vs. obtained clinical specimen was: 5/5 (2011), 9/11 (2012), 2/6 (2014), 2/4 (2015), 17/18 (2016) and 3/3 (2017). Isolation of viral RNA, amplification of the F and HN genes and Sanger sequencing of amplified segments. All procedures have been previously described in [8].

Phylogenetic analyses of the F and HN genes

Accession numbers of sequences obtained from the GenBank which were included in analysis are DQ072586-DQ 072589, JF912194-JF912196, NC_003443, AF533011, D00865, AF533012, AB176531, AB189948-AB189953, KT315655-KT315667 and KT595232-KT595234.

Alignments were performed by CloneManager Suite software (Scientific&Educational Software, USA).

For phylogenetic analyses, the coding region of the HN gene sequence was used. Multiple sequence alignments were made using the CLSUTAL W program integrated in software MEGA6.0 (http://www.megasoftware.net). The best evolutionary model for phylogenetic inference was determined by using Model selection implemented in MEGA6.0. The phylogenetic tree was generated in MEGA6.0 using maximum likelihood (ML) analysis with the Kimura 2-parametar substitution model and 1000 bootstrap replicates.

Immunisation of guinea pigs with HPIV2

Two strains of HPIV2 were used for immunisation of guinea pigs: HPIV2i/Zagreb.HR/42.14(1472) (acc.nr. KT315666) and HPIV2i/Zagreb.HR/47.14(1688) (acc.nr. KT315667). HPIV2i/Zagreb.HR/42.14(1472) has been previously shown to belong to genotype G3 and HPIV2i/Zagreb.HR/47.14(1688) belongs to genotype G1a [8].

All animal work was in accordance to Croatian Low on Animal Welfare (2017) which strictly complies with EC Directive (2010/63/EU). Female guinea pigs weighing 300–400 g (Institute of Immunology Inc., Zagreb, Croatia) were housed in a 12-h light/12-h dark cycle and a constant temperature of 22 °C. A standard guinea pig diets (Muchedola srl, Milano, Italy) and water were supplied ad libitum during all immunisation period.

Antigen for immunisation was prepared by growing the two HPIV2 strains on Vero cells, centrifugation of the cell culture supernatants at low speed for 10 min and then by ultracentrifugation at 28 000 rpm in SW28 rotor (Beckman Coulter, Indianapolis, USA) for 2 h. Finally, the pellet was resuspended in PBS, pH 7.4. Six animals were immunised with each HPIV2 strain while the mock group consisted of three animals received PBS only. Animals were given three s.c. immunisations (days 0, 18 and 36) each containing 4.3 logPFU. Fifty days following prime immunisation guinea pigs’ blood was collected by heart puncture and sera were prepared for further serological analyses by centrifugation at 1200g at 4 °C for 30 min and decomplementation by heating at 56 °C for 30 min.

Plaque reduction neutralisation test (PRNT)

PRNT was performed as in [20], with some modifications. Test sera were prepared in six fourfold dilutions starting from at least one dilution completely neutralising added virus and ending with at least one dilution having no neutralising capacity. The virus concentration was adjusted to 35 ± 10 PFU/0.1 ml with MEM-H with 2.5% FCS. The test was performed in 24-well plates in which Vero cells were seeded at concentration 3.5 × 10⁵ cells/well 24 h prior to the addition of mixture serum/virus. Each dilution of the serum was mixed with an equal volume of virus and incubated at 35 °C for 1 h. Then the medium was removed from the plate and 0.2 ml of the mixture was transferred in four wells in the plate. Each assay run included virus as a positive control and medium used for diluting the sera and viruses as a negative control. The plates were incubated at 35 °C for 1.5 h, the mixture serum/virus was removed and the cells were washed twice with medium without serum. Finally, cells were overlayed with 0.75 ml of a semisolid medium consisting of 1 v/v 2 × MEM-H without phenol red, 1 v/v 1.5% Noble agar (Sigma) and 0.75 µg/ml acetylated trypsin (Sigma). Plates were incubated at 35 °C in a humidified atmosphere of 5% CO₂. After six days the monolayers were stained with 0.5 ml of 0.05% neutral red (Sigma) and plaques were counted. Fifty percent end point titres were calculated using the Kärber formula: log₁₀ ND₅₀ = −max−Δ(∑p−0.5), where max is the log₁₀ of the highest dilution of serum in the test, Δ is log₁₀ of the dilution factor (for fourfold dilutions is −0.6) and ∑p is total number of plaques divided by the number of plaques in positive control sample.

A two-tailed unpaired Students’ t test was used to calculate statistical significance.

In silico B cell epitope prediction

Protein structures of the fusion protein and hemagglutinin neuraminidase of HPIV2i/Zagreb.HR/42.14(1472) and HPIV2i/Zagreb.HR/47.14(1688) was modelled using a homology SWISS-model platform [21–23] (http://swissmodel.expasy.org) in order to create PDB structures.

Three different servers were used to predict conformational epitopes in the fusion protein and hemagglutinin neuraminidase: (a) DiscoTope 2.0 server [24], (b) ElliPro [25] and (c) BepiPred-2.0 [26]. DiscoTope 2.0 and ElliPro servers were used with default settings and BepiPred-2.0 was used at treshold 0.55.

Only epitopic regions which were predicted by all three servers were selected.

Analysis of glycosylation sites within the fusion protein and hemagglutinin neuraminidase

Complete fusion protein and hemagglutinin neuraminidase sequence were submitted for prediction of N-linked and O-linked glycosylation sites by using glycosylation predictor program which identifies pairwise patterns surrounding glycosylation sites and use an odds ratio to weight their propensity of association with modified residues [27] (http://comp.chem.nottingham.ac.uk/glyco/).

Natural selection analysis

To evaluate the selection pressure effect on each codon/amino site in the predicted epitopic regions of the hemagglutinin neuraminidase the HN ORFs of all strains involved in the phylogenetic analyses were analysed by computing ratio between non-synonymous and synonymous substitutions. The selection analyses consisted of four independent maximum likelihood methods from the Datamonkey server [28] (http://www.datamonkey.org/): a single-likelihood ancestor counting (SLAC), a fixed-effects likelihood (FEL), an internal branch fixed-effects likelihood (IFEL) and a random effects likelihood (REL), which estimate synonymous and nonsynonymous substitution rates at every codon. Sites were considered positively selected if they met the cut-off criteria of P-value of <0.1 for SLAC, FEL and IFEL, or Bayes factors >50 for REL.

Results

Distribution of HPIV2 genotypes in Croatia in the period 2011–2017

A total of 38 HPIV2 specimens were collected from hospitalised patients in Zagreb, Croatia, in period 2011–2017: 16 isolates from period 2011 to 2014 were previously described in [8] and the new 22 isolates have been collected and isolated in the period 2015–2017. The age distribution for each season (Fig. 1) shows that age groups were evenly distributed throughout the seasons with the exception of the season 2016 where the number of patients in the age group 1–6 years is larger than in the age group <1 year. Unfortunately, the data about possible reinfection of patients from the age group 1–6 years are unavailable so it was not possible to further explain if these patients had a primary infection or they were re-infected.

Fig. 1. Age distribution of the patients from whom the samples were collected. n.c., none collected.

Our previous study [8] showed that in 2014 one isolate of genotype G1a appeared while the other isolate from 2014 and earlier isolates belong to genotype G3. This urged further analysis in order to confirm the replacement of the G3 genotype with G1a genotype in this area. The phylogenetic analysis for all 38 samples was performed and the phylogenetic tree is presented in Figure 2. The analysis revealed that the vast majority of isolates collected after 2014 are placed in genotype G1a (19/22, Figs 2 and 3). However, the G3 genotype has been still sporadically circulating in the population but with much lower intensity (3/22, Figs 2 and 3) indicating that G1a genotype is been replacing G3 genotype.

Fig. 2. Maximum likelihood phylogenetic tree of the HN gene ORF of HPIV2 strains isolated in Croatia 2011–2017 and sequences accessible in the GenBank constructed by the MEGA6.0 with 1000 bootstrap replicates. Croatian isolates are presented with their isolate numbers and sequences retrieved from the GenBank are presented with their accession numbers. The year of isolation is indicated for each sequence. The genotypes are indicated by the brackets on the right side. Only bootstrap values over 80% are displayed at the branch nodes. The scale bar indicates 0.01 nt differences per site, over the indicated region.

Fig. 3. Number of isolates from Croatia 2011–2017 belonging to G1a or G3 genotype. n.s., no specimen obtained.

Genotype specific neutralisation potential of guinea pigs’ immunised with different HPIV2 strains

A possible reason of G1a genotype becoming dominant over the G3 genotype is an incomplete antigenic match between the two genotypes. To analyse this, an immunisation of guinea pigs with two HPIV2 strains was performed in order to obtain neutralising sera specific for different genotypes. Both viruses used for immunisation were isolated from a specimen collected in 2014 but belonged to different genotype: HPIV2i/Zagreb.HR/42.14(1472) (G3 genotype) and HPIV2i/Zagreb.HR/47.14(1688) (G1a genotype). The results of the PRNT with guinea pigs’ sera and the two viruses are shown in Figure 4. Sera of animals immunised with either virus had strong neutralising activity. However, sera specific for the genotype G3 1472 virus had significantly lower neutralisation capacity for the genotype G1a 1688 virus (PRTN₅₀ titres 11 935.5 ± 1 046.5 for 1472 virus vs. 8 654.7 ± 972.1 for 1688 virus, P < 0.01, Fig. 4) and vice versa, sera specific for the genotype G1a had significantly lower capacity to neutralise virus of the genotype G3 (PRTN₅₀ titres 8566.1 ± 925.3 for 1472 virus vs. 12 407.7 ± 1389.7 for 1688 virus, P < 0.01, Fig. 4). To confirm that this finding is not specific for these two viruses only but can be applied at the genotype level, two more viruses from both, G1a and G3, genotypes were neutralised with these sera. Figure 4 shows that the neutralisation capacity of the sera is comparable as long as the genotype of the virus used for immunisation matches the virus used in neutralisation test.

Fig. 4. Neutralisation of HPIV2i/Zagreb.HR/47.14(1688), HPIV2i/Zagreb.HR/03.17(208), HPIV2i/Zagreb.HR/42.16(1995), HPIV2i/Zagreb.HR/42.14(1472), HPIV2i/Zagreb.HR/38.12(2858) and HPIV2i/Zagreb.HR/47.11(16 475) with guinea pigs’ sera immunised with either HPIV2i/Zagreb.HR/47.14(1688) or HPIV2i/Zagreb.HR/42.14(1472) (n = 6 for each virus). Values are mean ± STD. *P < 0.001.

These findings indicate that there are some antigenic differences between the two genotypes which are the cause of the lower neutralisation specificity and possibly the cause of waning of the genotype G3 from the population and establishing the dominance of the genotype G1a.

Analysis of in silico predicted conformational B cell epitopes in the hemagglutinin neuraminidase and fusion protein

Although the attachment protein is the major target of neutralising antibodies, antibodies raised against F proteins are also able to neutralise infections [1, 2]. To further analyse the mechanism of how one HPIV2 genotype is replacing another, conformational B cell epitopes for two viral surface glycoproteins, fusion protein and hemagglutinin neuraminidase, were predicted in silico. To increase the relevance of the predicted epitopic regions three servers which predict conformational B cell epitopes were used: DiscoTope2.0, ElliPro and BepiPred-2.0. Only those regions predicted to contain epitopes which were predicted by all three tools were selected.

Analysis of the hemagglutinin neuraminidase

The putative hemagglutinin neuraminidase sequences were subjected to homology modelling what yielded a homo-tetramer model with >48% of homology. A PDB structure was submitted to in silico prediction of potential epitopes. Each search resulted in several predicted epitopes of different size but only two regions were identified as potential epitopes by all three servers: a region between amino acids 314 and 361 and a region between amino acids 474 and 490 (Table 2). Although there are 24 unique sites throughout the hemagglutinin neuraminidase (24/571 amino acid residues (Supplementary Fig. S1)) alignments of the predicted epitope sequences show that these two regions contain amino acid residues unique for the genotype (Fig. 5a). Additionally, these two potential epitopic regions were searched for N-linked and O-linked glycosylation sites. Five sites with different glycosylation status for viruses 1472 and 1688 were identified: 338Thr, 244Ser and 481Thr were glycosylated in 1688 and 351Ser and 360Asn were glycosylated in 1472 (Table 3). Three of the differentially glycosylated sites are the genotype unique amino acid residues (Ser351Gly, Asn360Tyr, Met481Thr), while the other two differentially glycosylated sites are not but might be influenced by the neighbouring unique residues (Table 3).

Fig. 5. Predicted conformational epitopes of (a) hemagglutinin neuraminidase and (b) fusion protein of Croatian isolates belonging to genotypes G1a and G3. The genotypes are indicated by the brackets on the right side. Non-matching residues are shown. Potentially glycosylated residues are marked with asterisk (*) and positively selected codons are marked with an arrow (↓).

Table 2. In silico predicted B cell epitopes for hemagglutinin neuraminidase and fusion protein.

Different residues in HPIV2 1472 and HPIV2 1688 are bold and underlined.

Table 3. Glycosylated and positively selected sites in predicted epitopes in hemagglutinin neuraminidase

^a Virus with glycosylated site is in parenthesis.

Positive selection analysis of epitopic sites is advantageous because the mutations in these sites are highly favourable for adaptive evolution. To identify the existence of positive selection pressure at individual codon sites, four likelihood procedures (SLAC, FEL, IFEL, REL) were used for the analysis of the in silico identified epitope regions 314–361 and 474–490. The analysis was performed with the complete HN gene sequences of all HPIV2 used for phylogenetic analysis in section Distribution of HPIV2 genotypes in Croatia in the period 2011–2017. SLAC and FEL were not able to identify any positively selected site. In contrast, IFEL and REL identified six positively selected codons (Table 3, Fig. 5a), four of which encoded for amino acids located in the epitope region 314–361 and differ between HPIV2 1472 and HPIV2 1688 (Thr333Lys, Arg346Gln, Ala349Ile, Gly351Ser). The other two positively selected codons encoded for amino acid residues in the epitope region 474–490, but are identical in both HPIV2 1472 and HPIV2 1688.

Analysis of the fusion protein

The homology modelling of the putative fusion protein sequence resulted in a homo-tetramer model with >51% of homology. A PDB structure was submitted to in silico prediction of potential epitopes. A search of potential epitopes in the fusion protein resulted in a single region between amino acids 440 and 484 which were identified by all three epitope predictor servers (Fig. 5b, Table 2). The fusion protein is more conserved than hemagglutinin neuraminidase and there are only 11 unique sites throughout the fusion protein (11/551 amino acid residues (Supplementary Fig. S2)). None of them is located in the predicted epitopic region 440–484 although there is a difference Asp469Glu between HPIV2 1472 and HPIV2 1688 (Table 2) but this substitution is also present in several G1a isolates (Fig. 5b). This region was subjected to inspection for potential glycosylation sites but no difference was found for the two viruses.

Since it seems that the difference in the neutralisation capacity between sera of guinea pigs immunised with the two HPIV2 had no origin in the fusion protein, evolutional pressure selection was not analysed for this gene/protein.

Discussion

Human parainfluenza viruses are important causes of respiratory tract diseases affecting both upper and lower respiratory tract. HPIV2 infection sometimes necessitates hospitalisations, especially in children under 2 years of age or immunocompromised patients and therefore cannot be marginalised. Nevertheless, there are very few reports on the epidemiology and evolutionary dynamics of this virus. Phylogenetic analysis [7] identified two clusters within HPIV2 isolates but the low number of isolates gave poorly informative phylogenetic tree. Later, phylogenetic tree with more different worldwide HPIV2 isolates was constructed and suggested the existence of four clusters (G1–4) [6]. We have added 38 new HPIV2 isolates ([8] and this study) and thus have generated adequate number of HPIV2 sequences to create a solid phylogenetic framework (Fig. 2) which confirms the distribution of HPIV2 isolates into four genotypes identified by [6]. Our previous study of genetic diversity of HPIV2 in Croatia between 2011 and 2014 betoken that the G1a genotype would take over the dominance of the G3 genotype in this period [8]. This has been proven in the study presented here which expanded the period of analysis up to 2017. Although the G1a genotype is now dominant, the G3 genotype can still be sporadically found among the isolates (Fig. 3). A limited number of HPIV2 isolates in this study could be misleading. However, the number of collected isolates was relatively high in 3 years: 5 isolates in 2011, 9 isolates in 2012 and 17 isolates in 2016, what reassures that conclusions drawn out from this study are relevant. Although we analysed only hospitalised cases, they are a projection of the composition of the viruses circulating in the population. Furthermore, the severity of the disease demanding hospitalisation of otherwise healthy children indicates that there could be a shift in the virus population. The overall age distribution in each season was uniform except in 2016. The high number of patients from 1 to 6 years of age in this season can be explained as either coincidentally increased number of patients which later then usually encountered HPIV2 for the first time or as secondary infection which caused a severe disease due to the shift from G3 genotype to G1a genotype which resulted in an incomplete matching of antibodies developed during primary infection (G3) and the antigen causing the disease in the secondary infection (G1a). From one of the patients in 2016 an HPIV2 belonging to genotype G3 was isolated. This patient was hospitalised at the age of almost 16. In this case the G3 infection might have been a secondary infection while primary infection could have occurred more than a decade ago, at the time when unknown HPIV2 genotype was circulating in the region. Although the number of samples is relatively low, we observed two seasons with increased number of HPIV2 infections, 2012 and 2016. Such seasonal waves may occur at the time when a new genotype enters the population, wedges out the genotype which dominated until then and cause more severe disease due to incomplete antigenic matching.

While the evolutionary processes shaping HPIV2 virus diversity within one season (autumn and early winter) may be a stochastic process, here our hypothesis was that the replacement of one genotype with another across epidemic seasons is an immune-driven process with specific neutralising antibodies as a key determinant. Our hypothesis was based on the finding that there is an antigenic diversity of the hemagglutinin neuraminidase of different HPIV2 isolates [29, 30]. Also, [31] found that monoclonal antibodies which effectively neutralised the activity of hemagglutinin neuraminidase of the prototype Greer strain isolated in 1955 failed to show any detectable binding or neutralising activity against a number of natural isolates collected from 1986 and 1987. To test our hypothesis, guinea pigs were immunised with HPIV2 isolates belonging to either G1a or G3 genotype. Although epitopes recognised by sera of different species can differ, it was unfeasible to obtain undisputed human anti-HPIV2 sera specific for a single genotype since many HPIV2 infections are treated as a common cold and it is impossible to distinguish a pre-existing antibodies from those of acute infection. Also, we analysed mainly the effect of IgG present in sera while IgA may have a more important role since mucosal surfaces are primary sites of HPIV2 infection. A cross-neutralisation with guinea pigs’ anti-HPIV2 sera revealed the presence of genotype-specific antigenic determinants and confirmed our hypothesis. Similarly, it was suggested previously that vaccine raised specific antibodies in the population may lead to the replacement of the endemic genotypes with new genotype variants for Japanese encephalitis virus and measles virus [32–34].

To analyse in depths the mechanisms of differential recognition of HPIV2 by different sera, potential B cell epitopes were predicted in silico. Since attachment protein hemagglutinin neuraminidase is the major target of neutralising antibodies and antibodies to fusion protein can have neutralising activity [2] these two proteins were analysed in depth. Hemagglutinin neuraminidase is assembled as a tetramer and fusion protein as trimer on the surface of the viral particle. Such oligomeric states form quaternary structures which shape specific antigenic sites often not recognised in the monomeric form. As a result, antiviral antibodies against the extracellular domain of surface proteins predominantly recognise conformational epitopes [1]. To increase the relevance of predicted epitopes, three different methods for prediction of conformational B cell epitopes were used and only those epitopes which were predicted by all three methods were selected for further analysis. Two and one epitopic regions were identified for hemagglutinin neuraminidase and fusion protein, respectively. To analyse the origin of the antigenic difference between the two HPIV2 genotypes as the cause of the genotype replacement, a glycosylation and positive selection of these regions were analysed in more detail. Profound analyses of the predicted epitopes in the hemagglutinin neuraminidase indicate that indeed these regions (especially epitopic region 314–361) might be responsible for the decreased match between G1a and G3 genotypes. Five differentially glycosylated sites were found: four and one for epitopic regions 314–361 and 474–490, respectively. Differential glycosylation may be associated with feeble recognition of specific epitopes by changing or shielding the contact surface between specific antibody and antigen. A strong positive selection in this region (all positively selected codons identified in the hemagglutinin neuraminidase are located either in the region 314–316 or in the region 474–490) reveals that these regions are susceptible to evolutionary pressure possibly caused by antibodies what gives a strong verification to our hypothesis that neutralising antibodies are a key determinant in the inherently complex adaptive evolution of HPIV2 in the region. Hemagglutinin neuraminidase and fusion protein are major targets of neutralisation antibodies which give a long-term protective immunity against human and animal paramyxoviruses [35–39]. Although fusion protein can be a target of selection process [40] here it seems that fusion protein has a minor role in the selection process since a unique residue for the epitopic region of this protein was not identified. This finding additionally indicates that hemagglutin neuraminidase is the main target of neutralising antibodies while fusion protein has probably a very minor role in this process in HPIV2 infection.

Although the work presented here identified two epitopic regions in hemagglutinin neuraminidase and one in fusion protein there are most certainly more epitopes in these two proteins which are recognised by the neutralising antibodies and which may differ between viruses of different genotypes. There are numerous genotype unique residues in hemagglutinin neuraminidase (Supplementary Fig. S1) which may cause the viruses of the two genotypes to further diverge from each other.

The current progress in the development of human HPIV vaccine is reviewed by [41]. It is not very likely that vaccine against HPIV2 will be developed as a monovalent vaccine, rather it would be combined with another virus or viruses like HPIV1 or 3, RSV or similar. Whatever the approach (attenuated vaccine, T-cell inducing vaccines, vaccines targeting less variable proteins like nucleoprotein or matrix protein etc.), an extreme care should be taken considering the antigenic composition of the vaccine, the measures how to control it and readapt according to continuous virus surveillance of the virus circulating in the population.

It will be of interest to monitor the circulation of HPIV2 in the population to see if these two genotypes are exchanging each other or a new, third genotype will replace the current dominance of G1a.