The role of human Metapneumovirus genetic diversity and nasopharyngeal viral load on symptom severity in adults

Sample collection, symptom severity assessment, and HMPV genotyping

A total of 3706 consenting adult outpatients who presented with symptoms of acute URTI for not more than two weeks were recruited at the Primary Care Clinic of University of Malaya Medical Centre (UMMC) in Kuala Lumpur, Malaysia between February 2012 and May 2014. During enrollment, participants were interviewed to determine their demographics (age, gender and ethnicity), estimated number of days elapsed between symptom onset and clinic visit or enrollment date, and the presence and severity of common cold symptoms [16, 17]. The common cold symptoms assessed were sneezing, nasal discharge, nasal congestion, cough, sore throat, hoarseness of voice, headache and muscle ache. The severity of each symptom was then rated by a standardized four-category ordinal scale previously reported [16–19]: none (0), mild (1), moderate (2) and severe (3). The symptom ratings were summed to create a total symptom severity score (TSSS) for each participant with a maximum points of 24 from eight symptoms [20], in which greater symptom severity was indicated by a higher score.

Nasopharyngeal swabs were collected from the patients and transferred to the laboratory in universal transport media (Copan Diagnostics, California, USA). Total nucleic acid purification was performed using the NucliSENS easyMAG automated nucleic acid extraction system (bioMérieux, Marcy I’Etoile, France) according to manufacturer’s protocol [21]. The xTAG Respiratory Virus Panel (RVP) FAST multiplex RT-PCR assay (Luminex Molecular, Toronto, Canada) and Luminex’s proprietary Universal Tag sorting system on Luminex 200 IS platform (Luminex Corp., Austin, Texas, USA) were used to detect HMPV in the samples [22]. The genotype of HMPV-positive samples was first determined by performing amplification and sequencing of the fusion (F) and attachment (G) genes as previously described [15]. This was followed by phylogenetic tree reconstruction using the maximum-likelihood (ML) method which was heuristically inferred using subtree pruning and regrafting and nearest neighbor interchange algorithms with a general time-reversible (GTR) nucleotide substitution model, a proportion of invariant sites (+I) and four categories of gamma rate heterogeneity (+Γ₄), which were implemented in PAUP version 4.0 [23]. Kimura’s two-parameter model with a reliability of branching order analyzed by bootstrap replicates of 1000 was used.

HMPV viral load quantification

For improved quantification of HMPV RNA in the nasopharyngeal specimens, a comprehensive and updated list of reference genomes was used to design a quantitative one-step reverse transcriptase-quantitative polymerase chain reaction (RT-qPCR) assay. Newly designed primer pair and a fluorescent probe targeting the highly conserved M2 gene [24] of HMPV were developed based on a set of global HMPV complete genomes representing all genotypes A1, A2a, A2b, B1 and B2 (n = 135) available in GenBank (retrieved on 31 January 2016). The sequences were codon-aligned using the web-based multiple sequence alignment program MAFFT [25] to look for conserved regions of the complete genome. Highly conserved forward primer, reverse primer and the probe with a coverage of 99.3, 100, and 99.3%, respectively, based on the alignment of the global reference sequence (Additional file 1), were designed using the Primer Express Software v2.0 (Applied Biosystems, California, USA). With reference to the nucleotide numbering of NC_004148 (HMPV reference strain), the forward primer (designated as 4730f), reverse primer (4919r) and probe (4796fp) spanned the genetic regions corresponding to nucleotide positions 4730–4754, 4919–4893 and 4796–4814 nt, respectively. The probe was labeled with the 6-fluorescein amidite (FAM) at the 5′ end, and the non-fluorescent quencher (NFQ) and minor groove binder (MGB) at the 3′-end. A synthetic single-stranded RNA oligonucleotide with a randomly generated sequence was designed and used as internal control (IC) to check for potential PCR inhibition. The unique IC sequence (5′-ACATCGTAAGGCTCCATGCAAATATGAAGATAGAATGCTTAGGACCATCAGCGAAACTCTACAATAATATCAGGCGCAGGCAGAGAAGTA-3′) showed < 10% similarity with any published sequence in the GenBank (data not shown). The IC was flanked by sequences similar to the primer binding site of the newly designed HMPV primer set, with a unique non-HMPV probe designed for the IC (VIC-5′-TTAGGACCATCAGCGAAAC-3′-NFQ-MGB). In a single reaction, 0.2 μl of reverse transcriptase (40 U), 10 μl of 2× One-Step SensiFAST Probe Lo-ROX mix (Bioline, London, UK), 0.8 μl of each primer (20 μM), 1 μl of each probe (10 μM), and 6.0 × 10³ RNA copies/μl of IC in a final volume of 20 μl. The optimized thermal cycling profile used was as follows: reverse transcription at 48 °C for 8 min, initial denaturation at 95 °C for 2 min, followed by 40 cycles of 97 °C for 2 s, and 60 °C for 20s. The thermal cycling period of the assay was short, in which viral load quantification can be accomplished within approximately 50 min in a single run. A synthetic DNA oligonucleotide containing the M2 gene sequence was used to generate a 10-fold dilution series of standard concentrations, ranging from 2.0 × 10¹ to 2.0 × 10⁶ RNA copies/μl. HMPV quantification was performed in the ABI ViiA7 Real-time PCR System (Applied Biosystems, California, USA).

