The birth order, defined as the ordinal position of a child among siblings, has fascinated researchers for more than a century^1,2. Theorists in the early years proposed that first-borns receive greater parental investment and face higher expectations, while later-borns develop in the immunological and social wake of their older siblings^3,4. These ideas have generated a rich but fragmented empirical literature, with individual studies examining birth order in relation to specific diseases or developmental outcomes.

The most influential disease-specific finding concerns allergic and atopic conditions. Strachan’s original household-size observation proposed that younger siblings, owing to greater microbial exposure from older siblings during early life, develop stronger immune tolerance and lower allergy risk⁵. Subsequent studies confirmed the protective effects of subsequent birth order for hay fever, eczema and asthma^6,7,8,9. The mechanistic understanding has since evolved: the ‘old-friends’ hypothesis emphasizes that early exposure to diverse commensal and environmental microorganisms, rather than pathogenic infections per se, is critical for immune regulatory development^10,11,12.

A separate literature has examined birth order and neurodevelopmental conditions. Large Scandinavian registry studies observed that first-borns show slightly higher educational attainment and IQ, potentially reflecting differential parental investment^13,14,15. For autism, the relationship with birth order is complicated by reproductive stoppage or the tendency of parents to reduce the rate of childbearing after a diagnosis. This can create artifactual first-born enrichment^16,17. Interpretation is also complicated by parental-age effects, particularly the association between advanced parental age and autism risk¹⁸. Birth order has also been associated with psychiatric conditions and childhood mental health outcomes¹⁹, as well as metabolic diseases²⁰.

Despite this extensive literature, three limitations have impeded progress. First, nearly all studies examine a single disease or a small cluster, precluding systematic comparison of effect sizes across organ systems. Second, many existing designs have limited ability to separate birth-order associations from factors that covary with birth order, including parental age, family size, socioeconomic status and secular diagnostic trends^21,22. Third, sample sizes have often been too small for precise estimation, particularly for rarer conditions.

In this study, we address these limitations using a cohort of over 10 million individuals from 5.1 million two-child families in the Merative MarketScan commercial claims dataset. We use two complementary analytical designs: a between-family matched cohort used as a high-powered phenome-wide scan after adjustment for measured demographic, geographical, parental-age, follow-up and clinical covariates; and a within-family sibling comparison using conditional logistic regression as an internally controlled contrast that mitigates confounding by factors shared within families^23,24. We tested 569 diseases, applying sensitivity analyses, including state fixed effects, full-sibling restriction and a stricter clinically enriched rematching analysis, and validate our approach with prespecified positive and negative control diseases.

We screened 27,975,854 Merative MarketScan 2003–2024 families with at least two age-window candidate members and identified 5,135,006 two-child families (10,270,012 individuals) meeting our eligibility criteria: at least one inferred parent, exactly two non-parent children, each with 365 or more days of enrollment visibility and age at last observation of 12 years and older (Fig. 1a and Table 1). For family-size context before parent inference and individual eligibility filtering, 66,054,423 family identifiers contained at least one age-window candidate member; 57.6% contained one, 24.7% contained two, 10.8% contained three, 4.6% contained four and 2.2% contained five or more. Among the 27,975,854 family identifiers with at least two valid-sex and birth-year candidate children, 58.3% had two candidate children and 41.7% had three or more.

For the between-family analysis, we identified 1,616,881 matched sibling pairs by pairing a first-born from one family with a second-born from a different family, matched exactly on sex, birth year and urbanization tertile, with calipers on follow-up duration (± 50 days), paternal age (± 10 years), maternal age (± 10 years), and sibling age gap (± 2 years). After matching, standardized mean differences (SMDs) improved for all covariates (Fig. 1c). The strict clinical rematch further reduced the imbalance for parental-age and clinical baseline variables, while a separate within-family cohort of 5.1 million families was used for sibling comparisons. The underlying cohort was predominantly from high-urbanicity areas, drawn from across all four census regions, with birth years spanning 1978–2013 (Fig. 1b).

For the within-family analysis, we used all 5,135,006 cohort families directly, comparing the first-born to the second-born within each family using conditional logistic regression stratified on family identifier. This design reduces confounding by factors shared between siblings (for example, parental genetics, household socioeconomic status, geographical exposures, family health attitudes), at the cost of being powered only by disease-discordant sibling pairs²⁵.

Signal yields across the four analytical designs are summarized in Fig. 1d: the primary between-family design detected 150 Bonferroni-significant and 242 nominally significant diseases; the state fixed-effects specification, the stricter clinical rematch and the within-family design each produced broadly consistent counts.

Throughout, odds ratios (ORs) compare second-borns with first-borns: an OR below 1 denotes lower risk in second-borns (equivalently, increased first-born risk; ‘first-born excess’), whereas an OR above 1 denotes increased risk in second-borns (‘second-born excess’). Birth-order associations run in both directions across the phenome.

Phenome-wide birth-order atlas

To display the Bonferroni-significant birth-order associations that were also Bonferroni-significant and directionally concordant in the within-family analysis, we constructed a disease atlas organized according to clinical domain (Fig. 2 and Extended Data Fig. 1). In this atlas, each disease tile is colored according to the direction and magnitude of the birth-order effect, with blue indicating first-born excess and red indicating second-born excess. Tile color intensity is proportional to the absolute effect size (\(| {\mathrm{log}}_{2}(\mathrm{OR})|\)); only diseases reaching Bonferroni significance in both the primary between-family and within-family analyses with concordant direction are displayed.

The atlas across five key domains (Fig. 2) reveals that first-born excess is concentrated in the neuropsychiatric domain, whereas second-born excess is prominent in musculoskeletal, infectious and neurological diseases. Within the neuropsychiatric domain, first-born excess is observed across a broad diagnostic spectrum from the other/unspecified PDD code group and tics/Tourette syndrome (strongest effects, \(| {\mathrm{log}}_{2}(\mathrm{OR})| > 0.5\)) through autism, obsessive-compulsive disorder (OCD) and attention-deficit/hyperactivity disorder (ADHD) to milder effects such as anxiety, eating disorders and depression, while substance abuse is a notable exception with second-born excess.

