IMPACC study design and participant recruitment
IMPACC is a prospective, longitudinal COVID-19 study that enrolled more than 1,000 hospitalized patients, as previously described27,32,55,56,57,58. Adults aged 18 years or older were recruited from 20 hospitals at 15 academic institutions across the United States (Fig. 1a). All participants had SARS-CoV-2 infection confirmed by reverse transcription PCR (RT–PCR), and none had received a SARS-CoV-2 vaccine at enrollment.
Nasal swabs, blood samples and, when applicable, endotracheal aspirates from ventilated patients were collected within 72 h of hospital admission (visit 1), followed by collections on days 4, 7, 14, 21 and 28. Convalescent samples were collected at 3, 6, 9 and 12 months. Additional nasal, blood and endotracheal aspirate samples were collected, when possible, within 24 h and 96 h of escalation to intensive care or hospital readmission more than 48 h after discharge. These were classified as “escalation visits” in the IMPACC dataset and were used only for longitudinal generalized additive mixed-model (GAMM) and CyTOF cellular association analyses to limit the effects of more frequent sampling in critically ill participants.
Participants were assigned to one of five acute COVID-19 trajectory groups using latent class mixed modelling of a seven-point respiratory illness severity scale24. They were also assigned to one of four post-acute sequelae of COVID-19 (PASC) patient-reported outcome (PRO) groups using latent class mixed modelling of surveys completed at 3, 6, 9 and 12 months32. The models incorporated EQ-5D-5L responses59, a health recovery score ranging from 1 to 100, and the following PROMIS measures: Physical Function, Cognitive Function60, Global Health Mental 2a61, Psychosocial Illness Impact-Positive—Short Form 8a61 and Dyspnea Time Extension62. PROMIS measures were scored and standardized according to official PROMIS instructions. Cohort demographics and clinical characteristics are provided in Supplementary Table 1.
COVID-19 study ethics and informed consent
The Department of Health and Human Services Office for Human Research Protections and the National Institute of Allergy and Infectious Diseases determined that IMPACC qualified for a public health surveillance exemption. The study protocol was reviewed by the institutional review board (IRB) at each participating site. Twelve sites conducted IMPACC as a public health surveillance study, while three sites incorporated it into existing IRB-approved protocols: The University of Texas at Austin (IRB 2020-04-0117), University of California San Francisco (IRB 20-30497) and Case Western Reserve University (IRB STUDY20200573). Participants at these sites provided informed consent.
Individuals enrolled through public health surveillance sites received information sheets describing the study, sample collection procedures, planned analyses and data de-identification procedures. Participants who declined after reviewing the study information were not enrolled. Hospitalized participants were not compensated during their admission but received compensation for subsequent outpatient visits and surveys. The trial was registered at ClinicalTrials.gov (NCT04378777) and conducted according to the Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) guidelines (Extended Data Fig. 1a).
Biological sample processing and multi-omics assays
Samples were processed according to previously published IMPACC methods27,32,55,56,57,58 and the IMPACC study protocol55. At each visit, 10 ml of blood and a nasal swab were collected. Blood was processed within 6 h. A 2.5 ml Greiner Vacuette CAT serum-separating tube (SST; 454243P) was used for serum, while a 7.5 ml Sarstedt EDTA monovette (NC9453456) was used for whole blood, PBMCs and plasma. SST tubes remained upright at room temperature for at least 30 min before centrifugation at 1,000g for 10 min. Serum was aliquoted in 100 μl volumes for downstream testing.
After gentle mixing, 270 μl of EDTA-treated whole blood was aliquoted for CyTOF analysis. The aliquot was added directly to a Maxpar Direct Immune Profiling Assay tube containing the antibodies listed in Supplementary Table 9 and incubated for 30 min at room temperature. Smart Tube Prot1 Stabilizer (410 μl) was then added, followed by a 10 min room-temperature incubation. Samples were stored at –80 °C until shipment to the processing core. The remaining blood was centrifuged at 1,000g for 10 min at room temperature. Plasma was aliquoted in 500 μl volumes and stored at –80 °C for proteomic and metabolomic analysis. PBMCs were isolated from the remaining sample using the SepMate and Lymphoprep system (StemCell) according to the manufacturer’s instructions and previously described procedures55. PBMCs were stored at 2.5 × 105 cells in 200 μl of Qiagen RLT buffer containing β-mercaptoethanol at –80 °C.
Interior nasal turbinate swabs were stored in 1 ml of Zymo DNA/RNA Shield reagent. RNA was extracted in duplicate from 250 μl aliquots and purified with the KingFisher Flex system (Thermo Fisher) and Quick-DNA/RNA MagBead kit (Zymo Research). Duplicate extracts were pooled and aliquoted in 20 μl volumes for SARS-CoV-2 RT–qPCR and RNA sequencing.
For ventilated participants, endotracheal aspirates were collected in a 40 cm3 Argyle specimen trap and processed within 2 h. A 500 μl aliquot of aspirate diluted 1:1 with calcium- and magnesium-free Maxpar PBS was combined with 500 μl of DNA/RNA Shield in a Zymo tube containing lysis beads. Samples were stored at –80 °C for bulk RNA sequencing.
Samples were shipped to designated processing cores for nasal, PBMC and endotracheal aspirate RNA sequencing; plasma proteomics; serum cytokine proximity extension assay; EBV and CMV antibody measurements; whole-blood CyTOF; and plasma metabolomics27,32,55,56,57,58. Detailed protocols are available in the cited publications. Analytes measured by PBMC RNA-seq, nasal RNA-seq, serum cytokine and chemokine testing, and plasma metabolomics are listed in Supplementary Table 8.
Nasal, endotracheal aspirate and PBMC RNA sequencing
Nasal transcriptomic data were processed using a Galaxy-based workflow. Extracted RNA was DNase-treated and depleted of human ribosomal RNA before random-hexamer amplification with the Ribo-Zero Plus kit. Convalescent samples collected at 3, 6, 9 and 12 months instead underwent poly-A amplification using SMART-seq V4. Libraries were quantified with the Quant-iT dsDNA High Sensitivity Assay and Fragment Analyzer. Libraries containing adapter dimers representing more than 4% of the electropherogram area were excluded. K562 technical controls were compared between batches to monitor batch variation.
Libraries were normalized to 10 nM and sequenced on a NovaSeq6000 using 100 bp paired-end reads. Demultiplexed, unaligned BAM files were generated with Picard tools (https://broadinstitute.github.io/picard/) and converted to FASTQ format with Samtools bam2fq (v1.4)63. Adapter removal and quality trimming were performed with Trimmomatic (v0.36.5)64. Reads were trimmed at the 3′ end and filtered to a minimum base quality of Q30. Trimmed reads were aligned to the GRCh38 human reference genome65 with STAR (v2.4.2a)66, using Ensembl release 91 gene annotations67. Gene-level counts were generated with HTSeq-count (v0.4.1)68. Picard, FASTQC (v0.11.3) and Samtools (v1.2) were used for quality control. Samples were excluded when the median coefficient of variation in gene coverage exceeded 0.8 or aligned counts were below 1 million.
Endotracheal aspirate RNA was DNase-treated and depleted of human ribosomal RNA. cDNA was synthesized with random hexamers to capture coding and noncoding transcripts. Libraries were sequenced on a NovaSeq6000 using NovaSeq S4 flow cells and 100 bp paired-end reads, with a target of 50 million reads per sample. Human reads were aligned to GRCh38 and quality-controlled. Raw counts were normalized with the trimmed mean of M-values (TMM) method in edgeR69. Participant and control samples were co-sequenced within batches to reduce batch effects.
For PBMC transcriptomics, RNA was extracted from 2.5 × 105 stored cells using the Quick-RNA MagBead Kit (Zymo) with DNase treatment. RNA quality was assessed using a Qubit HS RNA assay and Fragment Analyzer. cDNA was generated from 10 ng RNA with the SMART-Seq v4 Ultra Low Input RNA Kit (Takara Bio). Nextera XT was used for library preparation after bead-based cleanup. Equimolar libraries were sequenced on an Illumina NovaSeq6000 using 100 bp paired-end reads and a target of 25 million reads per sample. Cutadapt (v1.14) was used for adapter trimming and quality filtering.
PBMC reads were aligned with STAR to a composite reference containing GRCh38, Ensembl release 9165, and SARS-CoV-2 strain MN908947.370. Gene counts were generated with HTSeq-count within the STAR workflow. Quality control used FASTQC (v0.11.5), Picard tools (v2.22) and STAR alignment logs. Metrics included average read quality above Q30, the percentage and number of uniquely mapped reads, and other alignment statistics.
Human-infecting virus taxonomic alignments from PBMC, nasal and endotracheal aspirate RNA-seq were obtained with CZID71. CZID removes host reads and aligns remaining sequences against NCBI nucleotide and non-redundant databases. A virus was considered detected when at least one read mapped to both databases. Water controls were included on each transcriptomics plate or batch; none contained detectable reads for the viruses evaluated, supporting the positivity threshold.
One nasal transcriptomics plate with unusually high HSV1 positivity, exceeding 60% compared with an approximately 10% batch average, was excluded from HSV1 detection and host transcriptomic analyses because of possible cross-contamination. When samples were sequenced in multiple batches, the mean reads per million (RPM) for each virus across the participant event was used. Additional RNA-seq data were obtained from the Mount Sinai COVID-19 Biobank (syn35874390)25,26 and processed through CZID to identify human-infecting viral reads.
Nasal SARS-CoV-2 RT–qPCR and viral load measurement
SARS-CoV-2 viral load was measured in nasal swabs by RT–qPCR targeting the N1 and N2 regions of the nucleocapsid gene, following the CDC protocol (https://www.cdc.gov/flu/php/laboratories/influenza-sars-cov-2-multiplex-assay.html). Reactions used Quantabio One-Step RT–qPCR ToughMix and were run on a QuantStudio 5 instrument. N1 and N2 cycle threshold (Ct) values were the primary viral-load measurements.
Whole-blood CyTOF immune-cell profiling
CyTOF-prepared whole-blood samples were thawed using the SmartTube Prot1 thaw and erythrocyte-lysis protocol. Samples were barcoded with the Fluidigm Cell-ID 20-Plex Palladium Barcoding Kit, pooled and stained with surface antibodies listed in Supplementary Table 9. Antibodies were added at 1 μl per 10 million cells in 100 μl per 10 million cells and incubated for 30 min on ice.
Samples were washed twice in CyFACS buffer at 800g for 3 min, fixed and permeabilized with BD Biosciences solution, and stained intracellularly for GZMB. Following paraformaldehyde fixation and iridium/osmium labeling, samples were acquired on a Fluidigm Helios mass cytometer. Data were normalized and concatenated with Fluidigm software. An internal Mount Sinai pipeline removed acquisition outliers, EQ beads and events with low DNA signal. Palladium barcoding and cosine similarity were used for demultiplexing, while low signal-to-noise cells and acquisition multiplets were removed.
Each sample was clustered into 1,000 k-means groups. A subset of clusters was manually annotated with Clustergrammer2 (https://github.com/ismms-himc/clustergrammer2) to create a reference matrix. Cell-type assignments were based on the highest or consensus similarity to reference types and then mapped to individual cells for downstream quantification. Antibodies are listed in the supplemental reporting summary.
Plasma proteomic analysis
Plasma proteins were depleted with perchloric acid to remove highly abundant proteins and improve detection of lower-abundance proteins72,73,74. Prepared samples were loaded onto Evotips and analyzed on an EVOSEP One system using the 60-samples-per-day method and a 21 min gradient. The system was coupled to a timsTOF Pro mass spectrometer operating in data-dependent acquisition parallel accumulation–serial fragmentation (DDA-PASEF) mode. HSV1 analysis included Swiss-Prot-reviewed proteins assigned to human herpesvirus 1 strain 17 (UniProt Proteome ID: UP000009294).
Plasma metabolomics profiling
Metabolomic profiling was performed by Metabolon using in-house standards75,76. Randomized samples underwent methanol precipitation and were separated into fractions for UPLC–MS/MS analysis in positive and negative ion modes using reverse-phase and hydrophilic-interaction chromatography. Metabolites were analyzed on a Thermo Q-Exactive mass spectrometer and identified using retention index, accurate mass within ±10 ppm and MS/MS spectral matching to the Metabolon reference library. Procedures followed Metabolomics Standards Initiative recommendations76.
Serum cytokine and chemokine measurements
Serum inflammation biomarkers were measured with the Olink Inflammation panel, which quantifies 92 inflammation-related proteins using proximity extension assay (PEA). Oligonucleotide-labeled antibody pairs bind to target proteins and generate PCR targets when brought into proximity. Following overnight incubation at 4 °C, targets were amplified and quantified using a Biomark microfluidic qPCR system (Fluidigm) with integrated quality-control measures.
EBV and CMV antibody titre testing
Serum antibodies against Epstein–Barr virus (EBV gp350) and cytomegalovirus (CMV gB) were measured using a Luminex assay as previously described58. Antigens were coupled to barcoded beads, and serum was diluted 1:400 in PBS containing 0.5% Triton X-100. Samples were incubated with antigen-coupled beads and Assay Chex Control beads, washed and detected with fluorescently labeled anti-human IgG or IgA antibodies. Plates were analyzed on a Luminex Flex 3D instrument using a minimum of 50 beads per antigen. The assay was performed in approximately half of the cohort (n = 479). EBV titres were normalized by regression against the four Assay Chex control beads, batch and plate.
CMV serostatus at visit 1 was also measured for 497 participants using a CMV IgG ELISA (Aviva, GWB-BQK12C). Samples were tested in duplicate. Average antibody index values above 1.1 were considered positive, values below 0.9 negative, and values from 0.9 to 1.1 equivocal according to the manufacturer’s instructions. Equivocal results were classified as seropositive. Results from both CMV assays were combined to compare CMV serostatus at hospital admission between patients with and without CMV transcripts (Fig. 3b).
Statistical analysis and study reproducibility
All IMPACC sites used a standardized biological-sample collection protocol55. Samples were randomized across processing batches, with longitudinal samples from the same participant assigned to the same batch whenever possible. Randomization was stratified by disease severity and age, and race, ethnicity, sex and enrollment site were balanced across batches. Because IMPACC was an observational study, biological replicates were prioritized over technical replicates except when technical replication was required to assess batch effects. Unless otherwise specified, analyses used a single measurement per participant and time point.
All P values were adjusted using the Benjamini–Hochberg procedure and are reported as adjusted P values unless adjustment was unnecessary. For multi-omics analyses, corrections were applied across comparisons or models conducted for each virus. Specific statistical methods are described below.
Clinical characteristics, demographics and outcomes
Associations between chronic viral transcript detection, trajectory group and age were tested with cumulative link mixed modelling using the ordinal R package (v2019.12-10)77. This approach accounted for ordinal COVID-19 severity and age categories while including enrollment site as a random effect. Because trajectory groups 4 and 5 contributed most endotracheal aspirate samples, viral comparisons in this tissue used Fisher’s exact test restricted to these two groups.
$${\rm{Trajectory}}\_{\rm{group}} \sim {\rm{virus}}\_{\rm{status}},{\rm{random}}={\rm{enrollment}}\_{\rm{site}}$$
$$\begin{array}{l}{\rm{Admit}}\_{\rm{age}}\_{\rm{quintile}} \sim {\rm{virus}}\_{\rm{status}}+{\rm{trajectory}}\_{\rm{group}},{\rm{random}}\\ \,=\,{\rm{enrollment}}\_{\rm{site}}\end{array}$$
Associations between virus status and ethnicity, sex, steroid use or remdesivir treatment were evaluated using chi-squared tests of independence with the rstatix R package (v0.7.2)78.
Within trajectory group 4, long-term mortality was analyzed using a Cox proportional hazards model incorporating steroid treatment and remdesivir treatment as fixed effects and enrollment site as a random effect. Clinical outcomes were recorded by clinical staff from electronic medical records and verified by the IMPACC Clinical and Data Coordinating Center. Participants without confirmed death were right-censored at the date of last follow-up or study completion.
$$\begin{array}{l}{\rm{Surv}}({\rm{time}}\_{\rm{to}}\_{\rm{event}},{\rm{right}}\_{\rm{censoring}}) \sim {\rm{ever}}\_{\rm{steroids}}\\ \,+\,{\rm{ever}}\_{\rm{remdesivir}}+{\rm{virus}}\_{\rm{status}}+(1|{\rm{enrollment}}\_{\rm{site}})\end{array}$$
Other clinical features, including complications, comorbidities and medication use, were analyzed with logistic mixed-effects models that included sex, age quintile and trajectory group as fixed effects and enrollment site as a random effect. Virus-positive participants were compared only with participants who had the relevant tissue sequenced at least once and remained negative at all assessed time points.
$$\begin{array}{l}{\rm{Clinical}}\_{\rm{feature}} \sim {\rm{virus}}\_{\rm{status}}+{\rm{sex}}+{\rm{age}}\_{\rm{quintile}}\\ \,+\,{\rm{trajectory}}\_{\rm{group}}+(1|{\rm{enrollment}}\_{\rm{site}})\end{array}$$
Additional logistic mixed-effects modelling evaluated whether Anelloviridae detection was independently associated with immunosuppressive medication use, COVID-19 severity, sex and age.
$$\begin{array}{l}{\rm{Anelloviridae}}\_{\rm{status}} \sim {\rm{ever}}\_{\rm{immsupp}}+{\rm{trajectory}}\_{\rm{group}}\\ \,+\,{\rm{sex}}+{\rm{discretized}}\_{\rm{admit}}\_{\rm{age}}\_{\rm{quantile}}+(1{\rm{|enrollment}}\_{\rm{site}})\end{array}$$
Associations between viral detection and PASC PRO groups were tested using a model that included virus status, immunosuppression, trajectory group, sex, age quintile and enrollment site.
$$\begin{array}{l}{\rm{PRO}}\_{\rm{group}} \sim {\rm{virus}}\_{\rm{status}}+{\rm{immunosuppressed}}\_{\rm{status}}\\ \,+\,{\rm{trajectory}}\_{\rm{group}}+{\rm{sex}}+{\rm{a}}{\rm{ge}}\_{\rm{quintile}}+(1{\rm{|enrollment\_site}})\end{array}$$
CyTOF cellular immunophenotyping analysis
Associations between chronic virus detection and circulating immune-cell phenotypes were tested with linear mixed-effects models. Models included trajectory group, sex, age quintile and visit number as fixed effects, with enrollment site and participant as nested random effects. Cell-type abundance was calculated relative to the non-granulocyte population, followed by log1p transformation and scaling.
$$\begin{array}{l}{\rm{Normalized}}\_{\rm{celltype}}\_{\rm{abundance}} \sim {\rm{virus}}\_{\rm{status}}+{\rm{trajectory}}\_{\rm{group}}\\ \,+\,{\rm{sex}}+{\rm{discretized}}\_{\rm{admit}}\_{\rm{age}}\_{\rm{quantile}}+{\rm{event}}\_{\rm{type}},{\rm{random}}\\ \,= \sim 1|{\rm{enrollment}}\_{\rm{site}}/{\rm{participant}}\_{\rm{id}}\end{array}$$
Cytokine, chemokine and plasma metabolomics analysis
Generalized additive mixed models from the gamm4 R package (v0.2-6)82 were used to compare Olink cytokine measurements and plasma metabolites between participants with chronic viral transcripts and those without detected human-infecting viral transcripts other than SARS-CoV-2. Analytes were modeled against days from admission using cubic regression splines. Models included interactions between chronic virus status and trajectory group, as well as sex, age quintile, trajectory group and nasal SARS-CoV-2 viral RPM. Only visits with both PBMC and nasal transcriptomic data were included. Because virus status was included as both a main effect and interaction term, an adjusted P value threshold of 0.01 was used.
$$\begin{array}{l}\mathrm{Analyte} \sim {\rm{s}}(\mathrm{days},\mathrm{bs}= \mbox{`} \mathrm{cr}\mbox{‘})+{\rm{s}}(\mathrm{days},\mathrm{bs}= \mbox{`} \mathrm{cr}\mbox{‘},\mathrm{by}= \mbox{`} \mathrm{virus}\_\mathrm{status}\mbox{‘})\\ \,+\,{\rm{s}}(\mathrm{days},\mathrm{bs}= \mbox{`} \mathrm{cr}\mbox{‘},\mathrm{by}= \mbox{`} \mathrm{trajectory}\_\mathrm{group}\mbox{‘})+\mathrm{virus}\_\mathrm{status}\\ \,+\,\mathrm{sex}+\mathrm{trajectory}\_\mathrm{group}+\mathrm{age}\_\mathrm{quintile}\\ \,+\,\mathrm{sarscovs}2\_\mathrm{nasal}\_\log \_\mathrm{rpm},\mathrm{random}\\ \,= \sim (1|\mathrm{enrollment}\_\mathrm{site}/\mathrm{participant}\_\mathrm{id})\end{array}$$
Differential gene expression in nasal and PBMC RNA-seq
The limma R package (v3.46.0)83 was used to identify differentially expressed genes associated with chronic viral transcript detection in nasal and PBMC RNA-seq data. Models included age quintile, sex, visit number, trajectory group, nasal SARS-CoV-2 viral RPM and virus status.
$$\begin{array}{l} \sim {\rm{age}}\_{\rm{quintile}}+{\rm{sex}}+{\rm{visit}}\_{\rm{number}}+{\rm{trajectory}}\_{\rm{group}}\\ \,+\,{\rm{sarscov}}2\_{\rm{nasal}}\_\log \_{\rm{rpm}}+{\rm{virus}}\_{\rm{status}}\end{array}$$
Upregulated and downregulated genes were analyzed separately using hypergeometric pathway enrichment against the Reactome database (v95)84 with the clusterProfiler R package (v3.18.1)85.
Research reporting summary
Additional information about the study design, methods and reporting standards is available in the Nature Portfolio Reporting Summary linked to this article.
Source: www.nature.com