Dynamic transcriptional responses during natural infections

To understand how pathogen and host responses change through time in the course of natural infections, we selected 410 samples representative of longitudinal viral episodes (see Methods), among which 102 samples were from children (age <10), 52 were from teenagers (age 14 to 18), and 256 from adults; 157 samples were from male and 253 were from female participants. Samples were classified as PRE samples (n = 98, the last negative samples collected less than 21 days prior to an infection event), DURING samples (n = 166, the first sample of an infection event), DURING2 (n = 30, the second sample of an infection event), POST1 samples (n = 76, the first negative sample collected less than 11 days after the last positive), and POST2 (n = 42, the second negative sample collected 11 to 20 days after the last positive).

Our results show a strong immune signal during the acute phase of infection, followed by a rapid decay of transcriptomic signal after infection clearance: 6,594 genes were differentially expressed (padj <0.05) between the PRE and the DURING1 phase. Among the differentially expressed genes, those that were up-regulated during infection (n = 1,422) had a much stronger significance (Fig 3A and S2 Data) with 210 genes more than 2-fold up-regulated and 23 genes more than 4-fold up-regulated. Gene set enrichment analyses revealed that the most significant biological processes associated with viral infection were those related to adaptive immune response, cytokine production, acute inflammatory response, lymphocyte- and neutrophil-mediated immunity (Gene Set Enrichment Analysis (GSEA) adj, p-value <0.05; all enrichment scores >2.0; see S3 Data). The 10 most significantly overexpressed genes (7.5*10⁻⁵³ <padj<1.8*10⁻³²) were known markers of activation for immune response against viruses, such as CXCL9 and CCL8, CXCL10, CXCL11 and BATF2 (associated with proinflammatory M1 macrophages), RSAD2, GBP1 and IFIT3, IFI6 (interferon-regulated genes), SLAMF7 (marker for plasma cells); see S2 Data for all differentially expressed transcripts.

Open in a separate window

Fig 3

Dynamic transcriptional responses during natural infections.

Panel (a): Volcano plots of differentially expressed genes between DURING1 and PRE, DURING1 and POST1, POST1 and PRE. The cutoff for fold change is >|2|; the default cutoff for p-value is 10^-6. (b) Relative abundance of immune cells, as predicted with Cibersort, across the PRE-DURING1-DURING2-POST1-POST2 samples. Significant pairwise comparisons are reported in S2 Data. (c) Heatmap of normalized expression of selected genes, known to be markers for immune cells, across PRE-DURING1-DURING2-POST1-POST2 samples. (d) Scores F₁₀ of longitudinal samples, calculated using the top 10 genes differentially expressed in the DURING1 vs. PRE comparison, normalized by the respective PRE sample in each episode. Red boxplot lines are the median of distributions, and blue boxes are the interquartile. Whiskers extend to all non-outliers of the distribution. Outliers are represented with red + symbol. For each group, we reported p-values of a Wilcoxon signed rank test for the hypothesis that the dataset comes from a distribution with median = 1. (e) Medians (straight lines) and interquartile (shaded areas) of longitudinal symptomatic (red) and asymptomatic (blue) episodes normalized by the respective PRE sample. Underlying data for panels (b–e) can be found in S5 Data.

Conversely, the difference between the PRE and POST1 groups was minor with no genes differentially expressed after FDR correction (Fig 3A and S2 Data). Furthermore, we did not see evidence of overexpression of recovery genes in POST1 versus DURING1 samples and despite some residual effects in the postinfection samples (i.e., the slower recovery of genes related to CD4 Th1 activation (MAF, ICOS)), the profile of POST samples closely resembles that of PRE samples (Fig 3C). We did not find obvious functional annotation for down-regulated genes during infection, although many transcripts are consistent with previous reports [15].

