Molybdenum is not a risk factor for changes in serum testosterone

Abstract Molybdenum is a bio-essential trace nutrient, and elevated urinary molybdenum has been implicated as correlating with lower testosterone levels in men. To properly address the contention that molybdenum constitutes a potential health risk factor, we reviewed statistical research findings on molybdenum exposure and somatic levels in humans in a biological plausibility context, which is an indispensable component of reliable risk factor determination. Our analysis applied advanced nonparametric statistical modeling that expands linear regressions to include more potential variables and common confounders. The analysis shows that previously published negative molybdenum-testosterone associations no longer retain statistical relevance. In place of molybdenum as a stressor responsible for lowered testosterone levels, body mass index, age, related hormone variability, and underlying global downward trends in testosterone are shown to be the factors responsible for this effect.


Introduction
Humans require molybdenum (Mo) as an essential nutrient. It acts via the molybdenum cofactor (MoCo) biosynthesis pathway in various enzymes to catalyze the carbon, nitrogen and sulfur cycles (Mayr et al. 2021), which are critical for multiple beneficial biological processes. However, whether and how excess exposure to molybdenum might create adverse effects in humans remains a complex area of study. This paper focuses on the biological plausibility and weight of evidence, both scientific and statistical, related to effects of Mo on serum testosterone (T) levels in men. Exploratory statistical analyses, using Mo concentrations in spot samples of urine as a surrogate for exposure and presented without consideration of the full biological picture, have been used to suggest that Mo could reduce serum T concentrations and potentially pose a risk to male fertility (Lewis and Meeker 2015;Meeker et al. 2010). Regardless of the statistical approach taken, past analyses have not assessed the weight of the biological evidence for or against the hypothesis that Mo causes (as distinct from merely being associated with) changes in T, and notably when research findings clearly conflict. For example, Alves et al. (2020), Rotter et al. (2016) and others have found positive correlations between Mo and T in blood. Other investigators found no relation between Mo and T. Using multiple linear regression, Zeng et al. (2013) found significant negative associations between urinary manganese and zinc, and serum levels of T among patients of a fertility clinic, but no relationship between urinary Mo and T. Building on Cox (2021), we further expanded previous analysis using non-parametric methods (including Classification and Regression Tree [CART] and Bayesian Networks [BN]) that confirm body mass index (BMI), age, and fasting time as robust predictors of T, without confirming the correlation with Mo identified by Lewis and Meeker (2015).
Biological plausibility has been proven as an indispensable component of reliable risk factor determination. Because the literature relying on such statistical evidence appears to be mixed and inconclusive with regard to T serum levels related to Mo, a review of the biological evidence was warranted and appears in the Discussion section of this paper, with key foundational points as follows.
Pinpointing causes of low T is often difficult because of multifactorial influences stemming from the environment, comorbidities, exposure, and other contributing factors. Also, an overall global downward trend toward lower T has been observed in the last several decades on a worldwide scale. The statistical analyses and visualizations herein suggest that BMI, along with uric acid, age, and fasting time, are mathematically robust predictors of T in both parametric and nonparametric models. By contrast, Mo (or its logarithm), income, cotinine, and cesium are not robust predictors, as they can be included or excluded from statistical models without significantly affecting explanatory power (Cox 2021).
We supplement our review of the biological evidence with additional analyses to clarify whether and to what extent statistical methods provide reliable evidence about a Mo-T relationship. The implications of the statistical analyses are discussed within a core context of biological plausibility, and further statistical detail, references, tables, and graphs are available in the Supplementary Material.

