Simultaneous Viral Whole-Genome Sequencing and Differential Expression Profiling in Respiratory Syncytial Virus Infection of Infants

Abstract

Targeted metagenomics using strand-specific libraries with target enrichment is a sensitive, generalized approach to pathogen sequencing and transcriptome profiling. Using this method, we recovered 13 (76%) complete human respiratory syncytial virus (RSV) genomes from 17 clinical respiratory samples, reconstructed the phylogeny of the infecting viruses, and detected differential gene expression between 2 RSV subgroups, specifically, a lower expression of the P gene and a higher expression of the M2 gene in RSV-A than in RSV-B. This methodology can help to relate viral genetics to clinical phenotype and facilitate ongoing population-level RSV surveillance and vaccine development.

Clinical Trials Registration. NCT03627572 and NCT03756766.

respiratory syncytial virus, transcriptome profiling, RNA sequencing, whole-genome sequencing, differential gene expression

Methods to provide high-throughput, whole-genome sequences for specific viruses without a bespoke assay would have substantial value in epidemiological and evolutionary studies. For example, considering the global impact of human respiratory syncytial virus (RSV) infection, the lack of a simple, cost-effective assay means that relatively few complete sequences have been determined, and little is known about potential links between viral genetics and pathogenicity. In this study we demonstrate that a generalized approach to pathogen sequencing, which we term targeted metagenomics, can provide a quantitative sequencing assay for a chosen virus that, because it relies only on library construction and enrichment (2 processes that do not depend on the particular target pathogen), can be applied to new viruses without the requirement for optimization. In addition, the methodology enables virus transcriptome profiling to detect genotype-specific patterns of transcription of potential phenotypic significance.

RSV, an unsegmented, single-stranded, negative-sense RNA virus with a 15.2-kb genome, consisting of 10 genes, is the leading pathogen causing lower respiratory tract infection among young children, globally responsible for an estimated 33.1 million episodes of lower respiratory tract infection, 3.2 million hospitalizations, and up to 149 400 deaths in children < 5 years of age annually [1]. RSV is classified into 2 antigenically and phylogenetically distinct subgroups: RSV-A and RSV-B. A detailed understanding of RSV clinical infection at the molecular level could aid evaluation of viral responses to immune pressure from new vaccines or monoclonal antibodies that are under development. In this study we adopted a targeted sequencing approach, ve-SEQ [2], designed for the detection of pathogen sequences, and adapted existing computational pipelines for use in RSV genomics and transcriptomics.

METHODS

Detailed methods are provided in the Supplementary Methods

Clinical Studies and Samples

All available nasopharyngeal swabs collected from a cohort of RSV-infected infants < 1 year of age enrolled in the clinical studies of the Respiratory Syncytial Virus Consortium in Europe project (ClinicalTrials.gov identifier: NCT03627572 and NCT03756766) from Santiago de Compostela, Spain, and Oxford, United Kingdom, during the 2017–2018 RSV season, were included in this study. The protocol was approved by the local ethics committee at each center. The parents or guardians of all participants provided written, informed consent.

Targeted Metagenomic Sequencing

NucliSENS easyMAG (BioMérieux) was used to extract 25 μL of total nucleic acids from 500 μL of each swab sample, following the manufacturer’s instructions. Sequencing libraries were constructed from 8 μL of each nucleic acid sample, following the ve-SEQ protocol applied to human immunodeficiency virus (HIV) [3] with modifications. Briefly, dual indexed Illumina libraries were constructed using the SMARTer Stranded Total RNA-Seq Kit v2—Pico Input Mammalian (Takara Bio USA) without RNA fragmentation. A 2-μL aliquot of library from each sample was pooled for sequence enrichment by a probe set targeting over 100 bacteria and viruses, including RSV (Supplementary Table 1) [4]. Polymerase chain reaction (PCR) of 16 cycles was carried out for postcapture amplification and the final product was purified by Agencourt AMPure XP (Beckman Coulter).

