Kashin-Beck disease (KBD) is an endemic, deformative osteoarthropathy characterized by distinct geographic clusters^[1]. Historically, 13 provinces and autonomous regions in China have been designated as KBD-affected areas. In severely endemic areas, the radiographic detection rate had reached up to 80%^[1]. Although comprehensive prevention and control measures, including dietary restructuring, grain substitution, and relocation of educational facilities, have successfully halted new pediatric cases, the persistent threat posed by environmental risk factors suggests that the relaxation of these measures could lead to the resurgence of KBD.

The most widely accepted etiological hypothesis for KBD is the combined “environmental selenium deficiency–fungal toxin exposure” model. Selenium deficiency creates a critical nutritional environment, while exposure to the T-2 toxin, produced by Fusarium spp. in cereals, directly triggers pathogenesis. Natural environmental factors associated with KBD-affected areas include cold and humid climates, significant diurnal temperature variations, and mountainous topography^[2]. Soil characteristics, such as water-holding capacity, organic matter, nitrogen content, and salinization-alkalization processes, directly influence soil microecology in these regions. These conditions favor fungal proliferation and toxin production while simultaneously regulating the biogeochemical cycling of selenium^[3]. Additionally, previous studies suggest that KBD predominantly affects economically disadvantaged areas with poorly developed transportation infrastructure^[4].

The current understanding of the environmental etiology of KBD primarily relies on descriptive summaries, with limited quantitative analysis of synergistic environmental factors. Furthermore, traditional statistical methods are inadequate for analyzing hierarchical data structures such as villages nested within different selenium zones. Ignoring within-cluster correlations may severely underestimate parameter standard errors, artificially narrow confidence intervals, and consequently, increase type-I error rates^[5]. To address these gaps, this study aimed to (1) quantify the impact of key environmental factors on KBD severity in endemic villages and (2) identify environmental factor combinations that present high risk, thereby providing a scientific basis for targeted prevention.

We used data from the National KBD Surveillance System of the Chinese Center for Disease Control and Prevention, which covers the 13 historically endemic provinces (including autonomous regions) in China. The analysis included 16,512 endemic villages (Supplementary Figure S1A; administrative map sourced from the National Platform for Common Geospatial Information Services; Review No. GS(2024)0650). KBD-endemic villages were categorized as mild, moderate, or severe according to the national criteria for designating and classifying endemic areas of KBD (GB 16395-2011; Supplementary Tables S1 and S2). Stratified random sampling was employed to select control villages without a history of KBD from each endemic province. In the present study, the outcome variable was ordinal. “Risk” refers to the probability or odds ratio of a village being classified within a more severe KBD category, serving as a measure of association between environmental factors and increased disease severity.

The environmental and socioeconomic predictive variables used in this study were obtained from publicly available national databases. Detailed descriptions and coding of the variables are presented in Supplementary Table S3. The topographic factors included altitude and slope. The climatic variables included mean annual temperature (Bio1), annual precipitation (Bio12), relative humidity (RH), and isothermality (Bio3). The soil physical properties included gravel content (T_GRAVEL), sand content (T_SAND), silt content (T_SILT), and porosity. The soil chemical properties included the exchangeable sodium percentage (T_ESP), electrical conductivity (T_ECE), pH, alkali-hydrolysable nitrogen (AN), and organic matter. Socioeconomic indicators included gross domestic product (GDP) and distance to the nearest railway (rail). Selenium zones were classified into three categories: low, marginal, and moderately-high selenium (Supplementary Table S4).

Univariate screening was performed before multilevel modeling. Continuous variables are expressed as mean ± standard deviation (x±s). Differences among severity groups were tested using one-way analysis of variance (ANOVA). Categorical variables were summarized as frequencies (percentages), and differences between groups were assessed using the chi-square test. Statistical significance was set at α = 0.01. All continuous variables were standardized using Z-scores before inclusion in the model. Multicollinearity was assessed using Pearson’s correlation coefficients (with |r| > 0.8 indicating strong correlation) and variance inflation factors (VIF; a VIF > 5 indicated severe multicollinearity).

A Bayesian two-level ordinal logistic regression model was employed to address the hierarchical structure of villages (level-1) nested within the selenium zones (level-2). Village-level environmental variables and selenium zone categories were modeled as fixed effects. Random intercepts captured between-zone variability, allowing a systematic analysis of the associations between each variable and KBD severity in endemic villages.

Let y_ij represent the KBD severity grade (ranging from 0 to 3) of the i-th village within the j-th selenium zone. Then,

To evaluate the significance of the random effects associated with selenium zones, a null model containing only an intercept and random effects was initially fitted. The model components are defined as follows:

In the above equations, y_ij denotes the severity grade of the i-th village within the j-th selenium zone; κ_k represents the k-th threshold parameter for ordinal classification; x_ij represents village-level predictors; z_ij indicates the selenium zone; u_j is the random intercept for the j-th selenium zone, assumed to follow u_j ~ N(0, σ²); and H(·) denotes the cumulative logistic distribution function, used to compute cumulative probabilities.

