This study was approved by the institutional review board of the Veterans Affairs St. Louis Health Care System, which granted a waiver of informed consent (protocol no. 1606333). All participants who were eligible for this study were enrolled; no a priori sample size analyses were conducted to guide enrollment. All analyses were observational, and investigators were aware of participant exposure and outcome status.

Participants were selected from the US Veterans Health Administration (VHA) electronic health database. The VHA delivers healthcare to discharged Veterans of the US armed forces in a network of nationally integrated healthcare systems including more than 1,415 healthcare facilities. Veterans enrolled for care in the VHA have access to extensive medical benefits, such as inpatient and outpatient services, preventative, primary and specialty care, mental health services, geriatric care, long-term and home healthcare, medications and medical equation and prosthetics. The VHA electronic health database is updated daily.

A flowchart of cohort construction is provided in Extended Data Fig. 1. We first identified users of the VHA with at least one positive SARS-CoV-2 test between 1 March 2020 and 6 April 2022 (n = 519,767), enrolling these participants at the date of first positive test (set as T0). Use of the VHA was defined as having record of use of outpatient or inpatient service, receipt of medication or use of laboratory service with the VHA healthcare system in the 2 years before enrollment. We selected those still alive 90 d after their first positive SARS-CoV-2 test (n = 489,779). We then further selected participants who experienced reinfection, defined as a positive SARS-CoV-2 test 90 d or more after the first infection, where reinfection could occur between 1 June 2020 and 25 June 2022, which spans the time frame in the US in which pre-Delta, Delta and Omicron variants predominated16,17,18,19. The 90-d minimum time frame to define reinfection was specified to minimize inclusion of repeat positive tests that may be related to the first infection16,17,18,19. There were 40,947 participants who had a reinfection, where the time of reinfection was set as T1. To ensure a similar distribution of follow-up time in the no reinfection and reinfection groups, participants in the no reinfection group were randomly assigned a T1 based on the distribution of T1 of those in the reinfection group who shared the same calendar month as the date of first infection, resulting in a group of 443,588 participants with no reinfection that were alive at their assigned T1.

We then constructed a noninfected control group. We first identified 5,760,792 VHA users between 1 March 2020 and 6 April 2022 with no record of a positive SARS-CoV-2 test. We then randomly assigned a T0 to each participant in the group on the basis of the distribution of the T0 dates in those with at least one positive SARS-CoV-2 test, selecting the 5,458,815 who were alive at their assigned T0. We selected those who were alive 90 d after their T0 (n = 5,408,880). After randomly assigning a T1, there were 5,334,729 in the noninfected control cohort. All cohort participants were followed until 25 June 2022.

Participant data were obtained from the VHA Corporate Data Warehouse. The patient and vital status domains provided data on demographic characteristics. VHA mortality information contains both inhospital and nonhospital deaths collected from the Veterans Affairs and non-Veterans Affairs sources including the VHA’s Beneficiary Identification Record Locator System and medical inpatient datasets, as well as Medicare Vital Status File, Social Security Administration’s Master File and information from death certificates and the National Cemetery Administration20. The outpatient and inpatient encounter domains provided information on health characteristics including details on date and place of encounter with the healthcare system and diagnostic and procedural information. The Pharmacy and Bar Code Medication Administration domains provided medication records, while the laboratory results domain provided laboratory test results for tests conducted in both inpatient and outpatient settings7,21. Information about SARS-CoV-2 tests and vaccinations were obtained from the COVID-19 Shared Data Resource. Positive SARS-CoV-2 tests consisted of results from PCR or antigen tests conducted in the Veterans Affairs or reported to the Veterans Affairs. The 2019 Area Deprivation Index at the residential address of each cohort participant was used as a contextual measure of socioeconomic disadvantage22. Information from the US Center for Disease Control and Prevention provided the proportion of SARS-CoV-2 variant by week in each Health and Human Services (HHS) region.

