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.