Viruses use diverse strategies to enhance their persistence and dissemination, including establishing chronic infection^1,8,9. This strategy is exemplified by human-infecting viruses, particularly members of Herpesviridae and Anelloviridae families, which establish lifelong infections in a significant portion of the human population^10,11,12. Although primary infection typically remains asymptomatic in immunocompetent individuals, some chronic viruses contribute to the development of autoimmune disorders and cancers, among other adverse health outcomes^1,13,14,15. These viruses typically remain dormant, but can reactivate during periods of stress, sleep deprivation, surgery, hormonal imbalances or in critical illness^16,17,18,19. The full range of immunological consequences from these viral reactivations remain largely unknown.

Since its emergence, SARS-CoV-2 has resulted in more than 774 million cases of COVID-19 and 7 million deaths^20,21. Owing to the physiological stress introduced by SARS-CoV-2, underlying viral infections may reactivate and potentially contribute to the immunological consequences of COVID-19. For example, reactivation of Herpesviridae, including Epstein–Barr virus (EBV), cytomegalovirus (CMV), human herpesvirus 6 (HHV6) and human herpesvirus 8 (HHV8), is associated with worse clinical outcomes in patients with COVID-19 (refs. ^2,3,4,5,6). Reactivation of CMV and EBV in particular have been linked to more severe outcomes, including increased mortality^2,7. Additionally, patients with long COVID develop increased EBV antibody titres, raising the possibility that reactivation of these viruses may contribute to long COVID^4,22.

Many of the foundational COVID-19 viral reactivation studies have been limited by small sample sizes, focused on a subset of Herpesviridae or relied solely on antibody responses to assess viral reactivation^{2,3,4,5,6,7,22,23}, as opposed to measuring transcripts of actively replicating viruses. Thus important gaps remain in our understanding of the dynamics and biology of viral reactivation during acute COVID-19 and their role in long COVID.

To address the knowledge gap in viral reactivation in COVID-19, we leveraged IMPACC, a longitudinal prospective observational study of 1,154 patients who were hospitalized for COVID-19, which evaluated patients during acute hospitalization and for 12 months post-hospitalization. We carried out longitudinal, multi-omic analyses of nasal swabs, peripheral blood mononuclear cells (PBMCs) and endotracheal aspirates and found significant reactivation of chronic viruses, particularly from the Herpesviridae and Anelloviridae families, associated with COVID-19 severity. By integrating host and viral transcriptomics, cytokine profiling, cellular immunophenotyping, metabolomics and proteomics, we observed distinct viral reactivation dynamics and striking associations between viral reactivation, clinical outcomes, immunologic features and patient demographics, both during acute COVID-19 and during long COVID. Our results provide insights into the endogenous virological landscape of patients with COVID-19, highlighting the complex interplay between SARS-CoV-2 infection, chronic viral reactivation, host immune responses and clinical outcomes.

The Immunophenotyping Assessment in a COVID-19 Cohort (IMPACC) consortium enrolled 1,154 patients who were hospitalized for COVID-19 across 20 US hospitals between May 2020 and March 2021 (Fig. 1a, Extended Data Fig. 1a and Supplementary Table 1). All participants were COVID-19 vaccine-naive at the time of enrolment. To assess COVID-19 severity, participants were assigned to one of five trajectory groups (TG1–TG5) using latent class mixed modelling of respiratory status over the first 28 days²⁴. Groups were classified as mild (TG1, length of stay (LOS) approximately 3–5 days, n = 228), moderate (TG2, LOS approximately 7–14 days, n = 263), severe (TG3, LOS approximately 10–14 days and discharged with limitations, n = 333), critical (TG4, critically ill with LOS greater than 28 days, n = 222) or fatal within 28 days (TG5, n = 108). From each participant, bulk RNA sequencing (RNA-seq) was performed on PBMCs, nasal swabs and for mechanically ventilated patients, endotracheal aspirates, at up to ten visits during one year post-hospital admission (Fig. 1b). Additionally, we assessed whole-blood immune cell populations by cytometry by time of flight (CyTOF), serum EBV and CMV antibody titres, serum cytokine levels by proximity extension assay (PEA), and the plasma proteome and metabolome by mass spectrometry.

Upon screening RNA-seq data for any human-infecting virus transcripts, we identified viral RNA in nasal, endotracheal aspirate and PBMC samples for a diverse range of human-infecting viruses beyond SARS-CoV-2, including EBV, CMV, HHV6, human alphaherpesvirus 1 and 2 (HSV1, HSV2) and several Anelloviridae and enterovirus species (Fig. 1c,d and Extended Data Fig. 1b–d). Unsurprisingly, SARS-CoV-2 was the most prevalent virus, and was primarily found in nasal and endotracheal aspirate samples (Fig. 1c). We confirmed that SARS-CoV-2 abundance measured by RNA-seq highly correlated with results from reverse transcription with quantitative PCR (RT–qPCR) (Extended Data Fig. 2a,b).

Beyond SARS-CoV-2, most of the other detected viruses were chronically infecting (for example, Herpesviridae and Anelloviridae); we focused on these owing to their higher prevalence in our cohort (Extended Data Fig. 1b,c). Among Herpesviridae, transcripts of HSV1, EBV and CMV were commonly detected across compartments during acute COVID-19 (first 40 days after admission) and were detected less frequently during the convalescent period (2 months or more post-admission) (Fig. 1c). Additionally, we detected a wide range of Anelloviridae and acute-infecting enterovirus species (Extended Data Fig. 2c,d), and analysed their collective viral load at the family and genus level, respectively, owing to the sparsity of individual strains. Reactivated viruses were common across recruitment sites (Extended Data Fig. 2e) and viral detection rates were comparable across sequencing cores (Extended Data Fig. 2f), suggesting minimal-to-no enrolment or sequencing site effects. Furthermore, read duplication within batches was low (Extended Data Fig. 2g), with detected duplication likely reflecting specific viral gene expression (Extended Data Fig. 2h). Finally, viral alignment E values were exceedingly low (Extended Data Fig. 2i), providing confidence in the taxonomic alignment.

Notably, each viral species displayed unique temporal dynamics of reactivation relative to hospital admission (Fig. 1d). For example, EBV reactivated early in the disease course, with 24% of participants having detectable transcripts near time of admission (days 1–8, 260 out of 1,080 participants), followed by a gradual decline over time. Unlike EBV, Anelloviridae transcript frequency remained constant until day 20 post-admission, followed by a slow decline. By contrast, HSV1 and CMV reactivated later in disease course, with HSV1 detected in 43% of participants with endotracheal aspirate samples (13 out of 30) and 15% of nasal samples (14 out of 90) and CMV detected in 8% of PBMC samples (9 out of 110) at 19–23 days post-admission (Fig. 1d and Extended Data Fig. 3b–e). These temporal dynamics were also reflected in viral reads per million (RPM) across samples (Extended Data Fig. 3a). Of note, as sample composition changed over time owing to participant death, drop-out and discharge, we also assessed viral detection rates by time from symptom onset, which showed similar patterns to the days-from-hospitalization analysis (Extended Data Fig. 3b–i).

