Influenza B Viruses Exhibit Lower Within-Host Diversity than Influenza A Viruses in Human Hosts

Within-host genetic diversity of IBV in natural infections.

All samples exhibited low genetic diversity. The vast majority had no iSNV above the 2% cutoff. Of the 99 samples with high-quality NGS data, 70 had no minority iSNV, 17 had one iSNV, 7 had two iSNV, and 3 samples had 3 iSNV (median, 0; interquartile range [IQR], 0 to 1) (Table 2). Two outliers had a large number of iSNV, with 8 and 20 iSNV. These two samples came from the same individual, with one collected at home and the second at the study clinic 2 days later. Most of the iSNV in these two samples were present at similar frequencies, 3 to 5% in the home sample and 17 to 23% in the clinic sample (Table 3), both of which were sequenced in duplicate on separate Illumina runs. The high number of mutations present at similar frequencies is suggestive of a mixed infection with distinct haplotypes or strains as opposed to de novo mutations arising on a single genetic background. The iSNV in the home-collected sample are all found in the subsequent clinic-collected sample, each with a similar change in frequency across the two samples. This further supports the conclusion that these mutations are on the same genome in a mixed infection with two distinct strains.

We examined how within-host diversity changes by day of sampling during IBV infections, as the virus population rapidly expands and contracts. As we have previously shown that specimen viral load can affect the sensitivity and specificity of variant identification (24), we sought to control for this variable in our analysis. Although viral load generally decreased with time after symptom onset (Fig. 1A), we found that within-host diversity as measured by number of identified minority iSNV did not vary with viral load (Fig. 2A) (P = 0.2996, adjusted r² = 0.0009) or with day of infection (Fig. 2B) (P = 0.62, analysis of variance [ANOVA]). The frequencies of the identified iSNV were consistent across replicate libraries from the same samples, indicating that our measurements of iSNV frequency are precise (Fig. 2C).

We detected minority iSNV across all eight genome segments (Fig. 3). The ratio of nonsynonymous to synonymous iSNV was 0.74, which given the excess of nonsynonymous sites across the genome suggests significant purifying selection. There was only one minority iSNV present in more than one individual; we identified a variant encoding a synonymous mutation in PB1 in two individuals from separate households infected with B/Yamagata in the 2016–2017 season. We did not identify any nonsynonymous minority iSNV in the known antigenic sites of IBV hemagglutinin, which suggests that positive selective pressure for variants that escape antibody-mediated immunity is not particularly strong within hosts. We found that there is no difference in distribution of the number of iSNV per sample between vaccinated and nonvaccinated individuals (Fig. 4A). During the first few seasons of the study, some individuals received trivalent vaccines, which contain only one of the two IBV lineages. We therefore repeated this analysis, excluding three individuals for whom we had no information about specific vaccine product and reclassifying six individuals who received trivalent vaccines and were infected with a lineage not included in that season’s trivalent formulation as “unvaccinated.” We again found no difference in the number of iSNV between groups (Mann-Whitney U [MWU] test, P = 0.9103). Together, these data indicate that vaccine-induced immunity is not a major diversifying force for IBV within hosts in our study population. This is consistent with our previous work on IAV in the HIVE as well as a randomized-controlled trial of vaccine efficacy (FLU-VACS), both of which showed no difference in intrahost diversity based on same-season vaccination status (19, 20). Intrahost diversity was similar between B/Victoria and B/Yamagata virus populations (Fig. 4B), consistent with our previous comparison of subtype A/H3N2 and A/H1N1 viruses (19).