Material and methods
As a base dataset, we compiled 2,271 cases of National Health and Nutrition Examination Survey (NHANES) data on urinary Mo and serum T from three consecutive two-year cycles (years 2011-2016), to produce a dataset nearly five times larger than that of Lewis and Meeker (2015) who used a single NHANES cycle (2011)(2012). We conducted further analyses in addition to those previously published by Cox (2021) to characterize the dependence of T on Mo for varying levels of other variables. These analyses included multiple linear regression analyses, as well as non-parametric ('modelfree') analyses of dependencies among variables, using CART analyses and BN learning. CART analyses iteratively partition data from matched sets, consisting of a paired dependent (response) variable and one or more corresponding independent (explanatory) variables, into progressively smaller groups. CARTs are part of the U.S. Environmental Protection Agency's Causal Analysis/Diagnosis Decision Information System "toolkit" (De'ath and Fabricius 2000;EPA 2022). The specific benefits of CART analysis for this purpose were reviewed in Cox (2021). The full methodology for BN learning is quite detailed; readers may wish to consult Nagarajan et al. (2013) for details of the approach.

Linear and ordinal correlations between individual variables
Linear correlations between individual variables (Supplementary Material, Table S1) suggest that, in addition to potential confounders such as age, income, body mass index, and serum cotinine (a marker of cigarette smoking status) considered by Lewis and Meeker (2015), fasting time and time of day (e.g., when urine and blood sampling occurs) are two additional time-related potential confounders that should be considered. These two additional confounders are variables easily overlooked without biological consideration of molybdenum toxicokinetics, including quick urinary excretion, which is understood to be generally biphasic in humans (ATSDR 2020). Urinary cesium, iron in refrigerated blood serum, sex hormone-binding globulin (SHBG) concentration in blood, and number of people living in the same house ("household size") were also identified as significant correlates of both Mo and T in the linear regressions (Supplemental Material).
We complemented these linear correlations by displaying ordinal correlates of serum testosterone and urinary molybdenum with all other variables in the dataset, which do not assume a linear relation between pairs of variables (Suppl., Figure S1). This exploratory analysis aims to automatically highlight possibly important correlations, without manual selection or deselection of variables, to avoid human error or bias. Urinary Mo is only a weak negative ordinal correlate of T. Several variables (e.g., blood lead, serum cotinine, blood cadmium, and urinary creatinine) are correlated with both T and Mo, suggesting that these should be considered in further statistical analyses.
In both linear and ordinal correlations between individual variables, SHBG was the top variable that is most related to (positively correlated with and influences) serum T.

Multiple linear correlations
Although the lognormal transformation of the dataset may not have been strictly necessary in Lewis and Meeker (2015), a key finding of their analysis was the identification of a significant negative linear regression coefficient between log(Mo) and log(T). Logtransformation was used by us only to facilitate reproduction of this finding in the 2011-2012 data and was likewise applied to confirm the observation in the analysis of the expanded dataset across three NHANES cycles (2011)(2012)(2013)(2014)(2015)(2016). The regression coefficient for log(Mo) remains significantly negative in the expanded dataset, with p ¼ 0.01 (Cox 2021).
However, including additional predictors in a larger multiple linear regression model (Suppl., Table S2), makes log(Mo) no longer a significant predictor of log(T) at the conventional 0.05 significance level (p ¼ 0.16). Instead, in the example model, backward variable selection identifies age, BMI, household size, fasting time, time of day, and uric acid as independent predictors of log(T). If uric acid is excluded, no other variables enter into the model; thus, a hypothesis that Mo affects T only through uric acid is not supported. All-subsets regression modeling showed that body mass index, uric acid, age, and fasting time, but not log(Mo), income, cotinine, or cesium, are included in all of the best-fitting statistical models as predictors of log(T).
Whether or not log(T) significantly depends on log(Mo) in multiple linear regressions, appears to be a model-dependent finding. Thus, log(T) depending on log(Mo) is driven by the choice of predictors used in the model.