These longitudinal data were further subjected to in silico cell sorting with CibersortX, in order to compare the relative abundance of different leukocytes. Our results suggest that onset of infection leads to an increase of proinflammatory macrophages M1 (Fold-Change = 4.4), activated CD4 T cells (FC = 1.8), and activated NK cells (FC = 2.9) (all Padj <0.05, Mann–Whitney tests) in nasopharyngeal samples. Soon after an infection is cleared, most changes are restored (S4 Data, Fig 3B, and corresponding S5 Data). We compared these findings with the changes in expression of gene markers for different immune cell types. In agreement with the results from Cibersort, samples during infection had higher levels of expression for genes such as Maf, Icos (CD4 T cell activation), Cd69, FCGR3 (NK cell activation), Cd80, CD86 (macrophages M), as well as markers of cytotoxicity (GZMB, GZMK) and inflammation (IFNG, TNF) (Fig 3C and corresponding S5 Data).

The signature scores F_N, based on the top 10, top 25, and top 100 differentially expressed genes between DURING1 and PRE (see Methods) were able to distinguish infection and postinfection samples from the paired preinfection samples (see Fig 3D for F₁₀ and S1 Fig for F₂₅ and F₁₀₀). To pair samples, we restricted to the (n = 98) longitudinal episodes that included a PRE sample, and we normalized the DURING and POST scores, simply dividing them by the score of their respective PRE. Results showed that in addition to being markedly distinguishable from DURING samples (median DUR1 normalized score = 1.53, p < 10⁻¹³ and median DUR2 normalized score = 1.33, p = 0.004), PRE samples were still weakly distinguishable from POST1 (median POST1 normalized score = 1.09, p = 0.03) but not from POST2 (median POST2 normalized score = 1.07, p = 0.12, significance assessed with Wilcoxon signed rank tests).

The effect of respiratory viral infections was not significantly different across sex and age group, even though some differences in the transcriptomic profile of immune genes across groups existed at baseline (i.e., in no-infection condition). Specifically, in the absence of viral infections, children displayed a larger expression of immune genes compared to adults that was significant for 94 of the 100 genes of signature F₁₀₀. S7 Data and S2 Fig show the infection effect in male versus female and in children versus teenagers versus adults, together with pairwise comparisons between different groups at baseline.

We also compared the transcriptomic signal, quantified by F₁₀ and F₁₀₀, across normalized symptomatic and asymptomatic episodes (an episode was symptomatic if the weekly symptom score associated with any DURING sample was above 9 and asymptomatic otherwise, see Methods). The transcriptomic signal of the acute phase increased more markedly in symptomatic than in asymptomatic episodes (p = 0.01 and p = 0.003, Wilcoxon rank sum test); however, the recovery phase of POST1 and POST2 did not differ significantly across the 2 (Figs ​Figs3E3E and S3). In the next section, we investigate the association between symptom severity and differential expression of immune-related genes in greater depth.

In summary, our longitudinal analyses reveal that viral infection generally leads to a strong activation of immune response, followed by a rapid decay of the expression of immune-related genes once the viral infection has cleared.

Specific transcriptional signatures are associated with symptom severity

To understand how asymptomatic infections differ from symptomatic ones, we compared the transcriptomic profiles of symptomatic and asymptomatic samples. We used a set of 699 samples (see sample selection in Methods), of which 209 were positive for at least 1 respiratory virus. This set of samples, featuring a median of 8 samples per individual, allowed comparison of the signal across symptom level and virus positivity status while controlling for the differences across the 81 contributing individuals. The samples were divided into 4 categories based on diagnosis (positive versus negative for respiratory viruses) and symptoms (symptomatic versus asymptomatic). Overall, there were 63 symptomatic positive (SP) samples, 146 asymptomatic positive (AP), 67 symptomatic negative (SN), and 423 asymptomatic negative (AN).

We found 516 genes differentially expressed (419 overexpressed) in SP compared with AP (Fig 4A and S9 Data). The 81 genes that were at least 2-fold up-regulated at significance level padj <0.05 are well-known modulators of host immune response and included components of interferon pathways (MX1, IFITs, and IFIs groups), cytokines of the CXCL and CCL groups, genes with known antiviral activity (MX1, ISG15, RSAD2, and the OASs and DDXs groups), and genes associated with inflammation, such as IL6. These results were consistent when different thresholds were considered for the symptomatic classification: 506 (342 overexpressed) when using a symptom score threshold of >19 in SP versus AP (S11 Data). In this latter design, there were 43 SP samples, 166 AP, 26 SN, and 464 AN.

Open in a separate window

Fig 4

