Association of Serum Bile Acids Profile and Pathway Dysregulation With the Risk of Developing Diabetes Among Normoglycemic Chinese Adults: Findings From the 4C Study Free

Better understanding of the early metabolic events underlying the pathogenesis of diabetes may facilitate the development of novel prevention strategies to reduce the incidence of type 2 diabetes mellitus (T2DM). Although multiple metabolic abnormalities, including amino acid, gut microflora metabolites, and lipid metabolism, among others, have been shown to associate with the risk of T2DM (1–3), comprehensive assessment of serum bile acids (BAs) aberrations before T2DM onset remains less extensively studied.

BAs, synthesized from cholesterol in hepatocytes, metabolized by endogenous, as well as intestinal bacteria-derived enzymes, are amphipathic steroid-signaling molecules that coordinately regulate metabolism and inflammation via the nuclear farnesoid X receptor and the Takeda G protein-coupled receptor 5 (TGR5) (4–6). Impaired BA metabolism has been associated with the development and progression of T2DM (7). In healthy subjects, 12α-hydroxylated (12αOH) BAs and ratios of 12αOH to non-12αOH BAs are associated with insulin resistance (8). Glycodeoxycholic acid (GDCA) positively correlates with impaired insulin secretion, while GDCA, glycocholic acid (GCA), and deoxycholic acid (DCA)-glycine conjugate were positively associated with insulin resistance (9). Fasting taurine-conjugated BA concentrations are associated with fasting and postload glucose (10), as seen in cross-sectional studies. However, whether serum/plasma BA profiles could predict the T2DM has been rarely investigated in prospective studies, and the results have been inconclusive (9–12). Using a nontargeted metabolomics approach in the Finnish Diabetes Prevention Study (DPS), researchers found GCA, taurochenodeoxycholic acid (TCDCA), glycochenodeoxycholic acid (GCDCA), GDCA, DCA, and cholic acid (CA) were associated with increased T2DM in individuals with impaired glucose tolerance who developed T2DM (n = 96) or did not convert to T2DM within the 15-year follow-up (n = 104) (12). In contrast, a high-quality meta-analysis comprising four European prospective cohorts reported only a marginal increased risk of diabetes with per SD-unit of DCA increment (hazard ratio [HR] 1.13, 95% CI 1.00–1.27) (11). Using an enzymatic assay, we reported that serum total BAs associate with the risk of incident T2DM in a Chinese population (13).

Most previous studies on human blood (serum/plasma) samples are based on a nontargeted metabolomics approach and do not achieve a satisfactory coverage to render investigation of BA profiles and BA pathway dysregulation in the development and progression of diabetes (9–12). These studies all had limited sample sizes. In the current study, we tested 23 BA species in human serum using ultraperformance liquid chromatography coupled to tandem mass spectrometry (UPLC-MS/MS), which confers both accurate quantitation and extensive coverage of essential BA species.

More importantly, the preceding BA studies were predominantly conducted among mixed populations with normal glucose regulation (NGR) and impaired glucose regulation at baseline (11,12), making it difficult to dissect whether the identified BAs changes are harbingers of early dysglycemia or whether dysglycemia preceded these BAs changes. A comprehensive evaluation of BA profile in individuals with NGR could therefore identify new markers and better understand the possible molecular pathways that define early T2DM pathogenesis.

The China Cardiometabolic Disease and Cancer Cohort (4C) Study, a nationwide, population-based, prospective cohort with up to 5 years of follow-up, provides an opportunity to identify metabolomics profiles associated with T2DM and to discern the metabolomics pathway regulation before T2DM development (14–17). The aims of the current study were to 1) characterize the absolute concentration and composition of circulating BAs in normoglycemic Chinese adults, 2) assess the association of BAs patterns with subsequent risk of incident T2DM, and 3) systematically investigate BAs pathway dysregulation in the prodromal stage of diabetes.

We conducted a nested case-control study of 1,707 matched case subject-control subject pairs within the 4C study. Details of 4C study have previously been described (14–17). The 4C study is a multicenter, population-based, prospective cohort study. During 2011–2012, 193,846 individuals (age ≥40 years) were enrolled from local resident registration systems of 20 communities from various geographic regions in China. During 2014–2016, all participants were invited to attend an in-person visit, and 170,240 participants (87.8%) were followed up.