The expanded atlas across all displayed clinical domains (Extended Data Fig. 1) additionally reveals dermatological first-born excess (acne, hirsutism, seborrheic dermatitis), respiratory associations (first-born excess for asthma and allergic rhinitis), endocrine and metabolic associations (first-born excess for lipid metabolism disorders and pubertal dysfunction, second-born excess for electrolyte/acid–base disorders), digestive second-born excess (gastritis and duodenitis, irritable bowel syndrome, appendiceal and esophageal disease) and musculoskeletal second-born excess concentrated in joint connective tissue conditions.

Domain-level summary

The distribution of significant birth-order effects across the 15 noncongenital/non-injury clinical domains displayed in the domain-level visualization is summarized in Extended Data Fig. 2. The atlas in Extended Data Fig. 1 displays 75 concordant Bonferroni-significant diseases across 14 clinical domains. The neuropsychiatric domain contributed the largest number of significant associations, with a striking predominance of first-born excess (Extended Data Fig. 2a). Dermatological and sense organ categories also showed predominantly first-born excess. In contrast, digestive, musculoskeletal, genitourinary, circulatory and infectious disease domains were enriched for second-born excess. Several domains, including respiratory and endocrine and metabolic, showed mixed directionality.

Within-domain effect size distributions (Extended Data Fig. 2b) reveal that the neuropsychiatric domain shows the widest spread of effect sizes, with median effects shifted toward first-born excess. The digestive and musculoskeletal domains show median effects that have shifted toward second-born excess. Most domains have median effects close to null, reflecting a mixture of excess diseases from the first-born and second-born within each category. To quantify domain-level clustering rather than relying only on visual inspection, we fitted an empirical-Bayes partial-pooling model to disease-level log-ORs and standard errors within each design (Extended Data Fig. 3). The strongest domain-level first-born shift was in the neuropsychiatric and behavioral domain, with pooled ORs of 0.902 (95% confidence interval (CI) 0.875–0.930) in the primary between-family scan, 0.914 (0.883–0.945) in the strict rematch and 0.937 (0.919–0.956) in the within-family design. Dermatological associations also showed consistent first-born shifts across between-family and within-family analyses. By contrast, digestive, genitourinary and reproductive, and musculoskeletal domains showed second-born shifts in the between-family analyses that were weaker in the within-family design.

Phenome-wide landscape

We defined 569 diseases using established International Classification of Diseases, Ninth Revision, Clinical Modification (ICD-9-CM) and International Classification of Diseases, Tenth Revision, Clinical Modification (ICD-10-CM) code groupings. Of these, 418 had 500 or more cases in the matched cohort and were included in the between-family analysis. Logistic regression adjusted for sibling age spacing, age at last observation, sex, parental ages, parental psychiatric history, urbanization, county (via clustered standard errors), ICD coding era, enrollment time and a full-sibling consistency flag.

Of these 418 diseases, 150 (35.9%) reached Bonferroni significance (P < 1.20 × 10⁻⁴) and 226 (54.1%) reached significance after Benjamini–Hochberg false discovery rate correction at q < 0.05. The landscape across the phenome (Fig. 3) displays all Bonferroni-significant diseases positioned according to prevalence and effect size, with the point size proportional to the number of cases and the colors denoting the category of the disease. Among the 150 Bonferroni-significant diseases, 79 showed first-born excess (OR < 1 for second-born) and 71 showed second-born excess (OR > 1); this deviation from a 50:50 split was not significant (exact binomial P = 0.568). The high rate of significant associations, together with the near-symmetric split between first-born and second-born excess, argues against a systematic bias inflating associations in one direction.

The landscape reveals that the largest effect sizes arise among rarer conditions: the other/unspecified PDD code group (prevalence < 1%, \({\mathrm{log}}_{2}(\mathrm{OR})\approx -0.81\)), tics/Tourette syndrome (\({\mathrm{log}}_{2}(\mathrm{OR})\approx -0.53\)), autism (\({\mathrm{log}}_{2}(\mathrm{OR})\approx -0.44\)) and OCD show pronounced first-born excess, while herpes zoster shows the strongest second-born excess (\({\mathrm{log}}_{2}(\mathrm{OR})\approx +0.43\)). Among highly prevalent conditions, ADHD, allergic rhinitis, asthma and acne show modest but precisely estimated first-born excess, while substance abuse and migraine show second-born excess (Fig. 3).

Strongest birth-order associations

The strongest and most clinically recognizable associations, together with their stability across alternative models, are summarized in Fig. 4, Table 2 and Supplementary 1.

First-born disease risk excess was most pronounced for neurodevelopmental conditions: the other/unspecified PDD code group (OR = 0.569, 95% CI 0.552–0.587, P = 2.7 × 10⁻²⁷⁸), tics/Tourette syndrome (OR = 0.693, 95% CI 0.672–0.715, P = 4.0 × 10⁻¹¹⁶) and autism (OR = 0.737, 95% CI 0.714–0.760, P = 1.2 × 10⁻⁸¹). First-born excesses were also observed for food allergy (OR = 0.797, P = 9.7 × 10⁻⁷³), acne (OR = 0.866, P < 10⁻³⁰⁰), anxiety/phobic disorder (OR = 0.889, P = 1.4 × 10⁻¹⁵¹) and allergic rhinitis (OR = 0.910, P = 5.4 × 10⁻⁹⁷).

Because this leading phenotype carried the historical source label ‘unspecified childhood psychoses,’ we audited its ICD definition before interpreting it biologically. The implemented phenotype consists of ICD-9 299.8x/299.9x and ICD-10 F84.8/F84.9 codes, corresponding to other or unspecified PDD codes rather than schizophrenia-spectrum psychosis codes. It contributed 18,557 cases to the primary matched analysis and 5,002 cases to the strict clinical rematch. Across the broader cohort, this code group included 46,640 cases, of whom 20,603 (44.2%) also had an autism-spectrum-disorder phenotype in the current disease map, indicating substantial but incomplete phenotypic overlap with the autism signal. Therefore, we refer to this association as an other/unspecified PDD code-group signal; this relabeling clarifies phenotype interpretation but does not change the estimated birth-order association.