Viral detection varied across compartments, with EBV transcripts found to be more common in PBMCs and HSV1 in nasal and endotracheal aspirate samples. Nevertheless, viral transcripts often correlated between compartments (Fig. 1e). Furthermore, among participants with viral reactivation (47.9%, 550 out of 1,148), most had only one virus detected in the acute period (67.6%, 372 out of 550), and co-detection of multiple viruses simultaneously was infrequent (Extended Data Fig. 4a–d).

Finally, we validated the detection of viral transcripts in an external cohort of COVID-19 whole-blood RNA-seq^25,26, observing strikingly similar temporal dynamics for EBV, CMV and Anelloviridae over the first 40 days post-symptom onset (Extended Data Fig. 5a–c).

Next, we evaluated how viral transcript detection within 40 days of hospital admission associated with COVID-19 severity, using IMPACC trajectory groups²⁴ (Fig. 2a and Extended Data Figs. 5d and 6a). Cumulative link modelling of trajectory groups demonstrated significant associations between COVID-19 severity and Herpesviridae and Anelloviridae transcripts (Supplementary Table 2). Specifically, we found associations between severity and transcripts from Anelloviridae (PBMC adjusted P = 5.02 × 10⁻⁵), CMV (nasal adjusted P = 4.83 × 10⁻⁴, PBMC adjusted P = 1.66 × 10⁻³), EBV (nasal adjusted P = 4.95 × 10⁻⁶, PBMC adjusted P = 2.55 × 10⁻⁹) and HSV1 (nasal adjusted P = 5.29 × 10⁻⁵). When limiting to TG4 participants, those with CMV transcripts in any respiratory compartment were more likely to die within 1 year (nasal adjusted P = 0.039, endotracheal aspirate adjusted P = 0.007). This was also true for TG4 participants with detectable nasal EBV (adjusted P = 0.025) and HSV1 (adjusted P = 0.015) (Fig. 2a). Notably, prevalence of chronic viruses was not significantly different between TG4 and TG5 (Supplementary Table 2).

We then evaluated whether viral transcripts varied with age, adjusting for COVID-19 severity (trajectory groups), and found a significant positive association between Anelloviridae in PBMCs and increasing age (Fig. 2b and Extended Data Fig. 5e, adjusted P = 0.04). Of note, Hispanic ethnicity was significantly associated with both CMV (adjusted P = 0.02) and EBV (adjusted P = 0.02; Extended Data Fig. 6b). However, no significant association were found between viral transcripts and biological sex, or treatment with remdesivir or steroids (Extended Data Fig. 6c–e).

We next evaluated the association of viral transcripts with comorbidities, medications and clinical complications (Fig. 2c and Supplementary Table 2). Anelloviridae transcripts in PBMCs were significantly associated with history of solid organ transplantation, immunosuppressive medications, shock, intensive care unit admission and myocardial infarction. Notably, although medication-mediated immunosuppression was strongly associated with Anelloviridae, these participants only comprised 17.4% of Anelloviridae-positive participants, demonstrating that presence of Anelloviridae is not exclusive to long-term medication-associated immunosuppression (Extended Data Fig. 6f). Furthermore, owing to the correlation between COVID-19 severity, age and immunosuppressive medications, we evaluated all three factors simultaneously and identified that Anelloviridae was independently associated with each (P = 8.6 × 10⁻⁵, 5.0 × 10⁻³ and 2.2 × 10⁻⁸, respectively; Supplementary Table 2).

For Herpesviridae (Fig. 2c), CMV in the nasal compartment was linked to pneumothorax, whereas in the PBMCs, CMV was associated with bacteraemia, pulmonary vascular disease, renal complications, shock and stroke. EBV transcripts in the nasal compartment were associated with intensive care unit-level care and shock, whereas EBV in PBMCs correlated with liver failure, concurrent infections, shock and the overall number of complications. HSV1 in the nasal compartment was significantly associated with acute venous thromboembolism and shock, and inversely associated with liver disease.

We next leveraged our multi-omic data to validate viral reactivation in COVID-19 and characterize the host immune responses, incorporating serum EBV and CMV antibody levels, immune cell frequencies, serum cytokines, plasma metabolomics and host transcriptomics.

First, we observed that participants with EBV transcripts in PBMCs had persistently elevated EBV IgG and IgA antibody titres, further validating the use of viral transcripts as a meaningful measure of viral activity (Fig. 3a and Supplementary Table 3). Similarly, participants with CMV transcripts had significantly higher CMV seropositivity rates at baseline (Fig. 3b, P = 0.0018). Furthermore, using mass spectrometry proteomics, we found that plasma HSV1 proteins were significantly more common in participants with HSV1 transcripts in nasal swabs (Fig. 3c, P = 0.016) and in TG4 participants (Extended Data Fig. 7a, P = 0.0003).

Using CyTOF, we observed significant associations between EBV transcription in PBMCs and increased proportions of B-cell plasmablasts, the primary host cells of the virus (Fig. 3d). Furthermore, we found that CMV transcripts were associated with a significant increase in CD4 and CD8 central memory T cell frequency, and a reduction in CD27low effector memory CD4 T cells. Of note, detection of both EBV (nasal) and CMV (PBMC) was associated with a higher frequency of activated CD4⁺ and CD8⁺ T cells.

To assess whether changes in circulating cell frequencies may explain the association between Herpesviridae and Anelloviridae transcripts and COVID-19 severity, we repeated our analysis across trajectory groups while controlling for immune cell frequencies, which vary with disease severity²⁷. Notably, even after this adjustment, viral transcripts remained significantly associated with COVID-19 severity (Extended Data Fig. 7b).

Next, we investigated whether Herpesviridae and Anelloviridae transcripts were associated with inflammatory protein changes. Using generalized additive mixed modelling (GAMM), we compared longitudinal cytokine dynamics in patients with or without viral reactivation, controlling for severity (trajectory groups), sex and age (Fig. 4a and Supplementary Table 4, adjusted P ≤ 0.01). This approach enabled identification of severity-independent cytokine changes associated with reactivation of different viruses.

We found that different viruses had distinct cytokine profiles (Fig. 4a), although several cytokines overlapped between viruses, including CXCL10 (Fig. 4b), CXCL11 (Fig. 4c) and IL-18 (Fig. 4d).