Sequencing was performed on the Illumina MiSeq platform (Illumina), generating 265-bp paired-end reads.

Viral Load Measurement

Viral load was measured by quantitative reverse transcription PCR (RT-PCR) from the same nucleic acid samples used for sequencing. PCR assays were prepared in triplicate using the genesig RSV Standard Kit (Primerdesign), targeting the N gene, and performed on the Applied Biosystems StepOnePlus Real-Time PCR System (Thermo Fisher Scientific) with amplification conditions of 55°C for 10 minutes, 95°C for 2 minutes, followed by 50 cycles of 95°C for 10 seconds and 60°C for 60 seconds.

Viral Genome and Transcriptome Reconstruction

shiver [5] was used to reconstruct RSV genomes and transcriptomes from the de novo assemblies of sequencing reads generated by either IVA (v. 1.0.8) [6] or SPAdes (v. 3.13.0) [7] with the deduplication option enabled. A minimum of 2 unique reads were required to support consensus base-calls.

Each prededuplicated binary alignment/map (BAM) file was divided into 2 files: one containing read pairs mapped in the antisense orientation (genomic reads, SAM flags: 83 or 163), and the other in the sense orientation (positive-strand reads, SAM flags: 99 or 147). The 2 files were then deduplicated and processed using shiver to generate consensus sequences and per-base coverage and composition.

Phylogenetic Analysis

Maximum-likelihood phylogenies were reconstructed from complete RSV sequences generated in this study with another 60 RSV whole-genome sequences downloaded from GenBank (Supplementary Table 2) using RAxML (v. 8.2.12) [8] with the general time reversible model and gamma-distributed rate heterogeneity among sites.

Differential Expression Analysis

Samples with complete RSV genome assembly were used for differential expression analysis, following the RNAseq123 workflow [9]. Briefly, HTSeq (v. 0.11.2) [10] was used to calculate the number of positive-strand reads mapped to each RSV gene from prededuplicated BAM files. Linear models were constructed to compare gene expression levels between RSV subgroups, severity, and host sex groups from the normalized and transformed count data. Severity was defined using the ReSVinet scale [11]; a score of 0–7 was defined as mild, a score of 8–13, moderate, and a score of 14–20, severe. Finally, empirical Bayes moderated t statistics were applied to compute the significance of differential expression with a minimum log₂ fold change of 1.2.

Statistical Analysis

All statistical analyses were performed using R. P values or adjusted P values < .05 were considered significant.

RESULTS

We constructed RNA-seq libraries from 17 RSV-positive samples, including 6 from Oxford, United Kingdom and 11 from Santiago de Compostela, Spain (Supplementary Table 3). The enriched libraries yielded a median of 30% (0.02%–75%) RSV reads, originating from both the positive (transcriptome and replicative intermediate) and negative (genomic) strands at a ratio of 5–11:1 (Supplementary Table 4 and Supplementary Data). Complete RSV genomes were recovered from 13 (76%) of the samples, including all samples with viral load > 2000 copies/mL or that generated > 10 000 unique (deduplicated) RSV reads (Figure 1A), with a sensitivity, achieved without further optimization, that makes the methodology suitable for studying a wide variety of infections.

Figure 1.

A, Sequencing yield and virus load. Yield of unique, respiratory syncytial virus (RSV)-mapped reads is associated with viral load and proportion of genome assembly. Shaded area around the regression line represents the 95% confidence interval. B and C, Maximum-likelihood phylogenies of RSV subgroups A and B. Phylogenies were reconstructed using RAxML from consensus sequences of the 13 complete genomes, augmented with 30 reference genomes representing sampled diversity of each of RSV-A and RSV-B. Sequences from this study are formatted as GB (United Kingdom) or ES (Spain) followed by a 2-digit number and genotypes are colored according to the legend. Reference taxa are labelled with country and year of sampling. Scale bar shows nucleotide substitutions per site. The phylogenies were rooted to the oldest sampled RSV sequences, dating from 1977. Abbreviations: AU, Australia; BE, Belgium; BR, Brazil; CN, China; DE, Germany; ES, Spain; GB, United Kingdom; IT, Italy; JO, Jordan; KE, Kenya; NL, Netherlands; NZ, New Zealand; PE, Peru; PH, Philippines; TH, Thailand; US, United States.