Second-born excess was strongest for herpes zoster (OR = 1.348, P = 4.7 × 10⁻¹⁰⁰), substance abuse (OR = 1.192, P = 3.8 × 10⁻²²⁷), biliary tract disease (OR = 1.179, P = 6.2 × 10⁻⁷⁰), gastritis and duodenitis (OR = 1.142, P = 4.3 × 10⁻⁸⁵) and migraine (OR = 1.128, P = 2.3 × 10⁻¹⁰⁷).

Within-family sibling comparison

We used conditional logistic regression (clogit) for the analysis of the within-family cohort. We adjusted for family-specific effects, birth cohort, age at last observation in the data, sex, length of enrollment, ICD coding era and calendar period of observation (5-year bins based on the midpoint of each child’s observation window). Of 569 diseases, 541 had 100 or more disease-discordant sibling pairs and were analyzed further. We fitted four model specifications per disease to assess robustness to age–period–cohort (APC) parametrization including (1) cohort-adjusted, (2) cohort plus calendar period, (3) period-only and (4) cohort with gap × birth-order interaction. The cohort-plus-period specification served as our primary within-family model (Supplementary Table 3).

The within-family results broadly corroborated the between-family findings. Among diseases significant in both designs, the direction and magnitude of birth-order effects were consistent: autism within-family OR = 0.804 (95% CI 0.786–0.823, P = 1.1 × 10⁻⁷⁴), ADHD within-family OR = 0.936, food allergy within-family OR = 0.901 and substance abuse within-family OR = 1.141.

Robustness across designs

We assessed the robustness of the key birth-order associations across seven analytical specifications spanning both between-family and within-family designs (Fig. 4). The robustness heatmap organizes diseases according to the direction and consistency of their effects, with four between-family columns (primary match, strict clinical rematch, state fixed effects and full-sibling restriction) and three within-family columns (primary, period-only and full-sibling models). The period-only within-family specification, which drops birth cohort and retains only calendar period, showed somewhat attenuated effects for several diseases, probably reflecting residual cohort confounding that is partially absorbed when the birth-year adjustment is omitted.

Among the diseases showing first-born excess, the other/unspecified PDD code group, tics/Tourette syndrome, autism, OCD, food allergy, acne, allergic rhinitis, ADHD and asthma reached Bonferroni significance across all or nearly all specifications, supporting strong robustness. Atopic dermatitis showed directional consistency but reached only nominal significance in most specifications.

Among the diseases showing second-born excess, migraine, gastritis and duodenitis, biliary tract disease, substance abuse, kidney infection and herpes zoster were consistently significant. Herpes zoster showed the strongest and most consistent second-born excess effect across all seven designs.

Cross-design concordance

Across all diseases informative in both designs, between-family and within-family ORs were positively correlated (Pearson r = 0.66, 74.2% directionally concordant; Extended Data Fig. 4a). Among the 150 Bonferroni-significant diseases in the primary between-family analysis, 127 (84.7%) showed the same direction of effect in the within-family analysis; 79 were also Bonferroni-significant in the within-family analysis. Of these 79 dual-significant diseases, 75 (94.9%) were directionally concordant. Among these 150 between-family hits, 110 (73.3%) were attenuated toward the null hypothesis in the within-family estimate. Therefore, we interpret the between-family design as a high-powered scan that can still contain residual between-family confounding, and the within-family design as the internally controlled complement that anchors interpretation when both designs agree.

The stricter clinical rematch produced highly concordant estimates with the primary between-family analysis (r = 0.93, 91% concordant, 92 Bonferroni-significant diseases overlapping with the primary analysis using the primary between-family Bonferroni threshold; Extended Data Fig. 4b and Supplementary Tables 6, 11 and 12), confirming that the primary results are not driven by residual clinical imbalance between matched first-borns and second-borns.

The state fixed-effects specification showed near-perfect agreement with the primary between-family model (r > 0.99, 98% concordant, 143 of 150 Bonferroni-significant diseases also significant; Extended Data Fig. 4c), indicating that geographical confounding has a negligible influence on the birth-order estimates.

Restricting the between-family analysis to the ‘full-sibling’ subset (families where parental-age differences are internally consistent with the sibling spacing, a heuristic for biological full siblings) produced highly concordant results (r = 0.998).

A small number of diseases showed directional discordance between designs. The most notable was obesity (between-family OR = 1.052 indicating second-born excess; within-family OR = 0.938 indicating first-born excess). Such discordances may reflect confounders that vary within families (such as differential parental feeding practices for first versus second children) or differential period effects on diagnosis.

Validation with positive and negative controls

We prespecified five positive controls and seven negative controls to validate the between-family design (Extended Data Fig. 5a). Positive controls were diseases with established birth-order associations, including allergic rhinitis, food allergy and asthma (predicted first-born excess per the sibling-exposure literature), acne (predicted first-born excess based on prior dermatological studies of sebaceous gland activity and healthcare-seeking patterns in first-borns) and substance abuse (predicted second-born excess per the behavioral literature). All five positive controls showed effects in the expected direction in both between-family and within-family designs: allergic rhinitis, food allergy, asthma and acne showed first-born excess (OR < 1), while substance abuse showed second-born excess (OR > 1). Between-family and within-family estimates were concordant in direction for all positive controls, with within-family estimates generally attenuated relative to between-family estimates (Extended Data Fig. 5a, left).

Negative controls included five diseases with primarily genetic or structural etiologies (type 1 diabetes mellitus, cystic fibrosis, Addison disease, Ehlers–Danlos syndrome, Turner syndrome) and two common acute diagnoses chosen as empirical null comparators (acute sinusitis and acute upper respiratory infection). We interpret these negative controls as specificity checks rather than uniformly powered falsification tests. The rare genetic or structural controls had limited power to exclude small birth-order effects at a Bonferroni threshold, whereas the two common acute controls provided higher-powered empirical null comparators. In the primary between-family design, the control estimates were close to the null hypothesis, with acute sinusitis (OR = 0.999) and acute upper respiratory infection (OR = 1.000) showing no meaningful between-family signal (Extended Data Fig. 5a, right). Acute sinusitis showed a small but Bonferroni-significant within-family deviation (OR = 0.965, P = 2.6 × 10⁻⁴²), so it is best viewed as a near-null specificity check rather than a globally clean negative control. Corresponding within-family estimates for the same control set are reported in Supplementary Table 5.

