The study was approved by the Swedish Ethical Review Authority (reference number 2020-06540). The requirement for informed consent was waived because the study is register-based, and all individual data were de-identified.

The Swedish nationwide health and administrative registers were linked via the unique personal identification number [25], with data available from inception dates of the registers to 31 December 2020 (Note S1). To ensure that all study participants had a similar opportunity to utilize Swedish healthcare services under the study period, the study population was identified from the Prescribed Drug Register (PDR) [26] and consisted of individuals who (1) had at least one dispensation record of any prescribed medication in 2007–2019 (n = 10,703,848) and (2) were BZDR naïve by 1 January 2007 (i.e., had no dispensation records of any BZDRs during a 1.5-year washout period from the PDR inception on 1 July 2005 and throughout 2006) and did not receive their first BZDR dispensation after 31 December 2019 (to allow for a 12-month follow-up before the study end on 31 December 2020) (n = 9,754,184). The study population was restricted to individuals without a lifetime diagnosis of epilepsy because of a complex relationship between epilepsy (one of indications for BZDs) and infections, particularly central nervous system infections [27]. Also, to enable the distinction between singletons and multiple births, we further restricted the study population to individuals born in Sweden with available linkages to biological mothers. Following these restrictions, 7,362,979 individuals remained available for constructing the study cohorts.

We compared incident BZDR recipients and individuals unexposed to BZDR using the triangulating design approach [23, 24]: (1) a demographically matched cohort, (2) a co-twin control cohort, and (3) an active comparator cohort. The demographically matched and active comparator cohorts included only singletons.

For all three cohorts, incident BZDR recipients were defined as those who had their first dispensation of any BZDR in 2007–2019, while the unexposed individuals were selected separately for each cohort, as described below. BZDR dispensation records were collected from the PDR using the Anatomical Therapeutic Chemical (ATC) codes for benzodiazepine derivatives in anxiolytics (N05BA), hypnotics/sedatives (N05CD), antiepileptics (N03AE), and Z-drugs (N05CF) (Table S1).

The construction of the three cohorts, together with the number of included and excluded individuals, is visualized in Fig. 1 and reported in Table S2. First, in the demographically matched cohort, incident BZDR recipients were 1:1 randomly matched to individuals from the study population who had no BZDR dispensations before or on the index date (i.e., when the exposed individuals were dispensed their first BZDR). Matching was performed by year and month of birth, sex, and county of residence in Sweden at the index date’s year. Second, to control for unmeasured familial confounding [28], we constructed a co-twin cohort of discordantly exposed twins or other multiple birth individuals, e.g., triplets by means of the Multi-Generation Register [29]. Multiple birth individuals were considered exposed if they were either the only or the first in the family to dispense their first BZDR prescription. Third, we created an active comparator cohort [30] where selective serotonin reuptake inhibitors (SSRIs) were chosen as a comparator, in line with prior literature [31, 32]. This approach helped reduce the potential for confounding by indication due to a partial overlap in psychiatric indications between BZDRs and SSRIs (e.g., anxiety disorders). Other measured covariates were balanced between BZDR recipients and SSRI recipients via the inverse probability of treatment weighting (IPTW) [33] method. Incident SSRI recipients were selected using the PDR records with N06AB ATC-codes (Table S1) if any SSRIs were dispensed for the first time in 2007–2019. Individuals with simultaneous first BZDR and first SSRI dispensations were excluded. Also, incident BZDR recipients were restricted to those with no SSRIs dispensed during the washout period (while incident SSRI recipients by default had no BZDR dispensations in the washout period, according to our study population definition). The index dates for BZDR recipients and SSRI recipients were defined separately as the dates of the first dispensation of the corresponding drugs.

Study population and the selection of the exposed individuals (BZDR recipients) and controls for each design. Detailed information on inclusion and exclusion criteria with the numbers of individuals included/excluded at the stage when each criterion was applied is reported in Table S2 (of the EMS). *Washout period refers to 1.5-year time period between 1 July 2005 (i.e., inception of the Prescribed Drug Register) and 31 December 2006. †At least one co-sibling from a multiple birth family did not have BZDR dispensation before and on the date when their exposed co-sibling(s) were dispensed the first BZDR (i.e., at least one co-sibling is unexposed and at least one is exposed). ‡Incident BZDR-recipients without precedent incident dispensation of SSRI (i.e., “BZDR first”). §Incident SSRI-recipients without precedent incident dispensation of BZDR (i.e., “SSRI first”). **The additional exclusion criteria are the same for all three cohorts and refer to (i) age ≥ 65 years at the index date or death before the index date, (ii) in-hospital stay for > 90 days between 1 July 2005 and the index date, (iii) having emigrated/re-immigrated between 2005 and the index date, and (iv) having records of any infections in the National Patient Register during 2 years prior to the index date. ††Refers to the incident BZDR recipients who were dispensed any SSRI for the first time during a 1.5-year washout period or after 31 December 2019 (i.e., in 2020). ‡‡If the same person was simultaneously (co-incidentally) dispensed the first BZDR and the first SSRIs at the same date. §§ BZDR recipients and SSRI recipients were included only if all information on the covariates used for the inverse probability of treatment weighting (IPTW) were available (i.e., a “complete case scenario”). BZDR benzodiazepines and related Z-drugs; IPTW inverse probability of treatment weighting; PDR Prescribed Drug Register; SSRI selective serotonin reuptake inhibitors