Using phylogenetic analysis, we identified 4 RSV-A (genotype ON1) and 9 RSV-B (genotype BA9) samples among the 13 completely sequenced genomes (Figure 1B and 1C). In all but 1 case, the nearest neighbor sequence for each new genome was another sample from our study, although not necessarily from the same country, reflecting the ongoing evolution and limited number of cocirculating RSV lineages.

Read coverage plots revealed characteristic, strand-specific patterns of read depth along the genome and antigenome. In contrast to the genomic strand (Figure 2A, left), positive-strand read coverage was marked by read depths in intergenic regions approximately 1%–10% of those in coding regions (Figure 2A, right), recapitulating the known mechanism of RSV transcription, with a single initiation point at the 3′ end of the genome and polyadenylation/termination at each gene-end signal before reinitiating at the next gene-start signal [12]. Positive-strand reads spanning gene boundaries presumably represent antigenomes and read-through transcripts. A general pattern of decreasing coverage of successively more-5′ genes reflects the dissociation of a fraction of polymerase complexes before reinitiation [12], leading to progressively weaker transcription. Under this model, due to an overlap between the M2 and L genes, expression of the 5′-most gene, L, requires polymerase to terminate at the gene-end signal of the preceding M2 gene, and then to move backwards to the L gene-start signal before reinitiation [12]. This inference is supported by a jump in read depth at the start of the L gene (Figure 2B). However, in light of this basic model, a more detailed examination of relative expression levels also provides indirect support for more complex, posttranscriptional control of transcript abundance: rather than uniformly decreasing transcript levels, transcripts from the NS2, P, and G genes were more abundant than those from the genes preceding them (Figure 2C), in contrast to patterns of transcription in cell-free systems [12].

Figure 2.

Coverage plots and differential gene expression. A, Read coverage of complete RSV genomes (left) and transcriptomes (right). Plots are normalized read depth per million mapped reads in each category. Legend below the plots represents the position and length of each RSV gene. B, Positive-strand read coverage across the M2 gene. There is little evidence for read depth dropping at the end of the M2 gene; instead coverage drops within the ORF M2-2, which is inconsistent with a single transcript covering both ORFs M2-1 and M2-2. Shaded areas represent gene or ORF overlaps. C, Transcription by gene and sample. Read counts are normalized by total read count per sample and gene length. D, Differential expression of M2 and P genes and ORFs M2-1 and M2-2. RSV-A has a higher expression of the M2 gene (fold change = 2.38; 95% CI, 1.44–3.95; adjusted P = .032) and a lower expression of the P gene (fold change = 0.57; 95% CI, .47–.70; adjusted P = .009) than RSV-B. Within the M2 gene, RSV-A has a higher expression of M2-1 than RSV-B (fold change = 2.29; 95% CI, 1.49–3.53; adjusted P = .017), while the expression of M2-2 is comparable between RSV-A and RSV-B (adjusted P = .703). Center line of each box denotes the median; box limits, first and third quartiles; whiskers, the highest and lowest values within 1.5 times the interquartile range from the box limits; and outlying points, outliers. Abbreviations: CI, confidence interval; CPM, counts per million; ORF, open reading frame.

