Inference of population genetic parameters from an irregular time series of seasonal influenza virus sequences

Abstract

Basic summary statistics that quantify the population genetic structure of influenza virus are important for understanding and inferring the evolutionary and epidemiological processes. However, the sampling dates of global virus sequences in the last several decades are scattered nonuniformly throughout the calendar. Such temporal structure of samples and the small effective size of viral population hampers the use of conventional methods to calculate summary statistics. Here, we define statistics that overcome this problem by correcting for the sampling-time difference in quantifying a pairwise sequence difference. A simple linear regression method jointly estimates the mutation rate and the level of sequence polymorphism, thus providing an estimate of the effective population size. It also leads to the definition of Wright's FST for arbitrary time-series data. Furthermore, as an alternative to Tajima's D statistic or the site-frequency spectrum, a mismatch distribution corrected for sampling-time differences can be obtained and compared between actual and simulated data. Application of these methods to seasonal influenza A/H3N2 viruses sampled between 1980 and 2017 and sequences simulated under the model of recurrent positive selection with metapopulation dynamics allowed us to estimate the synonymous mutation rate and find parameter values for selection and demographic structure that fit the observation. We found that the mutation rates of HA and PB1 segments before 2007 were particularly high and that including recurrent positive selection in our model was essential for the genealogical structure of the HA segment. Methods developed here can be generally applied to population genetic inferences using serially sampled genetic data.

Keywords: influenza virus; mismatch distribution; serial sample; summary statistics.

Figure 1

Coalescent tree of two viral sequences that are sampled at times different by τ. Assuming a constant rate µ of neutral (synonymous) mutation along a lineage, the expected neutral sequence difference is given by $(2E\left[T\right]+\tau )\mu$, where T is time to the coalescence of two contemporaneous sequences. Therefore, the expectation of synonymous sequence difference is greater than the scaled mutation rate, 2E[T]µ = 2 N_eµ, and the difference is τµ.

Figure 2

Pairwise nucleotide difference (d) per site of segments HA, NA, PB1, PB2, PA, and NP plotted against sampling time difference (τ, in days) for H3N2 data sequences. Data points are from 27-year data (1980–2006; black dots) and from the 10-year data (2007–2017; gray dots). Regression lines for 27- and 10-year data are shown in red and blue, respectively. The proportions of bootstrap samples in the tests for the statistical difference between 27- and 10-year periods of $\widehat{\mu }$, $\stackrel{⌣}{\pi }$, and $\widehat{N_e}$ are shown below each regression plot.

Figure 3

The TCMDs of six influenza virus segments in the 27-year (A) and 10- year (B) H3N2 data sets. To obtain $\stackrel{⌣}{d}$ and make histograms ${\tau }_{\max }=300$ and bin size w = 0.002 were used.

Figure 4

The average TCMD (black curve) for simulated data under neutrality (s = 0) with m = 0.004 and $K_{\max }=110$ that produce the best-fitting $F_{\text{ST}}$ and $\stackrel{⌣}{\pi }$ values to the observed data (red curve; the TCMD of HA segments in the 10-year data set). TCMDs of individual simulation replicates are shown in gray curves. To obtain $\stackrel{⌣}{d}$ and make histograms ${\tau }_{\max }=300$ and w = 0.002 were used.

Figure 5

The average TCMD (black curve) for simulated data under positive selection (s = 0.1 and ε = 10) with m = 0.00025 and $K_{\max }=6700$, which is congruent to the TCMD (red curve) of HA segment in the 10-year H3N2 data (KS test, p = 0.058). TCMDs of individual simulation replicates are shown in gray curves. To obtain $\stackrel{⌣}{d}$ and make histograms ${\tau }_{\max }=300$ and w = 0.002 were used.