The linear dynamic range of the HMPV qPCR assay was assessed using a 6-log₁₀ dilution series of HMPV M2 synthetic oligonucleotide. A standard curve plotted against quantification cycle (Cq) was built using the serial concentration. Linear regression analysis was performed to calculate the PCR efficiency and correlation coefficient based on the standard curve. In order to assess the intra- and inter-assay variability for HMPV viral load quantification, triplicate reactions were performed using a low (2.0 × 10¹ RNA copies/μl) and moderate (2.0 × 10³ RNA copies/μl) standard to determine the mean, standard deviation (SD) and coefficient of variance (CV).

Statistical analysis

Demographic (sex, age and ethnicity) and clinical (presence of the eight common cold symptoms and estimated number of days elapsed between symptom onset and enrollment date) characteristics of patients infected by different HMPV genotypes (A and B) and sub-lineages were first assessed using the bivariate analyses (Pearson’s chi-square for categorical variables, Independent Samples t-test and One-way ANOVA for continuous variables), similar to statistical techniques previously reported [12]. The overall severity of the symptoms was measured through the summation of eight individual symptom severity score (TSSS), which was modeled as a continuous variable. Association of symptom severity with virological factors (infection by different HMPV genotypes and sub-lineages, and viral load in log₁₀ RNA copies/μl) and demographic factors was performed using bivariate analyses (Independent Samples t-Test, One-way ANOVA, and Pearson’s bivariate correlation) and linear regression. Lastly, comparisons of viral load between different periods of enrollment after symptom onset (taking into account the infection by different genotypes/sub-lineages and the differences in demographics and symptom severity profiles) were made using similar analyses mentioned above. A two-sided p value lower than 0.05 was considered statistically significant. In order to control false positives in multiple statistical tests such as One-way ANOVA and linear regression, the Bonferroni correction was used to lower the critical p value of significance (performed by dividing the critical value by the number of comparisons corresponding to the number of levels in a group) [26]. All analyses were performed using the Statistical Package for Social Sciences version 22.0 (SPSS Inc., New York, USA).

Detection and genetic diversity of HMPV

During the study period, a total of 81/3706 (2.2%) nasopharyngeal swab specimens collected were tested positive for HMPV. Among them, only 7/81 (8.6%) specimens were cases of coinfection with other viruses (adenovirus [n = 1], enterovirus/rhinovirus [n = 3], coronavirus 229E [n = 1], influenza A (H3) + enterovirus/rhinovirus [n = 1] and influenza B [n = 1]). Thus, as coinfection of HMPV viruses was not common in the adult population, we included patients with viral codetection in subsequent analyses and was not considered as another variable in this study. Phylogenetic analysis of the F and G genes showed that 40/81 (49.4%) of detected HMPV viruses belonged to genotype A whereas 41/81 (50.6%) belonged to genotype B (Fig. 1, Additional file 2). Within genotype A, 25/40 (62.5%) were classified as sub-lineage A2b, whereas 15/40 (37.5%) belonged to a recently described sub-lineage of A2 (designated as unique A2 sub-lineage) [15]. Within genotype B, 25/41 (61.0%) and 16/41 (39.0%) were classified as sub-lineages B1 and B2, respectively.

Demographic and clinical association with HMPV genetic diversity

The presence of common cold symptoms assessed in this study (sneezing, nasal discharge, nasal congestion, cough, sore throat, hoarseness of voice, muscle ache and headache) was self-reported by the patients. The majority of the patients experienced cough (n = 80/81, 98.8%) followed by hoarseness of voice (n = 72/81, 88.9%), while muscle ache and sneezing (n = 56/81, 69.1%) were the least experienced by patients, respectively (Table 1). No significant differences were noted in the manifestation of all symptoms between different genotypes or between different sub-lineages. Though, there were noticeably more HMPV genotype B patients (n = 31/41, 75.6%) that experienced nasal congestion compared to genotype A patients (n = 23/40, 57.5%) (p = 0.084). The estimated number of days elapsed between symptom onset and enrollment date as reported by the patients ranged from 1 day to 2 weeks, where most patients were enrolled after a symptomatic period of 3–5 days (n = 36/81, 44.4%) (Table 1). While no significant association was found between different enrollment periods after symptom onset and HMPV genotypes or sub-lineages, it was observed that HMPV genotype A-infected patients enrolled slightly later after symptom onset (5.15 ± 3.63 days) compared to genotype B patients (4.12 ± 2.55 days) (p = 0.143).

Symptom severity of HMPV-infected patients