EBV in PBMCs was associated with increases in key cytokines, including IL-6 (Fig. 4e), CCL7, CCL2, IL-10 (Fig. 4f) and CXCL10, all of which are linked with COVID-19 severity^28,29 (Fig. 4a). Of these, only CXCL10 was also associated with EBV in the nasal transcriptomics. EBV in the nasal compartment was additionally correlated with increases in IL-18, CXCL11, CXCL10, IL-18R1, CD274, IL-15RA, IL22RA1, HGF, IFNγ, CCL8 and MMP10 and a decrease in KITLG, suggesting that host immune response to viruses may differ between compartments, with nasal EBV possibly reflecting a more severe state of viral reactivation. It has previously been observed that increased TGFβ preceded EBV reactivation in COVID-19-associated multisystem inflammatory syndrome in children³⁰, which we did not observe (Extended Data Fig. 7c,d), possibly suggesting that the dynamics of viral reactivation may differ across stages of COVID-19.

Similar to EBV, HSV1 (nasal) and CMV (PBMCs) were associated with a common set of pro-inflammatory serum cytokines (Fig. 4a). HSV1 and CMV both associated with increases in IL-18, CXCL11, CXCL10, IL-18R1, CD274, IL-15RA, CD40, CXCL9, CX3CL1, TNF and CDCP1, and HSV1 also correlated with increased CCL25, SLAMF1, CD5, FGF23, TNFRSF9 and IL-10RB, whereas CMV associated with increased HGF, IFNγ (Fig. 4g), CCL8, CCL7, LIFR, CD8A, ADA and CXCL8 and decreased MMP1 and IL-4. Finally, Anelloviridae were only associated with increases in CXCL11 and IL-18 and decreases in TNFSF11.

We next evaluated the association of viral reactivation with plasma metabolites (Extended Data Fig. 8 and Supplementary Table 5). Viral reactivation was associated with changes in metabolites, particularly those belonging to amino acid and lipid metabolism (Fig. 5a,b). Notably, CMV was associated with the greatest number of metabolome changes, including increased urea and TMAP (N,N,N-trimethyl-l-alanyl-l-proline betaine) levels (Fig. 5a,c and Extended Data Fig. 9a). Additionally, CMV and Anelloviridae were associated with increased long chain fatty acids (for example, erucate, arachidate and docosadienoate) (Fig. 5d and Extended Data Fig. 9b,c) and dimethylarginine (symmetric dimethylarginine and asymmetric dimethylarginine, regulators of nitric oxide synthesis) (Extended Data Fig. 9d). We also identified a shared metabolomic signature across multiple viruses, including Anelloviridae, HSV1, CMV and EBV, with notable reductions in S-methylcysteine sulfoxide and 6-bromotryptophan (Fig. 5e and Extended Data Figs. 8 and 9e).

To identify a signature for each chronic virus independent of COVID-19 severity and participant demographics, we evaluated nasal and PBMC transcriptomic data for differentially expressed host genes while controlling for COVID-19 severity, SARS-CoV-2 nasal viral load, sex, age and days from hospital admission (Supplementary Table 6).

In the PBMC transcriptomics, viral detection in either the PBMC or nasal compartments was associated with changes in PBMC gene expression (Extended Data Fig. 9f,g). Gene set enrichment analysis demonstrated that Anelloviridae and CMV in the PBMCs, and nasal HSV1 were associated with downregulation of diverse RNA processing and protein translation pathways. Similarly, EBV, CMV and HSV1 were associated with an upregulation of cellular replication pathways. Additionally, Anelloviridae and CMV in the PBMCs, and nasal HSV1 and CMV were significantly associated with signatures of neutrophil degranulation, potentially reflecting presence of low-density neutrophils in the PBMCs or involvement of immune genes that are non-specific to neutrophil degranulation. Finally, EBV in the PBMCs was uniquely associated with platelet activation and signalling.

By contrast, nasally reactivated viruses (EBV, CMV and HSV1) had the strongest associations with gene expression changes in the upper airway (Extended Data Fig. 9f,h). Furthermore, CMV, EBV and HSV1 transcripts were all associated with upregulation of interleukin signalling (including IL-10 signalling), lymphocyte immunoregulatory interactions and neutrophil degranulation (Extended Data Fig. 9h). Additionally, upper airway CMV and EBV transcripts were associated with increased anti-parasite and phagocytosis pathways. EBV was also associated with increased expression of genes related to T cell activation. Finally, Anelloviridae in PBMCs was also associated with nasal transcriptome changes, including keratinization and downregulation of proteasome-related proteins.

Given recent reports linking viral reactivation to long COVID^4,22,31, we assessed viral reactivation across patient-reported outcome (PRO) groups, generated from convalescent survey data³². The four PRO groups comprised participants reporting minimal to no-deficits (minimal), physical deficits (physical disability or fatigue), cognitive deficits (cognitive impairment) and global deficits (physical and cognitive).

Upon evaluating viral reactivation during the acute stage of COVID-19, no significant relationship was found with PRO groups (Fig. 6a and Extended Data Fig. 10a,b). However, the sample size was limited owing to the association of chronic viruses with mortality and participant drop-out in the convalescent stage.

We then investigated the correlation between viral detection during the convalescent period and PRO groups. Anelloviridae and enteroviruses were the most frequently detected viruses in convalescent samples (Figs. 1c and 6b and Extended Data Fig. 10c–e). Notably, Anelloviridae transcripts were significantly more prevalent in participants from the physical PRO group, which was characterized by high scores on the PRO Measurement Information System (PROMIS) physical function survey, even when controlling for sex, age, immunosuppressive medications and acute COVID-19 severity (adjusted P = 0.012). This finding suggests that detection of Anelloviridae transcripts may serve as a potential signature for persistent physical disability in patients with long COVID. We further observed that the Anelloviridae-associated gene expression signature was consistent between the acute and convalescent periods (Fig. 6c and Extended Data Fig. 7).

Chronically infecting viruses are prevalent in the general population, with individuals carrying 8–12 viruses at any given time¹. These viruses are generally innocuous and typically do not cause infectious symptoms. However, emerging evidence suggests that their reactivation may contribute to autoimmune diseases, malignancies and chronic fatigue syndrome^1,13,14,15. Additionally, Herpesviridae such as EBV and CMV commonly reactivate in severe infections (for example, malaria or pneumonia) and sepsis^18,33,34 and have more recently been implicated in long COVID^4,22,35. Thus, there is an urgent need to understand how chronically infecting viruses reactivate in COVID-19 and develop strategies to combat their reactivation.