Noting marked variation among isolates in the proportions of positive-strand reads assigned to each gene (Supplementary Figure 1), we undertook a formal analysis of differential gene expression. RSV subgroup, but not host sex or disease severity, was significantly associated with differences in specific transcript levels (Supplementary Figure 2). The mean level of M2 transcripts in RSV-A (samples GB-01, -03, -05, and -06) was 240% of that in RSV-B (GB-02, -04, and all Spanish isolates), while the expression of the P gene in RSV-A was 57% of that in RSV-B (Figure 2D). While intersample variation was marked, perhaps reflecting clinical heterogeneity between samples, other genes were not significantly differentially expressed.

The M2 gene is known to contain 2 overlapping open reading frames (ORFs): M2-1, whose product interacts with the P gene product (phosphoprotein) to induce antitermination (promoting transcription of full-length subgenomic mRNAs), and M2-2, whose protein conversely acts to favor production of antigenomes over full-length transcripts [12]. Under the accepted model of RSV molecular biology, both ORFs belong to a single transcriptional unit. Given their contrasting roles, we scrutinized our transcription data more carefully to see whether the change in M2 transcription affected both ORFs. In our data (Figure 2B), (1) read depth corresponding to transcription drops by an aggregate 0.5–1 log-fold downstream of M2-1 but within M2-2; (2) levels of transcription in the remainder of M2-2 are similar to intergenic levels in the rest of the genome; and (3) the signal for enhanced transcription of the M2 gene in RSV-A is confined to the M2-1 region (Figure 2D).

DISCUSSION

In this study, we demonstrate the simultaneous recovery of RSV genomes and transcriptomes directly from clinical samples, without a prerequisite for culture or purification of viral particles. Using strand-specific libraries with target enrichment, we were able to examine separately the RSV genome and transcriptome sequences, and to relate viral gene expression to the phylogenetic subgroup of the infecting virus. We identified intriguing differences between RSV-A and RSV-B transcriptional profiles, which could provide novel pathogenetic and therapeutic insights.

The nonconcordant expressions of the ORFs M2-1 and M2-2 lead us to hypothesize that the 2 ORFs may represent separate transcriptional units, and that the M2-2 protein may be translated from a low-abundance RNA of unknown length (including possibly read-through of M2-1 mRNA or genomic intermediate) and to conclude that the M2 differential expression signal relates specifically to M2-1. Because the M2-1 protein facilitates the synthesis of polycistronic read-through mRNAs [12], a testable hypothesis consistent with our transcriptomics result is that the RSV-A samples should have less marked declines in positive-strand read depth within intergenic regions than RSV-B due to more read-through mRNAs.

Differential expression of the P gene between subgroups is also intriguing because the direct interaction between the M2-1 and P proteins [13] means that the contrasting changes in gene expression could combine in complex ways. For example, a skew towards transcription in RSV-A could foreseeably affect infection at both the cellular and whole-host levels, and may lead to generally higher prevalence of RSV-A in annual epidemics and proposed greater transmissibility of RSV-A than RSV-B [14]. However, during the study period of 2017–2018, RSV-B was predominant over RSV-A, accounting for 72% of the RSV samples collected in our multicenter study. Although no study has found differences in transcriptional machinery that would potentially cause the differential expressions between subgroups, we speculate that sequence differences or polymerase dissociation in intergenic regions may play a role.

Limitations of this study include the modest sample size, but our study design was robust to technical and analytical bias (see Supplementary Methods). Another limitation is that we did not perform quantitative PCR to support the findings of differential expression. Nevertheless, RNA-seq has been shown to be highly accurate for quantifying transcriptomes, providing estimates of gene expression levels closely correlated with measurements using quantitative PCR [15].

More generally, our methodology facilitates population-level RSV surveillance in the context of current efforts in global RSV surveillance and vaccine development. With our approach, it is straightforward to relate viral load, sequence variation, and transcriptional profiles from sequencing data in order to investigate the epidemiological and phenotypic consequences (eg, host biomarkers) of virus evolution. Targeted metagenomics provides a powerful methodology for characterizing sets of pathogens that may co-occur in a single sample type.