The model estimation employs weakly informative priors. Fixed-effect parameters (β) were assigned N(0, 1.5) priors for regularization. The model was constructed using the brms package in R with the cmdstanr backend employing the No-U-Turn Sampler (NUTS) for Markov chain Monte Carlo (MCMC) sampling (three chains, 1,200 iterations each). Posterior odds ratios (ORs) and 95% credible intervals (CrIs) were calculated to evaluate predictor effects. Model convergence was assessed using the potential scale reduction factor (Rhat), with values approaching 1.0 indicating stable convergence. To investigate nonlinear dose–response relationships, restricted cubic splines (RCS, four degrees of freedom) were incorporated for significant predictors. The knot positions were automatically selected using splines::ns() based on variable quantiles. The P for overall assessed global significance and P for nonlinear tested deviations from linearity (significance level < 0.01).

The final analysis included 33,022 villages: control (50.0%), mildly endemic (28.0%), moderately endemic (7.3%), and severely endemic (14.7%). Significant differences (P < 0.01) among these groups were observed in altitude, slope, climatic factors (Bio1, Bio12, RH, and Bio3), soil properties (T_GRAVEL, T_SAND, T_SILT, porosity, T_ESP, T_ECE, pH, AN, and organic matter), socioeconomic factors (GDP and distance to railways), and selenium backgrounds (Supplementary Table 5S). Eighty-four villages were located in areas with undefined selenium landscapes; these were excluded from the analyses. However, sensitivity analyses remained statistically significant after their exclusion (Supplementary Table 6S). Further analysis revealed a strong negative correlation between T_SAND and T_SILT (r = −0.91), with VIF exceeding 5, indicating severe multicollinearity. T_SAND was retained in the final model, whereas the other variables showed acceptable multicollinearity (VIF < 5; |r| < 0.8) (Supplementary Figure 2S; Supplementary Table 7S).

The Bayesian multilevel ordinal logistic regression model showed strong convergence (Rhat ≤ 1.01). A significant random effect at the selenium zone level (ICC = 0.47; σ = 1.87; P ≤ 0.01) justified the use of hierarchical modeling (Supplementary Table 8S). The primary analysis indicated that altitude, slope, Bio3, Bio12, T_SAND, pH, and increased distance to railways were associated with greater disease severity in KBD-endemic villages. Conversely, Bio1, porosity, T_ESP, and organic matter exhibited protective effects. GDP, RH, T_GRAVEL, T_ECE, and AN were not significantly associated with the severity of KBD. Compared to low-selenium zones, marginal and moderately-high selenium zones tended toward reduced risk, although their credible intervals included zero (Table 1, Figure 1A). Sensitivity analyses further validated these findings: (1) combining marginal and moderately-high selenium zones into a non-low selenium group provided results consistent with the primary analysis (Supplementary Table 9S; Supplementary Figure 3S); (2) converting the outcome into a binary classification (endemic/non-endemic villages) also revealed significant clustering (ICC = 0.42), with T_ECE emerging as a risk factor, whereas the effects of other predictors remained consistent (Supplementary Tables 10S–11S, Supplementary Figure 4S).

Figure 1.  Association of environmental factors with KBD severity based on Bayesian multilevel ordinal logistic regression. (A) Forest plot of the fixed effects for primary environmental factors. (B) Forest plot of the interaction effects between selenium zones and soil properties. Green squares represent the posterior OR, with their size proportional to the estimation precision. The horizontal lines depict the 95% CrIs. The vertical dashed line indicates the null effect (OR = 1). Abbreviations: OR, odds ratio; CrI, credible interval; GDP, gross domestic product; rail, railway length; Bio1, mean annual temperature; Bio12, annual precipitation; RH, relative humidity; Bio3, isothermality; T_GRAVEL, gravel content; T_SAND, sand content; T_ESP, exchangeable sodium percentage; T_ECE, electrical conductivity; AN, alkali-hydrolysable nitrogen; Se, selenium.

RCS analysis identified significant nonlinear relationships (P < 0.001) between KBD severity and altitude, slope, Bio1, Bio12, porosity, and organic matter. Specifically, at altitudes ranging from 237 to 4111 m, the risk of KBD first increased, reaching a peak at 953 m, followed by a decline, and then a sharp rise above 4111 m. Slope and Bio12 displayed N-shaped associations with disease risk. The risk initially increased and then decreased for slope values between 2.7% and 26.2%, peaking at 12.3%. For Bio12, the highest risk occurred at approximately 816 mm, within the range of 578–1079 mm. Consistent with the linear analysis, Bio1, porosity, and organic matter negatively correlated with the severity of KBD (Figure 2). The spatial distribution of KBD-endemic villages overlaid with key environmental determinants is shown in Supplementary Figure 5S.