Specific transcriptional signatures are associated to symptom severity.

(a) Volcano plots of differentially expressed genes between SP samples and AP samples. The cutoff for fold change is >|2|; the default cutoff for P-value is 10e-6. (b) Score factor F₁₀ across the 699 samples calculated with the top 10 genes differentially expressed in the longitudinal analysis and applied to the 4 symptom-infection groups (AN, SN, AP, and SP). Samples expression values are normalized using the DESEQ2 normalization factor. Boxplots show the median (red line) and interquartiles (blue lines) of the distribution. Whiskers of boxplots extend to all non-outliers in the distribution. Outliers are shown as red “+”. (c) Log fold changes +/-SE for top DE genes in SP, AP, and SN when compared to AN samples assumed as baseline. The list of genes was selected as the top 25 DE genes in SP vs. AN. Panel (a) is based on data in S9 Data. Panels (b and c) are based on data in S10 Data. AN, asymptomatic negative; AP, asymptomatic positive; SN, symptomatic negative; SP, symptomatic positive.

The top processes enriched in symptomatic positive versus asymptomatic positive, all with FDR < 0.05, included both innate and adaptive immune response processes (S12 Data), including immune response related to inflammation, such as the production of Type 1 Interferon, IL1, IL12, and IL6.

We then applied the expression signature score (factor score F₁₀), calculated by combining the expression of the top 10 DE genes in the longitudinal analysis (CXCL11, CXCL10, CCL9, IFI6, IFIT3, CCL8, SLAMF7, BATF2, RSAD2, GBP1) to the 4 groups AN, SN, AP, and SP. The score F₁₀ was significantly different among the groups of negative samples, positive asymptomatic, and positive symptomatic samples. In particular, the mean score for SP and AP was statistically different (Wilcoxon rank sum test, p = 0.02) and the mean score for AP samples was statistically different from both AN and SN (Wilcoxon rank sum test, p < 10⁻³² and p < 10⁻¹³, respectively). Results are consistent when using the metrics F₂₅ and F₁₀₀, respectively, calculated with the top 25 and top 100 genes (see S4 Fig).

We also statistically compared (with DESEQ2) the expression profiles of SP, AP, and SN to AN. The log2FC from a selection of the top 25 overexpressed genes in SP (from a total of 1,456 overexpressed genes; see S9 Data and S10 Data) shows a clear hierarchical structure in the 4 groups, with the log2FC of symptomatic positives on average 2 times that of the log2FC for asymptomatic positives (Fig 4C).

Overall, the severity of infection correlates with increased specific transcriptomic levels (see also analysis with symptom severity as a continuous variable in S11 Data). Interestingly, the expression profile of asymptomatic positives is markedly distinguishable from the negative samples (Fig 4B and 4C).

Interactions between bacteria and respiratory viruses

We tested for potential viral and bacterial interactions during the course of natural infection. For this, we detected reads assigned to a number of bacterial species, common to the human respiratory tract, that can cause opportunistic diseases. We then assessed whether such bacterial abundance (measured as normalized fraction of reads) had a differential effect, at the host transcriptomic level, in individuals before and during a viral infection. This was done by performing differential gene expression analyses, using samples representative of longitudinal viral episodes (PRE and DUR1 samples).

We found significant viral/bacterial interactions, at the gene expression level, for Haemophilus influenzae and 2 species of the genus Streptococcus: S. pneumoniae and S. salivarus (padjs <0.05 in DESEQ2 analyses). According to the GSEA results, the abundance of these bacteria in people infected with respiratory viruses was associated with higher expression of genes related to activation of lymphocytes (B cells, NK cells, and granulocytes) and interleukin (IL1, IL6) production (Fig 5A and S14 Data). Although we observed significant interactions at the transcriptomic level, abundance of these bacteria was not significantly different between positive and negative samples (Fig 5B and S14 Data).

Open in a separate window

Fig 5

Interactions between bacteria and respiratory viruses.

Panel (a): Top 10 significantly enriched biological processes, as found with GSEA (FDR < 0.05), derived from the analysis of transcriptomic interactions between abundance of bacteria and viral infection. Panel (b): Comparison of bacteria abundance in PRE vs. DURING samples. Underlying data for Fig 5 can be found in S14 Data. GSEA, Gene Set Enrichment Analysis.

