Genetic drift and purifying selection shape within-host influenza A virus populations during natural swine infections

Two swine IAV subtypes co-circulated at the fair

The six-day agricultural fair took place in 2014 and involved the showing of a total of 423 exhibition swine. Nasal wipes were collected daily from all pigs present at the fair (Fig 1A). All 2159 of the collected nasal wipes were screened for IAV using a real-time reverse transcription-polymerase chain reaction (rRT-PCR) assay targeting the M segment. On day 1 of the fair, few pigs were found to be PCR-positive for IAV (33 of 419 samples; S1 Table). Over the next five days, PCR positivity rates increased from 7.9% (day 1) to 88.0% (day 6; 139 of 158 samples) (Fig 1A). A subset of the positive samples were tested using an IAV subtype-specific assay, which revealed that both H3N2 and H1N1 viruses were circulating at the fair (S1 Table). Samples with an M segment rRT-PCR cycle threshold (Ct) value of ≤31.0 were selected for next generation sequencing, resulting in sequenced samples from a total of 253 pigs (Fig 1B and S2 Table). These sequenced samples had a read depth across the genome of approximately 4000 reads per nucleotide (S1 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Overview of sampling and sequencing of pigs at the county fair.

(A) Sampling, testing, and sequencing effort over the course of the fair. The decrease in the number of samples collected on days 5 and 6 stems from departure of a subset of pigs. (B) Distribution of the number of sequenced samples per pig. Pigs with no sequenced samples either never tested positive for IAV or tested positive at some point but positive samples all had Ct values exceeding 31.0.

 

https://doi.org/10.1371/journal.ppat.1012131.g001

To characterize diversity of IAV genotypes present, we inferred consensus sequences for each sample for each of the eight IAV gene segments. Haplotype networks reconstructed from these consensus sequences indicated that two distinct lineages co-circulated at the fair (Fig 2), consistent with the subtyping results. We labeled these two lineages I and II, each comprising a set of eight gene segments that were grouped together based on the associations that were evident in the data. By comparing the outcomes of the subtype-specific tests against the assigned lineages of the consensus sequences, we were able to conclude that the hemagglutinin (HA) and neuraminidase (NA) gene segments of lineage I belonged to swine IAV subtype H1N1 and that the HA and NA gene segments of lineage II belonged to swine IAV subtype H3N2. Specifically, the lineage I HA and NA segments were from the gamma-c2 and classical H1N1 swine clades, respectively, and the lineage II HA and NA segments were from the IBV1 and 2002B H3N2 swine clades, respectively. Each of the two lineages comprised a single dominant consensus sequence and many consensus ‘singletons’. This indicates that the outbreak at the fair likely originated from a very small number of pigs (possibly only 2) that arrived at the fair already infected.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Haplotype networks of IAV segments support co-circulation of two distinct lineages at the county fair.

Consensus sequences from all successfully sequenced samples were inferred and used to construct haplotype networks. Each circle represents a unique haplotype (excluding differences caused by ambiguous bases), with circle size of a haplotype being proportional to the number of consensus sequences that belong to it. Numbers in parentheses refer to the number of differentiating single nucleotide polymorphisms (dSNPs) that distinguish the dominant haplotypes from each lineage in each internal gene segment. Because of the high genetic divergence between the HA of lineage I and the HA of lineage II, and between the NA of lineage I and the NA of lineage II, dotted lines rather than solid lines are shown between the lineages.

 

https://doi.org/10.1371/journal.ppat.1012131.g002

Despite the majority of samples having consensus sequences whose gene segments belonged either exclusively to lineage I or to lineage II, a small number of samples had consensus sequences with at least one gene segment belonging to a different lineage than the others (e.g., vial 14SW5957). These putatively reassortant samples indicate that at least some of the infected pigs at the fair were likely to have been heterosubtypically coinfected. To identify samples that were heterosubtypically coinfected, we looked at viral diversity below the consensus level. For this purpose, we selected representative consensus sequences of each segment from both lineages to serve as references for read mapping. We first assessed whether there was evidence for HA coinfection or NA coinfection. Because the consensus sequences of the HA and of the NA segments of lineages I and II are sufficiently distinct at the nucleotide level to support unambiguous mapping, we detected co-occurrence of these segments by mapping the samples’ reads to H1, N1, H3, and N2 reference genomes. Samples with reads mapping to only one HA and one NA subtype were classified as singly infected at each of these gene segments (Fig 3A and S2 Table). Samples that had an appreciable number of reads mapping to both H3 and H1 were classified as coinfected on the HA gene segment (Methods). Similarly, samples that had an appreciable number of reads mapping to both N2 and N1 were classified as coinfected on the NA gene segment (Methods).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Classification of gene segments and samples into lineages.

(A) Mapping of sample reads on to reference H1, N1, H3, and N2 gene segments for four selected samples. Based on this mapping, a sample’s HA gene segment was classified as belonging to either lineage I (with reads mapped to H1; vial#14SW6024), lineage II (with reads mapped to H3; vial#14SW6331 and vial#14SW5957), or coinfected (with reads mapped to both H1 and H3; vial#14SW6164). Similarly, a sample’s NA gene segment was classified as belonging to either lineage I (with reads mapped to N1; vial#14SW6024), lineage II (with reads mapped to N2; vial#14SW6331 and vial#14SW5957), or coinfected (with reads mapped to both N1 and N2; vial#14SW6164). (B) Detection of dSNPs characteristic of lineages I and II for each of the six internal gene segments of the four selected samples. (C) The proportion of samples classified as lineage I, II, or coinfected, by gene segment. Segments were labeled as unknown when there was insufficient coverage across the focal gene segment to determine nucleotide identities at the dSNP sites or to robustly map to the HA and NA gene segments. Samples (‘All’) were classified as belonging to lineage I, lineage II, coinfected, or reassortant based on their constituent segments. A sample was defined as a lineage I or II sample if at least five of its gene segments were successfully classified and found to be either singly infected lineage I or singly infected lineage II, with no evidence of coinfection. A sample was considered a reassortant if at least five segments were successfully classified and a combination of lineage I singly infected and lineage II singly infected segments were present in the same sample. A sample was considered coinfected if at least one of its gene segments was classified as coinfected. In panels (A) and (B), each row corresponds to a different sample classification. Vial#14SW6024 corresponds to a sample with a lineage I virus, vial#14SW6331 corresponds to a sample with a lineage II virus, vial#14SW6164 corresponds to a coinfected sample, and vial#14SW5957 corresponds to a reassortant virus.

 

https://doi.org/10.1371/journal.ppat.1012131.g003

To classify a sample as coinfected based on one of the six internal gene segments (PB2, PB1, PA, NP, M, and NS), a different approach was needed because the internal segments of lineages I and II share 94–98% average nucleotide identity. As such, most internal gene segment reads would map to the reference genomes for both lineage I and lineage II, and coinfection would therefore be concluded erroneously. To assess coinfection at these internal gene segments, we identified the set of single nucleotide polymorphisms (SNPs) that differentiated lineage I from lineage II in each of the six internal gene segments (S3 Table). Each of these differentiating SNPs (dSNPs) was characteristic of and unique to one of the two lineages. For each sample, we counted the number of lineage I and the number of lineage II dSNPs that were supported by the mapped reads (Fig 3B). Based on these results, we classified a given internal gene segment as either singly infected with lineage I, singly infected with lineage II, or coinfected (S2 Table). Once every gene segment of a sample was classified as either singly infected or coinfected, we summarized these results at the sample level (Fig 3C, ‘All’). The majority of samples contained only lineage II gene segments (subtyped as H3N2). Samples containing only lineage I gene segments (subtyped as H1N1) were also apparent. We further found a non-negligible number of coinfected samples (n = 70; S2 Table), and a small number of reassortant samples (n = 9; S2 Table).

To improve our ability to call variants in our downstream analyses, we selected 384 samples to re-sequence at even greater depth using two NovaSeq 6000 lanes (S1 and S2 Tables). Re-sequencing of these samples resulted in a read depth of approximately 20,000 reads per nucleotide, substantially higher than the read depth in our initial sequencing runs (S1 Fig).

Within-host IAV diversity is low and largely synonymous in singly infected samples

To evaluate viral evolutionary dynamics without the complications introduced by co-infection, we decided to focus on within-host IAV diversity in pigs with samples that contained only lineage I or lineage II gene segments throughout their duration of infection. Samples classified as ‘unknown’ were presumed to derive from the same lineage as prior or later samples from the same pig and therefore did not affect categorization of pigs as persistently infected with lineage I or II. However, samples from pigs with any positive evidence of coinfection or reassortment were excluded. Among the re-sequenced samples, there were 260 samples from 168 pigs that met these criteria. We called intrahost Single Nucleotide Variants (iSNVs) in these 260 samples using a variant calling threshold of 3% at sites that exceeded 500x coverage. The majority of these samples contained only low levels of IAV diversity, with 90% of the samples containing fewer than 15 iSNVs (Fig 4A). Most identified iSNVs were present at low frequencies of <10% (Fig 4B), paralleling findings in human IAV infections [19,20]. Overall, for both lineage I and lineage II samples, nonsynonymous iSNVs were present at lower frequencies than synonymous iSNVs, consistent with patterns one would expect in the presence of purifying selection. We further found that the ratio of the number of nonsynonymous to synonymous iSNVs generally was significantly lower than the neutral expectation, given by the ratio of nonsynonymous to synonymous sites (Fig 4C). This was the case for both lineages and across genes, and again points towards the role of purifying selection in shaping IAV diversity within swine. In two instances (lineage I M gene and lineage I NS gene) the calculated NS/S ratio was higher than the neutral expectation. However, in these cases, the 95% confidence intervals were exceptionally broad, with lower bounds that were lower than the neutral expectation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Singly infected samples harbor relatively low levels of genetic diversity and show evidence of purifying selection.

(A) Distribution of the number of iSNVs identified per sample, across samples. (B) The proportion of detected iSNVs that fall below a given frequency, as specified on the x-axis. Approximately 70% iSNVs are detected at frequencies below 10%. Results are shown by sample lineage (I or II) and separately for nonsynonymous and synonymous iSNVs. (C) The ratio of the number of nonsynonymous to synonymous iSNVs, by gene segment and lineage. These ratios are shown alongside the neutral expectation, given by the ratio of nonsynonymous to synonymous sites. Black whiskers show the 95% confidence interval of the ratio.

 

https://doi.org/10.1371/journal.ppat.1012131.g004

The presence of spurious iSNVs in our dataset would result in NS/S ratios that are closer to the neutral expectation than might otherwise be the case. Because of this, we plotted the number of iSNVs called in a sample against the Ct value of the sample (S2 Fig). A positive association between Ct value and the number of iSNVs could indicate that PCR amplification during sequencing may have generated spurious iSNVs. Despite our conservative variant calling threshold and our high coverage requirement, we did observe a positive association between Ct value and the number of iSNVs, indicating that some of the iSNVs that we detected are likely spurious. As such, the observation that nonsynonymous to synonymous ratios generally were lower than the neutral expectation despite the likely presence of spurious iSNVs reinforces our conclusion of purifying selection at play. Finally, the presence of spurious iSNVs in our dataset would contribute to the number of iSNVs observed, and as such, levels of IAV diversity are likely lower than those shown in Fig 4A.

A strong role for genetic drift in shaping within-host IAV evolution

To further characterize the evolutionary processes acting on swine IAV populations within singly infected pigs, we analyzed patterns of viral diversity in pigs with two or more longitudinal samples available. S3 Fig shows iSNV dynamics, by segment, for all singly infected pigs with at least two re-sequenced time points (n = 82 pigs). To summarize these dynamics, we plotted the number of iSNVs in a sample against the timepoint during the pig’s infection at which the sample was taken (Fig 5A), as given by the number of days following the pig’s first PCR positive sample. We did not find a temporal pattern in iSNV richness over the course of infection, mirroring the finding in human IAV infections, where the number of iSNVs does not appear to change over the course of infection [19].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Longitudinal dynamics of iSNVs in singly infected pigs.

(A) The number of iSNVs detected in each sample, plotted according to the number of days since that animal had its first PCR positive sample. Results are stratified by IAV lineage. (B) Changes in the frequency of synonymous (left) and non-synonymous (right) iSNVs, plotted according to the number of days between samples. iSNVs that were detected at both time points are shown in blue. iSNVs that were detected only in the first sample are shown in red. iSNVs that were detected only in the second sample are shown in yellow. (C) Cumulative distribution showing the proportion of iSNVs that persist between one day and the next with a frequency change that is less than or equal to the frequency change shown on the x-axis. Red line shows results from the swine data; the blue line shows results from humans IAV infections, as calculated from the data provided in [19]. For both datasets, we called iSNVs as present if found at frequencies of ≥3% and ≤97%.

 

https://doi.org/10.1371/journal.ppat.1012131.g005

To further quantify patterns of within-host viral evolution, we calculated, as in [19], changes in iSNV frequencies from one sequenced sample to the next sequenced sample. We plotted these changes as a function of the time between sample collection, stratifying by whether the iSNV was a nonsynonymous or a synonymous variant (Fig 5B). The apparent similarity in pattern between the nonsynonymous and synonymous iSNV frequency changes points towards the role of genetic drift playing a large role in within-pig IAV evolution. Considering all iSNVs together, we further found that a large proportion of iSNVs are only observed once in a pair of samples sequenced one day apart (Fig 5B). Specifically, in the subset of cases where there was sufficient depth and sequence quality in the later time-point to make an iSNV theoretically detectible, we calculated that 34% of iSNVs identified on the first day were not observed on the following day. Similarly, 30% of iSNVs identified in the second sample of a pair were not observed in the sample taken one day earlier. In comparison, using the data provided in [19] with an analogous variant-calling threshold of 3%, we calculated the percentage of iSNV loss and gain in human IAV infections to be 24% and 50%, respectively, between samples taken a day apart. In both swine and humans, the frequencies of iSNVs that persisted between samples taken one day apart often times changed dramatically (Fig 5B), with >50% of iSNVs that persisted from one day to the next exhibiting greater than a 12% change in frequency (Fig 5C). Further, the frequency change patterns in our swine data are remarkably similar to those calculated from the human IAV data provided in [19], suggesting that the strength of genetic drift acting on IAV populations is similar in pigs and humans.

Ethics statement

Sampling of swine was reviewed and approved by The Ohio State University Institutional Animal Care and Use Committee under protocol number 2009A0134.

Sample collection

Nasal wipes [32] were collected upon arrival at the agricultural fair during the swine exhibition weigh-in procedure, and then subsequently every night of the fair at 24-hour intervals as previously described in detail [33]. Some animals departed the fair prematurely and were thus not sampled for the full six days. In brief, individual pig identification tags were recorded with samples, which were preserved on dry ice in the field after collection and for transportation to laboratory and long-term storage at -80°C. Animals were not systematically assessed for clinical signs of disease due to the large number of samples that needed to be collected over a short timeframe. IAV infection was assessed using rRT-PCR targeting the M segment: National Veterinary Services Laboratory PCR primer protocol (no. SOP-BPA-9034.04) with SuperScript One-Step RT-PCR (Invitrogen). Samples were defined as IAV positive using a Ct cutoff of <45.0. Virus isolation on MDCK cells was attempted on a subset of samples and HA/NA subtypes of the viral isolates were evaluated using subtyping rRT-PCR (VetMAX-Gold SIV Subtyping Kit; Life Technologies, Austin, TX USA).

RNA extraction and sequencing

Short read processing and alignment

Raw reads were processed using bbtools to remove sequencing adapters and trim low-quality bases [35]. We then used IRMA (Iterative Refinement Meta-Assembler) [36] to generate consensus sequence for each sample. These consensus sequences indicated that there were two IAV subtypes co-circulating at the fair: one from a swine H1N1 lineage and one from a swine H3N2 lineage (referred to as lineages I and II in this manuscript, respectively). For any given segment, the within-lineage variation in consensus sequence was low (average pairwise distance was within one nucleotide difference). Accordingly, the major haplotype of each segment from each lineage was selected to be used as a reference sequence (16 total reference sequences: 8 from lineage I and 8 from lineage II). Reads from each sample were split according to segment and subtype to which they shared the greatest nucleotide similarity. To do this, all reads were aligned to the 16 reference sequences using BLAT [37]. Each subtype and segment specific bin of short reads was then mapped onto their respective reference sequence using bbmap with the following settings: 90% minimum percent identity, toss ambiguously mapped reads, and global alignment mode. The coverage for all samples, by segment, is shown in S4 Fig.

Reconstruction of haplotype networks and identification of differentiating SNPs (dSNPs)

Using the Run 1 sequencing data, consensus sequences were generated for all segments in each sample by identifying the major allele at each site with ≥10x coverage. Sites with <10x coverage were recorded as gaps. Since the internal gene segments (PB2, PB1, PA, NP, M and NS) share on average 93–97% nucleotide identity between lineages, reads were mapped separately to lineage I and lineage II references, but total coverage and allele frequencies were calculated using alignments to both references. This method was adopted because, although the similarity of internal segments was sufficiently high to permit mapping of all reads to the same reference, this practice made it difficult to evaluate the significance of low-frequency variants. Low quality samples and potentially coinfected samples were excluded from alignments by removing all samples that contained ≥1% polymorphic sites. Using these alignments containing only sequences with high coverage from samples unlikely to be coinfected, we inferred haplotype networks using a custom python script (GitHub: https://github.com/Lowen-Lab/swineIAV). SNPs that differentiate lineage I from lineage II viruses (dSNPs) were identified as those SNPs that distinguish the dominant lineage I and lineage II haplotypes. The number of detected dSNPs for all samples is shown in S4 Fig.

Classifying sequenced segments and samples into lineage designations

Because lineage I and II HA and NA segments share little nucleotide similarity, short reads map reliably to only one reference or the other. Mapping coverage is therefore a reliable method of distinguishing H1 from H3, and N1 from N2. Accordingly, for HA and NA segments, positive evidence of lineage I or II genotypes required that reads map to at least 10% of all sites in a segment with 50x coverage or higher in Run 1 and 500x coverage in Run 2.

Internal segments (PB2, PB1, PA, NP, M and NS) were classified into lineages based on the alleles present at the dSNPs sites. The base quality and mapping quality thresholds to count a putative dSNP allele were set empirically based on the distribution of quality and mapping statistics observed in our sequence data and using our read alignment protocol. Specifically, only sites with ≥50x coverage (or 500x for re-sequenced samples), an average mapping score of the major allele ≥43, and average mismatch and indel counts of reads containing the major allele ≤1.5 and ≤0.5, respectively. At those sites, only those alleles present at ≥3% frequency, with an average phred score of ≥37, average mapping quality score of ≥40, and average location of that allele in each mapped read ≥30 bases from the nearest end of a read were included. A sample’s gene segment was classified as supporting infection with a lineage I gene segment if more than 10% of the total dSNP sites for a segment contained the alleles defining lineage I. Similarly, a sample’s gene segment was classified as supporting infection with a lineage II gene segment if more than 10% of the total dSNP sites for a segment contained the alleles defining lineage II. If a sample’s gene segment supported infection with both a lineage I and a lineage II gene segment, it was classified as coinfected. Gene segments that were not classified as either lineage I or lineage II were classified as “unknown”. As such, gene segments with low sequencing coverage were generally classified as “unknown”.

At the sample level, a sample was classified as coinfected if at least one of the eight segments was classified as coinfected. Of the remaining samples, those with five or more segments classified as lineage I or lineage II were classified as lineage I or II, respectively, at the sample level. Samples with five or more segments classified as “unknown” were classified as “unknown”. Samples with some segments classified as lineage I and others classified as lineage II were classified as reassortant at the sample level, but only if the number of “unknown” segments was less than five. We manually inspected the classification of all samples to verify the accuracy and consistency of the genotyping process.

Identification of intrahost single nucleotide variants (iSNVs)

Run 2 data were used for calling iSNVs. Similar to the identification of dSNPs, the cutoffs for including iSNVs were set empirically. First, sites were evaluated based on their total coverage and the average quality and mapping statistics of reads containing the major allele. In particular, only sites with ≥500x coverage, an average mapping score of the major allele ≥43, and average mismatch and indel counts of reads containing the major allele ≤1.5 and ≤0.5, respectively, were included. These major allele focused site exclusion criteria were introduced out of necessity and functioned primarily to exclude false positives driven by mapping errors, such as at the borders of defective viral genomes (DVGs). For minor variants at these sites to be included in subsequent analyses, they were required to be present at ≥3% frequency, have an average phred score of ≥37, and the reads that contained the minor allele at any given site also had to have sufficient mapping quality to justify inclusion. Specifically, reads containing the minor allele needed an average mapping quality score of ≥40, the average location of the minor allele needed to be ≥30 bases from the nearest end of a read, and the reads overall needed to have ≤2.0 and ≤1.0 average mismatch and indel counts, respectively, relative to the reference sequence.