Figure 2.  Nonlinear dose-response relationships of environmental factors with the risk of being classified as KBD-endemic areas (RCS analysis). Note: The red solid line represents the estimated ORs of KBD severity across the range of each environmental variable, and the blue shaded area indicates the 95% CrI. The P for overall assesses the overall significance of the variable, while the P for nonlinear evaluates the degree of deviation from linearity. Abbreviations: OR, odds ratio; CrI, credible interval; Bio1, annual mean temperature; Bio12, annual precipitation.

Overall, this study systematically evaluated the environmental determinants influencing the severity of KBD-endemic villages. Among the socioeconomic factors, regional GDP was not significantly associated with KBD risk, suggesting that natural environmental influences may outweigh regional economic status. In agreement with previous findings, villages farther from railways exhibited an increased risk, potentially because of poor grain storage conditions, limited access to external food sources, and inadequate coverage of nutritional intervention in remote areas^[4].

The findings of this study revealed that cold and humid conditions increase the risk for villages being classified as KBD-endemic. This aligns with previous experimental evidence that cold and damp storage significantly increases T-2 toxin accumulation in grains^[5], suggesting a critical pathogenic mechanism. Disease risk increased notably when annual precipitation ranged from 578 to 1079 mm and slope gradients varied between 2.7% and 26.2%. This might be because steep slopes form cold, humid, and poorly drained valley microenvironments, favoring mold growth and exacerbating soil erosion, thereby reducing trace elements such as selenium^[6]. Contrary to conventional assumptions, isothermality exhibited a weak positive correlation with KBD risk, suggesting the need for further examination using more precise climate–temperature variability indicators^[6]. Additionally, precipitation was significantly correlated with an increased KBD risk, whereas relative humidity showed no significant association. This discrepancy may arise because rainfall directly elevates grain moisture, promoting Fusarium infestation and toxin production, whereas macroscale relative humidity measures may inadequately reflect the humidity within grain storage microenvironments.

This study identified soil sand content as a risk factor for KBD, which is consistent with previous findings that sandy soils promote Fusarium survival, decrease selenium retention, and limit selenium uptake by plants^[7]. Conversely, soil porosity (defined as the pore space relative to soil volume, influencing water retention and aeration) acted as a protective factor, likely by reducing waterlogging conditions that facilitate Fusarium toxin production^[7].

Interaction analysis indicated that the selenium zone classification served as both an essential background context and a significant modifier, substantially influencing local environmental risk (Figure 1B, Supplementary Table 12S). Exchangeable sodium (reflecting soil alkalization) and electrical conductivity (indicating salinization) emerged as significant risk factors within moderately high selenium zones. This effect is likely mediated by enhancing Fusarium proliferation, pathogenicity, and toxin production, thereby overwhelming the antioxidant defense mechanisms in the host^[8]. In contrast, alkali-hydrolysable nitrogen, soil pH, and organic matter exhibited protective roles, potentially because adequate selenium improves plant disease resistance, alkaline conditions increase selenium bioavailability, and organic matter facilitates selenium chelation and availability^[9]. In marginal selenium zones, alkali-hydrolysable nitrogen was a risk factor, possibly because elevated nitrogen levels dilute the already limited selenium availability and provide additional nutrients that the pathogens can exploit^[10].

Methodologically, this study utilized a Bayesian multilevel model to effectively manage hierarchical data structures and address the biases inherent in traditional methods. Theoretically, these findings further substantiate the “environmental selenium deficiency–mycotoxin poisoning” hypothesis. These results provide a foundation for targeted zoning of KBD risk. However, this study also has some limitations. First, the use of categorical selenium zones rather than quantitative selenium measurements may compromise the accurate characterization of exposure–response relationships. Second, this ecological study, conducted at the village level, did not include individual-level factors (e.g., age and dietary habits), restricting the assessment of personal risk determinants. Third, although the multilevel approach accounted for clustering within selenium zones, it did not adjust for spatial autocorrelations. Finally, although historical environmental data largely coincided with the main KBD epidemic period and the chosen factors were generally stable across time and space, precise temporal alignment with individual disease onset was not achieved. Future studies should incorporate quantitative selenium measurements and integrate critical individual-level covariates, account for spatial effects, and acquire temporally aligned exposure data. Such advancements would enable more precise predictive modeling and enhance early warning in high-risk regions.

The risk of villages for being KBD-endemic results from the synergistic influence of multiple environmental factors, notably cold and humid climates, steep terrain, and soil with low porosity and limited organic matter. Considerable heterogeneity was observed among the different selenium zones. Specifically, high soil nitrogen content emerged as a significant risk factor in marginal selenium zones, whereas soil salinization–alkalization presented as the primary threat in moderately-high selenium areas. These results offer a scientific basis for effective risk zoning and targeted interventions, including improved grain mold prevention strategies and enhanced selenium nutrition monitoring.