During a median follow-up of 3.03 (interquartile range 2.87–3.24) years, among 54,807 subjects defined as NGR based on a 75-g oral glucose tolerance test (OGTT) at baseline, 1,864 developed diabetes. A total of 124 serum samples were not available at baseline, thus a final number of 1,707 incident case subjects with diabetes were included in the current study. The control group of 1,707 NGR individuals at baseline was selected using propensity score matching (PSM) (18) with a logistic model that included age, sex, BMI, and fasting plasma glucose (FPG) (Supplementary Fig. 1). The median (interquartile range) follow-up duration was 3.03 (2.88–3.25) years in the case group and 3.03 (2.87–3.24) years in the control group (P = 0.2945).

The study protocol was approved by the Institutional Review Board of Ruijin Hospital affiliated to the Shanghai Jiao-Tong University School of Medicine. Informed consent was obtained from study participants.

Demographic characteristics, lifestyle, dietary factors, and medical history were collected by using standard questionnaires at baseline and the follow-up visit. Physical activity was assessed using the International Physical Activity Questionnaire (19). Moderate and vigorous physical activity was defined as ≥150 min/week of moderate-intensity physical activity, 75 min/week of vigorous aerobic activity, or an equivalent combination of moderate-intensity and vigorous aerobic activities. Healthy dietary score was calculated according to the recommendation of the American Heart Association with replacement of whole grains with bean consumption (20). Anthropometric measurements and blood samples were obtained by trained study nurses following a standard protocol. Blood specimens were processed within 2 h of blood collection at the field center. Sera were shipped by air on dry ice to the central laboratory located at Shanghai Institute of Endocrine and Metabolic Diseases, which is certified by the College of American Pathologists. Aliquots of serum were stored at −80°C in the Ruijin biobank (21). Serum insulin, triglycerides (TGs), total cholesterol (TC), HDL cholesterol, LDL cholesterol, AST, and ALT were measured at the central laboratory by using an ARCHITECT ci16200 autoanalyzer (Abbott Laboratories, Abbott Park, IL).

At both baseline and follow-up visits, all participants underwent an OGTT, and plasma glucose was obtained at 0 and 2 h during the test. Plasma glucose concentrations were analyzed locally by using the glucose oxidase or the hexokinase method within 2 h after collecting the blood sample under a stringent quality control program. Incident diabetes was defined as the following: FPG ≥126 mg/dL, 2-h postload plasma glucose (2h-PG) ≥200 mg/dL, or self-reported previous diagnosis of diabetes by physicians and the current use of antidiabetic medications.

Serum samples collected at baseline were analyzed for BAs using the UPLC-MS/MS with multiple reactions monitoring methods in an ACQUITY UPLC system (Waters) coupled to a triple quadrupole mass spectrometer (Waters), as described previously (22,23). Detailed information on BA measurement is shown in Supplementary Materials 1.0.

BAs were classified into four groups based on conjugation degree (24), including

1. unconjugated primary BAs (PBAs): CA and chenodeoxycholic acid (CDCA);

2. conjugated PBAs: GCA, GCDCA, taurocholic acid (TCA), TCDCA, GCDCA-glucuronid, sulfated glycochenodeoxycholic acid (GCDCS), and sulfated taurochenodeoxycholic acid (TCDCS);

3. unconjugated secondary BAs (SBAs): DCA, ursodeoxycholic acid (UDCA), and lithocholic acid (LCA); and

4. conjugated SBAs: GDCA, GDCA-glucuronide, taurodeoxycholic acid (TDCA), glycolithocholic acid (GLCA), taurolithocholic acid (TLCA), glycoursodeoxycholic acid (GUDCA), tauroursodeoxycholic acid (TUDCA), GDCA, taurodeoxycholic acid (TDCS), sulfated glycolithocholic acid (GLCAS), and sulfated taurolithocholic acid (TLCAS).

Four selected ratios reflective of enzymatic activities in the liver were calculated, including TCA-to-CA, GCA-to-CA, TCDCA-to-CDCA, and GCDCA-to-CDCA. Also calculated were 12αOH BA species (CA, DCA, and their conjugates), and 12αOH–to–non-12αOH ratios were calculated (8). The hydrophobicity index (HI) of the measured BAs was calculated by using the previously published method (25).

Variables with or without normal distribution are presented as mean ± SD or median (interquartile range), respectively. Comparison of the means and medians between groups used the Student t test and paired Wilcoxon rank sum test, respectively. Spearman correlation analysis was used to assess the association of BAs with the main clinical parameters, including BMI, systolic blood pressure (SBP), diastolic blood pressure (DBP), FPG, 2h-PG, HOMA of insulin resistance (HOMA-IR), liver enzymes, and lipids at baseline.