Introduction
Multiple linear regression modeling is limited by assuming linear relationships between its explanatory and dependent variables, as well as by other assumptions such as normally distributed error terms with constant variance. Exploratory visualizations of the data, however, indicate substantial departures from linearity. By way of example, the nonparametric smoothing regression curve (LOWESS) for serum T vs. body mass index falls well outside the 95% confidence bands for the linear model, indicating that the linear regression model may be too restrictive to accurately describe the data (Suppl., Figure S2). Similar conclusions hold true for multivariate linear regression models.
Classification and regression Tree (CART) analysis As a nonparametric alternative to multiple linear regression analysis, CART analysis results with T as the dependent (response) variable (see Suppl. for details), and each of the possible independent (explanatory) variables are presented in Figure 1. Urinary Mo, correcting for creatinine or not, was not identified as a useful independent predictor of T in this CART analysis, which is in contrast to the implications of Lewis and Meeker (2015). Instead, Figure 1 demonstrates that the SHBG in mena variable omitted by Lewis and Meekerwas the most important predictor for T. Since SHBG transports androgens and estrogens in blood and regulates their access to target tissues, its inclusion in the analysis is a valuable expansion to understanding the true correlates for T.
Bayesian network (BN) analysis Another way of looking at possible statistical dependences between variables is by learning a BN from the data. BNs indicate the presence or absence of a statistical dependence between variables (i.e., if knowing the value of one variable helps to predict the value of another). The fact that neither unadjusted urinary Mo nor creatinine-adjusted urinary Mo is directly linked to T in such BNs (Suppl., Figures S4 and S5) is consistent with the null hypothesis that serum T does not depend directly on urinary Mo.

Discussion
We assembled a robust NHANES dataset to explore exposure-response correlations and regression coefficients for urinary Mo and serum T. This analysis finds that in multiple linear regressions, whether or not T significantly depends on Mo is model-driven, and further, that outcomes are influenced by the choice of predictors employed in the models. In contrast to urinary Mo, other variables such as BMI, age, and fasting time are included in all the best-fitting multiple linear regression models. The narrative that follows considers various biological factors that enhance our statistical observations.
Body mass index and obesity All our statistical approaches identified BMI as a robust predictor of serum T. Biologically, obesity is a key factor in interpreting testosterone trends and correlations (El-Shehawi et al. 2020). Within the selected NHANES dataset, 744 of 2,258 respondents (33%) had a BMI of 30 or more (there was no BMI data for 13 out of 2,271 samples). Increasing understanding of how obesity adversely affects male reproductive health and induces low serum testosterone levels has improved our ability to have confidence in the statistical findings. The current work expands on the interpretation of Meeker et al. (2010), who acknowledged that BMI is inversely associated with testosterone (but who also found an inverse correlation with molybdenum in blood) among 219 patients at a Michigan fertility clinic. When considering the more recent literature and non-obese subjects, the opposite has been observed: Alves et al. (2020) studied 60 runners and found blood molybdenum in runners (range 0.1-3.33 mg/L, mean 0.62 ± 0.59 lg/L) to be positively correlated with serum T (range 4.13-9.94 ng/mL, mean 6.52 ± 1.14 ng/mL). Likewise, Maynar et al. (2018) found Mo content in serum from 80 professional athletes to be higher among sportsmen (2.083 ± 0.674 lg/L Mo for 24 anaerobic athletes; 1.622 ± 0.711 lg/L for mixed athletes; and 0.880 ± 0.822 lg/L for 28 aerobic athletes) than in a 31-person control group (0.339 ± 0.115 lg/L Mo; statistically significant with p < 0.001 for each of the three sports types). A key difference between the control and athlete groups was BMI (controls had higher BMI of 23.81 ± 3.14 versus athlete BMI 20.25 ± 1.73, p < 0.001). Maynar et al. (2018) explain this trend: Mo participates as an integral part of oxide-reduction processes via an enzyme (xanthine dehydrogenase, XDH) which catalyzes the hypoxanthine transformation of xanthine to uric acid (an antioxidant), and this biochemical process is augmented especially among athletes under anaerobic exercises. Differences in Mo values were attributed to a higher production of free radicals by ischemia-reperfusion in the athlete group: Mo presence eases uric acid formation and decreases the damage caused by superoxide anions generated by xanthine oxidase in the ischemia-reperfusion processes, "[a] situation induced by high intensity muscular activities" (Maynar et al. 2018). While statistical correlates between T and uric acid were reported (Suppl., Table S1), this was not the case with Mo and uric acid, which conforms to the biological evidence noted here from recent studies in athletes. Figure S1 confirms the classic negative correlation of T with increasing age. Typically, males will show a decline in testosterone and related hormone levels with increasing age. Not all datasets exploring serum T and Mo evaluate the same age ranges, but controlling for age is essential for proper conclusions regarding T. For example, Rotter et al. (2016) evaluated 313 men aged 50-75 years, compared to younger men aged 18-55 in Lewis and Meeker (2015), Cox (2021) and this paper. Rotter et al. concluded that molybdenum and serum T are positively correlated for age 50-75 (in contrast to Lewis and Meeker (2015) and their reported inverse correlation for men aged 18-55). Thus, data-driven Mo-T associations may vary by age. This calls for the careful inclusion of age ranges in all statements of causality or correlation.

