Study design and participants

A previously published medical center-based case–control study was conducted from 2011 to 2019 in central Taiwan (Chang et al. 2016; Chung et al. 2014). Every UC case was diagnosed with histological confirmation and further to distinguish the tumor location. The controls participants were recruited randomly from receiving health checks at the Department of Family Medicine covering nearly all the population living in the corresponding area. All of the control participants had no UC history. After excluding no collections of blood samples, incomplete data of metals in blood as well as LINE-1 methylation, 478 UC patients (266 urinary bladder; UB as well as 204 upper tract UC; UTUC) and 569 controls were analyzed in the present study. This study was approved by the Research Ethics Committee of China Medical University Hospital, Taichung, Taiwan. Also, all study participants signed an informed consent before recruitment.

Clinical information on participants and biochemistry examinations

Data about demographic characteristics and clinical variables, such as cigarette smoking, education, and receiving Chinese herbal medicine as well as individual history were collected through face-to-face interviews with a structured questionnaire and medical records. Meanwhile, general biochemistry examinations were examined, including baseline anthropometric, blood pressure, and plasma levels of triglycerides, total cholesterol, low-density lipoprotein cholesterol (LDL-C), high-density lipoprotein cholesterol (HDL), fasting plasma glucose, insulin, uric acid, and blood creatinine after an 8-h fasting period. Homeostatic Model Assessment for Insulin Resistance (HOMA-IR) was used to measure insulin resistance (IR) and calculated as fasting plasma glucose (mg/dL) × insulin (mU/L)/405 (Matthews et al. 1985). Estimated glomerular filtration rate (eGFR) was based on a single moderate or established kidney function of eGFR level < 60 mL/min/1.73 m2, calculated using the Chronic Kidney Disease Epidemiology Collaboration (CKD-EPI) equation (Ebong et al. 2013). We defined smoking status as nonsmokers, former smokers who had quit smoking at the time of recruitment, and current smokers. The status of clinical comorbidities were verified by clinicians and also based on medical history. All comorbidities included hypertension (yes, no), diabetes mellitus (yes, no), urolithiasis (yes, no), and chronic kidney diseases (CKD) (yes, no). Study subjects with HDL-C ≤ 40 or LDL_C > 190 or triglycerides > 500 were categorized as hyperlipidemia (Expert Panel on Detection and Treatment of High Blood Cholesterol in 2001; Kreisberg and Oberman 2003). Cutoff value ≥ 0.3 mg/dL of high-sensitivity C-reactive protein (hs-CRP) levels was dichotomized into low and high groups (Simanek et al. 2011).

Pyrosequencing to measure the LINE-1 methylation

The blood samples were centrifuged at 3000 g for 10 min to separate the buffy coat and the genomic DNA was extracted from buffy coat. The detailed protocol was published in a previous study (Chung et al. 2017). In addition, genomic DNA from tumor tissues was extracted using the commercially available the ChargeSwitch® gDNA Micro Tissue Kit (Invitrogen). Approximately 500 ng of leucocyte and tissue DNA were sodium bisulfite modified using the Epitect Plus DNA Bisulfite Kit (Qiagen) according to manufacturer’s protocol. The bisulfite modification could convert unmethylated cytosines to uracil and each batch performed with positive methylated controls (EpiTect Control DNA, Qiagen, Lot No. 151022455) and unmethylated controls (EpiTect Control DNA, Qiagen, Lot No. 151022699). The converted DNA containing three CpG sites for LINE-1 (positions 819, 826, and 829) was then amplified using the polymerase chain reaction with forward primer: 5′-TTTTGAGTTAGGTGTGGGATATA-3′, reverse primer 5’-Biotin-AAAATCAAAAAATTCCCTTTC- 3′. Total volume of 25 μL contained PyroMark PCR Master Mix, 1.5 mM MgCl2, Q solution, 0.2 μM each primer, RNase-free water and 1 μL of converted DNA with the following cycling profile: 95 °C for 15 min followed by 40 cycles of denaturation at 95 °C for 30 s, annealing at 58 °C for 1 min, extension at 72 °C for 1 min and final extension at 72 °C for 10 min (146 bp). Total volume of 80 μL contained 15 μL PCR product, binding buffer, RNase-free water and streptavidin-coated beads were prepared to measure the methylation levels of LINE-1 through pyrosequencing reaction with sequencing primer 5′-AGTTAGGTGTGGGATATAGT-3′ (PyroMark Q24 platform, Qiagen). Methylation quantification was calculated for individual CpG site of LINE-1 as percentage of methylated cytosines divided by the sum of methylated and unmethylated cytosines in the automatic software. Finally, the average values (percentage) of LINE-1 methylation from above 3 CpG sites were acquired (Supplementary Fig. 1). To ensure quality control, a random 5% of the samples were confirmed repeatedly.