The data of serum BA concentrations were not normally distributed, thus log transformation was applied before analysis using continuous BA variables. To explore the dose-effect relationship, quartile analysis (using the quartile values in control subjects as the cutoff points) was conducted as well. Odds ratios (ORs) and 95% CIs of developing T2DM represented as both quartiles and per SD increase in each BA species were calculated by conditional (matched pairs) logistic regression analysis. We adjusted for potential confounders, including age, sex, BMI, smoking status, alcohol intake, physical activity, education attainment, family history of diabetes, systolic blood pressure, FPG, TGs, and TC. The P value was corrected for multiple testing via false discovery rate (FDR) using the Benjamini-Hochberg method. In the sensitivity analyses, we analyzed the effects of adjustment for diet habits, liver enzymes, including ALT and AST, 2h-PG, and HOMA-IR. In ad hoc analysis, baseline metabolic traits were better matched (see Supplementary Materials 2.1).

The R package MEGENA was used to build correlation network from differentially correlated BA pairs. Differential correlation was calculated using the R package DGCA. Only BA pairs with differential correlation (P < 0.05) were included for analyses (see Supplementary Materials 2.2 for details).

All statistical comparisons were two sided, and a P value < 0.05 was considered statistically significant. Statistical analyses were performed using SAS 9.3 software (SAS Institute, Cary, NC) and R 3.3.1 software.

In addition to age, sex, BMI, and FPG matched under PSM, other baseline characteristics, including smoking status, alcohol intake, physical activity, diet habit, and education status, were also well-matched between case subjects and control subjects. Case subjects showed higher levels of 2h-PG, TGs, TC, AST, ALT, and HOMA-IR than control subjects. No difference was observed for LDL cholesterol between the two groups (Table 1).

PBAs comprised 65.8% of the serum total BAs in the control subjects, the greatest contributors being GCDCA and CDCA, whereas TDCS, TUDCA, and TLCA were present in relatively low concentrations (Fig. 1A). Case subjects had higher level of GCDCA, GCA, GUDCA, GCDCAS, TCDCA, TCA, TUDCA than control subjects (all P < 0.05), whereas the levels of DCA, CA, and GCDCA-glucuronide were lower (all P < 0.05) (Fig. 1A). The composition of individual BAs was statistically different in case subjects compared with control subjects. Case subjects had a higher percentage of conjugated PBAs (P < 0.001) (Fig. 1B). Serum 12αΟΗ BA species did not differ between case subjects and control subjects (P = 0.9090). Serum total taurine- and glycine-conjugated BAs were significantly higher in case subjects compared with control subjects, with the median (interquartile range) of 136.81 (74.39–252.30) vs. 123.05 (66.84–221.85) nmol/L (P = 0.0006), and 2,073.81 (1,279.31–3,575.20) vs. 1,863.79 (1,154.35–3,158.72) nmol/L (P < 0.0001), respectively.

Spearman correlation analysis revealed that individual and subgroup of baseline serum BAs composition was correlated with a range of biochemical measurements and metabolic parameters in normoglycemic individuals (Fig. 1C). At baseline, the sum of all measured serum BAs was positively associated with ALT, AST, FPG, SBP, DBP, and waist-to-hip ratio (WHR), and negatively associated with 2h-PG and LDL. The 12αOH–to–non-12αOH BA ratio was negatively associated with BMI, WHR, body weight, waist, hip, TGs, DBP, SBP, AST, and ALT, but not associated with FPG, 2h-PG, and lipids. Additionally, TCA-to-CA and GCA-to-CA ratios were inversely associated with FPG and SBP. TCDCA-to-CDCA and GCDCA-to-CDCA ratios were negatively associated with BMI, TGs, and SBP. The TCA-to-CA, GCA-to-CA, and TCDCA-to-CDCA ratios correlated with HOMA-IR and 2h-PG positively.