Sibling age spacing modulates birth-order effects

We examined whether the magnitude of birth-order effects varied with sibling age gap using a gap × birth-order interaction model, stratified into age gap categories (< 4, 4–6, 7–10, > 10 years) (Extended Data Fig. 5b).

For autism, the first-born excess was strongest at gaps of 4–6 years (stratum OR ≈ 0.60) and attenuated at very short (< 4 years) and very long (> 10 years) gaps, producing a U-shaped pattern. ADHD showed a similar pattern, with first-born excess increasing from short to medium gaps and remaining stable at longer gaps. Allergic rhinitis showed progressive attenuation of the first-born protective effect with increasing gap, which is consistent with the microbial diversity framework (closer spacing provides more microbial sharing from the older sibling). Food allergy showed pronounced first-born excess at short gaps that weakened substantially at wider spacing. Substance abuse showed a decreasing second-born excess with greater spacing, suggesting that peer-influence effects of older siblings weaken when the age difference grows. Anxiety and phobia and depression showed stable first-born excess across gap categories (Extended Data Fig. 5b). Gap × birth-order interactions were tested for all 418 diseases; 127 (30.4%) showed significant interactions at the Bonferroni threshold (P < 1.20 × 10⁻⁴). Stratified ORs for selected diseases are shown in Supplementary Table 7. We also tested targeted effect modification for seven high-priority diseases. Autism showed significant birth-order interactions with sibling gap (omnibus P = 1.9 × 10⁻⁶), paternal age (P = 1.7 × 10⁻⁷), maternal age (P = 6.4 × 10⁻¹²) and calendar period (P = 3.2 × 10⁻⁶), but not sex (P = 0.97). Allergic rhinitis and food allergy also showed gap-dependent effects; substance abuse showed heterogeneity according to gap, sex, parental age and calendar period. In a supplementary full-cohort-adjusted T-learner analysis for these same diseases, models included sibling gap, sibling sex composition, parental age, calendar period and baseline comorbidity, alongside individual age, follow-up, sex, birth-year and ICD-era covariates. The standardized second-born − first-born risk differences were consistent with first-born excess for autism, food allergy, allergic rhinitis, tics/Tourette syndrome and the other/unspecified PDD code group; with second-born excess for substance abuse; and with a near-null average contrast for ADHD, whose predicted contrasts spanned both directions (Supplementary Figs. 1 and 2 and Supplementary Table 15). These analyses are exploratory and descriptive; they identify where the observed association is strongest but do not convert the birth-order contrast into an individualized clinical prediction model.

Healthcare use sensitivity analysis

To assess whether differential healthcare contact according to birth order could inflate diagnosis rates, we computed the total number of distinct claim days per individual in the matched cohort from the diagnostic claims database. First-borns had a mean of 24.0 distinct claim days compared with 22.9 for second-borns (median 11 versus 10), a clinically modest 4.5% difference. We then reestimated all between-family models in the strict clinical cohort with log-transformed visit count as an additional covariate. Visit-adjusted birth-order ORs were highly concordant with the unadjusted strict-cohort estimates (r = 0.99), indicating that the observed birth-order associations are not driven by differential healthcare use. ORs attenuated modestly toward the null hypothesis after adjustment (for example, the autism OR shifted from 0.799 to 0.880), which is consistent with the visit count acting partly as a mediator rather than as a pure confounder (Supplementary Table 13).

Reproductive stoppage

To assess whether reproductive stoppage, or the tendency of parents to curtail childbearing after a child is diagnosed with a serious condition, could explain the first-born enrichment observed for neurodevelopmental conditions, we conducted family-level logistic regression analyses (Supplementary Table 4). The adjusted model used 10,016,101 complete-case families from the full MarketScan database.

A first-born autism diagnosis was associated with a modest reduction in the probability of having a second child (OR = 0.870, 95% CI 0.856–0.885, P = 2.6 × 10⁻⁶²), corresponding to approximately a 13% relative reduction. Tics/Tourette syndrome showed a marginal 4% reduction (OR = 0.960, P = 5.5 × 10⁻⁴). Critically, ADHD showed negligible stoppage (OR = 1.011, P = 1.7 × 10⁻⁴); the other/unspecified PDD code group, the strongest first-born excess finding (primary OR = 0.569), showed a stoppage OR of 1.056 (P = 3.3 × 10⁻⁷), indicating that parents of children with diagnoses in this code group were more likely to have a second child, the opposite direction expected under reproductive stoppage.

Sex-stratified analysis

To examine whether birth-order effects differ by sex, we reestimated all between-family models separately in males and females from the strict clinical cohort, dropping sex from the covariate set within each stratum (Supplementary Table 14). Across 310 diseases analyzable in both sexes, male and female \({\mathrm{log}}_{2}(\mathrm{OR})\) estimates were moderately correlated (r = 0.65, P = 6.6 × 10⁻³⁹), indicating broadly consistent directionality with some sex-specific modulation. For several neurodevelopmental conditions, the first-born excess was more pronounced in males (for example, ADHD: male OR = 0.880, female OR = 0.984, P_diff = 3.1 × 10⁻¹²; developmental delay: male OR = 0.799, female OR = 0.961, P_diff = 6.0 × 10⁻³), which is consistent with the known male predominance in these disorders. A similar male-predominant pattern was observed for immune-allergic and dermatological conditions: allergic rhinitis (male OR = 0.859, female OR = 0.933, P_diff = 6.8 × 10⁻¹²), asthma (male OR = 0.945, female OR = 0.996, P_diff = 2.8 × 10⁻⁴), acne (male OR = 0.801, female OR = 0.841, P_diff = 2.4 × 10⁻⁷) and ear infection (male OR = 0.958, female OR = 0.992, P_diff = 9.8 × 10⁻⁴). By contrast, second-born excess conditions, such as substance abuse (male OR = 1.158, female OR = 1.256, P_diff = 9.3 × 10⁻⁵) and migraine (male OR = 1.079, female OR = 1.162, P_diff = 5.6 × 10⁻⁴) showed stronger effects in females.

