Abstract
At the end of 2019 Wuhan witnessed an outbreak of “atypical pneumonia” that later developed into a global pandemic. Metagenomic sequencing rapidly revealed the causative agent of this outbreak to be a novel coronavirus - SARS-CoV-2. Herein, to provide a snapshot of the pathogens in pneumonia-associated respiratory samples from Wuhan prior to the emergence of SARS-CoV-2, we collected bronchoalveolar lavage fluid samples from 408 patients presenting with pneumonia and acute respiratory infections at the Central Hospital of Wuhan between 2016 and 2017. Unbiased total RNA sequencing was performed to reveal their “total infectome”, including viruses, bacteria and fungi. Consequently, we identified 37 pathogen species, comprising 15 RNA viruses, 3 DNA viruses, 16 bacteria and 3 fungi, often at high abundance and including multiple co-infections (12.8%). However, SARS-CoV-2 was not present. These data depict a stable core infectome comprising common respiratory pathogens such as rhinoviruses and influenza viruses, an atypical respiratory virus (EV-D68), and a single case of a sporadic zoonotic pathogen – Chlamydia psittaci. Samples from patients experiencing respiratory disease on average had higher pathogen abundance than healthy controls. Phylogenetic analyses of individual pathogens revealed multiple origins and global transmission histories, highlighting the connectedness of the Wuhan population. This study provides a comprehensive overview of the pathogens associated with acute respiratory infections and pneumonia, which were more diverse and complex than obtained using targeted PCR or qPCR approaches. These data also suggest that SARS-CoV-2 or closely related viruses were absent from Wuhan in 2016-2017.
Introduction
The emergence of COVID-19 at the end of 2019 has had a profound impact on the world. The causative agent, a betacoronavirus (Coronaviridae) termed SARS-CoV-2, has high transmissibility, causes mild to severe respiratory symptoms in humans1-3. The first documented case of SARS-CoV-2 was reported in Wuhan, Hubei province, China3. Despite intensive research into SARS-CoV-2 and COVID-19, important questions remain unclear regarding the emergence of this virus, including for when and where it appeared and how long it was circulating in human populations prior to its initial detection in December 2019.
Meta-transcriptomics has several advantages over traditional diagnostic approaches based on serology or PCR4: it targets all types of micro-organisms simultaneously, identifies potential pathogens without a priori knowledge of what micro-organisms might be present, reveals the information (RNA) expressed by the pathogen during infection that is central to agent identification and studies of disease association. This method has been proven highly successful in revealing the entire virome and microbiome in a diverse range of species5-8, including the initial identification of SARS-CoV-2 from patients with severe pneumonia1.
Acute respiratory infections and pneumonia are a significant public health concern on a global scale. However, far less is known about the total “infectomes” associated with respiratory infections and pneumonia. Herein, we report total infectome surveillance of 408 patients presenting with pneumonia and acute respiratory infections at Wuhan Central Hospital prior to the SARS-CoV-2 epidemic. The purpose of this study was to use an un-biased meta-transcriptomics tool to characterize the total infectome within these patients. Nevertheless, since the sampling period occurred before the outbreak of COVID-19, this represents the first opportunity to characterize the entire range of pathogens simultaneously within a cohort and determine the microbial composition of the population in which SARS-CoV-2 was initially reported.
Results
Patient context
We considered 408 patients clinically diagnosed with pneumonia or acute respiratory infection at the Central hospital of Wuhan in Wuhan, China. The sampling lasted for 20 months and covered the period between May 2016 to December 2017, two years before the onset of the COVID-19 pandemic (Fig. 1A). The male-to-female ratio among the patients sampled was 1.4 (Fig. 1B), with age ranging from 16 to 90 years (medium, 62). Pre-existing medical conditions present in these patients included hypertension (n=108), diabetes (n=46), bronchiectasia (n=31), chronic obstructive pulmonary disease (COPD, n=23), cancer (n=12), and heart disease (n=10). Based on evaluations made by clinicians at the hospital, 27 patients were described as severely ill, with 381 presenting with non-severe syndromes (Fig. 1E). The mortality for the entire cohort was 0.74% (n=3) and the average duration of hospitalization was 8 days (range 2-322, medium 9).
(A) Sampling frequency and intensity during the study period. (B) Male-to-female ratio of the recruited patients. (C) Age structure of all patients involved in this study depicted using violin plots. (D) Type and frequencies of pre-existing conditions. (E) The number of severe, non-severe and control cases.
Total infectome
Meta-transcriptomic analysis of the BALF samples identified a wide range of RNA viruses, DNA viruses, bacteria and fungi. For the purposes of this study, we only characterized those likely associated with human disease (i.e., pathogens). This included (i) existing species that are known to be associated with human disease, and (ii) potentially novel pathogens that have not been previously characterized. For the latter, we only considered DNA and RNA viruses that are related to a virus genus or family that have previously been shown to infect mammals and are at relatively high abundance level (i.e., >0.1% of total RNA, or 1000 RPM). Other than new pathogens, the abundance threshold for pathogen positives was set as 1 RPM. Furthermore, commensal bacteria population was not considered here.
Based on these criteria we did not identify any potential novel viral pathogens. All the microbes identified belonged to those previously characterized as human pathogens, comprising 15 RNA viruses, 3 DNA viruses, 16 bacteria and 3 fungal pathogens (Fig. 2). The case positive rate for all pathogens was 71.25% (n=249, Fig. 2A), many of which were only associated with RNA viruses (27.3%, n=151) or bacteria (28.9%, n=160). Co-infection with two different pathogens was also commonplace, comprising a total of 71 (12.8%) cases (Fig. 2A). Among the pathogens identified, most were common respiratory pathogens such as influenza viruses, rhinoviruses, Pseudomonas aeruginosa and Haemophilus influenzae (Fig. 2B). In addition, we identified a number of unconventional respiratory pathogens that are often not included in respiratory pathogen screening panels but known to cause severe infections in respiratory tract or lungs, including enterovirus D68 and Chlamydia psittaci (see below).
(A) Proportion of cases infected with RNA virus, DNA virus, bacteria, fungus and with mixed infections. (B) Prevalence of each pathogen, ordered by the number of cases. (C) Heat map showing the prevalence and abundance of pathogens in diseased and control samples. Escherichia coli is not regarded as pathogen but is shown here as an example of non-pathogen contamination. The samples (x-axis) are divided into “Patient” and “Control” groups, each ordered chronically. The pathogens (y-axis) are divided into four categories: RNA viruses, DNA viruses, bacteria, and fungi.
Finally, none of the pathogens described here appeared in the blank controls. Since the blank control samples were generated using the same procedures for RNA extraction, library preparation and sequencing as the experimental groups, these results effectively exclude the possibility that the pathogens described above were of contaminant origin.
Viruses
RNA viral pathogens exhibited both great diversity (15 species) and abundance (up to 52% of total RNA) in the BALF samples examined here. The most frequently detected RNA viruses were human rhinoviruses A-C (HRV, n=55), followed by influenza A virus (IAV, n=29), human parainfluenza virus type 3 (HPIV3, n=20), influenza B virus (IBV, n=8), and human metapneumovirus (HPMV, n=7) (Fig. 2B). While the majority were found throughout the study period, a few had more specific time-scales (Fig. 2C). For example, influenza B viruses were mostly identified in 2017, whereas enterovirus D (ENV-D, n=6) was only detected in the summer of 2016. In addition, we identified all four types of common cold associated coronaviruses - OC43 (n=4), HKU1 (n=4), 229E (n=6) and NL63 (n=1) - all of which had a relatively low prevalence in our cohort. Importantly, none of the libraries contained any hit to SARS-CoV or SARS-CoV-2, these results were confirmed by both read mapping and blasting against corresponding virus genomes.
In comparison to RNA viruses, the DNA viruses identified were limited in diversity and abundance. All three major types of human herpesviruses were identified - HSV1 (n=3), CMV (n=3), and EBV (n=2) - although with low abundance in all cases (up to 30 RPM or 0.003%) (Fig. 2C). Another common DNA virus that causes respiratory disease – adenovirus - was also identified in several cases, although at abundance levels lower than the 1 RPM threshold such that it was considered a ‘negative’ result in this context.
Bacteria and fungi
The most common bacterial pathogens identified included Pseudomonas aeruginosa (n=31), Haemophilus influenzae (n=19), Staphylococcus aureus (n=18), and Mycoplasma pneumoniae (n=11), all of which are common respiratory pathogens. Acinetobacter bacteria were also prevalent in our cohort, including Acinetobacter baumannii (n=5) and A. pittii (n=1), both of which are commonly associated with hospital-acquired infections. Other important respiratory pathogens, such as Mycobacterium tuberculosis (n=2), Legionella peumophila (n=2), Streptococcus pneumoniae (n=2), Klebsiella pneumoniae (n=3) and Moraxella catarrhalis (n=3), were also detected, although at a relatively low prevalence. Of particular interest was the identification of a single case of Chlamydia psittaci – a potentially bird-associated zoonotic pathogen – present in the BLAF at relatively high abundance (6396 RPM). Conversely, all the fungal pathogens identified here – Candida albican (n=3), C. tropicalis (n=1), and Aspergillus spp. (n=1) – were common pathogens known to cause respiratory infections.
qPCR confirmation of pathogen presence and abundance
For each of the pathogens identified here, we performed a qRT-PCR assay to confirm their presence and validate their abundance level as measured using our meta-transcriptomic approach. Strikingly, strong correlations were observed between the abundance measured by qPCR (i.e., CT value) and those estimated by read count after log 2 conversion (−0.8 < Pearson’s R < -1, Fig. 3). Hence, the quantification by the two methods is strongly comparable.
Abundance by qPCR methods is measured by cycle threshold, or CT value, while those by meta-transcriptomics are measured by “log2(number of reads)”. The comparisons are performed on four most abundant pathogens: influenza A virus, influenza B virus, human rhinoviruses, and human parainfluenza virus 3. Pearson’s correlation coefficient between CT values and log2(number of reads) is estimated for each pathogen.
In-depth phylogenetic characterization of pathogens
Although no novel viral pathogens were identified in this study, those viruses detected were characterized by substantial phylogenetic diversity, reflected in the presence of multiple viral lineages that highlight their complex epidemiological history in Wuhan (Fig. 4). We identified more than 14 genomic types of rhinovirus A, 7 of rhinovirus B, and 4 of rhinovirus C. A similar pattern of the co-circulation of multiple viral lineages was observed in other viruses. For example, the influenza A viruses identified in this study can be divided into the H1N1 and H3N2 subtypes, each containing multiple lineages that clustered with viruses sampled globally and reflecting the highly connected nature of Wuhan (Fig. 4). In this context it was striking that in the case of several coronaviruses - OC43, HKU1 and 229E - the Wuhan sequences grouped directly with those identified from the United States (Fig. 4), although this may reflect limited sampling.
Phylogenetic trees were estimated using the maximum likelihood method implemented in PhyML. Sequences identified from this study were marked with red solid circle. For larger trees, we only show lineages or sub-lineages which contain sequences identified in this study.
Pathogen presence and abundance in diseased and healthy individuals
Our meta-transcriptomic analysis revealed that many RNA viruses and bacteria detected were present at extremely high abundance levels (>1%, and up to 52% of total RNA) and hence likely indicative of acute disease. This was particularly true of eight species of RNA viruses – EV-D68, the influenza viruses, HRV, HPIV3, 229E – as well as two species of bacteria (Haemophilus influenzae and Pseudomonas aeruginosa) (Fig. 5). Together, these comprise a total of 54 cases (13.2% of total diseases cases).
Boxplots for patient and control groups for each pathogen identified in this study, including RNA virus (green), DNA virus (yellow), bacteria (orange) and fungus (blue). For clarity, only non-zero abundance levels are reported.
In marked comparison, high levels of abundance were never observed in the healthy control group (Fig. 5). The highest abundance was 1302 RPM (0.013% of total RNA) in this group was for Haemophilus influenzae. Indeed, for most of the pathogens, the abundance level in healthy group was either undetectable or well below that observed in the diseased group, with the exception of Human parainfluenza virus type 1 (HPIV) and Mycoplasma orale for which the abundance levels were higher in the control group, although the sample size for both pathogens was relatively small. Conversely, commensal microbes, particularly Escherichia coli, exhibited similar abundance levels in the disease and control groups.
Discussion
Our study provides a critical snapshot of the respiratory pathogens present in Wuhan prior to the emergence of SARS-CoV-2. Our unbiased metagenomic survey in patients presenting with pneumonia or acute respiratory infections provides strong evidence that SARS-CoV-2 or any related SARS-like viruses were absent in Wuhan approximately two years prior to the onset of pandemic, although a number of common cold coronaviruses (HKU1, OC43, 229E, and NL63) were commonly detected in our cohort. Indeed, the earliest COVID-19 case, identified by qRT-PCR or next-generation sequencing-based assays performed at designated authoritative laboratories, can currently only be traced back to early December 2019 in Wuhan9. In addition, a retrospective survey of 640 throat swabs from patients with influenza-like illness in Wuhan from the period between 6 October 2019 and 21 January 2020 did not find any evidence of SARS-CoV-2 infection prior to January 202010, such that the ultimate origin of SARS-CoV-2 remains elusive11.
The data presented provide a comprehensive overview of the infectome associated with pneumonia or acute respiratory infections in Wuhan, which is clearly more diverse and complex than described using previous surveys based on targeted PCR or qPCR approach alone12. In light of our observations, we can divide the respiratory infectome into three categories based on epidemiological characteristics: (i) the “core” infectome that is commonly found in patients with respiratory infection and is expected to occur each year globally, (ii) an “emerging” infectome that occurs during outbreaks but are not typically found in the geographic regions under investigation, and (iii) the sporadic infection or zoonotic infectome of new or rare pathogens.
The core infectome comprised of a wide range of common respiratory or systemic pathogens that are subject to frequent screening in hospitals. These include influenza viruses, HMPV, RSV, Moraxella catarrhalis, Acinetobacter spp., Klebsiella pneumoniae, Mycoplasma spp., Haemophilus influenzae, Pseudomonas aeruginosa, Staphylococcus aureus and Streptococcus pneumoniae. However, the remaining pathogens identified here, including rhinoviruses, parainfluenza viruses, coronaviruses, and herpesviruses, have often received far less attention from clinicians and are sometimes ignored entirely, most likely due to the lack of association with severe disease in adults13-15. Nevertheless, our results showed these “neglected” respiratory viruses had high diversity, abundance and prevalence in the cohort of pneumonia or acute respiratory patients studied here in comparison to healthy controls, such that their role as agents of disease should not be underestimated. One scenario is that they represent opportunistic pathogens that take advantage of weakened immunity, such as herpesviruses associated acute respiratory distress (ARDS)16. It is also possible that their pathogenic effects have yet to be identified and may extend to disease manifestations beyond respiratory infections. For example, deep sequencing of a brain biopsy sample suggesting that the OC43 coronavirus may be associated with fatal encephalitis in humans17.
We also identified a potential “emerging” infectome, in this case comprising a single virus - EV-D68 - that may represent a regional or national outbreak of an unconventional respiratory pathogens. The prevalence of EV-D68 remained low from its discovery in 1967 until 201418, when a major outbreak started in the United States and spread to more than 20 countries19, causing severe respiratory illness with potential neurological manifestations such as acute flaccid paralysis in children20, 21. In China, EV-D68 was only sporadically reported22, 23, although serological surveys suggest a much wider prevalence for both children and adults since 200924. We identified six EV-D68 cases, all adults that occurred within a relatively narrow time window between June and December 2016. Phylogenetic analysis revealed that the sequences of these viruses were closely related to each other and to other Chinese strains from the same period (Fig. 4), suggesting that it may be a part of a larger outbreak in China. The EV-D68 cases identified here showed moderate to severe respiratory symptoms, although viral abundance was generally very high, with four of six cases showing >106 RPM (i.e., >10% of total RNA) in the BALF sample. This highlights the active replication and massive proliferation of viruses within the respiratory system of these patients.
Finally, our zoonotic infectome also comprised a single pathogen, Chlamydia psittaci, that is associated with avian species but causes occasional outbreaks in domestic animals (i.e., pigs, cattles, and sheep) and humans25. In humans, C. psittaci infections often starts with influenza-like symptoms, but can develop into serious lung infections and even death26. The single case of C. psittaci identified here was at relatively high abundance level (6396 RPM) and caused a relatively severe disease, with the patient experiencing expiratory dyspnea, severe pneumonia and pleural effusion, and was subsequently transfer to intensive care unit (ICU) for further treatment. Since the patient had no travel history for a month prior to illness, this discovery underlines the risk of local exposure to this bacterium.
Our phylogenetic analysis of the metagenomic data generated revealed extensive intra-specific diversity in each virus species identified highlights the complex epidemiological history of these pathogens. Indeed, the influenza A viruses (H1N1 and H3N2), influenza B virus, HPIV3 and HCoV-OC43 discovered here all comprised multiple lineages (Fig. 4), suggesting these viruses were introduced from diverse sources. Since some of the viruses were closely related to those circulating in other countries it is possible that they represent overseas importations: this is not surprising given that Wuhan is a major domestic travel hub and well linked internationally. Indeed, the rapid and widespread transmission of SARS-CoV-2 between Wuhan and other major cities globally was key to seeding the global pandemic of COVID-19, with early cases in a number of localities all linked to travel from Wuhan26-29.
The methodology used here served as an unbiased investigation of potential emerging pathogens. Meta-transcriptomics has several advantages over traditional diagnostic approaches based on serology or PCR: it targets all types of micro-organisms simultaneously, identifies potential pathogens without a priori knowledge of what micro-organisms might be present, reveals the information (RNA) expressed by the pathogen during infection that is central to agent identification and studies of disease association. With the popularization of next-generation sequencing platform in major hospitals, the approach outlined here can be easily integrated into diagnostic practice with much greater speed and significantly more information output than traditional technologies, providing a broad-scale understanding of infectious disease in general.
Material and Methods
Ethics statement
The sampling and experimental procedures for this study were reviewed by the ethics committees of the Central Hospital of Wuhan and the National Institute for Communicable Disease Control and Prevention, Chinese Center for Disease Control and Prevention. Written informed consents were taken from all patients and volunteers recruited in this study. In addition, for child patients, written consents were also obtained from their parents or guardians. Physicians were informed of results of the pathogen discovery exercise as soon as the meta-transcriptomic results were obtained.
Sample collection from patients and controls
More than 1000 patients were recruited from the Central Hospital of Wuhan between 2016 and 2017. The target clinical conditions were community-acquired pneumonia and acute respiratory infections based on the initial diagnosis made by clinicians. All patients were hospitalized and subject to bronchoalveolar lavage fluid (BALF) collection required by the initial diagnosis for pneumonia or acute respiratory distress syndrome and independent of this study. The BALF sample was divided into two parts for the clinical laboratory test and this study, respectively. Of the BALF samples collected, 409 were subjected to meta-transcriptomic analysis based on their condition and the time between hospitalization and sample collection. No diagnostic information was provided prior to sample selection. To establish a healthy control group, 5ml-10ml of BALF samples were also taken from 32 volunteers without respiratory symptoms. We also included 10 blank controls where only RNase free water was used for nucleic acid extraction and library construction, although only four of these produced viable RNA sequencing results.
Meta-transcriptomic pathogen discovery pipeline
We followed a standard protocol for meta-transcriptomics analysis for each BALF sample. Total RNA was first extracted from 200–300ul of each sample using the RNeasy Plus Universal Kit (Qiagen, USA) according to the manufacturer’s instructions. From the extracted RNA, we performed human rRNA removal and low concentration library construction procedures with the Trio RNA-Seq kit (NuGEN Technologies, USA). The libraries were then subjected to 150bp pair-end sequencing on an Illumina Hiseq 4000 platform at Novagene (Beijing), with target output of 10G base pairs per library. For each of the sequencing results generated, we removed adaptor sequences, non-complex reads, as well as duplicated reads using the BBmap software package. Human and ribosomal RNA (rRNA) reads were subsequently removed by mapping the de-duplicated reads against the human reference genome (GRCh38/hg38) and the comprehensive rRNA sequence collection downloaded from SILVA database30.
The remaining sequencing reads were subject to a pathogen discovery pipeline. For virus identification, the reads were directly compared against reference virus databases using the blastn program and against the non-redundant protein (nr) database using diamond blastx31, with an e-value threshold set at 1E-10 and 1E-5 for blastn and diamond blastx analyses, respectively. Viral abundance was summarized from both analyses, calculated using the relation: total viral reads/total non-redundant reads* 1 million (i.e., reads per million, RPM). To identify highly divergent virus genomes, reads were assembled using megahit32 into contigs before comparison against the nt and nr databases. Those reads with <90% amino acid similarity to known viruses were treated as potential novel virus species. For bacterial and fungi identification, we first used MetaPhlAn233 to identify potential species in both groups. Relevant background bacterial and fungal genomes were subsequently downloaded from NCBI/GenBank and used as a template for read mapping. Based on the mapping results for each case, we generated relevant contigs for blastn analyses against the non-redundant nucleotide (nt) database to determine taxonomy to the species level. The abundance level of bacterial and fungi pathogens was also calculated in the form of RPM based on genome and mitochondrial genome read counts, respectively.
A microbe was considered as “positive” within a specific sample if its abundance level was greater than 1 RPM. To prevent false positives resulting from index hopping, we used a threshold of 0.1% for viruses present in the same sequencing lane: that is, if the libraries contain less than 0.1% of the most abundant library it is treated as “negative”.
Confirmatory testing by conventional methods
For pathogen positive samples, the same RNA used for meta-transcriptomics analysis was also subject to a qRT-PCR assay with primers sets designed for a specific or related group of pathogens (Table S1). RNA was first reverse transcribed by SuperScript™ III First-Strand Synthesis SuperMix for qRT-PCR (Invitrogen, California), and then amplified by TaqPath™ ProAmp™ Master Mix (Applied Biosystems, California). A cycle threshold (CT) value of 38 and above was treated as negative.
Pathogen genomic analyses
For viruses at high abundance levels (i.e., > 1000 RPM), complete genomes were assembled using Megahit and confirmed by mapping reads against the assembled contigs. These genomes were then aligned with related reference virus sequences downloaded from NCBI/GenBank using MAFFT program34. Ambiguously aligned regions were removed using the Trimal program35. Phylogenetic trees were estimated using the maximum likelihood approach implemented in PhyML36, employing the GTR model of nucleotide substitution and SPR branch swapping. The support for each node in the tree was estimated with an approximate likelihood ratio test (aLRT) with the Shimodaira-Hasegawa-like procedures.
Data Availability
All non-human reads have been deposited in the SRA databases under the project accession PRJNA699976. Relevant consensus virus genome sequences have been deposited in NCBI/ GenBank under the accessions MW567157- MW567162, MW570805- MW570808, MW571087- MW571107, MW587035- MW587095.
Data availability
All non-human reads have been deposited in the SRA databases under the project accession PRJNA699976. Relevant consensus virus genome sequences have been deposited in NCBI/ GenBank under the accessions MW567157-MW567162, MW570805-MW570808, MW571087-MW571107, MW587035-MW587095.
Author Contributions
Yong-Zhen Zhang conceived and designed the study. Su Zhao, Yi Hu, Wen Yin, Fang Ni, Hong-Ling Hu, Shuang Geng, and Li Tan performed the clinical work and sample collection. Bin Yu, Jun-Hua Tian and Ying Peng performed epidemiological investigation and sample collection. Mang Shi, Wei-Chen Wu, Yan-Mei Chen, Wen Wang, and Zhi-Gang Song performed the experimental works. Mand Shi, Edwards C Holmes and Yong-Zhen Zhang analysed the data. Mang Shi, Edward C Holmes and Yong-Zhen Zhang wrote the paper with input from all authors.
Competing Interests
All other authors declare no competing interests.
Acknowledgments
This study was supported by the National Natural Science Foundation of China 32041004 (YZZ), 31930001 (YZZ), 81861138003 (YZZ) and 81672057 (YZZ). E.C.H. is supported by an ARC Australian Laureate Fellowship (FL170100022).
References