Here we conducted a prospective, multi-omic analysis of 1,154 patients who were hospitalized due to COVID-19, which revealed widespread reactivation of chronic viruses from the Herpesviridae and Anelloviridae families. By integrating cellular and cytokine immunophenotyping, metabolomics and transcriptomics, our findings expand on the complex interplay between SARS-CoV-2 and chronic viral reactivation and provide a deeper understanding of the host immune response. Importantly, our findings demonstrate an association of viral reactivation with COVID-19 and long COVID outcomes.

Despite prior studies evaluating viral reactivation in COVID-19 (refs. ^{2,3,4,5,6,7,22,23}) and other infections^18,33,34, the exact timing and rates of reactivation for different viruses have not been clearly established, with prior studies limited to single time points in small cohorts or relying on antibody titres to identify viral reactivations. Here, leveraging a large longitudinally sampled multi-institution cohort, we define timing and duration of chronic viral reactivations. Specifically, we found that EBV transcripts in PBMCs peaked early in acute disease following hospital admission and then decreased over time. Similarly, Anelloviridae transcripts were most common early during hospitalization and declined several weeks later. By contrast, CMV and HSV1 reactivated later and were primarily in respiratory samples, peaking about 22 days post-hospitalization. Furthermore, we replicated the rates and dynamics of viral reactivations in an external cohort, demonstrating robustness of our findings^25,26. Of note, patients with EBV transcripts in either the nasal or PBMC compartments at admission had higher EBV IgG antibody titres, suggesting that EBV reactivation may happen prior to hospitalization for some patients. Furthermore, since viruses primarily reactivate in tissues (for example, HSV1 and HSV2 in sensory ganglia, CMV in myeloid or dendritic cells and EBV in lymphoid tissue), our results based on systemic detection (PBMCs) may underestimate tissue-specific reactivations. Collectively, our results provide novel insights into the dynamics of the diverse virological landscape of COVID-19 signifying ‘dysvirosis’ (a dysregulation of the human virome) and reveal that the time course of reactivation varies for different viruses.

Exactly how and whether viral reactivations contribute to COVID-19 clinical outcomes remains unknown. Although our data cannot causally determine the direct contribution of viral reactivations to acute COVID-19 disease severity, we identified several key observations that implicate viral reactivation in contributing to disease pathology.

First, multiple cytokines and chemokines, including IL-6, IL-10, CXCL10 and CXCL11, which have been previously associated with COVID-19 pathology, were elevated in participants with viral reactivations when controlling for COVID-19 severity. Notably, increased IL-10 associated with EBV and CMV, which are known to upregulate human IL-10 and produce viral IL-10, which may further displace human IL-10 from IL-10R^36,37. We also identified increased circulating activated CD4⁺ and CD8⁺ T cells and increased cellular replication gene expression in PBMCs, suggesting expanding lymphocytes potentially responding to these viral reactivations. Furthermore, we observed that for PBMC viruses, the associated host inflammatory changes were limited to the blood and did not generalize to other compartments where the virus was not detected (such as the nasal compartment). Although it is possible that chronic viral reactivation may be secondary to disease severity, reactivations correlated with expected inflammatory changes of a responding immune system. Additionally, we observed increased EBV IgG titres at time of hospital admission, suggesting that EBV reactivation likely occurs prior to severe COVID-19 symptoms, raising the suspicion of the potential role of EBV in disease pathophysiology. Thus, the integration of findings across our multi-omics datasets suggests that viral reactivation may not just be epiphenomena but may magnify ongoing inflammation and evoke additional immune responses.

Whereas viral reactivation has frequently been observed in immunosuppressed patients^{38,39,40,41,42,43}, we find that viral reactivation occurs frequently in immunocompetent patients with COVID-19. If viral reactivations were only indicative of a dysregulated immune system, we would expect reactivation to associate not only with severity, but also with immunosuppressive status, which we did not observe for the Herpesviridae in our cohort. Thus, viral reactivation may represent a broader phenomenon than previously appreciated and should be considered in a wider range of clinical contexts beyond immunosupression.

Further integrating our multi-omic dataset, we observed associations between chronic viral reactivations and plasma metabolites. In particular, levels of methylcysteine sulfoxide and 6-bromotryptophan, which are involved in protection against oxidative stress^44,45, were negatively associated with Herpesviridae, suggesting that viral reactivation may contribute to oxidative stress and cellular damage. Furthermore, 6-bromotryptophan has previously been associated with COVID-19 complications^27,46 and impaired kidney function⁴⁴, and its persistently low levels in the context of viral reactivation could imply additional risk to kidney health. This finding aligns with sustained increases of urea and TMAP in CMV-infected patients, which is indicative of renal impairment. Additionally, these data raise the possibility that monitoring for CMV reactivation may be beneficial in patients with COVID-19 presenting with acute kidney damage. The association between CMV and renal dysfunction is also interesting, as CMV infection is one of the most common viral complication in kidney transplant recipients⁴³, suggesting that future studies should evaluate the role of CMV reactivation across different renal disorders. Collectively, these integrated data across our multi-omic study suggest that viral reactivation is associated with specific metabolic perturbations, possibly reflecting strategies by pathogens to manipulate host cell metabolism to their advantage.

Our findings raise important questions about clinical applications of viral reactivation monitoring to guide patient care in acute infections such as COVID-19. Notably, there is already validated infrastructure to evaluate chronically infecting viruses, via quantitative PCR for different Herpesviridae and Anelloviridae (for example, Torque teno virus (TTV)). Furthermore, the increasing affordability of metagenomic sequencing could further broaden the range of detectable viruses. Additionally, for Herpesviridae in particular, there are existing antiviral therapies and several antivirals in clinical trials. Thus, these diagnostic and treatment tools could be readily adapted into broader viral reactivation monitoring and treatment strategies as future work develops specific protocols in patient care.

Importantly, COVID-19 can have ramifications that extend beyond the acute period into persistent sequelae, known as long COVID. Patients with long COVID develop heterogeneous symptoms that may persist for months to years after acute COVID-19 that profoundly affect individuals’ overall health, often resulting in disability and loss of income^32,47,48. Long COVID is estimated to affect 10–30% of individuals after COVID-19 (ref. ⁴⁹). One of the emerging factors associated with long COVID is chronic viral reactivation, particularly reactivation of EBV^4,22,31.

We demonstrate that EBV reactivation is extensive in severe acute COVID-19 and establish the timing of EBV reactivation relative to SARS-CoV-2 infection and hospitalization. Although we did not find a higher rate of EBV transcripts in the acute phase of COVID-19 for patients who developed long COVID, this result raises the question of why increased EBV antibody titres have been frequently observed in long COVID^4,22 and whether the decoupling of viral load in the acute period and antibody titre in the convalescent period reflects an underlying pathology in the immune response. Of note, the PROs in this study were designed early in the pandemic prior to the definition of long COVID and thus may not fully capture current long COVID phenotypes⁵⁰, which may have limited our ability to detect association between EBV reactivation and long COVID. Thus, future work evaluating the dynamics of antibody titres versus viral transcription may continue to reveal the role of EBV in long COVID.