This study provides a large-scale phenome-wide assessment of birth-order effects on disease risk in US commercial claims data, leveraging two complementary epidemiological designs in over 10 million siblings. We identified 150 diseases with Bonferroni-significant birth-order associations spanning neurodevelopmental, psychiatric, immune-allergic, dermatological, gastrointestinal and cardiovascular domains.

The disease atlas view (Fig. 2 and Extended Data Fig. 1) reveals that birth-order effects are not confined to a few well-studied conditions but instead pervade nearly every clinical domain tested, with a striking concentration of first-born excess in neuropsychiatric conditions and second-born excess in digestive and musculoskeletal diseases. The domain-level summary (Extended Data Fig. 2) and partial-pooling analysis (Extended Data Fig. 3) further demonstrate that the directionality of birth-order associations is domain-specific, suggesting that different biological and social mechanisms may contribute across organ systems.

The consistency of results between designs (Extended Data Fig. 4), the performance of validation controls (Extended Data Fig. 5) and the robustness to geographical confounding collectively support the conclusion that birth order is associated with widespread, albeit modest in size, differences in disease risk.

The pattern of associations is consistent with at least three broad interpretive pathways. First, the first-born excess for allergic and atopic conditions is consistent with the sibling-exposure pattern first identified in the allergy birth-order literature⁵ and now more precisely interpreted through the old-friends and microbial diversity framework, in which exposure to diverse commensal and environmental microorganisms promotes immune regulatory development^10,11. The attenuation of the allergic rhinitis effect with wider sibling spacing (Extended Data Fig. 5b) further supports this interpretation. Second, the first-born excess for neurodevelopmental conditions, strongest for the other/unspecified PDD code group (OR = 0.57) followed by tics/Tourette syndrome, autism and OCD, is biologically and phenotypically plausible but not mechanistically resolved by claims data. A recent Pregnancy and Childhood Epigenetics meta-analysis reported birth-order-associated differences in neonatal blood DNA methylation across 16 cohorts²⁶, providing independent evidence that birth order can be associated with measurable biology at birth without establishing a mechanism for the disease associations observed in this study. In our data, parental-age confounding remains important, but autism remained first-born-enriched in the strict parental-age rematch (OR = 0.799) and in the within-family design (OR = 0.804), suggesting that parental age alone does not explain the result. The most conservative interpretation is that neurodevelopmental associations probably reflect a mixture of pregnancy-order biology, parental surveillance, diagnostic timing and residual confounding rather than one uniform mechanism. Third, the second-born excess for substance abuse aligns with sociological theories of later-born risk-taking behavior and elder sibling modeling effects³.

Our findings are concordant with, and extend, prior disease-specific studies. The strongest first-born excess was for the other/unspecified PDD code group (OR = 0.57), whose implemented ICD definition points to broader neurodevelopmental diagnostic coding rather than schizophrenia-spectrum psychosis. The autism first-born excess (OR = 0.74) is consistent with prior birth-order studies but larger in magnitude, probably reflecting the combined contribution of biological birth-order effects and residual reproductive stoppage^16,27. Importantly, reproductive stoppage cannot explain the first-born excess for most neurodevelopmental conditions. Even for autism, where stoppage is detectable, the effect is substantially smaller than the 26% first-born excess in the primary analysis (OR = 0.737); the within-family analysis, which is robust to second-born reproductive stoppage by construction, confirmed the first-born excess (within-family OR = 0.804). The allergic rhinitis (OR = 0.91) and asthma (OR = 0.97) effects are consistent with the established sibling-exposure literature on allergic disease^6,8. The substance abuse second-born excess (OR = 1.19) is consistent with this broader later-born behavioral framework. Our study adds hundreds of previously unexamined diseases, including strong associations for acne (OR = 0.87), adjustment disorder (OR = 0.86), biliary tract disease (OR = 1.18) and herpes zoster (OR = 1.35).

Several limitations merit discussion. First, claims data capture diagnoses that lead to billable healthcare encounters, not true disease incidence. Healthcare-seeking behavior may also vary with birth order. For example, if parents bring first-borns to the doctor more readily, this could inflate first-born diagnosis rates independently of true disease risk. The persistence of effects in the within-family design partially mitigates this concern (within-family comparisons control for family-level healthcare-seeking tendencies), but within-family differences in birth-order-dependent parental attention could remain. Arguing against blanket ascertainment bias, among all 418 tested diseases, 230 (55%) showed point estimates in the direction of second-born excess; among the 150 Bonferroni-significant associations, the split between first-born and second-born excess was near-symmetric (79 versus 71; binomial P = 0.568). If first-born diagnoses were systematically inflated by greater parental attention, one would expect a strong skew toward first-born excess; the observed near-parity is inconsistent with this explanation. Moreover, 71 Bonferroni-significant diseases showed second-born excess, including clinically acute conditions (herpes zoster, kidney infection, acute renal failure, biliary tract disease) whose presentation is driven by objective pathology rather than differential parental surveillance.

Second, the APC problem is inherent in any birth-order study. Within sibling pairs, the older child was born earlier, is observed at different ages during any calendar period and experienced different diagnostic standards. We addressed this through cohort and period adjustment in regression, multiple-model specifications in the within-family analysis, and negative controls; however, residual APC confounding cannot be fully excluded.

Third, residual confounding according to parental age persists despite matching and regression adjustment. The structural correlation between birth order and parental age at birth (first-borns have younger parents by definition) means that parental-age effects cannot be fully disentangled from birth-order effects without strong parametric assumptions. For autism, where parental-age effects are strongest, tightening parental-age matching from a ± 10-year caliper to a ± 1-year caliper attenuated the between-family OR from 0.737 to 0.799. The within-family estimate (OR = 0.804) was similar to this tighter between-family estimate; however, within-family comparisons do not fully eliminate parental-age concerns because parental age at birth differs between first-born and second-born siblings. Thus, these analyses reduce but do not fully resolve parental-age confounding.

