NHANES is a complex, multistage sampling design, nationally representative study to assess the health and nutritional status of non-institutionalized civilians in the United States [14]. Since 1999, it has been a continuous project, with each cycle lasting two years [14]. NHANES was conducted by the Centers for Disease Control and Prevention (CDC) and and approved by the institutional review board of the National Center of Health Statistics. All participants provided written informed consent.

We used nine cycles of the NHANES from 2001 to 2018 (vitamin D data were not available in the NHANES 1999–2000). Inclusion criteria for the study were as follows: 1) age ≥ 18 years; 2) diagnosed with CKD [eGFR < 60 ml/min/1.73 m2 or urine albumin-to-creatinine ratio (UACR) ≥ 30 mg/g)] [13]; 3) did not start dialysis within 12 months before baseline; 4) eligible for follow up.

The serum phosphorus was measured using a Hitachi model 737 multichannel analyzer (Boehringer Mannheim Diagnostics, Indianapolis, IN), by observing changes in absorbance at 365 nm following phosphomolybdate formation from inorganic phosphorus in an acidic solution, which correlate directly with phosphorus levels [15].

Age, gender, race/ethnicity, co-medications, and comorbidities were collected from household interviews using standardized questionnaires [14]. Body weight, height were obtained when people participated in physical examinations at a mobile examination center [14]. BMI was calculated as weight in kilograms divided by height in meters squared. Race was classified as non-Hispanic White, non-Hispanic Black, Mexican American, or other. Additionally, serum calcium, creatinine, serum albumin, hemoglobin, urinary albumin, and creatinine were measured at baseline when the participants provided their blood and urine samples [14]. UACR was calculated by dividing urinary albumin by urinary creatinine [13]. eGFR values were calculated using the creatinine equation developed by the Chronic Kidney Disease Epidemiology Collaboration (CKD-EPI) [16], and CKD was staged based on eGFR according to the recommendations of the KDIGO guidelines [13], with stages 1 (eGFR ≥ 90 mL/min/1.73m2), 2 (eGFR 60–89 mL/min/1.73m2), 3a (eGFR 45–59 mL/min/1.73m2), 3b (eGFR 30–44 mL/min/1.73m2), 4 (eGFR 15–29 mL/min/1.73m2), or 5 (eGFR < 15 mL/min/1.73m2). Diabetes was defined as self-reported doctor-diagnosed diabetes, or use of insulin or oral hypoglycemic drugs, or fasting blood glucose ≥ 7.0 mmol/L, or random blood glucose ≥ 11.1 mmol/L, or glycated hemoglobin A1c (HbA1c) ≥ 6.5% [17]. Hypertension was defined as self-reported doctor-diagnosed hypertension, or systolic blood pressure ≥ 140 mmHg, or diastolic blood pressure ≥ 90 mmHg [17]. CVD included self-reported coronary heart disease, congestive heart failure, heart attack, stroke, and angina [17].

Mortality from any cause and CVD was ascertained by linkage to the National Death Index through 31 December 2019. The ICD-10 was used to determine disease-specific death. CVD mortality was defined as the primary cause of death being any disease of the circulatory system (ICD-10 codes I00-I09, I11, I13, I20-I51, or I60–I69).

Sample weights, clustering, and stratification were incorporated into all analyses because of the complex sampling design of the NHANES [14]. Participants were followed up to death or the date of 31 December 2019, whichever comes first. We take 2/9 of the 4-year weight for each individual sampled from 2001 to 2002, and 1/9 of the 2-year weight for each individual sampled from 2003 to 2018.

Serum phosphorus was analyzed as a continuous variable, or categorized into three groups: < 3.5 mg/dL, 3.5 to < 4.5 mg/dL, and ≥ 4.5 mg/dL. Baseline characteristics were summarized as means (standard errors) for continuous variables and numbers (percentages) for categorical variables. Time-to-event data were described using Kaplan–Meier curves, and the between-group difference was compared using log-rank test. The Cox proportional hazards model was used to analyze the association between serum phosphorus and all-cause or CVD mortality, with or without adjustment for confounders, including age, sex, race, hemoglobin, eGFR, serum albumin, serum calcium, 25 (OH) vitamin D, obesity, hypertension, diabetes and CVD. Additionally, a restricted cubic spline Cox regression, with the smallest Akaike information criterion (AIC) of three knots, was performed to test for linearity and explore the shape of the dose–response relation of serum phosphorus concentrations and mortality.

We further performed subgroup analysis stratified by baseline characteristics including eGFR (< 60, ≥ 60 ml/min/1.73 m2), sex, age (< 65, ≥ 65 years), presence/absence of CVD and obesity. We examined possible effect modification by introducing multiplicative interaction terms between phosphorus and the grouping factor into our primary Cox regression model.

Additional sensitivity analyses were conducted to verify the robustness of the results. Firstly, considering the interaction between vitamin D, blood calcium, parathyroid hormone (PTH), and serum phosphorus. We further adjusted for PTH (only provided in NHANES 2003–2006). Secondly, to reduce potential reverse causality bias, we included participants with a follow-up of more than 2 years. Thirdly, additional adjustments were conducted for co-medications (RAS inhibitors, other anti-hypertension drugs, lipid-lowering drugs, and hypoglycemic) and survey years. Finally, We made multiple imputations for missing values to avoid selection bias.

All analyses were carried out with R version 4.2.0, and a two-tailed P < 0.05 was considered statistically significant.