We also report a novel association between the long COVID physical PRO group and Anelloviridae transcripts post-hospitalization, even when controlling for age, chronic immunosuppression by medications and acute COVID-19 severity. Anelloviridae are a large family of negative-sense single-stranded DNA viruses, and some members of this family (such as TTV) are found in 80–90% of the population⁵¹. Anelloviridae have been previously linked with chronic conditions, such as chronic fatigue syndrome and multiple sclerosis^52,53. These conditions often present with physical symptoms similar to those reported by patients with long COVID, including fatigue, cognitive dysfunction and post-exertional malaise. More broadly, Anelloviridae have also been associated with immunosenescence states⁵¹, suggesting that increased viral transcription may underscore a potential state of immune dysfunction and immunosenescence in long COVID. Previously, the role of Anelloviridae in disease has been debated owing to the lack of an identified pathogenic role⁵¹. However, in our study we observed an association between Anelloviridae reactivation and enrichment of genes involved in neutrophil degranulation, suggesting that Anelloviridae may contribute to an inflammatory host response. Additionally, TTV has been shown to replicate preferentially in activated T cells, raising a concern for possible pathogenicity rather than just a passive ‘passenger’ role⁵⁴. Our data highlight the need for future research to dissect the role of Anelloviridae in patients with long COVID.

There were several limitations to our study. First, the use of transcripts to identify viral loads is less common than RT–qPCR. However, the extensive sequencing depth of our samples (25–50 million reads) provided sufficient viral identification. Additionally, we only had six time points during acute COVID-19, limiting the ability to identify viral reactivation between sample collection. Nonetheless, with more than 1,000 participants, we were able to calculate the global longitudinal reactivation patterns for viruses. Furthermore, we only sequenced three tissues, which might have missed reactivation in other compartments, and endotracheal aspirate samples were only collected in ventilated patients, limiting our ability to analyse reactivation in this compartment. Similarly, for enteroviruses and Anelloviridae, owing to the sparsity of individual virus species, we collapsed reads at the genus and family level, respectively, which may have failed to capture more specific relationships for individual species. Another limitation was participant drop-out during the convalescent stage, particularly in patients with acute viral reactivation, limiting the power of analyses testing the association of acute viral reactivation with long COVID. Finally, IMPACC participants were unvaccinated, hospitalized and primarily exposed to ancestral SARS-CoV-2 strains; thus, future studies should ascertain whether more recent strains similarly drive viral reactivation in acute COVID-19 and long COVID in populations with hybrid immunity or milder disease.

In this study, we integrated clinical, immunologic, virologic and multi-omic data from cellular and cytokine immunophenotyping, metabolomics, proteomics and transcriptomics in a longitudinal cohort of 1,154 patients with COVID-19 to investigate for reactivation of chronic viruses and their association with clinical outcomes in one of the largest COVID-19 studies to date. We found that multiple chronic viruses reactivate during acute COVID-19 infection, particularly from the Herpesviridae and Anelloviridae families. Of note, immunosuppression did not explain the high rates of viral reactivation, suggesting that viral reactivation may have a larger role in COVID-19, even in immunocompetent individuals. Furthermore, we delineated the dynamics and rates of reactivation for various chronic viruses and validated these findings in an external cohort. Additionally, we report association of viral reactivation with the host immune response and molecular pathways, as well as acute and chronic clinical sequelae of COVID-19. Notably, our results raise the possibility that viral reactivation may contribute to the development of severe acute COVID-19 and long COVID. Our findings underscore the emerging role of viral reactivations, highlighting the potential contribution to hyperinflammation, and raise the possibility that management of viral reactivations in conditions such as acute COVID-19 and long COVID could improve patient outcomes.