In terms of age, mean TSSS was significantly higher in the young and middle-age adults who were < 65 years old (12.86 ± 4.55) as compared to the elderly population who were ≥ 65 years old (9.50 ± 3.71) (p = 0.008) (Table 2). Moreover, when age was modeled as a continuous variable, it had a significant negative correlation with TSSS (r = − 0.335, p = 0.002). In a simple linear regression model, for the increase of every unit (years) in age, there was a significant decrease of 0.080 unit (score) of TSSS. Thus, our results indicated that age was the only significant predictor of TSSS, with the elderly patients experiencing less severe symptoms than the young and middle-age adults in an outpatient setting. Though, there was no significant difference observed in the estimated enrollment period after symptom onset between patients who were < 65 years old (4.58 ± 3.36 days) and ≥ 65 years old (4.81 ± 2.20 days) (p = 0.797).

The impact of viral load on symptom severity

The performance of the improved one-step RT-qPCR assay targeting the M2 gene of HMPV was first evaluated using a 6-log₁₀ dilution series of HMPV M2 synthetic oligonucleotide standard. The assay produced a typical amplification plot and a standard curve with a correlation coefficient of 0.999 and amplification efficiency of 96.78%, with a coefficient of determination (R ² ) of 0.996 (Additional file 3). The intra- and inter-assay variability was estimated within standard range (Additional file 4). Nasopharyngeal viral load in 78/81 HMPV-positive specimens were quantified, with a success rate of 96.5% (Additional file 2). Besides being able to detect all different genotypes and sub-lineages of HMPV, the lowest quantifiable concentration of the assay was estimated at 13 RNA copies/μl, while the highest viral load recorded in our specimen was 731,917 RNA copies/μl.

At different periods of enrollment after symptom onset (≤1–2 days, 3–5 days, and ≥ 6 days), we observed that patients who enrolled in 3–5 days had the highest viral load (3.77 ± 1.20 log₁₀ RNA copies/μl, n = 35) compared to those who enrolled in ≤1–2 days (3.55 ± 1.53 log₁₀ RNA copies/μl, n = 20) and in ≥6 days (2.93 ± 1.16 log₁₀ RNA copies/μl, n = 23) after symptom onset (Fig. 2a, Additional file 5). Even though no bivariate association between viral load and periods of enrollment after symptom onset was observed (p = 0.052), a post-hoc test using the Bonferroni procedure showed that those who enrolled in 3–5 days had higher viral load (p = 0.048, insignificant when Bonferroni-corrected p < 0.0083 instead of p < 0.05 was used) compared to those who enrolled in ≥6 days after symptom onset.

However, the impact of viral load on symptom severity was investigated by performing bivariate correlation between the estimated viral load and TSSS, in which no significant correlation was observed (r = − 0.118, p = 0.303) (Table 2). Taking into account the different periods of enrollment after symptom onset, there was also no significant correlation observed between viral load and TSSS in these periods (Additional file 5). When the patients’ TSSS was grouped into three categories with scores that range from 1 to 8, 9–16, and 17–24 to reflect the relative symptom severity of mild, moderate, and severe, respectively [20] (Fig. 2b), no significant differences in viral load between the three categories of symptom severity were observed. Though, among the group of patients with a TSSS of 1–8 (indicative of having milder symptoms) (n = 19), those who enrolled ≥6 days after symptom onset had a significantly lower viral load (2.58 ± 1.07 log₁₀ RNA copies/μl, n = 3) compared to those who enrolled 3–5 days after symptom onset (3.99 ± 0.81 log₁₀ RNA copies/μl, n = 9) (p = 0.043), but the difference was insignificant after Bonferroni correction (Fig. 2b, Additional file 5). Overall, our analyses showed that nasopharyngeal viral load could not predict the severity of symptoms caused by HMPV infection, as well as the TSSS scoring system could not predict the amount of viral load in the patient, hence both variables were found to be not correlated in this study.

Association of viral load with HMPV genetic diversity and demographic factors

By comparing the genetic diversity of HMPV, we observed that after 3–5 days of symptom onset, genotype A patients had a significantly higher viral load (4.44 ± 1.07 log₁₀ RNA copies/μl, n = 15) compared to genotype B patients (3.27 ± 1.07 log₁₀ RNA copies/μl, n = 20) (p = 0.003) (Fig. 2c, Additional file 5). However, viral load differences between sub-lineages were not significant after Bonferroni correction (p = 0.024), though the A2b sub-lineage had a higher viral load (4.55 ± 0.96 log₁₀ RNA copies/μl, n = 11) compared to B2 (3.07 ± 1.33 log₁₀ RNA copies/μl, n = 8, p = 0.037) sub-lineage during this 3–5 days period (Fig. 2d, Additional file 5). On the other hand, we also observed that genotype A patients who enrolled ≥6 days had a significantly lower viral load (after Bonferroni correction) compared to those who enrolled 3–5 days after symptom onset (p = 0.012, post-hoc test with Bonferroni procedure) (Fig. 2c). Lastly, no significant viral load differences were found between the different demographics (sex, ethnicity, and age) of patients at any of the three periods of enrollment after symptom onset (Additional file 5).