We compared the within-host genetic diversity of IBV to our previously published data on IAV from the HIVE cohort (19). Here, IBV exhibits lower within-host diversity than IAV (Table 4) (P < 0.001, MWU test). IBV samples had a lower median number of minority iSNV (median, 0; IQR, 0 to 1) than IAV samples (median, 2; IQR, 1 to 3); 71% of IBV samples contained no minority iSNV, compared to 20% of IAV samples. The difference in within-host diversity was robust to relaxation of the iSNV frequency cutoff to 1% and 0.5%. IBV also exhibited a lower within-host diversity, as measured by nucleotide diversity (π), which takes both the number and the frequency of iSNV into account (Fig. 4C and Table 4). To ensure that our results were not an artifact of overly stringent quality thresholds, we also identified minority iSNV with less conservative read mapping quality (MapQ) and base quality (Phred) scores. We identified the same set of minority iSNV with a MapQ cutoff of 20 as with the original cutoff of 30. Similarly, reduction of the Phred base-quality cutoff to >25 in addition to a MapQ score cutoff of >20 resulted in only 20 more minority iSNV, 8 of which were found in the individual with a mixed infection. The other additional 12 minority iSNV were dispersed across specimens and did not significantly change the overall distribution of within-host diversity. We also examined whether our results were biased by use of a single B/Yamagata and B/Victoria reference for alignment and variant calling, which were both drawn from the 2012-2013 season (see Materials and Methods). We realigned sequence data from 43 of the original 99 samples to season-specific reference genomes isolated in southeastern Michigan. We found that the overall alignment rate for any given specimen was similar between the original reference and the new season-matched reference. Variant identification based on the new references and the original quality thresholds resulted in the same distribution of within-host diversity, although the identities of some iSNV were different (Fig. 4D).

Together, these results indicate that our measurements of within-host diversity are robust to several technical aspects of variant identification, which are unlikely to account for the lower observed diversity of IBV. Because these data are from the same cohort and were generated using the same sequencing approach and analytic pipeline as our previous IAV data sets, the observed differences likely reflect true biological differences between IAV and IBV.

Identification of household transmission pairs.

We compared viral diversity across samples from individuals in the same household to investigate the genetic bottleneck that influenza B viruses experience during natural transmission. Over the seven influenza seasons, 39 households in the HIVE cohort had two or more individuals positive for the same IBV lineage within a 7-day interval (Table 1). This epidemiologic linkage is suggestive of transmission events but does not rule out coincident community-acquired infection (19). We identified 16 putative transmission pairs for which we sequenced at least one sample from each individual. In one of these pairs, the putative recipient was the individual with a mixed infection. The donor did not have evidence of a mixed infection based on number of iSNV, which would imply that the recipient may have been infected twice or that the second virus was lost from the donor by the time of sampling. This pair was excluded from the between-host analysis, leaving 15 putative transmission pairs for which we have high-quality sequencing data on both donor and recipient influenza populations.

We used our sequencing data to determine which of these epidemiologically linked household pairs were actual IBV transmission pairs. We generated maximum likelihood phylogenetic trees for samples from the two IBV lineages using the concatenated coding consensus sequences. Phylogenetic analysis provided genetic evidence that the 15 epidemiologically linked pairs were indeed true transmission pairs, as epidemiologically linked pairs were found nearest each other in each tree (Fig. 5A and B; vertical bars with household identifiers [ID]). We also validated these transmission pairs by analyzing the genetic distance across viral populations. True transmission pairs should have genetically similar populations exhibiting low genetic distance, while individuals with coincident community acquisition are more likely to have populations with a higher genetic distance. We compared the genetic distance between epidemiologically linked household pairs and random community pairs from the same season and infected with the same IBV lineage, using L1-norm as the measurement of genetic distance (Fig. 5C). The distribution of random community pairs functions as a null model of genetic distances among locally circulating strains. All of the 15 putative transmission pairs fell on the tail of this distribution, below the 5th percentile of the community pair L1-norm distribution, indicating that they are true transmission pairs (Fig. 5C). While the L1-norm is a function of both the consensus sequence and the iSNV, this signal was driven predominantly by consensus differences, as reflected in the phylogenetic analysis.

Comparison of viral diversity across transmission pairs.