In a similar way, we also assessed the interaction between symptomatology in viral infections and bacterial abundance. At gene expression level, we only found significant interactions between symptomatology and S. pneumoniae. According to GSEA, the abundance of these bacteria in symptomatic people infected with respiratory viruses was associated with higher expression of transcripts related to adaptive immune response (humoral immune response, B cell activate) as well as leukocyte proliferation and NK activation (S5 Fig and S15 Data). Within the participants infected with a respiratory virus, the abundance of bacteria was significantly different in symptomatic versus asymptomatic people only for S. salivarius and S. oralis, where there was a significant decrease of these bacteria in symptomatic samples (Mann–Whitney test: corrected p-value <0.05; S5 Fig).

Specimen and symptoms data collection and classification

Duration of enrollment was heterogeneous for the participants: Some individuals remained enrolled for the entire study period (19 months), whereas others only for a few weeks (the shortest enrollment was 1 week). Demographic and enrollment information for all individuals in the Virome Project are summarized in [22,24] and S22 Data. Two nasopharyngeal samples were collected once a week from all available participants by the study coordinators using minitip flocked swabs, irrespective of the participants’ symptoms. Individuals were not always available for sampling collection, which resulted in incomplete series as described in [22,24]. The 2 samples were stored jointly in 2 mL DNA/RNA Shield (Zymo Research, Irvine, California, United States of America) at 4 to 25°C for up to 30 days and then stored at −80°C in 2 aliquots. Nucleic acids were extracted from 200 μl of sample and 10 μl of internal control using the EasyMAG NucliSENS system (bioMerieux, Durham, North Carolina, USA). Samples were then screened for viruses with the multiplex PCR assay eSensor XT-8 respiratory viral panel (RVP; GenMark Dx, Carlsbad, California, USA) [44]. The multiplex PCR assay tested for 18 commonly circulating respiratory viruses: rhinovirus, influenza A (H1N1, H1N1pdm2009, H3N2, and any subtype), influenza B, RSV (A and B) parainfluenza (1, 2, 3, 4), coronavirus 229E, NL63, OC43, HKU1, adenovirus (B/E and C), and metapneumovirus. Per manufacturer specification, samples PCR-positive for a particular virus were identified by an electrical signal intensity of ≥2 nA/mm², except for coronavirus OC43, for which positive results were identified by an intensity of ≥25 nA/mm². Demographic and enrollment information for all individuals in the Virome Project and features of the 4,215 samples collected are summarized in [22,24].

Collection of RNA-seq data

Among the 4,215 samples collected, we selected a subset of samples for RNA sequencing. Selection criteria were: samples testing positive by the RVP for one of 18 respiratory pathogens along with corresponding preinfection and postinfection samples available for the same participant where at least 40% of the samples had RNA integrity score (RIN) above 4 as assessed with Agilent Bioanalyzer (Santa Clara, California, USA). For the selected samples, the remaining quantity of eluted RNA was sequenced with Illumina at the Columbia Genome Center with the Ribo-Zero rRNA Removal Kit, target 30 M single-end 100 bp reads for the first batch (Protocol single in model design) and 20 M paired-end 100 bp reads for the second batch (Protocol paired in model design).

Specifically, for single-end sequenced samples libraries were sequenced with Illumina HiSeq2500/HiSeq4000 at Columbia Genome Center. Samples were multiplexed in each lane, which yields targeted number of single-end/paired-end 100 bp reads for each sample, as a fraction of 280 to 400 million reads for the whole lane. RTA (Illumina) was used for base calling and bcl2fastq2 (version 2.17) for converting BCL to fastq format, coupled with adaptor trimming. The reads were mapped to human reference genome NCBI/build37.2 using STAR(2.5.2b) [45] and featureCounts (v1.5.0-p3) [46]. For paired-end sequenced samples libraries were sequenced with Illumina NovaSeq 6000 at Columbia Genome Center. Samples were multiplexed in each lane, which yields targeted number of paired-end 100 bp reads for each sample. RTA (Illumina) was used for base calling and bcl2fastq2 (version 2.20) for converting BCL to fastq format, coupled with adaptor trimming. A pseudoalignment was performed to a kallisto index created from transcriptomes (Human: GRCh38) using kallisto (0.44.0) [47].