Outcomes were prespecified on the basis of previous evidence1,2,3,4,5,6,7,8,21,23,24,25,26,27,28,29. Outcomes included all-cause mortality, hospitalization, having at least one sequela and organ system disorders including cardiovascular disorders, coagulation and hematological disorders, diabetes, fatigue, gastrointestinal disorders, kidney disorders, mental health disorders, musculoskeletal disorders, neurological disorders and pulmonary disorders. Organ system disorders were defined as a composite outcome of a set of prespecified individual sequelae in that system at the date of first incident sequela in that system during follow-up. Organ system disorders were defined on the basis of inpatient or outpatient diagnostic codes, medication prescriptions or laboratory values. A list of the individual sequelae by organ system are provided in Supplementary Table 13. The outcome of ‘at least one sequela’ was defined at the time of occurrence of first incident sequela among all individual sequelae. For a participant, for a given outcome, each individual sequela was included in the assessed outcome only when there was no record of that health condition in the 2 years before T0. Participants were excluded from the analysis of an outcome if they had previous history of all the individual sequelae that contributed to the outcome being examined. Hospitalization was defined as first inpatient admittance during follow-up. In analyses of kidney disorders, participants with a previous history of end-stage kidney disease were excluded and follow-up was censored at the time of end-stage kidney disease (Supplementary Table 13).

Covariates included a set of variables that were predefined based on previous knowledge4,5,6,7,21,23,25,26,27,30,31,32,33 and a set of variables that were selected algorithmically. Predefined covariates included demographic information (age, race and sex), contextual information (Area Deprivation Index) and measures of healthcare use in the 2 years before T0, which included the number of outpatient visits, inpatient visits, unique medication prescriptions, routine laboratory blood panels and use of Medicare services, as well as a previous history of receiving an influenza vaccination. Smoking status was also included as a covariate. Characteristics of the participants’ health history included record of anxiety, cancer, cardiovascular disease, cerebrovascular disease, chronic kidney disease, chronic obstructive pulmonary disease, dementia, depression, type 2 diabetes mellitus, estimated glomerular filtration rate, immunocompromised status, peripheral artery disease, systolic and diastolic blood pressure and body mass index on the basis of inpatient or outpatient diagnostic codes, medication prescriptions, laboratory values and vital signs. Immunocompromised status was defined according to the US Center for Disease Control and Prevention definitions by a history of organ transplantation, advanced kidney disease (an estimated glomerular filtration rate <15 ml min −11.73 m−2 or end-stage renal disease), cancer, HIV or conditions with prescriptions of more than 30-d use of corticosteroids or immunosuppressants including systemic lupus erythematosus and rheumatoid arthritis.

We also included a set of covariates related to the acute phase of the first infection: severity of the acute phase of the disease, defined in mutually exclusive groups of nonhospitalized, hospitalized and admitted to the intensive care unit during the acute phase and whether the participant received SARS-CoV-2 treatment of antivirals, antibodies and immunomodulators including corticosteroids, interleukin-6 inhibitors and kinase inhibitors. We also included—as measures of spatiotemporal differences—the calendar week of enrollment and geographical region of receipt of care defined by the Veterans Integrated Services Networks. We also adjusted for vaccination status, which was defined as receiving no, one, two and three or more Janssen (Ad26.COV2.S; Johnson & Johnson), Pfizer-BioNTech (BNT162b2) or Moderna (mRNA-1273) vaccination shots. In consideration of the dynamicity of the pandemic, additional covariates included hospital system capacity (the total number of inpatient hospital beds), inpatient bed occupancy rates (the percentage of hospital beds that were occupied) and a measure of the proportions of SARS-CoV-2 variants by HHS region33. These measures were ascertained for each participant in the week of cohort enrollment at the location of the healthcare system they received care at.

In addition to the predefined covariates, we leveraged the high dimensionality of Veterans Affairs electronic health records by employing a high-dimensional variable selection algorithm to identify additional covariates that may potentially confound the examined associations34. We used the diagnostic classifications system from the Clinical Classifications Software Refined v.2021.1, available from the Healthcare Cost and Utilization Project sponsored by the Agency for Healthcare Research and Quality, to classify more than 70,000 International Classification of Diseases, 10th revision diagnosis codes in the 2 years before T0 for each participant into 540 diagnostic categories35,36,37. Using the Veterans Affairs national drug classification system, we also classified 3,425 different medications into 543 medication classes38,39. Finally, on the basis of Logical Observation Identifiers Names and Codes, we classified laboratory results from 38 different laboratory measurements into 62 laboratory test abnormalities, defined by being above or below the corresponding reference ranges. Of the high-dimensional variables that occurred at least 100 times in participants in each group, we selected the top 100 variables with the highest relative risk for differences in group membership in first infection or reinfection.