Supplementary Data

Supplementary materials are available at The Journal of Infectious Diseases online. Consisting of data provided by the authors to benefit the reader, the posted materials are not copyedited and are the sole responsibility of the authors, so questions or comments should be addressed to the corresponding author.

Notes

RESCEU Investigators. Harish Nair, Harry Campbell, Steve Cunningham (University of Edinburgh); Philippe Beutels (Universiteit Antwerpen); Louis Bont, Joanne Wildenbeest (University Medical Centre Utrecht); Andrew Pollard, Christopher Butler, Matthew Snape, Simon Drysdale, Gu-Lung Lin, Daniel O’Connor, Elizabeth Clutterbuck, Kimberley Jefferies, Joseph McGinley (University of Oxford); Peter Openshaw, Ryan Thwaites, Dexter Wiseman (Imperial College London); Federico Martinon-Torres, Alberto Gómez-Carballa (Servicio Galego de Saude); Terho Heikkinen (University of Turku and Turku University Hospital); Adam Meijer (National Institute for Public Health and the Environment); Thea Kølsen Fischer (Statens Serum Institut); Maarten van den Berge (University of Groningen); Carlo Giaquinto (PENTA Foundation); Michael Abram (AstraZeneca); Kena Swanson (Pfizer); Amanda Leach (GlaxoSmithKline); Clarisse Demont, Scott Gallichan (Sanofi Pasteur); Jeroen Aerssens, Deniz Öner (Janssen); Brian Rosen (Novavax); Eva Molero (Team-It Research).

Author Contributions. G.-L. L., T. G., D. B., R. B., and A. J. P. conceived and designed the work. G.-L. L, S. D., K. J., J. A., L. B., P. O., F. M.-T., and A. J. P. conducted and supervised the clinical studies. M. A. A. designed the probe set that was used for capture. G.-L. L., A. B., and M. d. C. performed the experiments. G.-L. L., T. G., D. O’C., and R. B. analyzed and interpreted the data. G.-L. L. drafted the manuscript and T. G., R. B., and A. J. P. substantively revised it. All authors have approved the submitted version and agreed to submit the manuscript.

Disclaimer. The views expressed in this article are those of the authors and may not be understood or quoted as being made on behalf of or reflecting the position of the organizations with which the authors are employed/affiliated.

Financial support. This work was supported by the National Institute for Health Research (NIHR) Oxford Biomedical Research Centre; the NIHR Thames Valley and South Midlands Clinical Research Network,; the British Research Council; and the Respiratory Syncytial Virus Consortium in Europe (RESCEU) project. RESCEU has received funding from the Innovative Medicines Initiative 2 Joint Undertaking (grant number 116019). This Joint Undertaking receives support from the European Union Horizon 2020 Research and Innovation Program and European Federation of Pharmaceutical Industries and Associations.

Supplement sponsorship. This supplement is sponsored by RESCEU (REspiratory Syncytial Virus Consortium in EUrope).

Potential conflicts of interest. J. A. is an employee of Janssen Pharmaceutica NV. All other authors report no potential conflicts. All authors have submitted the ICMJE Form for Disclosure of Potential Conflicts of Interest. Conflicts that the editors consider relevant to the content of the manuscript have been disclosed.

Presented in part: 7th Applied Bioinformatics and Public Health Microbiology Conference, Hinxton, UK, 5–7 June 2019.

Data availability. All raw sequencing reads are accessible at the European Nucleotide Archive under the study accession PRJEB34042 (https://www.ebi.ac.uk/ena/data/view/PRJEB34042). The complete RSV sequences generated in this study have been deposited in GenBank with the accession numbers LR699315, LR699726, LR699734, and LR699736–44.

Author notes

G.-L. L. and T. G. contributed equally to this work.

R. B. and A. J. P. contributed equally to this work.

Members of the study group are listed at the end of the text.