After eliminating failed and duplicate samples, there were 847 unique samples sequenced from 104 different individuals, with a median number of 7 samples per participant (see S1 Data for details on number of samples per individual and date of collection).

Quality control and pathogen detection

Quality control of the raw reads was assessed with FASTQC (see sequence quality plots in S9 Fig; underlying data for this figure available at GEO accession number GSE223679) [48]. In order to differentiate between human and nonhuman reads, and quantify reads related to respiratory viruses, we ran the pipeline Pandora: an open-source pipeline developed by the Rabadan lab that has been validated and used previously [49]. The Pandora workflow takes RNA-seq data as input and produces the annotated microbial reads present in the sample as outputs. The pipeline runs through 5 main modules: (1) mapping to the human host genome using STAR [45] and bowtie2 [50]; (2) processing of the host-mapped reads with featureCounts [46] to obtain the host gene expression profile; (3) de novo assembly of host-subtracted short (2 × 100 bp) reads using Trinity [51] to create long contiguous full-length transcripts (contigs); (4) contigs alignment to curated sets of NCBI sequences from viruses, bacteria, fungi, and other eukaryotic parasites with BLAST [52]; and (5) filtering and parsing the BLAST results into a report specifying the detected microbial abundances in the sample.

To assess positivity for respiratory viruses in the 847 sequenced samples we first isolated pathogen reads attributed to human viruses falling in the picornaviruses, adenoviruses, polyomaviruses, parvoviruses, coronaviruses, paramoxyviruses, and pneumoviruses families. We then set a reads threshold for assessing positivity to correct for (1) short reads fragments typically present in metagenomics samples; and (2) index hopping, a phenomenon occurring in multiplex sequencing at different scales for single and paired-end protocols (typically 2 orders of magnitude different [53]). The threshold was set different for samples testing positive and negative by qPCR (qPCR+ and qPCR-). Samples qPCR+ for a specific family of viruses were considered positive if containing at least 15 reads of a virus in that family and negative otherwise. Samples qPCR- were considered negative for a respiratory virus if containing less than 500 reads of that respiratory virus. Conversely, qPCR- samples with more than 500 reads of a virus were considered positive for that virus, with an additional correction for index hopping that became relevant from samples with over 10^5 reads of a specific virus. More specifically, the threshold of positivity for a virus was max(500, IH*10^-3) for single-end protocols and max(500, IH*10^-5) for paired-end protocols, with IH being the largest reads count of that virus within samples the same batch. Samples qPCR- containing more than 500 reads but less than max(500, IH*10^-3) for single-end or max(500, IH*10^-5) for paired-end were considered Undetermined and excluded from subsequent analysis.

In S8 Fig, we evaluate the agreement of qPCR and RNAseq for different RNAseq thresholds using only viruses included in the RVP multiplex.

We defined an infection event as a group of consecutive weekly specimens from a given individual that were positive for the same virus. Samples from the same individual that were positive for the same virus were considered part of the same infection event if separated by a week gap (i.e., a missing test). Coinfections (samples positive for multiple respiratory viruses) were considered one infection event if all samples of the event were positive for the same viruses, or if later samples were positive for a subset of viruses (e.g., coronavirus+rhinovirus on week 1 and rhinovirus only on week 2). Conversely, if a new virus was detected in addition to a previous one (e.g., rhinovirus on week 1 and coronavirus+rhinovirus on week 2), we considered the infections as separate events.

Nonhuman reads were further used to identify and quantify bacterial species commonly found in the respiratory tract (H. influenzae, H. Parainfluenziae, S. pneumoniae, S. Pseudopneumoniae, S. Oralis, S. Salivarius). These reads were analyzed using Kracken [54], a taxonomic sequence classifier that analyzes short reads, in combination with Bracken, a companion program to Kracken that accurately classifies and quantifies reads assigned at species resolution using a Bayesian approach [55].

Fig 1 describes the structure of the complete Virome dataset leveraged in this analysis and available at the web-resource The Virome of Manhattan Project: Virome Data Explorer.