The associations between per SD difference of individual and grouped BAs and the risk of incident T2DM are shown in Table 2. In the multivariable model with adjustment of age, sex, BMI, smoking status, alcohol intake, physical activity, education status, diabetes family history, SBP, FPG, TGs and TC, we observed 7 of 9 PBAs and 2 of 14 SBAs were associated with incident diabetes (P < 0.05). Per SD difference of unconjugated PBAs (CA and CDCA) and SBA (DCA) were inversely associated with diabetes risk. The OR (95% CI) was 0.89 (0.83–0.96) for CA, 0.90 (0.84–0.97) for CDCA, and 0.90 (0.83–0.96) for DCA. By contrast, conjugated PBAs (GCA, TCA, GCDCA, TCDCA, and GCDCS) and SBA (TUDCA) were positively associated with T2DM. The ORs (95% CIs) per SD difference were GCA, 1.18 (1.08–1.26); TCA, 1.19 (1.10–1.28); GCDCA, 1.13 (1.05–1.21); TCDCA, 1.18 (1.09–1.26); GCDCS, 1.11 (1.03–1.20); and TUDCA, 1.13 (1.05–1.21). The OR (95% CI) per SD difference was 1.14 (1.06–1.23) for conjugated PBAs and 1.05 (0.97–1.12) for conjugated SBAs, respectively. Total taurine- and glycine-conjugated BAs were positively associated with diabetes risk. The OR (95% CI) per SD difference was 1.16 (1.08–1.25) and 1.14 (1.06–1.22), respectively. Four selected ratios of primary conjugated to unconjugated BAs were positively associated with T2DM. The ORs (95% CIs) per SD difference were TCA-to-CA, 1.19 (1.09–1.31); GCA-to-CA, 1.16 (1.05–1.29); TCDCA-to-CDCA, 1.22 (1.13–1.31); and GCDCA-to-CDCA, 1.19 (1.11–1.28).

The above-mentioned BAs remained statistically significant after adjusting for multiple comparisons (FDR <0.05). No significant association was detected between 12αOH–to–/non-12αOH BA ratios, HI, and the risk of incident T2DM. Sensitivity analyses with adjustment for additional baseline risk factors, including diet score, liver enzyme levels, 2h-PG, and HOMA-IR, showed the similar results. In the fully adjusted model including all the confounding factors, including liver enzymes, HDL, 2h-PG, HOMA-IR, and waist circumference, the ORs did not show significant change (Table 2). Likewise, adjusting for baseline HbA_1c did not substantially change estimates (Supplementary Table 1).

The risk for developing T2DM across BA quartiles is shown in Supplementary Table 2. In multivariable adjusted logistic regression model, 7 of 9 PBAs and 2 of 14 SBAs were associated with diabetes risk (P for trend < 0.05). Using the lowest quartile as reference, the OR (95% CIs) for T2DM in the highest quartile of CA and CDCA were 0.74 (0.60–0.90) and 0.76 (0.62–0.93), respectively, whereas the OR (95% CI) in the highest was 1.50 (1.23–1.84) for GCA, 1.34 (1.09–1.64) for GCDCA, 1.61 (1.31–1.97) for TCA, 1.51 (1.23–1.84) for TCDCA, 1.35 (1.10–1.67) for GCDCS, and 1.44 (1.18–1.77) for TUDCA. The TCA-to-CA, GCA-to-CA, TCDCA-to-CDCA, and GCDCA-to-CDCA ratios were positively associated with the risk of T2DM (P for trend <0.0001). In the fully adjusted model including liver enzymes, HDL, 2h-PG, HOMA-IR, and waist circumference, CA, GCA, TCA, TCDCA, and TUDCA, and four ratios were associated with the risk of T2DM (P for trend <0.05).

Results from ad hoc analysis on those matched for baseline metabolic parameters, including blood pressure, serum lipids, liver enzymes, 2h-PG, and HOMA-IR, are shown in Supplementary Tables 3 and 4. The risk estimates were comparable to those in the fully adjusted model.

To evaluate for perturbed BA coregulation underlying incident diabetes, we conducted multiscale embedded correlation analysis (26) to visualize differential correlation between various BAs in T2DM relative to control subjects. Perturbations in BA coregulation that preexisted before diabetes onset were identified (Fig. 2). First, CA was positively correlated with GCA, TCA, and GLCAS in control subjects but not in subjects with incident diabetes (brown line, 0/++), and the positive relationship between CA and GCDCA was stronger in control subjects than in case subjects (dark blue line, +/++). Second, CDCA was positively correlated with LCA in control subjects, but not in case subjects (red line, 0/+), whereas CDCA became negatively associated with GLCAS in case subjects (green line, –/0). Third, LCA and most of the conjugated PBAs (green dots) or SBAs (blue dots) were connected by the pink lines (++/0), indicating that LCA became positively associated with the conjugated BAs in case subjects. In addition, the positive correlation between unconjugated SBAs (LCA and UDCA) became negatively significant in case subjects (light blue, −/++). Taken together, differential correlation networks revealed that even before diabetes onset, perturbations in intraclass (i.e., primary and secondary) and interclass (i.e., unconjugated and conjugated) BA coregulation have already been altered (Fig. 2).