Fourth, the stricter clinically matched sensitivity analysis should be interpreted as a robustness analysis rather than a primary causal estimate. By matching on early baseline comorbidity and medication burden, it reduces residual clinical imbalance between first-born and second-born children from different families, but it may also partially control for early manifestations that lie on the pathway between birth order and later diagnosis. Reassuringly, this stricter rematch mainly pruned the original between-family signal rather than reversing it (Extended Data Fig. 4b).

Fifth, MarketScan captures employer-insured individuals, who are predominantly working-age, higher-income and disproportionately White. Our findings may not generalize to uninsured, Medicaid-covered or non-European ancestry populations.

Finally, the restriction to two-child families introduces selection. This design choice ensures clean identification of first-born versus second-born status and eliminates confounding by completed family size, but families with two children may differ systematically from larger families (for example, in socioeconomic resources or reproductive preferences). Future work should evaluate whether the same phenome-wide patterns replicate in larger sibships and in independent cohorts (for example, Nordic registry data or Medicaid claims). However, we note that internal triangulation in two distinct analytical designs, the concordance of our estimates with published disease-specific studies and the validation of positive and negative controls collectively argue against a purely spurious pattern.

These results have several implications. Although individual effect sizes are modest (median OR = 0.89 for first-born excess, 1.10 for second-born excess), birth order is a universal exposure. Modest relative risks applied to the entire pediatric population can be translated into nontrivial population-level burden; the consistency of these effects between independent designs argues against noise. For clinicians, they highlight birth order as a modestly informative risk marker across multiple disease domains, which may be relevant for family counseling and screening prioritization. For researchers, the phenome-wide catalog of birth-order effects provides a resource for hypothesis generation and for benchmarking future studies. For epidemiologists, the dual-design approach demonstrated in this study offers a template for studying other nonrandomizable familial exposures.

In conclusion, birth order is associated with disease risk more broadly than previously appreciated. The convergence of evidence from between-family and within-family designs, combined with validation controls and robustness analyses, is consistent with a mixture of physiological, immunological and social mechanisms operating throughout the human disease phenome.

This study used de-identified secondary administrative claims data from the Merative MarketScan research databases, which are statistically de-identified to meet Health Insurance Portability and Accountability Act privacy requirements. The analyses used existing de-identified records, involved no direct contact with human participants and were conducted without access to direct identifiers. The University of Chicago Institutional Review Board determined this study to be exempt from human participant review. Informed consent was not applicable.

Data source

We used Merative MarketScan Commercial Claims and Encounters data from the 2003–2024 annual releases, which capture inpatient, outpatient and pharmacy claims for approximately 200 million unique covered lives in employer-sponsored health insurance plans across the United States. The raw yearly releases are distributed as SAS7BDAT files. From these annual releases, we constructed extracted analysis databases containing patient demographics (sex, birth year, enrollment family identifier), enrollment histories (start and end dates for coverage intervals) and diagnostic codes (ICD-9-CM and ICD-10-CM) recorded at each encounter²⁸. We used the MarketScan enrollee identifier as the individual identifier and the enrollment family identifier as the family identifier. In the submitted analytical cohort, all 10,270,012 sibling rows had distinct individual identifiers, all 5,135,006 families contained exactly two analytical siblings and no cohort individual identifiers mapped to more than one family identifier in the extracted demographics table. These checks were used to verify internal identifier consistency in the extracted data, but they cannot rule out unobserved reenrollment under a new identifier after employer or insurance-plan changes without a vendor crosswalk.

All Merative MarketScan data were accessed and analyzed under an institutional license agreement with Merative, and were used in compliance with the terms of use of that license agreement.

Cohort definition

We identified enrollment families in extracted enrollment databases derived from the annual Merative MarketScan releases, containing at least two members born between 1978 and 2013, aged 12–60 years in the current data year. Within each family, we inferred parents as the youngest male and female members who were aged 15 years or older than the oldest candidate child and whose age at the youngest child’s birth fell within 18–69 years. At least one inferred parent was required (eliminating spouse–pair misclassification).

After excluding inferred parents, we required exactly two remaining children (a ‘true two-child family’ restriction applied before individual eligibility filters to prevent families with a third ineligible child from being misclassified as two-child families). Both children were required to have 365 or more days of enrollment visibility, age at last observation of 12 years or older and valid first and last observation years. The older child was designated sib_order = 1 (first-born) and the younger sib_order = 2 (second-born). The sibling age gap was computed as the absolute difference in birth years. Parental psychiatric history was ascertained by searching all diagnostic codes for the relevant parent against a curated set of 302 ICD-9-CM and 252 ICD-10-CM psychiatric diagnostic codes spanning schizophrenia-spectrum disorders (ICD-9: 295.xx; ICD-10: F20–F29), mood and bipolar disorders (296.xx; F30–F39), anxiety, stress and somatoform disorders (300.xx, 308–309.xx; F40–F48), personality disorders (301.xx; F60–F69) and sleep disorders (327.xx; G47.xx), among others. The complete code list is provided in Supplementary Table 10. A parent was flagged as having a psychiatric history if any qualifying code appeared in their claims record during the enrollment window.

Geographical annotation

Each individual was assigned a three-digit ZIP code (ZIP3) in our extracted demographics database derived from the annual Merative MarketScan files. ZIP3 was mapped to a dominant county (FIPS code) using the HUD ZIP–TRACT crosswalk, weighted according to residential ratio²⁹. County population estimates from the U.S. Census Bureau Vintage 2023 county file (co-est2023-alldata.csv) were then used to assign both census region and county population^30,31. Specifically, we used the census REGION code (1 = Northeast, 2 = Midwest, 3 = South, 4 = West) and recoded it as our analysis variable direction (E, N, S, W), where E corresponds to the census Northeast and N corresponds to the census Midwest. County population was also used to define the urbanization tertile (low, medium, high).

Between-family matching