Transmission bottlenecks restrict the genetic diversity that is passed between hosts. With a loose transmission bottleneck, many unique genomes will be passed from donor to recipient. Because this will allow two variants at a given site to be transmitted, sites that are polymorphic in the donor are more likely to be polymorphic in the recipient. However, in the case of a tight or stringent bottleneck, sites that are polymorphic in the donor will likely be either fixed or absent in the recipient. We have previously demonstrated that influenza A virus experiences a tight transmission bottleneck of 1 to 2 unique genomes (19). Across our 15 IBV transmission pairs, we found no sites that were polymorphic in the donor and the recipient (Fig. 6). iSNV present in the donor were either fixed (100%) or absent (0%) in the recipient. These data suggest a stringent transmission bottleneck for influenza B virus, similar to that of influenza A virus. As there were fewer samples, transmission pairs, and iSNV in our IBV data set, we were unable to obtain a robust and precise estimate of bottleneck size.

Description of the HIVE cohort.

The HIVE study is a prospective, community-based household cohort in southeastern Michigan based at the University of Michigan School of Public Health (34–39). The cohort was initiated in 2010, with enrollment of households with children occurring on an annual basis and an active surveillance period lasting from October through May. In 2014, active surveillance was expanded to take place year-round. Participating adults provided informed consent for themselves and their children, and children ages 7 to 17 years provided oral assent. Individuals in each household were followed prospectively for acute respiratory illness, defined as two or more of the following: cough, fever or feverishness, nasal congestion, chills, headache, body aches, or sore throat. Study participants meeting the criteria for acute respiratory illness attended a study research clinic at the University of Michigan School of Public Health, where a combined throat and nasal swab, or a nasal swab only for children less than 3 years old, was collected by the study team. Beginning in the 2014-2015 season, study participants with acute respiratory illnesses took an additional nasal swab at home at the time of illness onset, collected either by themselves or by a parent. The study was approved by the Institutional Review Board of the University of Michigan Medical School.

Viral detection, lineage typing, and viral load quantification.

We processed upper respiratory specimens (combined nasal and throat swab or nasal swab) for confirmation of influenza virus infection by reverse transcription-PCR (RT-PCR). We extracted viral RNA with either QIAamp viral RNA mini kits (Qiagen) or PureLink Pro 96 viral RNA/DNA purification kits (Invitrogen) and tested samples using the SuperScript III Platinum one-step quantitative RT-PCR system with ROX (Invitrogen) and primers and probes for universal detection of influenza A and B viruses (CDC protocol, 28 April 2009). Specimens positive for influenza virus were tested using subtype/lineage primer and probe sets, which are designed to detect influenza A (H3N2), A (H1N1)pdm09, B (Yamagata), and B (Victoria) viruses. An RNAseP primer/probe set was run for each specimen to confirm specimen quality and successful RNA extraction.

We quantified the viral load in each sample by RT-qPCR using primers specific for the open reading frame of segment 8 (NS1/NEP): forward primer 5′-TCCTCAACTCACTCTTCGAGCG-3′, reverse primer 5′-CGGTGCTCTTGACCAAATTGG-3′, and probe 5′-(FAM)-CCAATTCGAGCAGCTGAAACTGCGGTG-(BHQ1)-3′. Each reaction mixture contained 5.4 μl of nuclease-free water, 0.5 μl of each primer at 50 μM, 0.1 μl of ROX dye, 0.5 μl SuperScript III RT/Platinum Taq enzyme mix, 0.5 μl of 10 μM probe, 12.5 μl of 2× PCR buffer master mix, and 5 μl of extracted viral RNA. To relate genome copy number to threshold cycle (C_T) value, we used a standard curve based on serial dilutions of a plasmid control, run in duplicate on the same plate.

Amplification, library preparation, and sequencing.

We amplified viral cDNA from all eight genomic segments using the SuperScript III one-step RT-PCR Platinum Taq high-fidelity kit (Invitrogen). Each reaction mixture contained 5 μl of extracted viral RNA, 12.5 μl of 2× PCR buffer, 2 μl of primer cocktail, 0.5 μl of enzyme mix, and 5 μl of nuclease-free water. The primer cocktail was a mixture of B-PBs-UniF, B-PBs-UniR, B-PA-UniF, B-PA-UniR, B-HANA-UniF, B-HANA-UniR, B-NP-UniF, B-NP-UniR, B-M-Uni3F, B-Mg-Uni3F, B-M-Uni3R, B-NS-Uni3F, and B-NS-Uni3R (sequences and proportions are listed in reference 40). The thermocycler protocol was as follows: 45°C for 60 min, 55°C for 30 min, and 94°C for 2 min, followed by 5 cycles of 94°C for 20 s, 40°C for 30 s, and 68°C for 3 min 30 s, followed by 40 cycles of 94°C for 20 s, 58°C for 30 s, and 68°C for 3 min 30 s, and a final extension of 68°C for 10 min. We confirmed IBV genome amplification by gel electrophoresis. We sheared amplified cDNA (100 to 500 ng) on a Covaris ultrasonicator with the following settings: time, 80 s; duty cycle, 10%; intensity, 4; cycles per burst, 200. We prepared sequencing libraries with NEBNext Ultra DNA library prep kits (NEB) and sequenced them on an Illumina NextSeq with 2 × 150 paired end reads (mid-output run, v2 chemistry). To increase the specificity of variant identification, samples with a viral load between 10³ and 10⁵ genome copies/μl of transport medium were amplified and sequenced in duplicate. Samples amplified from B/Victoria and B/Yamagata plasmid clones were included on each sequencing run to account for sequencing errors. The plasmids used in the control reactions were generated by segment-specific RT-PCR from clinical samples of B/Victoria and B/Yamagata strains from the 2012-2013 season, followed by gel extraction and TOPO-TA cloning (Invitrogen). The sequence of each plasmid was determined by Sanger sequencing. We generated the plasmid control amplicons included on each Illumina sequencing run by using the same multiplex amplification protocol but with cloned plasmid DNA as the template.

Identification of iSNV.

Intrahost single-nucleotide variants (iSNV) were identified using a previously described analytic pipeline (24). We identified iSNV in samples that had an average genome coverage greater than 1,000× and a viral load greater than 10³ genome copies/μl of transport medium in the original sample. Sequencing adapters were removed with Cutadapt (41), and reads were aligned to the sequences derived from the B/Victoria and B/Yamagata plasmid controls with Bowtie2 (42). Duplicate reads were marked and removed with Picard and samtools (43). Putative variants were identified with the R package deepSNV using data from the clonal plasmid controls of each sequencing run (44). Minority iSNV (<50% frequency) were identified using the following empirically derived criteria: deepSNV P value, <0.01; average mapping quality, >30; average Phred score, >35; and average read position, in the middle 50% (positions 37 and 113 for 150-bp reads). For samples processed in duplicate, we used only variants that were present in both replicates; the frequency of the variant in the replicate with greater coverage at that site was used. Lastly, variants with a frequency of <2%, which have a higher false-positive rate from RT-PCR and/or sequencing errors, were not included in downstream analyses.

In our previous work on IAV, we found that there were multiple sites with mutations that were essentially fixed (>0.95) relative to those of the plasmid control and in which the base in the plasmid control was therefore identified as a minority variant in the sample (19). At these sites, deepSNV is unable to estimate the base-specific error rate and cannot distinguish true minority iSNV; however, we found that we could accurately identify minority variants at these sites at a frequency of 2% or above (19). This frequency threshold was incorporated into the pipeline for iSNV identification at these sites. Therefore, we identify intrahost variants with frequencies between 2 and 98%. Minority iSNV are the subset of these variants with a frequency between 2 and 50% relative to the sample consensus, which we use as a metric of within-host diversity. For each minority iSNV, we identify the majority iSNV present at a frequency of 50 to 98%. Any sites that were monomorphic after application of quality filters were assigned a frequency of 100%. Nucleotide diversity (π) was calculated using identified iSNV in each sample by the formula described by Zhao and Illingworth (45).

Data and code availability.

Raw sequence data, with human content filtered out, are available at the NCBI Sequence Read Archive under BioProject accession number PRJNA561158. Code for the variant identification pipeline is available at http://github.com/lauringlab/variant_pipeline. Analysis code is available at http://github.com/lauringlab/Host_level_IBV_evolution.