In the current study, using a UPLC-MS/MS approach, we evaluated circulating BAs profiles in a nested case-control study drawn from a nationwide, population-based, prospective cohort of normoglycemic Chinese adults. This is the largest and most comprehensive BAs study investigating the association between serum BAs profiles and risk of developing diabetes in individuals with NGR at baseline. Unconjugated PBAs are negatively associated with incident T2DM. By contrast, conjugated PBAs are positively associated with incident T2DM. Moreover, differential correlation network analysis revealed that perturbations in intraclass (i.e., primary and secondary) and interclass (i.e., unconjugated and conjugated) BAs coregulation preexisted before diabetes onset.

Our study contributes a systematic evaluation of serum BAs profile changes predictive of incident diabetes in individuals with NGR. We tested and identified far more BAs associated with risk of T2DM than other prospective studies, including DPS (12) and the Swedish meta-analysis (11) (Supplementary Table 5). Our results confirmed the previous observations that GCA, GDCA, TCDCA, and GCDCA were associated with increased risk of diabetes in European populations (11,12). DCA was, however, found to be associated with increased risk of diabetes in the DPS cohort (12) and the Swedish meta-analysis (11), which was inconsistent with the finding of our study. This discrepancy might be partially explained by the differences in the sample size, study design (including only NGR at baseline in the current study vs. both NGR and DM in Swedish meta-analysis study vs. impaired glucose regulation in the DPS cohort), and other cohort epidemiological features, including ethnicity-related genetic and lifestyle variations.

Interestingly, our study revealed for the first time that all unconjugated PBAs (CA and CDCA) identified were associated with decreased risk of diabetes. These results are consistent with the findings in experimental studies. CA has been shown to have antidiabetic effects by increasing production of insulin in pancreatic β-cells in vivo (27) and in vitro (28). CDCA treatment also induces beneficial effects on glucose metabolism by upregulating farnesoid X receptor expression in liver tissue and increasing GLP-1 and glucagon secretion (29,30). By contrast, we noticed that most of the conjugated PBAs (GCA, GCDCA, TCA, TCDCA, and GCDCS) identified were related with increased T2DM risk. Moreover, ratios of conjugated to unconjugated PBAs (TCA-to-CA, GCA-to-CA, TCDCA-to-CDCA, and GCDCA-to-CDCA) were positively associated with increased T2DM risk, suggesting that the reduced production of unconjugated PBAs and excess transition from unconjugated PBAs to the conjugated form may play essential roles in the pathogenesis of T2DM development. These observations point to a possibility of altered activity of enzymes related to synthesis and metabolism of BAs in incident diabetes. The conjugation of the unconjugated PBAs to taurine or glycine is catalyzed by BACS and BAAT (24). Data from the European Prospective Investigation into Cancer and Nutrition (EPIC) study indicated that the variants in gene of BACS may be etiological for T2DM (31). Our study provides supporting evidence that BACS and BAAT may play critical roles in maintaining whole-body glucose homeostasis. Further genetic analyses were needed to elucidate the precise mechanisms.

Furthermore, our correlation network analyses revealed the association between serum individual BAs and the perturbed coregulation between conjugated BAs versus unconjugated BAs, and PBAs versus SBAs, in subjects with T2DM relative to control subjects, which confer useful BAs-centric insights on the etiology of diabetes. Owing to the diversity of structure and hydrophilic/hydrophobic properties in PBAs versus SBAs and in conjugated BAs versus unconjugated BAs, they exert distinct biological effects on T2DM (24,29,30,32,33). Our observations therefore point to a possibility that altered BA coregulation in incident diabetes might account for their positive correlation with diabetes risk, which could subsequently skew the compositional profiles of serum unconjugated BAs toward conjugated forms. The precise causal links between these observations await further mechanistic elucidation, and whether increases in conjugated PBAs and altered composition of the BAs pool denote cumulating or compensatory events of diabetes onset remains an interesting open question.