We constructed a matched cohort by pairing one first-born from family A with one second-born from family B, subject to exact match on sex, birth year, census direction and urbanization tertile and caliper match on follow-up duration (± 50 days), paternal age at birth (± 10 years), maternal age at birth (± 10 years) and sibling age gap (± 2 years). We also had the constraint that the two individuals came from different families.

Within each exact-match stratum, we applied a greedy nearest-neighbor algorithm with randomized order and randomized tie-breaking^32,33. The distance metric was a weighted sum of 2.0 × ∣Δdays_visible∣ + 1.0 × ∣Δfather_age∣ + 1.0 × ∣Δmother_age∣ + 0.5 × ∣Δgap∣. Matching used 45 parallel workers with a fixed random seed (2025) for reproducibility. Post-match SMDs for paternal (0.30) and maternal (0.36) age exceeded the conventional 0.1 balance threshold because, by construction, second-borns are drawn from families with structurally older parents at their birth; therefore, parental ages are included as categorical regression covariates rather than relying on matching alone to control for these differences.

Post-match quality control confirmed zero caliper violations, correct sib_order composition for all 1,616,881 pairs and improved SMDs for all covariates³⁴ (Supplementary Table 2). This procedure was standard greedy nearest-neighbor caliper matching; we did not derive a new matching estimator or doubly robust estimator. Because matching was performed at the individual level, the same family could contribute to more than one matched pair. A dependence audit found reciprocal unordered family-pair matches to be rare (ten of 1,616,881 primary matched pairs and 24 of 529,760 strict matched pairs); covariance estimator sensitivities are reported in the supplementary source data.

Stricter clinical matching sensitivity analysis

As an additional robustness analysis, we rebuilt the between-family matched cohort using substantially tighter parental-age and baseline clinical matching. To preserve overlap, we did not require exact state matching in this sensitivity analysis; instead, geography remained controlled by exact census direction and urbanization tertile in the match and by a separate state fixed-effects analysis in the regression. Pairs were required to come from different families and opposite birth order; they were matched exactly on sex, birth year, census direction, urbanization tertile, age at start of observation, sibling age gap and baseline Charlson comorbidity score. We imposed a caliper of ± 50 days on enrollment visibility, ± 1 year on paternal and maternal age at birth, and imposed a ± 0.1 caliper on the baseline Medication-Based Disease Burden Index, a pharmacy-claim-based comorbidity measure³⁵. The Charlson comorbidity index was computed from diagnostic codes using the Quan adaptation³⁶. The baseline clinical indices were derived from the first 365 observed enrollment days before outcome ascertainment. This stricter rematch yielded 529,760 pairs (1,059,520 individuals); we refitted the primary between-family logistic models in this cohort to quantify concordance with the primary between-family estimates.

Disease phenotyping

We defined 569 diseases using ICD-9-CM and ICD-10-CM diagnostic code groupings adapted from an established phenotyping system previously applied to longitudinal claim data³⁷. The complete mapping from disease names to ICD codes is provided in Supplementary Table 9. A disease was considered present for an individual if any qualifying diagnostic code appeared in their claims record during the entire observation window. We imposed a minimum of 500 cases in the matched cohort (between-family analysis) and 100 discordant sibling pairs (within-family analysis) for a disease to be included.

Clinical domain classification

Each disease was assigned to one of 17 clinical domains based on primary organ system or clinical category: neuropsychiatric/behavioral, neurological, infectious, musculoskeletal, sense organs, dermatological, respiratory, endocrine/metabolic, congenital, pregnancy-related, digestive, circulatory, genitourinary, general symptoms, neoplasms, injury/toxicology and hematological/immune. Domain assignments were made by the study team based on established clinical groupings of ICD codes and were fixed before the analysis. Diseases that could plausibly belong to multiple domains were assigned to the most clinically conventional category (for example, migraine to neurological rather than general symptoms).

For the revised domain-level analysis, we fitted a normal-normal empirical-Bayes partial-pooling meta-analysis to disease-level log(ORs) and CI-derived standard errors, separately for the primary between-family, within-family and strict between-family analyses. Parameters were estimated using maximum likelihood with scipy.optimize. The model estimates a global log-OR mean, between-domain heterogeneity, within-domain residual heterogeneity and domain-specific posterior means. This analysis was used to quantify domain-level clustering and shrinkage; it did not replace the prespecified disease-level Bonferroni inference.

Between-family logistic regression

For computational efficiency, individual-level data were collapsed into frequency tables (one row per unique covariate pattern within each disease), with case counts as weights. For each disease, we fitted a weighted logistic regression model with disease status (0/1) as the outcome and birth order (sib_order: 1 = first-born, 2 = second-born) as the primary exposure. Covariates included sibling age gap (categorical: < 4, 4–6, 7–10, > 10 years), age at last observation (categorical: ≤ 6, 7–10, 11–13, 14–16, 17–18, > 18 years), sex, paternal and maternal age at birth (categorical, 5-year bins), parental psychiatric history (father, mother), urbanization group, ICD coding era (ICD-9 versus ICD-10 based on first observation year), calendar period (5-year bins from observation mid-year) and a full-sibling consistency flag. Standard errors were clustered at the county level using HC1 robust variance³². Follow-up duration entered the model as a log-transformed covariate. We used a logistic model because each phenotype was analyzed as an ever-versus-never diagnosis indicator during the observed enrollment window, not as a recurrent event count. Therefore, follow-up duration was included as an adjustment covariate rather than as a person-time offset. As model-form sensitivity analyses, we refitted the 418 primary between-family disease models using log-link Poisson regression with robust standard errors and complementary-log–log regression on the same aggregated covariate tables; both alternative links were direction-concordant with the submitted logistic estimate for all 418 diseases. These link-function sensitivities used the same prespecified aggregate covariate tables; therefore, they do not assess continuous or spline parameterizations of binned covariates. County-level clustering was chosen to allow for local correlation in coding practice, provider availability and claim ascertainment after ZIP3-to-county geographical annotation. For seven high-salience diseases, inference was also compared across model-based, HC1, county-clustered, state-clustered, matched-pair-clustered and family-clustered covariance estimators.