In all three cohorts, BZDR recipients and their counterparts were considered to be eligible for analyses if they were (1) alive and aged below 65 years at the index date, (2) residing in Sweden (i.e., not emigrated/re-immigrated) between 2005 and the index date, (3) not hospitalized for more than 90 days between the start of washout period and the index date (to further clarify the exposure status since the PDR does not cover in-hospital medication use), and (4) free from any serious infections during 2 years before the index date, according to the National Patient Register. We focused on individuals aged below 65 years to reduce age-related differences in pharmacokinetics and pharmacodynamics of BZDRs (e.g., in older patients, BZDRs have prolonged elimination half‐lives, which could lead to accumulation of the drugs) [34, 35]. As a result, the final analytical sample of (1) the demographically matched cohort included 713,896 pairs of BZDR recipients and matched unexposed individuals, for whom information on matching variables was available; (2) the co-twin cohort comprised 9197 BZDR recipients and their 9298 unexposed co-twins (of which 98.8% were twin pairs and the rest were other multiple-birth siblings) from 9189 discordantly exposed families; and (3) the active comparator cohort incorporated 434,900 BZDR recipients and 428,074 SSRI recipients with complete information on all covariates used for weighting in IPTW (Fig. 1 and Table S2).

We defined outcomes as the first record of any infection which required an inpatient or specialist outpatient care and was recorded as the main or secondary diagnosis in the National Patient Register, or as death due to any infection which was considered the underlying death cause in the Cause of Death Register. The diagnoses in question and the corresponding International Classification of Diseases, Tenth Edition (ICD-10) codes are reported in Table S3. To assess the risk of infections in proximity to initiation of BZDR treatment and in line with prior studies [16], the outcomes were collected if they occurred within 12 months after the index date. We constructed outcome variables as “any serious infections,” “viral infections,” “bacterial infections,” and “mixed and others,” with the latter referring to simultaneously recorded viral and bacterial infections or to infections with unknown origin. The outcome definition did not include cases of COVID-19 infection.

We collected demographic data (birth year and month, sex, county of residence at the index date), socioeconomic characteristics (household disposable income the year before the index date), calendar year of the index date, history of psychiatric and somatic diagnoses before the index date (with ICD-10 codes for the diagnoses listed in Table S4), and records of other medications dispensed within 3 months prior to the index date (with ATC-codes listed in Table S1). Also, to rule out a general susceptibility to infectious diseases, we collected data on history of inpatient and specialist outpatient diagnoses of any infections recorded before 2-year infection-free period that preceded the index date (using the same ICD-10 codes from the National Patient Register as for the outcomes).

We used Cox proportional hazards regression models with the days since the index date as the underlying time scale. In all three cohorts, BZDR recipients and their unexposed counterparts were followed-up from the index date until the date of outcome event, emigration, death other than by infection, end of a 12-month follow-up period, or end of study on 31 December 2020, whichever came first. Unexposed individuals in the demographically matched and co-twin cohorts as well as BZDR recipients and SSRI recipients in the active comparator cohort were additionally censored at the date they changed the exposure status, if such change occurred during the follow-up. The hazard ratio and 95% confidence intervals were estimated for any serious infections and, separately, for viral, bacterial, and mixed/other infections. For each cohort, baseline crude models were followed by models adjusted for all study covariates. For the demographically matched and co-twin cohorts, the covariates included in the fully adjusted model are listed in Tables S5 and S6, respectively. For the active comparator cohort, individuals were weighted by the inverse of the predicted probability of their observed exposure to incident BZDR or SSRI use conditional on all potential confounders listed in Table S7. We truncated stabilized weights at the 1st and 99th percentiles to minimize impact of extreme values [33, 36]. The use of the inverse probability of treatment weighting resulted in balanced covariates between BZDR recipients and SSRI recipients (all standardized differences were < 0.05 [36]; Table S7). In all three cohorts, the main analyses were followed by an exploratory analysis where models were stratified by sex.

A series of additional analyses were performed in the demographically matched cohort with “any serious infections” as the outcome of interest. First, we re-ran the main analysis with the exposure subdivided by the type of the initial BZDR (as anxiolytics, hypnotics/sedatives, antiepileptics, Z-drugs, based on their ATC codes, or as poly-BZDRs, if more than one BZDR was dispensed at the treatment initiation). Second, we further subdivided the exposure by the specific initial BZDR drug. Third, exposure to incident BZDR use was categorized according to the cumulative dosage [measured by the defined daily dose (DDD)], for which we additionally collected information on all BZDR dispensations under the follow-up for each BZDR recipient. Cumulative dosage was assessed as a time-varying exposure and categorized based on tertiles into ≤ 20 DDDs, > 20 DDDs ≤ 65, and > 65 DDDs (Note S2 provides the details). Lastly, we assessed the effect modification by age at the index date. We modeled age using restricted cubic splines with 5 knots placed at 5, 27.5, 50, 72.5, and 95 quantiles and added an interaction term between exposure and the splined age to the fully adjusted model.

For all tests, we employed two-tailed significance set at p < 0.05. Data management and analyses were performed using SAS, version 9.4 (SAS Institute Inc).

No patients were involved in setting the research question or the outcome measures, nor were they involved in the design and implementation of the study. As this is a register-based study based on de-identified data, there are no plans to disseminate the results of the research directly to study participants. The study is co-authored by researchers affiliated to the Stockholm Health Care Services, Region Stockholm, where results will be widely disseminated. Dissemination to the broader Swedish population will be achieved through media outreach on publication of this study.