Mean (s.d.) and frequency (percentage) of characteristics are reported for those with no SARS-CoV-2 reinfection, SARS-CoV-2 reinfection and the noninfected control group, where appropriate. We provide information on the distribution of frequency of reinfections, time between infections and variant of reinfection (defined by predominant variant given the calendar week and HHS region of the residential location of cohort participants when reinfection occurred).

All associations were estimated based on weighting approaches combined with survival analyses. We conducted a primary analysis to evaluate the risk and burden of reinfection compared to no reinfection (Supplementary Fig. 2). Logistic regressions were constructed to estimate the propensity score of group membership; regressions included predefined covariates, high-dimensional covariates and time from T0 to T1 as a means to adjust for residual differences in duration of follow-up. A reference cohort of the overall infected cohort at T0 was used as the target population. Inverse probability weighting was then used to balance the covariates. A weighted Cox survival model was then used to estimate the average risk and event rate difference between those with a reinfection and those with no reinfection. Standard errors were estimated by applying the robust sandwich variance estimator method. Covariate balance among all predefined and high-dimensional variables were assessed through the standardized mean difference, where a difference <0.1 was taken as evidence of balance. We estimated the incidence rate difference (referred to as excess burden) between groups per 1,000 participants at 6 months after the start of the follow-up based on the difference in survival probability between the groups. These analyses were repeated by subgroup on the basis of the number of vaccination shots received (0, 1 or 2+) before reinfection using an overlap weighting approach. To test whether the risk on the multiplicative scale differed between participants with different duration between T0 and T1, a model with a linear interaction term between reinfection status and duration was constructed and the corresponding P value is reported for the outcomes of all-cause mortality, hospitalization and having at least one sequela.

To examine whether risks associated with a reinfection were present in the acute and postacute phases of reinfection, we conducted analyses to examine risks in 30-d time intervals starting at the time of reinfection up to 180 d after reinfection. HRs and 30-d burdens were estimated independently for each 30-d time interval. During each 30-d interval, outcomes were defined at the time of first occurrence within this interval in those who did not have that outcome in the 2 years before the first infection.

We then used a doubly robust approach to examine the risk and cumulative burden per 1,000 persons at one-year after first infection of sequelae associated with one, two and three or more infections versus a noninfected control (Supplementary Fig. 3). A third or more infection was defined as a positive test at least 90 d after the second infection. The number of infections and outcomes were assessed in the 360 d after T0 + 30 d. Since those with three or more infections predominantly had their first infection before the Omicron variant was present, to enhance comparison across groups, we restricted this analysis to those with a T0 period before Omicron became the predominant variant in at least one HHS region (11 December 2022). Because participants with three or more infections must have not died during the follow-up period to have that third (or more) infection, we did not examine the outcome of all-cause mortality due to immortal time bias.

We examined, as positive outcome controls, the risk of fatigue in those with a SARS-CoV-2 infection compared to the noninfected control group as a means of testing whether our approach would reproduce established knowledge4,5,25,26,27.

The application of a negative outcome control may help detect both suspected and unsuspected sources of spurious biases. Therefore, we examined the difference in risks of atopic dermatitis and neoplasms between those who had reinfection and the first infection, where no previous knowledge suggested that an association should be expected. The testing of positive and negative outcome controls may lessen, although not eliminate, concerns about biases related to study design, covariate selection, analytical approach, outcome ascertainment, unmeasured confounding and other potential sources of latent biases40,41.

All analyses were two-sided. In all analyses, a 95% CI that excluded unity was considered evidence of statistical significance. All analyses were conducted in SAS Enterprise Guide v.8.2, and all figures were generated in R v.4.0.4.

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