IMPACC is a prospective longitudinal study that enrolled >1,000 hospitalized patients with COVID-19, as previously described^{27,32,55,56,57,58}. Participants 18 years and older were recruited from 20 hospitals across 15 academic institutes within the USA (Fig. 1a). All participants were confirmed to be SARS-CoV-2-positive by reverse transcription PCR (RT–PCR) testing and no participants were vaccinated for SARS-CoV-2 at time of enrolment. Nasal swabs, blood, and endotracheal aspirate (for ventilated patients) were collected within 72 h of hospital admission (visit 1) and on days 4, 7, 14, 21 and 28 post-hospital admission in addition to convalescent samples at 3, 6, 9 and 12 months. Additionally, nasal swabs, blood and endotracheal aspirate were also collected (when possible) within 24 h and 96 h of escalation to intensive care unit-level care or when a participant was readmitted to the hospital >48 h after discharge. In the IMPACC dataset, these samples were indicated as ‘escalation visits’, and in this Article, the escalation visit samples were only used for the longitudinal GAMM and the CyTOF cellular association analyses to minimize effects of more frequent sampling for critically ill participants. Participants were characterized into one of five trajectory groups based on latent class mixed modelling of a seven-point ordinal scale that characterized degree of respiratory illness and reflected acute COVID-19 severity²⁴. Similarly, participants were clustered into 4 PRO groups from latent class mixed modelling trajectories of convalescent survey responses collected at 3, 6, 9 and 12 months³². The modelling used participant responses from the EQ-5D-5L⁵⁹, health recovery score (visual analogue scale of 1–100 to indicate overall physical and mental function compared to pre-COVID function), and the following PROMIS forms (https://commonfund.nih.gov/promis): PROMIS Item Bank v2.0—Physical Function, PROMIS Item Bank v2.0—Cognitive Function⁶⁰, PROMIS Scale v1.2—Global Health Mental 2a⁶¹, PROMIS Item Bank v1.0—Psychosocial Illness Impact-Positive—Short Form 8a⁶¹, and PROMIS Pool v1.0—Dyspnea Time Extension⁶². All PROMIS measures were scored and standardized following PROMIS standardized instructions. Clinical characteristics and demographics for the entire cohort are reported in Supplementary Table 1.

Ethics statement

The Department of Health and Human Services Office for Human Research Protections (OHRP) and the National Institute of Allergy and Infectious Diseases (NIAID) concurred that the IMPACC study qualified for public health surveillance exemption. The study protocol was sent for review to each site’s institutional review board (IRB), with twelve sites conducting as a public health surveillance study, and three sites integrating the IMPACC study into IRB-approved protocols (The University of Texas at Austin, IRB 2020-04-0117; University of California San Francisco, IRB 20-30497; Case Western Reserve University, IRB STUDY20200573) with participants providing informed consent. Participants enrolled at sites operating as a public health surveillance study were provided information sheets describing the study including the samples to be collected and plans for analysis and data de-identification. Participants who requested not to participate after review of the study plan and information were not enrolled. Participants were not compensated while hospitalized but were subsequently compensated for outpatient visits and surveys. This study was registered at clinicaltrials.gov (NCT04378777) and followed the strengthening the reporting of observational studies in epidemiology (STROBE) guidelines (Extended Data Fig. 1a).

Sample processing and assays

Samples were processed as previously described^{27,32,55,56,57,58}, with the sample protocol extensively documented in the IMPACC study design and protocol paper⁵⁵. In brief, 10 ml of blood and nasal swabs were collected at each visit, with blood processed within 6 h of collection. Blood was collected in both a 2.5 ml Greiner Vacuette CAT Serum Separating Tube (SST) (454243P) for serum and a 7.5 ml Sarstedt Venous blood collection monovette EDTA (NC9453456) for whole blood, PBMCs and plasma. The SST was kept vertical at room temperature for at least 30 min before centrifuging at room temperature for 10 min at 1,000g. Serum was then aliquoted at 100 μl for downstream assays.

From the EDTA tube, it was briefly inverted to mix before aliquoting 270 μl of whole blood for CyTOF. The 270 μl CyTOF aliquot was added directly to the Maxpar Direct Immune Profiling Assay tube (MDIPA antibodies listed in Supplementary Table 9) and incubated for 30 min at room temperature. After incubation, 410 μl of Smart Tube Prot1 Stabilizer (SmartTube) was added with a 10 min incubation at room temperature before storage at –80 °C until shipment to the respective processing core. The remaining blood was centrifuged at room temperature for 10 min at 1,000g before aliquoting and storing 500 μl of plasma at –80 °C for proteomic and metabolomics. PBMCs were then isolated from the remaining sample using the SepMate and Lymphoprep system (StemCell) following manufacturer protocol and as previously described⁵⁵. PBMCs were then stored at 2.5 × 10⁵ cells in 200 μl of RLT Buffer (Qiagen) and β-mercaptoethanol at –80 °C.

Interior nasal turbinate swabs (herein referred to as nasal swabs) were collected and stored in 1 ml of Zymo-DNA/RNA shield reagent (Zymo Research), before RNA was extracted twice in parallel from 250 μl of sample and purified with the KingFisher Flex sample purification system (ThermoFisher) and the quick DNA-RNA MagBead kit (Zymo Research). The duplicated RNA was pooled and aliquoted at 20 μl for the downstream assays (SARS-CoV-2 RT–qPCR and RNA-seq).

When participants were ventilated, an endotracheal aspirate was also collected in a 40 cm³ Argyle specimen trap and was processed within 2 h of collection. First, 500 μl of 1:1 diluted endotracheal aspirate with Maxpar PBS (Ca²⁺ and Mg²⁺ free) was mixed with 500 μl of DNA/RNA shield in a Zymo tube with lysis beads and subsequently stored at –80 °C for bulk RNA-seq.

Collected samples were then shipped and processed for nasal, PBMC, and endotracheal aspirate RNA-seq, plasma proteomics, serum cytokine PEA, serum EBV and CMV antibody titres, whole-blood CyTOF, and plasma metabolomics at their respective processing cores as previously described^{27,32,55,56,57,58}. Each assay is described in brief below with additional technical details in the prior publications^{27,32,55,56,57,58}, and a list of analytes evaluated from the PBMC RNA-seq, nasal RNA-seq, serum cytokines and chemokines, and plasma metabolomics assays is reported in Supplementary Table 8.

Nasal, endotracheal aspirate and PBMC RNA-seq

Nasal transcriptomics data were processed using a workflow managed on the Galaxy platform. RNA extracted from the nasal swabs was DNase-treated and depleted of human ribosomal RNA before amplification with random hexamers (Ribo-Zero Plus kit), except for the convalescent samples (3, 6, 9 and 12 months) which instead used poly-A amplification for library prep (SMART-seq V4). Before sequencing, libraries were quantified on both a Quant-it dsDNA High Sensitivity Assay and Fragment Analyzer (Advaced Analytical; kit ID DNF474), with samples containing adapter dimers as more than 4% of electropherogram area failed before sequencing. Technical controls (K562, Thermo Fisher Scientific, AM7832) were compared across batches to ensure minimal batch variability. Libraries were normalized to 10 nM before base calls were generated on the NovaSeq6000 instrument (RTA v3.1.5) at 100 bp paired-end read length, and demultiplexed unaligned BAM files were produced using Picard’s ExtractIlluminaBarcodes and IlluminaBasecallsToSam tools (https://broadinstitute.github.io/picard/). These BAM files were converted to FASTQ format using Samtools bam2fq (v1.4)⁶³. Adapter trimming and quality filtering were performed using Trimmomatic (v0.36.5)⁶⁴, with reads trimmed by one base at the 3′ end and further trimmed from both ends to ensure a minimum base quality score of Q30. Adapter sequences were also removed. Trimmed reads were aligned to the GRCh38 human reference genome⁶⁵ using STAR (v2.4.2a)⁶⁶ with gene annotations from Ensembl release 91 (ref. ⁶⁷). Gene-level counts were generated using HTSeq-count (v0.4.1)⁶⁸. Quality control metrics were compiled using Picard (v1.134), FASTQC (v0.11.3) (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and Samtools (v1.2). Samples with poor quality (defined as median coefficient of variation (CV) in gene coverage >0.8 or aligned counts <1 million) were excluded from downstream analysis.

RNA extracted from endotracheal aspirate specimens was DNase-treated and depleted of human ribosomal RNA. Complementary DNA (cDNA) was synthesized using random hexamers to capture both coding and noncoding transcripts. Libraries were sequenced with 100 bp paired-end reads on the NovaSeq6000 instrument with NovaSeq S4 flow cells (200 cycles) targeting 50 million reads per sample. Human reads were aligned to the GRCh38 reference genome⁶⁵ and quality-controlled. Raw counts were normalized across libraries using the TMM method implemented in the edgeR package⁶⁹. To control for batch effects, participant and control samples were co-sequenced within each batch.

For the PBMC transcriptomics, RNA was extracted from the 2.5 × 10⁵ cells stored in RLT Buffer (Qiagen) using the Quick-RNA MagBead Kit (Zymo) with DNase digest. RNA quality was then evaluated with both a Quibit HS RNA assay and Fragment Analyzer (Agilent) before cDNA was prepared with the SMART-Seq v4 Ultra Low Input RNA Kit (Takara Bio) from an input of 10 ng of RNA. After a bead-based clean-up, the Nextera XT DNA Library Preparation kit (Illumina) was used to prepare the sequencing libraries which were validated by capillary electrophoresis with a Fragment Analyzer (Agilent). The resulting libraries were pooled at equimolar concentrations and sequenced on an Illumina NovaSeq6000 at 100 bp paired-end read length with a target of 25 million reads per sample. Adapter trimming and quality filtering were performed using Cutadapt (v1.14). Reads were aligned using STAR (v2.4.2a) to a composite reference genome that included the human genome (GRCh38, Ensembl release 91)⁶⁵ and SARS-CoV-2 (NCBI strain MN908947.3)⁷⁰. Gene counts were computed using HTSeq-count internally within the STAR alignment step. Quality control metrics were assessed using FASTQC (v0.11.5), Picard tools (v2.22), and STAR log outputs. QC metrics included average base quality per read (>Q30), per cent and absolute counts of reads uniquely mapped to annotated transcripts, and other alignment-based quality statistics.

Taxonomic alignments for human-infecting viruses from the PBMC, nasal and endotracheal aspirate RNA-seq data were obtained from CZID⁷¹, which removes host reads before aligning remaining reads against the National Center for Biotechnology Information (NCBI) nucleotide and non-redundant databases. A sample was considered positive for a virus if it had at least one read that mapped to both the nucleotide and non-redundant database. At least one water control was included on each plate/batch from all transcriptomics, with no water control samples having detectible reads for the human-infecting viruses evaluated in this Article, supporting the low threshold for positivity. When evaluating for potential batch effects of viruses, we identified one plate from the nasal transcriptomics that had an above average rate of HSV1 positivity (>60% compared to the average batch approximately 10%), which may have been due to cross-contamination from a sample in the batch with an extremely high HSV1 viral load (about 250,000 RPM). As a result, we excluded this one plate from contributing to identifying HSV1 positivity and from all host transcriptomic analyses. Additionally, some samples were sequenced multiple times across batches, in the case of a sample sequenced multiple times, the mean RPM of each virus across the samples was used for that participant event.

In addition to the IMPACC PBMC RNA-seq, RNA-seq data were retrieved from the Mount Sinai COVID-19 Biobank (syn35874390)^25,26. The retrieved raw fastq were processed as outlined above through CZID to identify reads belonging to human-infecting viruses.

Nasal SARS-CoV-2 RT–qPCR

SARS-CoV-2 viral load was assessed from nasal swab samples using 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 were performed using Quantabio One-Step RT–qPCR ToughMix on a QuantStudio 5 instrument. Cycle threshold (Ct) values for N1 and N2 were used as the primary readout.

Whole-blood CyTOF

Whole-blood CyTOF prepared samples were thawed according to the SmartTube Prot1 Thaw/Erythocyte lysis protocol before samples were barcoded with the Fluidigm Cell-ID 20-Plex Palladium Barcoding Kit. After barcoding, samples were pooled and additional surface antibody staining was performed (antibodies listed in Supplementary Table 9) at a final concentration of 1 µl antibody per 10 million cells in a volume scaled to cell number (100 µl for every 10 million cells), with a 30 min incubation on ice. After surface staining, samples were washed twice (each 1 ml of CyFACS (1× PBS + 0.2% BSA + 0.05% NaN₃) at 800g × 3 min) and subsequently fixed and permeabilized with BD Biosciences Fixation/Permeabilization solution (100 µl per 5 million cells) with a 20 min incubation on ice. Samples were then washed twice with 1× BD Perm Wash buffer (800g × 3 min) before staining with an intracellular antibody for GZMB (Supplementary Table 9) at a final concentration of 1 µl antibody per 10 million cells in a volume scaled to cell number (100 µl for every 10 million cells), with a 30 min incubation on ice. Finally, samples were fixed with paraformaldehyde and simultaneous iridium/osmium cell labelling. CyTOF samples were then acquired using the Fluidigm Helios mass cytometer and normalized/concatenated with Fluidigm’s CyTOF software. Further cleaning was done using Mt Sinai’s internal pipeline, which removed acquisition outliers, EQ beads, and low DNA signal events. Demultiplexing was performed using Pd barcoding and cosine similarity, removing low signal-to-noise cells and acquisition multiplets. Each sample was clustered into 1,000 k-means groups. A subset was manually annotated via Clustergrammer2 (https://github.com/ismms-himc/clustergrammer2) to generate a reference matrix. For annotation, cluster similarity to reference cell types was computed, and assignments were made based on highest or consensus similarity. Cell-type labels were then mapped back to single cells for downstream quantification. Antibodies used are listed in the supplemental reporting summary.

Plasma proteomics

Plasma samples underwent protein depletion using perchloric acid to remove the most abundant proteins, enhancing the detection of lower abundance proteins^72,73,74. The prepared samples were loaded onto Evotips and analysed using the EVOSEP One system (EVOSEP). The system operated using the 60 samples per day method, which corresponds to a 21 min gradient, optimizing throughput without compromising data quality. The EVOSEP One was coupled to a timsTOF Pro mass spectrometer (Bruker Daltonics) operating in Data-Dependent Acquisition Parallel Accumulation–Serial Fragmentation (DDA-PASEF) mode. HSV1 proteins evaluated included the Swiss-Prot reviewed proteins attributed to human herpesvirus 1 (strain 17, UniProt Proteome ID: UP000009294).

Plasma metabolomics

Plasma metabolite profiling was performed by Metabolon using their in-house standards^75,76. Samples were randomized, extracted with methanol precipitation, and divided into fractions for analysis via UPLC–MS/MS under both positive and negative ion modes using RP and HILIC chromatography. Analyses were conducted with a Thermo Q-Exactive mass spectrometer. Metabolites were identified based on retention index, accurate mass (±10 ppm), and MS/MS spectral matching to Metabolon’s reference library, following Metabolomics Standards Initiative guidelines⁷⁶.

Serum cytokines and chemokines

All samples were analysed using the Olink Inflammation panel (Olink Bioscience), which measures 92 inflammation-related proteins via PEA. In brief, oligonucleotide-labelled antibody pairs bind to target proteins, enabling formation of PCR targets upon proximity. After overnight incubation at 4 °C, PCR amplification was performed, and protein levels were quantified using a microfluidic qPCR system (Biomark, Fluidigm), including built-in controls for quality assurance.

Serum EBV and CMV antibody titre

Serum anti-viral antibodies for EBV (gp350) and CMV (gB) were measured via a Luminex platform as previously described⁵⁸. The gp350 and gB antigens (Sino Biological) were conjugated to barcoded beads as recommended by Luminex. Serum samples were diluted 1:400 in PBS/0.5% Triton X-100 before 25 μl was added to an assay plate containing the antigen-coupled beads and incubated on an orbital shaker at 500–600 rpm for 2 h at room temperature. Each well also had Assay Chex Control beads (Radix Biosolutions) to measure non-specific binding. The plate was then washed with a Bio-Tak Magnetic washer (ELX-405, Bio-Tek) before a secondary antibody goat-anti human IgG (Fc fragment, NC9822979) or Anti-IgA Fc coupled to Phycoerythrin (501941614) were diluted and added to the plate and incubated on an orbital shaker at 500–600 rpm for 30 min at room temperature. The plate was then washed again with again with a Bio-Tak Magnetic washer (ELX-405, Bio-Tek) before 130 μl of wash buffer was added to the wells and read on a Luminex Flex 3D instrument (lower bound of 50 beads per target antigen). This assay was only performed for half of the cohort (n = 479). The data for the EBV titre values were normalized by regressing the values of the four different Assay Chex control beads as well as the batch and plate before analysis.

For an additional 497 participants, CMV serostatus at visit 1 was also measured using a CMV IgG ELISA assay (Aviva, GWB-BQK12C). Serum samples were assessed in duplicate, with an average antibody index value >1.1 considered positive, <0.9 considered negative, and values between 0.9 and 1.1 equivocal per manufacturer specifications. For this study, equivocal status was considered seropositive. The resulting calculated serostatus from both CMV assays were combined to evaluate rate of serostatus at hospital admission between CMV transcript positive and negative patients (Fig. 3b).

Statistics and reproducibility

All IMPACC sites followed a standard protocol⁵⁵ for biological sample collection. To mitigate batch effects, samples were randomized to batches while ensuring longitudinal samples from the same individual were run in the same batch. The randomization was stratified by disease severity (mild and moderate versus severe) and age (younger versus older) with representation across batches. Furthermore, race, ethnicity, gender, and enrolment site were then confirmed to be well represented across batches. As this was an observational study, we prioritized biological replicates over technical replicates, except where needed to evaluate for batch effects. Thus, all data used in this study were reflective of data collected as single measurements for each participant at each time point, unless otherwise specified in the methods.

All P values calculated in this manuscript were adjusted using the Benjamini–Hochberg procedure and are indicated as adjusted P values (circumstances where no adjustments were necessary report as only P values). Across all multi-omic analyses, adjustments were corrected for all comparisons or models ran on a per virus basis. The exact statistical approaches used for each analysis are detailed below.

Clinical features and demographics

For testing the association of detection of chronic viral transcripts with trajectory groups and age, we used cumulative link mixed modelling (clmm) from the ordinal (v2019.12-10) R package⁷⁷. Due to both COVID-19 severity and age quantiles being ordinal, cumulative link mixed modelling allowed for this ordinal relationship to be accounted for in addition to including enrolment site as a random effect. However, due to the TG4 and TG5 contributing a majority of the endotracheal aspirate samples, for comparison of the endotracheal aspirate viruses a Fisher’s exact test comparing only TG4 and TG5 was used instead of a clmm. The following R formulae were used with the clmm2 function from the ordinal (v 2019.12-10) R package⁷⁷:

For the testing of association of virus positivity status with ethnicity, sex, and steroid and remdesivir administration in Extended Data Fig. 6b–e, a chi-squared test of independence was calculated using the rstatix (v0.7.2) R package⁷⁸.

For the testing of virus positivity status with long-term mortality within trajectory group 4, a Cox proportional hazards model was used that included steroid and remdesivir administration as fixed effects and enrolment site as a random effect. All clinical outcomes, including mortality, were captured by clinical staff from the electronic medical record in accordance with protocol and then verified by the IMPACC Clinical and Data Coordinating Center. If death was not ascertained, the model was right censored with the date of participant lost to follow up or completion of the study. The following R formula was used with the coxme function from the coxme R package (v2.2-16)⁷⁹:

For testing of association with other clinical features including complications, comorbidities and medication usage, a logistic mixed effect model was used that also included sex, age quintile and trajectory group as fixed effects and enrolment site as a random effect. Of note, participants positive for a virus in a given tissue were only compared against participants who also had the respective tissue sequenced at least once and were negative at all time points checked. The following R formula was used with the glmer function from the lme4 R package (v 1.1-28)⁸⁰:

To further decouple association of Anelloviridae with COVID-19 severity, immunosuppressive medications and age, we conducted an additional logistic mixed effect model that included sex, age quintile, trajectory group and immunosuppressive medication status as fixed effects with enrolment site as a random effect. The following R formula was used with the glmer function from the lme4 R package (v 1.1-28)⁸⁰:

For testing of association of virus detection status with the PASC PRO groups, a linear model was used that also included sex, age quintile, trajectory group and immunosuppressed status by medication as fixed effects in the model. The following R formula was used with the glmer function from the lme4 R package (v 1.1-28)⁸⁰:

Cellular immunophenotyping

For testing the association of participant events positive for viruses with changes in immunophenotypes of circulating cells, a linear mixed effect model was used that also included trajectory group, sex, age quintile, and visit number as additional main effects in addition to enrolment site and participant as nested random effects in the model. The normalized abundances for cell types were computed by calculating the relative abundance after removing granulocytes from the total before a log1p transformation and scaling. The following R formula was used with the lme function from the nlme (v 3.1-149) R package⁸¹:

Serum cytokines and plasma metabolomics

To evaluate both the serum proximity extension array cytokine assay (Olink Inflammation panel) and the plasma metabolomics, GAMM from the gamm4 (v 0.2-6) R package⁸² was used to evaluate for differences in individual analytes between patients who had detected transcripts for a chronic virus compared to patients who never had human-infecting viral transcripts detected (other than for SARS-CoV-2). Analytes were modelled against days from admission using cubic regression splines with interactions of both status for a given chronic virus (that is, detected transcript at any collected time point) and trajectory group in addition to fixed effects of status for the chronic virus being evaluated, trajectory group, sex, age at time of admission sorted into quintiles, and SARS-CoV-2 nasal viral RPM at the time of that sample. As chronic viral status is both a main effect and interaction term in the model, we used a lower adjusted P value cutoff of 0.01 to account for the fact that a feature was significant if either term was significant. Only participant visits with both PBMC and nasal transcriptomics were used for this analysis. The following R formula was used with the gamm function from the gamm4 (v0.2-6) R package⁸²:

Nasal and PBMC RNA-seq

To analyse the signature of host gene expression associated with chronic viruses, the limma (v 3.46.0) R package⁸³ was used for both the nasal and PBMC RNA-seq to evaluate differential expressed genes associated with participants who had detectable chronic viral transcripts. The following R formula was used:

Upregulated and downregulated differentially expressed genes were then evaluated separately with hypergeometric pathway enrichment using the Reactome database (v95)⁸⁴ and the clusterProfiler R package (v 3.18.1)⁸⁵.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.