For sensitivity, we reestimated each disease model with state fixed effects (50 state indicators plus District of Columbia).

Within-family conditional logistic regression

For each disease, we fitted a conditional logistic regression (Cox proportional-hazards model with case–control sampling, implemented via clogit in R) stratified on family identifier, with birth order as the exposure. Covariates included birth-cohort year (centered), age at last observation (categorical: ≤ 6, 7–10, 11–13, 14–16, 17–18, 19–22, 23–30, > 30 years), sex, log(follow-up duration), ICD coding era and calendar period of observation (categorical 5-year bins based on observation-window midpoint).

We fitted four model specifications to assess sensitivity to APC parametrization: (1) cohort-adjusted (birth year + age + ICD era); (2) cohort plus calendar period (primary); (3) period-only (dropping birth year); and (4) cohort with gap × birth-order interaction. The cohort-plus-period specification served as the primary within-family model. The analyses were repeated in the full-sibling subset as a sensitivity analysis.

Gap × birth-order interaction

In the between-family design, we estimated gap-stratified birth-order ORs from a logistic regression that included a birth-order × gap-category interaction (< 4, 4–6, 7–10, > 10 years). In the within-family design, we included a gap × birth-order interaction term. As additional exploratory heterogeneity analyses, we first fitted targeted primary between-family interaction models for selected high-priority diseases using sibling age gap, sex, paternal age, maternal age and calendar period. We report stratum-specific birth-order ORs and omnibus likelihood-ratio tests for the interaction terms. We then fitted a full-cohort-adjusted T-learner for the same seven diseases. For each disease, separate gradient-boosted logistic outcome models were trained among first-born and second-born children and included individual covariates (age at last observation, log(follow-up), sex, birth year, calendar period and ICD era) together with family-level modifiers (sibling gap, sibling sex composition, paternal and maternal age, baseline Charlson and Medication-Based Disease Burden Index comorbidity summaries, urbanization, census region and direction, full-sibling-like flag and family observation period). For each two-child family, we predicted the second-born − first-born risk difference over both observed sibling covariate profiles and averaged the two contrasts, yielding a standardized family-level risk-difference contrast. We used a T-learner rather than causal forests or Bayesian additive regression trees because it retains an explicit first-born versus second-born outcome-model contrast while allowing flexible modifier discovery. Because birth order is fixed by family structure rather than assigned through a covariate-dependent propensity mechanism, we treated this analysis as exploratory modifier discovery rather than as individualized causal risk prediction. We did not use family size as a modifier because the primary cohort is restricted to exactly two-child families.

Reproductive stoppage analysis

We conducted family-level analyses using an extracted MarketScan demographics database (not restricted to two-child families). Among all families with at least one child meeting the basic eligibility criteria, we identified whether the first-born had a diagnosis for each of four key first-born excess conditions (autism, ADHD, tics/Tourette syndrome and the other/unspecified PDD code group) and whether the family had a second eligible child. For each disease, we fitted a logistic regression with ‘has second child’ (0/1) as the outcome and ‘first-born diagnosis’ as the exposure, adjusting for first-born sex, birth year, enrollment duration, parental ages and calendar period, with HC1 robust standard errors.

Healthcare use sensitivity analysis

To test whether differential healthcare contact according to birth order could confound disease ascertainment, we computed the total number of distinct claim days (unique service dates with any diagnostic code) for each individual in the stricter clinically matched cohort from the diagnostic claims database. We then reestimated all between-family logistic regression models in this cohort with the log-transformed visit count (\(\mathrm{log}({n}_{\mathrm{visits}}+1)\)) included as an additional covariate alongside all covariates in the primary model. We assessed concordance between unadjusted and visit-adjusted birth-order ORs using Pearson correlation of \({\mathrm{log}}_{2}(\mathrm{OR})\) estimates.

Sex-stratified analysis

To assess whether birth-order effects are modulated by sex, we reestimated all between-family logistic regression models separately in males and females from the strict clinical cohort, omitting sex from the covariate set as it is constant within each stratum. A minimum of 200 cases per stratum was required for model convergence. We assessed male–female concordance using Pearson correlation of \({\mathrm{log}}_{2}(\mathrm{OR})\) estimates across diseases analyzable in both sexes.

Cross-design concordance analysis

To quantify the agreement between alternative analytical designs, we computed three metrics for each pairwise comparison (for example, primary between-family versus within-family): (1) Pearson correlation of \({\mathrm{log}}_{2}(\mathrm{OR})\) estimates across all diseases analyzable in both designs; (2) directional concordance, defined as the proportion of diseases for which both designs yield an OR < 1 or both yielded an OR > 1; and (3) overlap of Bonferroni-significant diseases. We generated scatter plots of \({\mathrm{log}}_{2}(\mathrm{OR})\) from the primary between-family analysis against each alternative specification (Extended Data Fig. 4).

Positive and negative controls

We prespecified five positive controls (diseases expected to show birth-order effects) and seven negative controls (diseases expected to show no effect) before examining phenome-wide results. Positive controls were: allergic rhinitis, food allergy and asthma (predicted first-born excess per the sibling-exposure literature^5,6); acne (predicted first-born excess based on prior dermatological literature); and substance abuse (predicted second-born excess per the sociological literature³). Negative controls included five conditions with primarily genetic or structural etiologies unlikely to be influenced by birth order (type 1 diabetes, cystic fibrosis, Addison disease, Ehlers–Danlos syndrome, Turner syndrome) and two common acute conditions serving as empirical null comparators (acute sinusitis and acute upper respiratory infection).

Multiple-testing correction

We applied Bonferroni correction (α = 0.05/n_tests, where n_tests = 418 for between-family and 541 for within-family) and the Benjamini–Hochberg procedure³⁸ at a false discovery rate of q < 0.05.

Software

Cohort extraction and data preparation were performed in Python v.3.11 (pandas, sqlite3, multiprocessing). Statistical analyses were performed in R v.4.3 (survival, sandwich, lmtest packages). Figures were generated in Python using matplotlib v.3.8.

Reporting summary

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

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.