Measurement of heavy metals in blood

About 5–6 mL of blood was collected from each participant during recruitment for measurement of heavy metals. The detailed protocol for determination of heavy metals of including As, Cd, Cr, Ni, and Pb in blood was previously described (Yang et al. 2018b). In brief, multiple levels of the above metals in whole blood digested with nitric acid were quantified using inductively coupled plasma-mass spectrometry (ICP-MS; Agilent 7700c, Agilent Technologies, Inc., Palo Alto, CA, USA). The correlation coefficients for fitted calibration curves were ≥ 0.99, and each recovery rate was in the range of 85–115%. The standard reference materials (Seronorm™ Trace Elements Whole Blood) of blood heavy metals were used for quality control of the measurements. We divided the values by the square root of two, when the levels of blood heavy metals were less than the detection limit values: 0.03 ppb for As, 0.007 ppb for Cd, as well as 0.067 ppb for Pb. For experimental reliability, the repeated samples were randomly selected and analyzed. The coefficients of variation were within the range of 5–10% for each metal.

Statistical analysis

All analyzed data, including mean and standard deviation (SD) for continuous variables as well as numbers and frequencies for categorical variables, were compared for the subgroups of UC, UB, UTUC group, and control group. Mann–Whitney U tests or Chi-square tests were used to compare the differences of the demographic characteristics, UC-related risk factors and clinical biomarkers among comparison groups. Univariate and multiple logistic regression models were used to evaluate the associations of the levels of five metals in blood (As, Cd, Cr, Ni, and Pb) as well as LINE-1 DNA methylation with the risks of UC.

The multicollinearity among five metals in blood was managed through least absolute selection and shrinkage operator (LASSO) regression, which performed both variable selection and regularization by shrinking less significant coefficients towards zero. Through LASSO regression analysis, we separately examined UB and UTUC patients to determine which heavy metals had the most significant impact on individual subtypes of UC. Variable importance (VIMP) in Random forest was chosen for exploring the explainability of UB or UTUC. The methodology of VIMPs determines variable importance through the assessment of increased out-of-bag prediction error following the random permutation of specific predictor variables (Breiman 2001). In this study, we employed two measures of variable importance: Mean Decrease Accuracy and Mean Decrease Gini. The Mean Decrease Accuracy quantifies the decrease in model accuracy when a variable is randomly permuted, while the Mean Decrease Gini measures the average decrease in node impurity across all trees when a given variable is used for splitting (Efron 2020; Epifanio 2017).

Finally, we performed a causal mediation analysis to evaluate the direct effect of every kind of exposure variables on UC risk independent of a mediator (LINE-1 DNA methylation) (VanderWeele). The roles that the mediators played in the indirect effect as well as the percentages mediated (adjusted for confounding variables) were also calculated. Also, a proportion-mediated percentage was calculated from the direct (DE) and indirect effect (IE) odds ratios using the formula ORDE(ORIE − 1)/(ORDEORIE − 1)(VanderWeele).

All data were analyzed using the SAS statistical package (SAS, version 9.4, Cary, NC) or R software version 3.6.3 (The R Foundation for Statistical Computing; Vienna, Austria). A two-sided p value < 0.05 was considered statistically significant.