Age (and related overall global trend)
Studies in various regions and multiple countries are showing more age-specific decreases in total T and free T levels (Andersson et al. 2007;Chodick et al. 2020;Travison et al. 2007). Amongst men of the same age, multiple studies indicate that men born in later periods showed significantly reduced T. Andersson et al. (2007) examined over 5,000 men in Denmark, aged 30 to 70 in birth cohorts from 1921 to 1970, for serum T and SHBG levels. When age-matched, there were trends toward both lower T and SHBG levels amongst later birth cohorts and study periods. This finding was significant once adjustments for BMI were made. The authors concluded that factors beyond age and BMI are influencing global-scale decreases in T compared to earlier cohorts, with a distinct biological role for SHBG levels (which factor prominently into the findings shown in Figure 1, as the most important predictor for T). Even more recently, Chodick et al. (2020) observed testosterone levels in over 100,000 Israeli men between the ages of 13 to 80, studied from 2006 to 2019. These data show that age-independent global declines in testosterone are still being observed in relatively recent populations, and that the declines are not explained by BMI trends.
Population studies that examine T levels across birth cohorts, time periods and regions, as cited above, suggest a widespread rise of low T phenotypes, perhaps reflecting changes over time in factors such as disease incidences, fertility rates, lifestyle, industrialization, demographics, and so forth. Decreases in T levels over time have been confirmed in studies across several developed nations, so it is overly simplistic to posit a relationship caused by environmental agents only on statistical correlations, without first more fully understanding biological inputs to T levels and trends. In addition to the influence of BMI, Figure 1 showed NHANES cycle as a predictor of T (an indirect or secondary representation of "population age" overall, simply by virtue of the difference in years of analysis), and thus this unexplored hypothesis that the measured declines in T are merely mirroring the global decline in T deserves further consideration.
Fasting time: biological support for statistical predictor of serum testosterone This section seeks to clarify the weight of biological evidence surrounding fasting time (and related sampling time of day) in the context of a possible causal relationship between Mo and T. Similar to BMI, fasting time was a robust predictor of serum T in various statistical procedures, although it did not have the same weight in the Figure 1 CART analyses as BMI and hormone influences (SHBG and estradiol). An often-cited study suggests that fasting boosts serum T (statistically significantly) at 56 h, but not at 8 h (R€ ojdmark 1987). In addition, Long et al. (2015) recommend (since the natural diurnal variation in T levels tends to diminish with age) that it is acceptable to test men aged 45 and older before 2 pm, but variation is still present in younger men, whereby data taken morning and afternoon may not be comparable. Sample collection times for NHANES vary widely: the dataset includes both sampling time of day (morning, afternoon, or evening) and fasting times, with twenty (out of 2,271) subjects who had fasted longer than 24 h prior to serum collection, with a suspected but unknown impact of fasting status and collection time on T levels independent of Mo. When evaluating fasting time and sampling time of day in relation to Mo, another layer of complexity is added: An overall dietary Mo urinary half-life of 12 h (ATSDR 2020), paired with evidence that Mo elimination is approximately biphasic (with mean half-times of 30 min and 6.6 h (ATSDR 2020), suggests that fasting time could be important for Mo as well as T. Hays et al. (2016) found urine to be sufficiently reliable as a biomarker for Mo in the general public despite diurnal variation, noting the less homogeneous the intake (such as in an occupational setting), the more important it could be to make adjustments to achieve a representative sample. Overall, fasting time may affect conclusions of multiple regression analyses for serum T or urine Mo, or both.

Uncertainty and future directions: understanding genetic variability
The biological foundation for low T may include factors not included in the statistical analysis or discussion above. For example, genetic variability is an important confounder not directly accounted for in these analyses. Genome-wide data from the UK biobank studied by Fantus et al. (2021) revealed that there are up to 141 various loci across 41 cytobands that may influence variances seen in low T phenotypes. Despite previous studies linking genetic heritability to androgen levels, newer data suggest that genetic variability can account for up to 57% of differences seen across individuals, with multiple genes contributing to the low T phenotype (Fantus et al. 2021).
In another genome-wide study, Ohlsson et al. (2011) observed that genetic variations in the SHBG locus and portions of the X chromosome correlated with the varying levels of T seen in men with established low T levels. Their study identified key polymorphisms that alter binding affinity of SHBG, resulting in lower free T levels. This underlying genetic variability precedes chemical exposure to substances purported to cause lower T levels. As mentioned above, SHBG was the top variable, most-related to T, in our linear and ordinal pairwise comparisons (Suppl., Table S1 and Figure S1) and is another example of the importance of considering sources of inter-individual biological variability, and that risk factors for low T go far beyond overly simplistic correlations.
Studies of genetic variability highlight the limitations of exposure studies that only look to single timepoint measurements of T. Individuals in exposure studies were not evaluated for genetic susceptibility to lower T levels or for associations of genotypes with toxicokinetics and toxicodynamics of molybdenum.
Further, as noted in Cox (2021), men from Black and Hispanic subpopulations were oversampled in the NHANES dataset relative to the US population, so race-specific genetic variability remains a confounder. Meeker et al. (2010) found that reported "non-White" race was also associated with higher T levels (median 402 ng/dL vs. 381 ng/dL) in a Michigan fertility clinic study. This may or may not be linked to the SHBG locus variability and warrants further exploration.

Conclusions
In health and environmental risk assessment, statistically significant exposure-response correlations and regression coefficients alone should not necessarily be interpreted as providing evidence that reducing exposure will reduce risk (the dependent variable). Especially when model-dependent regression findings do not have a biological basis, as shown in this paper, a parallel analysis of plausible biological causality is particularly warranted when evaluating the statistical output. Numerous and varied risk factors and confounders are involved in pinpointing any given "cause" of an apparent or observed lowering of T, and this paper has demonstrated that both statistical and biological viewpoints are needed to confirm one another.
Conclusions differ from multiple linear regression analyses when non-parametric analyses (which are not limited by the assumptions of a multiple linear regression modeling framework) are paired with biology. Unlike the simplistic findings of Lewis and Meeker (2015), this appropriately expanded (biological and nonparametric statistical) reanalysis shows that urinary Mo does not remain as an informative predictor of serum T. Particularly in the light of global downward trends in T over the last several decades, without any accompanying worldwide increase in urinary Mo, the plausibility of Mo as a significant cause of reduced T is not well supported by the weight of biological evidence and recent public health research.
Future research endeavors should include robust reviews of genetic factors, prominent phenotypic profiles and lifestyle adaptations as possible contributors to androgenic health effects which lead to population level decline in T. Future statistical analyses in correlation studies with population-level biomonitoring databases (such as NHANES) would benefit from the robust nonparametric methods used here, carefully paired with biological mode of action understanding, to ensure that statistical model choices do not impede a fuller understanding of health risk factors.