The strengths of our study include the large sample size, prospective study design, and hence, the ability to study the association with T2DM, the inclusion of populations from 20 communities from various geographic regions in China, population being normoglycemic at baseline, adjustment for many potential confounders, and a series of sensitivity analyses. In the sensitivity analysis with additional adjustment of 2h-PG, we provided strong evidence that BA species can predict incident diabetes independent of baseline glucose. In terms of the BAs analytical approaches, our work extensively covers 23 key BA classes that are essential in the overall homeostatic balance of BA metabolism.

However, several limitations need to be addressed. First, the 3-year follow-up duration limits the predictive potential of our identified BA panel to a relatively narrow time window. Further studies with longer follow-up duration might be additionally informative. The diagnosis of diabetes, based on the follow-up OGTT was not repeated/confirmed. This could limit the accuracy of the diagnoses of diabetes.

Second, serum BAs might be indicative of both dietary and metabolic influences, affected by interplay of complex exogenous and genetic factors, which we were unable to tease out because detailed diet information and gene analysis from participants is lacking.

Third, BAs were measured only on baseline fasting samples, and data on postprandial BAs were not available, which hampers further analysis about the role of BAs in glucose metabolism. It would also be interesting to see whether BA composition continued to change with incident disease.

Fourth, the PSM-matched case-control design hampers the use of matching factors (age, sex, BMI, FPG) for risk prediction in the case-control samples. Although the current study was conducted in normoglycemic individuals at baseline, substantial metabolic differences were observed between case subjects and control subjects. In network analysis, the potential impact of the metabolic risk factors might exist. Thus, we cannot exclude that these changes may be secondary to the ongoing progressive disturbances in this population. Nevertheless, in the ad hoc analysis in a better matched population, the results were similar.

Fifth, gut flora was not assessed in the current study. Given that the BA changes have marked primary and secondary distributions, this might be an important element to fully understand the metabolic pathways in the future.

Finally, all of the participants in this study were Chinese, and further work is needed to determine whether our findings can be extrapolated to other races and ethnicities.

In conclusion, our findings indicate that serum BAs are differentially associated with risk of T2DM, which supports the importance of recognizing that unconjugated and conjugated BAs exert heterogeneous effects and the potential value of BA metabolism early in the pathogenesis of diabetes. Further studies are warranted to test whether serum BA measurements might help identify candidates for interventions to reduce diabetes risk and to elucidate the biological mechanisms of BAs in the onset and progression of T2DM. These studies may provide insights for the diagnosis as well as novel therapeutic targets for treatment of T2DM.

Acknowledgments. The authors would like to thank the study participants and investigators of the 4C study in China States, without whom this research would not be possible.

Funding. This work was supported by the Ministry of Science and Technology of China under award numbers 2016YFC1305600, 2016YFC1305202, 2016YFC1304904, 2016YFC0901200, 2017ZX09304007, and 2018YFC1311800, the National Natural Science Foundation of China under award numbers 81930021, 81970728, 81970691, 81670795, 81730023, and 81621061, Shanghai Outstanding Academic Leaders Plan under award numbers 18XD1402500 and 20XD1422800, Clinical Research Plan of Shanghai Hospital Development Center under award number SHDC2020CR3064B, Shanghai Science and Technology Committee under award number 20Y11905100, Shanghai Municipal Natural Science Foundation under award number 18ZR1433100, and Shanghai Medical and Health Development Foundation.

Duality of Interest. No potential conflicts of interest relevant to this article were reported.

Author Contributions. J.L., S.W., M.L., C.H., and R.D. did the statistical analysis. J.L., M.L., and Y.B. drafted the manuscript. J.L., G.N., Y.B., and W.W. conceived and designed the study. G.N. supervised the study. Z.G., Y.X., X.Z., Yi.Z., R.L., R.H., L.S., R.Z., Q.S., J.W., Y.C., X.Y., L.Y., T.W., Z.Z., X.W., Qi.L., G.Q., Q.W., G.C., M.X., M.D., D.Z., X.T., G.W., F.S., Z.L., Y.Q., Li.C., Y.H., Qia.L., Z.Y., Yin.Z., C.L., Y.W., S.W., T.Y., H.D., D.L., S.L., Y.M., Lu.C., J.Z., and G.X. contributed to acquisition, analysis, or interpretation of data. All authors revised the report and approved the final version before submission. J.L., G.N., Y.B., and W.W. are the guarantors of this work and, as such, had full access to all the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis.