Genomic analysis of respiratory syncytial virus infections in households and utility in inferring who infects the infant

Study location, design and samples

The study was undertaken within Kilifi County, which is located in coastal Kenya. A detailed description of the study location and study design was provided elsewhere²⁸. Briefly, 47 households scattered across an area of approximately 21 km² were followed up over a 6-month period beginning December 2009 and ending June 2010 coinciding with the RSV peak activity months in the area²⁹. Households were defined as a group of people sharing a compound and eating from the same kitchen²⁰. The selected households (abbreviated HHs) were given designated identifiers from 1 to 57 (HH01 to HH57). Twice weekly throughout the study period, a nasopharyngeal-flocked swab (NPS) was obtained from each member regardless of symptom status. The NPS samples were screened for RSV using a multiplex real-time RT-PCR method which subtyped RSV positives into RSVA and RSVB³⁰. For whole genome sequencing (WGS), we targeted 20 select households that documented RSV infection of ≥2 members. A geographical map showing the distribution of the study households is provided in the Additional File: Fig. S1.

RNA extraction, amplification and whole nucleotide sequencing

Viral RNAs from the positive samples of selected households were obtained using the QIAamp viral RNA extraction Kit (QIAGEN, Hilden, Germany) following the manufacturer’s instructions. Complementary DNA (cDNA) synthesis and RSV whole genome amplification was achieved using a six-overlapping PCR fragments strategy (each ~2.5 kb) as previously described²¹. Sequencing libraries were prepared using Nextera DNA Library Prep kits and nucleotide sequencing performed using Illumina MiSeq platform multiplexing 15–20 samples per run to generate approximately 1 million paired-end reads (150 bp × 2) per sample²¹.

Whole genome sequence assembly and multiple sequence alignments

The short reads from the MiSeq instrument were de-multiplexed, quality checked (median read Phred score of ≥35) and trimmed using QUASR v6.08³¹. Reads passing quality checks were de novo assembled into longer contigs using the SPAdes v3.5.0³². RSV contigs were identified by matching to a database of RSV sequences using USEARCH program³³, examined for completeness of the expected open reading frames (ORF) using Geneious v8.1.6 (https://www.geneious.com) and, where necessary, partial contigs were further combined to longer ones using Sequencher v5.0.1³⁴. These were subsequently checked presence of intact ORFs, sorted by household, re-aligned and positions of nucleotide variation double-checked if these were supported by majority of the raw reads associated with that sample¹⁰. Multiple sequence alignments were prepared in MAFFT v7.220³⁵.

Phylogenetic analysis

Sequence phylogenies were inferred using Maximum Likelihood (ML) methods in MEGA7³⁶ and RAxML v8.2.12³⁷. The best-fitting models of nucleotide substitution for each alignment were in IQ-TREE v1.4.3³⁸. Best tree search was performed by Nearest Neighbor Interchange (NNI). Branch support was evaluated by bootstrapping with 1,000 replicates. Pairwise genetic distances were calculated in pairsnp 0.0.6³⁹. The phylogenetic relatedness of the household RSVA and RSVB genomes was assessed at three levels; (i) in combination with global sequences deposited in GenBank (RSVA, n = 657 collected between 1977–2015 while RSVB, n = 416 collected between 1978–2016), (ii) among the households viruses alone and (iii) among viruses collected from same households only. The potential transmission networks within and between households for each group were inferred in PopART package v1.7.2 using median joining tree (MJT) method with an epsilon of zero⁴⁰. Evolutionary analyses were determined in maximum-likelihood-based TreeTime program⁴¹. Phylogenetic trees were visualized and annotated in FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).

Identifying who infected the household infant (s)

Infants were defined as the participants aged <1-year-old during the study²⁰. We grouped the other participants into 5 age-groups: (i) toddlers (1 to <3 years), (ii) pre-schoolers (3 to <6 years), (iii) school-aged (6 to <12 years), (iv) adolescents (12 to <18 years) and adults (>18 years). We attempted to identify who among these other age-groups were the most likely infectors of the infants by examining the relatedness of the virus genome(s) obtained from the infant to all the viral genomes obtained from the other members in the same household. Information on the dates of the sampling of the sequenced samples was taken into account to position the infant in the transmission network/chain.

Sequence nomenclature and accession numbers

The sequence nomenclature of the household samples has four digits that include the household identifier (first two digits) and subject identifier (the last two digits). All the new 112 full or partial RSVB genome sequences from this study were deposited in GenBank under the accession numbers MH594350 – MH594461. The RSVA genomes are deposited in GenBank under accession numbers KX510136-KX510266.

Ethical approval

The samples were collected after obtaining informed written consent from each participant if aged ≥18 years or through a guardian or parent if aged <18 years. In addition, children aged above 5 years were asked for assent. The study protocol approved by both the Scientific and Ethics Review Unit (SERU) of the Kenya Medical Research Institute (KEMRI), Nairobi, and Coventry Research Ethics Committee of UK²⁰. All study procedures were performed in accordance with the approved protocol guidelines and in compliance with the relevant regulations.

RSV infections and whole genome sequencing

We targeted 20 households with a total of 226 occupants (range 4–37 persons per household) for WGS. Details of the demographic characteristics of the analysed households, total specimens collected, diagnostic results, genome sequencing success and the observed phylogenetic clades (defined later) are summarized in Table 1. Over the six-month period (December 2009 – June 2010), a total of 7,695 nasopharyngeal-flocked swabs (NPS) were collected from the 20 HHs, 415 (5.4%) of which were determined to be RSV real-time RT-PCR positive (cycle threshold (Ct) value of <35.0; 189 RSVA, 214 RSVB and 12 RSVA/B co-infections) these originating from 130 participants. Of the 415 positive specimens, 374 (90.1%) samples were processed for WGS²¹ with successful amplification and assembly of RSV contigs of >1000 nucleotides length in 246 samples (65.7%). Of these 191 samples (51.1%) yielded contigs >14000 nucleotides (nt) (103 RSVA and 88 RSVB i.e. >90% of RSV full-length genomes) hereafter referred to as genomes. In eight and 14 HHs, two or more RSVA or RSVB genomes were recovered, respectively, allowing our investigation into within-household RSV transmission and variation, Table 1. Genome sequencing success negatively correlated with increasing diagnostic RT-PCR Ct value. These results, together with details on the metadata of the sequenced RSVB viruses, GenBank and Sequence Read Archive accession numbers and assembly metrics are provided in the Additional File and Supplementary Dataset, respectively.

Diversity of the viruses isolated in the study

From G gene phylogeny, all RSVA and RSVB viruses sequenced were genotypes GA2 and BA, respectively (results not shown). The genome-based maximum likelihood (ML) phylogenetic trees are shown in Fig. 1. The RSVA genomes formed a single monophyletic cluster on the global phylogeny while household RSVB genomes formed 5 distinct phylogenetic clusters interspersed with sequences from other global locations, Fig. 1, panel a. On their own, both RSVA and RSVB household genomes formed multiple phylogenetic clusters (several apparently genetically distinct and supported by >60% bootstrap values and we later assigned these into clades and sub-clades – see below). On the household genomes only ML tree, these clusters appeared to be mostly household specific with a few exceptions, Fig. 1, panel b.

The time-resolved ML trees and temporal signal in nucleotide divergence of the household RSVA and RSVB viruses are shown in Fig. 2. The time to Most Recent Common Ancestor (tMRCA) estimate of all the sampled RSVB viruses was estimated to be December 2004 (Lower and Upper boundaries for 90% Highest Posterior Density (HPD) of October 1997 and September 2008), which is much earlier than the equivalent estimate for RSVA viruses of December 2008 (Lower and Upper boundaries for 90% HPD of January 1987 and December 2009), with both point estimates showing a wide uncertainty interval.

We quantified the genetic diversity observed within the two RSV groups by calculating the number of pairwise single nucleotide polymorphisms (SNPs) (pairwise distance) of viruses within the same group, Fig. 3, panel a. We found this value to range from 0–35 (median: 19, mean: 16.6) for RSVA and 0–177 (median: 134, mean: 99.7) for RSVB. Overall within-group pairwise distances among RSVB viruses were 6.5 times higher than those of RSVA (mean distance of 0.006094 vs 0.001065). The distribution of the number of pairwise SNPs within clusters of the household viruses observed on the global phylogeny are shown in Fig. 3, Panel b–f.

To facilitate further analysis, we assigned the household viruses into “clades” and “sub-clades” defined by both their clustering patterns on global phylogenies (Fig. 1, panel a), the inferred divergence dates of the strains (Fig. 2, panel a) and, the number of pairwise SNP (Fig. 3). We grouped viruses in the same clade if they occurred as a monophyletic group on the global phylogeny, had <60 pairwise SNPs across the genome with every other member of that clade and diverged more than a year prior to their date of collection. Viruses within the same clade were further assigned into sub-clades if they showed >10 pairwise SNPs differences across the genome and were estimated to have diverged more than six months prior to their date of collection (Figs 2 and 3). Using these criteria, we assigned all household RSVA strains into a single clade named RSVA/I while household RSVB strains were assigned into 5 clades named RSVB/I through RSVB/V. Viruses within clade RSVA/I were assigned into five sub-clades; RSVA/Ia through RSVA/Ie, viruses within RSVB/I clade were assigned into two sub-clades RSVB/Ia and RSV/Ib, and viruses within RSVB/II were assigned into two sub-clades RSVB/IIa and RSV/IIb.

Virus transmission within and between households

We investigated the genomics and temporal and spatial patterns of RSVA and RSVB virus clades observed within and between households. An analysis using minimum spanning network which depicts shared differences without regard to an evolutionary model was used to detect patterns in the RSVA and RSVB genomes and examine potential intra- and inter-household transmission patterns (Fig. 4, panel a). Similar to the ML phylogenies, the majority of the viruses clustered by household with the major clusters corresponding to the clades and sub-clades observed in the ML trees. Notably clades/sub-clades RSVA/Ia, RSVA/Ie, RSVB/Ia, RSVB/Ib, RSVB/IIa, RSVB/IIb, and RSVB/IV were observed in multiple HHs indicating potential transmission linkage of the involved HHs during the epidemic. In the timeline of viruses identified (Fig. 4, panel b), all except five households (HH06, HH26, HH38, HH41 and HH42) had a single RSV clade sequenced. The exceptional households had two virus clades infecting members but mostly one of the two clades predominated e.g. in HH06, HH41 and HH42. On the other hand, in the remaining two households distinct RSVA and RSVB outbreaks occurred: HH38 in which the first outbreak was RSVB/I and at a later date a second outbreak of RSVA/I, and HH26 with concurrent RSVA/I and RSVB/IV.

The relationship between the geographical distance between the households and the RSVA and RSVB clades that circulated in these households is shown in Fig. 4, panel c. Paradoxically, some of the households that were in very close proximities experienced infections with viruses from different clades or sub-clades e.g. HH41 and HH42 were <30 meters apart, yet none of the virus clades circulating in these 2 households were shared (Fig. 4, panel c). In contrast, HH35 and HH38, separated by a distance of ~3 kilometres, shared the same virus clade (RSVB/Ia). There was no apparent correlation between inter-HH distance and genetic relatedness or between sampling dates and virus transmission, i.e. no correlation between geo-temporal-spatial patterns of virus transmission within and between households.

Intra-host, inter-host and inter-house virus variation

The SNP abundance in samples collected from same host during repeat visits and in presumed single household outbreaks are shown in Figs 5 and 6. Overall, intra-host SNPs ranged from 0–17 (median: 1, mean: 1.75, Fig. 5) while intra-household SNPs ranged from 0 to 21 (median: 3, mean: 6.2, Fig. 6). Nucleotide changes were, in general, rare intra-host during the shedding period of a presumed single episode. When changes were evident, they were usually multiple SNPs occurring simultaneously and mostly affecting the last few positive samples collected from the subject. For nine subjects who remained virus positive for more than 21 days, we compared the recovered genome sequences to determine if these represented more than one infection (Fig. 7, panel a and b). Four of these individuals showed zero change despite the sequenced samples spanning a period of over a month. For the individuals that showed SNPs, these were few (<6 SNPs). In the intra-household analysis, it appeared that the households with a higher number of SNPs (>5 i.e. falling in the upper quartile) may have experienced multiple introductions of viruses from the same clade or sub-clade e.g. in HH26 for RSVA (see Additional File: Fig. S10 sample from 2605 collected on 26-Mar-2010), HH38 for RSVB (see Fig. 8, sample from 3803 collected on 19-Feb-2010).

To track independent viruses that were either introduced from elsewhere into the study area during the epidemic or were local but diverged outside the 2009/10 season we coined the word “epidemiological strain”. Genetically, viruses referred to as same epidemiological strain had <10 SNPs across their genomes and belonged to the same clade and sub-clade (where assigned). In total, we identified 12 epidemiological strains (five within RSVA and seven within RSVB) that occurred in the study area during the six-month surveillance, eight (66.7%) of which were observed in multiple households while four were found in a single household. For the epidemiological strains that occurred in multiple households, between 5–33 (median:12, mean 15.3) SNPs were observed across their genomes. A comparison of SNP abundance intra-host, inter-host and inter-household is provided in Fig. 7, Panel c. SNP abundance appeared to increase linearly across these three levels.

Who infected the infant(s) in the study households?

There were 22 infants from the 20 HHs. By our diagnostics, the infant in HH18 did not get RSV infected during our surveillance period. The household-by-household time-resolved infection patterns, genome alignments, phylogenies and minimum spanning sequence networks are provided in the Additional Files S3–21. We present the infection and genomic patterns of HH38 as an example in Fig. 8. Patterns of RSVA infection in HH14 and HH38 can be found in our previous publication¹⁰. Following examination of the patterns from all the 20 households, the summary of our deductions on who most likely infected the infant is provided in Table 2. Overall, we could infer the single most likely individual to infect the study infant for only 19% (4/21) of the infants, Table 2. For a further 19% we identified the top two individuals who most likely to have infected the infant. Note that in HH38, the infant was infected in both the RSVA and RSVB outbreaks that occurred in this household. All except one of the suspected infant infectors were aged <12 years-old.