Pleiotropic Effects of Functional MUC1 Variants on Cardiometabolic, Renal, and Hematological Traits in the Taiwanese Population

MUC1 is a transmembrane mucin involved in carcinogenesis and cell signaling. Functional MUC1 variants are associated with multiple metabolic and biochemical traits. This study investigated the association of functional MUC1 variants with MUC1 DNA methylation and various metabolic, biochemical, and hematological parameters. In total, 80,728 participants from the Taiwan Biobank were enrolled for association analysis using functional MUC1 variants and a nearby gene regional plot association study. A subgroup of 1686 participants was recruited for MUC1 DNA methylation analysis. After Bonferroni correction, we found that two MUC1 variants, rs4072037 and rs12411216, were significantly associated with waist circumference, systolic blood pressure, hemoglobin A1C, renal functional parameters (blood urea nitrogen, serum creatinine levels, and estimated glomerular filtration rate), albuminuria, hematocrit, hemoglobin, red blood cell count, serum uric acid level, and gout risk, with both favorable and unfavorable effects. Causal inference analysis revealed that the association between the variants and gout was partially dependent on the serum uric acid level. Both gene variants showed genome-wide significant associations with MUC1 gene-body methylation. Regional plot association analysis further revealed lead single-nucleotide polymorphisms situated at the nearby TRIM46–MUC1–THBS3–MTX1 gene region for the studied phenotypes. In conclusion, our data demonstrated the pleiotropic effects of MUC1 variants with novel associations for gout, red blood cell parameters, and MUC1 DNA methylation. These results provide further evidence in understanding the critical role of TRIM46–MUC1–THBS3–MTX1 gene region variants in the pathogenesis of cardiometabolic, renal, and hematological disorders.


Introduction
Mucin family members are classified into two types on the basis of their localization: secreted or membranous [1][2][3][4]. Mucin 1 (MUC1), also named polymorphic epithelial mucin and encoded as MUC1 in humans, is a single-pass type I transmembrane phosphoprotein with substantial O-linked glycosylation. In terms of its functions, MUC1 was initially assumed to only protect from the external environment and lubricate the epithelium, whereas MUC1 is currently considered to play a crucial role as a multifunctional protein in cell signal transduction and cell-cell interaction [1,3]. MUC1 is normally expressed in the glandular or luminal epithelial cells of the mammary gland, esophagus, stomach, duodenum, pancreas, uterus, prostate, and lungs, and to a lesser extent, in hematopoietic cells [5][6][7]. Furthermore, MUC1 is overexpressed in breast, ovarian, lung, pancreatic, and prostate cancers and is a marker of poor prognosis in gastric cancer [8,9]. A recent proteomic analysis of urine revealed that the urinary excretion of five proteins including MUC1 is associated with the risk of renal impairment in the general population. In addition, this study indicated that the diagnostic value of urinary MUC1 in predicting a decline in the estimated glomerular filtration rate (eGFR) was higher than that of microalbuminuria [10].
MUC1, located on chromosome 1q22, is a causative gene for autosomal dominant tubulointerstitial kidney disease (ADTKD), which is a hereditary disease characterized by progressive tubulointerstitial nephropathy that ultimately leads to end-stage renal disease (ESRD) [11][12][13]. Several genome-wide association studies (GWASs) have revealed a significant association of MUC1 polymorphisms with the risk of gastric cancer [14][15][16]. In addition, recent GWASs reported associations of two MUC1 single-nucleotide polymorphisms (SNPs), namely the functional variants rs4072037 and rs12411216, with serum magnesium and uric acid levels, blood pressure, impaired renal function, and albuminuria [17][18][19][20][21]. The Taiwan Biobank (TWB) conducted a large-scale population-based cohort study by recruiting volunteers aged between 30 and 70 years with no history of cancer [22,23]. In the present study, we hypothesized that both genetic and epigenetic effects of the MUC1 gene may play a crucial role in the pathogenesis of cardiometabolic, renal, and hematological disorders. We investigated the association of functional MUC1 variants with various metabolic, biochemical, and hematological parameters in more than 80,000 TWB participants and with MUC1 DNA methylation in 1686 TWB participants. In addition, we examined whether causal inference occurred between the associated phenotypes, identified variants that played more crucial roles in determining the association with various phenotypes, and evaluated the role of nearby gene loci variants by performing a regional plot association study. Table 1 lists the demographic data, clinical and biochemical data, lipid profiles, and hematological traits of TWB participants stratified by sex. Compared with women, men had a significantly higher BMI, waist circumference, waist to hip ratio, systolic, diastolic, and mean blood pressure, fasting plasma glucose, serum triglyceride, uric acid, creatinine, BUN, AST, ALT, γ-GT, albumin, total bilirubin levels, hematocrit, red blood cell count, and hemoglobin (all p < 0.0001). By contrast, HDL-cholesterol, total cholesterol levels, eGFR, and platelet counts (all p < 0.0001) were lower in men than in women.

Association of MUC1 Genotypes with Clinical, Metabolic, and Biochemical Phenotypes and Hematological Parameters
A total of over 80,000 volunteers participated in the genotype-phenotype association analysis. We analyzed the association of MUC1 rs12411216 and rs4072037 genotypes with various clinical phenotypes and laboratory parameters. By using an additive model, after adjustment for age, sex, BMI, and smoking status and after Bonferroni correction, we observed a significant association of both MUC1 variants with waist circumference, systolic blood pressure, HbA1C, renal functional parameters (BUN, serum creatinine levels, and eGFR), albuminuria, red blood cell parameters (hematocrit, hemoglobin, and red blood cell count), serum ALT and uric acid levels, and gout (Table 2). A near complete linkage between the rs12411216 and rs4072037 polymorphisms was found with an r 2 of 0.9832.  Abbreviations, adjusted conditions, and subjects recruited for analysis as in Table 1. Adjusted p value: with Bonferroni correction, n = 28. * were analyzed with the exclusion of participants with previous history of hypertension; ** were analyzed with the exclusion of participants with previous history of diabetes mellitus; *** were analyzed with the exclusion of participants with previous history of gout; **** were analyzed with the exclusion of participants with previous history of hyperlipidemia.

Association of MUC1 Genotypes with Risk Factors for Atherosclerosis
We examined the association of MUC1 rs12411216 and rs4072037 genotypes with risk factors for atherosclerosis. The results revealed that gout and microalbuminuria were significantly associated with MUC1 gene variants ( Figure 1 and Table S1).

Subgroup Analysis of the Association of MUC1 Genotypes with Studied Parameters and Atherosclerotic Risk Factors
We investigated whether sex affects the association between the two variants and the studied phenotypes (Tables S2 and S3). The results of the t test revealed that sex affected only the association between the MUC1 variants and albuminuria (t = −2.6977 and p < 0.01 for the rs12411216 genotype and t = −2.7698 and p < 0.01 for the rs4072037 genotype). Both the variants showed a more significant association with albuminuria in men.  . p values were adjusted for age, sex, body mass index, and current smoking. Significantly higher frequencies of gout and lower frequencies of microalbuminuria were noted for A allele of rs407203 and C allele of rs12411216.

Subgroup Analysis of the Association of MUC1 Genotypes with Studied Parameters and Atherosclerotic Risk Factors
We investigated whether sex affects the association between the two variants and the studied phenotypes (Tables S2 and S3). The results of the t test revealed that sex affected only the association between the MUC1 variants and albuminuria (t = −2.6977 and p < 0.01 for the rs12411216 genotype and t = −2.7698 and p < 0.01 for the rs4072037 genotype). Both the variants showed a more significant association with albuminuria in men.

Mediational Analysis: MUC1 Genotypes and Various Phenotypes
The two-stage least square instrumental variable (IV) analysis was performed to examine the causality between the studied phenotypes such as the eGFR, hematocrit, uric acid level, gout, and albuminuria with rs4072037 and rs12411216 genotypes used as IVs (Tables S4 and S5). Our data revealed an association between the MUC1 variants and the studied phenotypes, and this association remained highly significant after adjustment for other phenotypes, suggesting no evidence of a causal relationship between these study phenotypes. The only exception was that the association between each variant and gout disappeared after adjustment for the serum uric acid level (p value decreased from 2.94 × 10 −8 to 0.0004 and from 3.66 × 10 −8 to 0.0006 for rs4072037 and rs12411216 genotypes, respectively). These results suggested that the association between both variants and gout would be at least partially mediated by the serum uric acid level.

Association between the Functional Polymorphisms of the MUC1 Gene and the Nearby MUC1 DNA Methylation Status
We further examined whether the two MUC1 functional variants are associated with the DNA methylation status of the MUC1 gene. We found no significant association between the two variants and 12 MUC1 promoter methylation sites (Table S6). By contrast, genome-wide significant associations were observed between both variants and four MUC1 gene-body DNA methylation sites, namely cg06339768, cg19011149, cg15646096, and cg00126087, with a stronger association noted for the rs4072037 genotype, as shown in Figure 2 (p = 1.01 × 10 −28 , 4.23 × 10 −66 , 1.34 × 10 −121 , and 5.84 × 10 −27 , respectively, after adjustment for age, sex, BMI, and smoking with Bonferroni correction). . p values were adjusted for age, sex, body mass index, and current smoking. Significantly higher frequencies of gout and lower frequencies of microalbuminuria were noted for A allele of rs407203 and C allele of rs12411216.

Mediational Analysis: MUC1 Genotypes and Various Phenotypes
The two-stage least square instrumental variable (IV) analysis was performed to examine the causality between the studied phenotypes such as the eGFR, hematocrit, uric acid level, gout, and albuminuria with rs4072037 and rs12411216 genotypes used as IVs (Tables S4 and S5). Our data revealed an association between the MUC1 variants and the studied phenotypes, and this association remained highly significant after adjustment for other phenotypes, suggesting no evidence of a causal relationship between these study phenotypes. The only exception was that the association between each variant and gout disappeared after adjustment for the serum uric acid level (p value decreased from 2.94 × 10 −8 to 0.0004 and from 3.66 × 10 −8 to 0.0006 for rs4072037 and rs12411216 genotypes, respectively). These results suggested that the association between both variants and gout would be at least partially mediated by the serum uric acid level.

Association between the Functional Polymorphisms of the MUC1 Gene and the Nearby MUC1 DNA Methylation Status
We further examined whether the two MUC1 functional variants are associated with the DNA methylation status of the MUC1 gene. We found no significant association between the two variants and 12 MUC1 promoter methylation sites (Table S6). By contrast, genome-wide significant associations were observed between both variants and four MUC1 gene-body DNA methylation sites, namely cg06339768, cg19011149, cg15646096, and cg00126087, with a stronger association noted for the rs4072037 genotype, as shown in Figure 2 (p = 1.01 × 10 −28 , 4.23 × 10 −66 , 1.34 × 10 −121 , and 5.84 × 10 −27 , respectively, after adjustment for age, sex, BMI, and smoking with Bonferroni correction).
The C allele of rs12411216 and the A allele of rs4072037 were significantly associated with the hypermethylated status of the MUC1 gene body ( Figure 3A-H and Table S6). The C allele of rs12411216 and the A allele of rs4072037 were significantly associated with the hypermethylated status of the MUC1 gene body ( Figure 3A-H and Table S6).    Furthermore, we analyzed the association between cg15646096 and the studied clinical and laboratory parameters and found that the methylation status of cg15646096 significantly decreased with age (Tables S7 and S8). We performed a regional plot analysis to examine the association between genetic variants from position 155.0 to 155.5 Mb on chromosome 1q22 and the studied phenotypes. Genome-wide significant associations for the lead SNPs were noted for uric acid level (rs2070803, p = 3.63 , and red blood cell count (rs76872124, p = 1.18 × 10 −9 ). Our data revealed that the lead SNP for each phenotype was situated at the TRIM46-MUC1-THBS3-MTX1 gene region ( Figure 4A, Figures S1-S12, and Table 3). After conditional analysis, we observed a genome-wide significant association of only TRIM5 rs76872124 with BUN after adjustment for the lead SNP rs4971101 ( Figure S6). Most of the lead SNPs showed strong LD with rs4072037 and rs12411216 (with r 2 > 0.87), whereas the lead SNPs for red blood cell parameters, such as rs4971093 and rs76872124, exhibited a moderate LD with rs4072037 (r 2 = 0.705 and r 2 = 0.544, respectively). The results are summarized in Table 3 and Figure 4B.

Discussion
To the best of our knowledge, this is the first study to investigate the association of MUC1 functional SNPs with multiple clinical and laboratory parameters by using both candidate gene and regional plot association approaches. We observed that both MUC1 variants were associated with renal functional parameters (BUN, serum creatinine level, and eGFR), albuminuria, and serum uric acid levels; this finding is in accordance with that of previous studies. Furthermore, we identified novel genome-wide significant associations between MUC1 functional SNPs and gout, red blood cell parameters (hematocrit, hemoglobin, and red blood cell counts), and MUC1 gene-body DNA methylation. Our results suggested that both SNPs play a critical role in determining different phenotypes. We also provide the first evidence of the possibility that the previously reported association between rs4072037 and the total or alternative spliced form of gene expression may be secondary to MUC1 gene-body methylation status. The findings of causal inference analysis showed that the association of both variants with gout was partially dependent on the serum uric acid level. No cause-effect relationship was found in other phenotypes, indicating mostly pleiotropic effects. The results of regional plot association analysis revealed lead SNPs at the nearby TRIM46-MUC1-THBS3-MTX1 gene region for the studied phenotypes, providing further evidence for the critical role of MUC1 as a candidate gene locus for cardiometabolic, renal, and hematological disorders.

Association between MUC1 Polymorphisms and Renal Function and Albuminuria
In a GWAS analyzing molecular insights into chronic kidney disease-defining traits, Xu et al. [20] reported the causal effects of three genes (NAT8B, CASP9, and MUC1) on the eGFR, identified rs4072037 as a common alternative splice variant in MUC1, and observed increased renal expression of a specific MUC1 mRNA isoform as a plausible molecular mechanism underlying the GWAS association signal. Moreover, according to combined annotation-dependent depletion scores, rs4072037 and rs12411216 belonged to the top 6% and 10%, respectively, of most functionally significant SNPs in the human genome [20]. Of these, rs12411216 maps onto the CpG island within the promoter region for MUC1, whereas rs4072037 operates as an alternative splice site acceptor. Our data showed a genome-wide significant association of MUC1 variants with BUN and eGFR in TWB participants, thus supporting the critical role of MUC1 in renal functional impairment.
Higher albuminuria is associated with adverse clinical outcomes such as ESRD, cardiovascular disease, and mortality [24]. Teumer et al. [21] performed a GWAS to identify genetic loci associated with albuminuria in diabetes by analyzing 7877 participants of European descent and found MUC1 as one of the candidate loci. By performing trans-ethnic (n = 564,257) and European-ancestry-specific meta-analyses of the GWASs of the urine albumin-to-creatinine ratio (UACR), including ancestry-and diabetes-specific analyses, Teumer et al. [25] identified 68 UACR-associated loci and showed a genetic correlation among proteinuria, hyperlipidemia, gout, and hypertension. Furthermore, fine-mapping and trans-omics analyses with gene expression showed MUC1 as a gene potentially operating through differential expression in the kidney. With only the urine albumin level available for analysis, our data also revealed a genome-wide significant association between albuminuria and both MUC1 variants. Moreover, a lower frequency of microalbuminuria was found in the C allele of rs12411216 and the A allele of rs4072037, supporting MUC1 as a candidate gene for albuminuria in both European and Asian populations.

MUC1 Polymorphisms and Serum Uric Acid Levels and Gout
Two GWASs including 109,029 and 121,745 Japanese individuals have shown rs4072037 as the lead SNP for the uric acid level in the MUCI gene region [17,26]. In addition, our results showed a strong association of rs4072037 and rs12411216 with both the serum uric acid level and gout. Furthermore, Sandoval-Plata et al. [27] indicated that the nearby MTX1 SNP rs760077 showed a significant genome-wide association with gout compared with normal uric acid level controls in UK Biobank participants. All these findings suggest the importance of MUC1 in the development of hyperuricemia and gout.

MUC1 Polymorphisms and Blood Pressure Status and Waist Circumference
Our results demonstrated a borderline significant association of the MUC1 rs4072037 genotype with blood pressure status. Takeuchi et al. [19] analyzed 289,038 individuals of East Asian and European descent and found a significant genome-wide association between MUC1 rs4072037 and diastolic blood pressure. Moreover, ADTKD-MUC1 with the MUC1 mutation was found to be associated with ESRD and hypertension. For waist circumference, the lead SNP determined in our study was MTX1 rs760077 (p = 0.0001; Table 3), which showed a near complete linkage with rs4072037 (r 2 = 0.965). Recently, a cross-trait GWAS mapped the THBS3 SNP rs72704117 to waist circumference using participants from the UK Biobank [28]. Their results indicated that MUC1 and its nearby gene region may play a critical role in determining blood pressure status and central obesity.

Association between MUC1 Polymorphisms and Red Blood Cell Parameters
MUC1 can be expressed in nonepithelial cells including hematopoietic cells [7,29]. In bone marrow differentiating cells, MUC1 was strongly and selectively expressed during erythropoiesis but was absent in circulating erythrocytes. During erythropoiesis, MUC1 expression was transcriptionally regulated through phosphorylation. The sialylated MUC1 glycoforms selectively expressed on erythroid cells can act as a ligand for CD169, a macrophage-restricted adhesion molecule that is vital for erythropoiesis. These results suggested that MUC1 acts as a cross-talk molecule between erythroblasts and surrounding cells during erythropoiesis [7]. In addition, MUC1 is overexpressed and aberrantly glycosylated in hematological malignancies [1]. Our data revealed that both MUC1 variants were significantly associated with red blood cell parameters including hematocrit, hemoglobin, and red blood cell counts. A recent study examining variations in human blood cell traits reported MTX1 rs760077 as a candidate variant associated with all the three red cell parameters [30].

Previous Functional Studies for MUC1 rs4072037 and rs12411216 Polymorphisms
Both MUC1 rs4072037 and rs12411216 polymorphisms have been shown to be functional. MUC1 rs4072037 located in the 5' side of exon 2 operates as an alternative splice site acceptor. By conducting the RNA-ligase-mediated rapid amplification of the 5 complementary DNA end procedure, the rs4072037 variant was found completely associated with two major MUC1 transcripts in the gastric epithelium: variants 2 and 3 [15]. The allele A at rs4072037 causes a 9-amino acid deletion in exon 2 and consequently modifies both the signal peptide and N-terminal amino acid of the mature protein by changing the signal peptide cleavage site [31]. The effect of rs4072037 on the MUC1 mRNA alternative splicing isoform, the second most abundant isoform in the kidney, was examined by performing transcriptome analysis [20]. Carriers of one and two copies of the allele A of this splice variant exhibited the intermediate and highest expression levels of the alternatively spliced MUC1 isoform and the total renal expression of MUC1 gene when compared with the reference genotype, suggesting that the rs4072037-driven allelic effect on the expression of a specific MUC1 mRNA isoform can be the key biological mechanism underlying the genetic association observed in previous GWASs. In addition, the association of the rs4072037 genotypes with various phenotypes may be due to a strong LD with other functional SNPs. EFNA1 rs12904, an miR-200c binding site SNP that shows a strong LD with rs4072037 and is located in the 3'-untranslated region of the EFNA1 gene, can modulate EFNA1 expression and is associated with gastric cancer susceptibility [32]. The MUC1 rs12411216 polymorphism is located in the CpG island within the promoter region for MUC1 and is in near complete LD with rs4072037. Jiang et al. [33] revealed that the rs12411216 genotype was associated with mild cognitive impairment and glucocerebrosidase (GBA) expression in Parkinson's disease. By performing an electrophoretic mobility shift assay, the authors reported that the rs12411216 polymorphism could affect the binding efficiency of the transcription factor E2F4 and regulate GBA expression. Low GBA expression promoted the prion-like spread of α-Syn interpolymer complexes, facilitated the pathogenesis of Parkinson's disease, and increased cognitive damage [34,35]. These results characterized rs12411216 as a pathogenic variant of GBA in patients with Parkinson's disease with mild cognitive impairment.

Association between Functional MUC1 Variants and MUC1 Gene Methylation
DNA methylation is a crucial epigenetic marker, with the cytosine of CpG dinucleotide being the major methylation site. Hypermethylated CpG islands in promoters may silence the transcription of genes [36]. In gene-body methylation, exons were more highly methylated than introns, and sharp transitions of methylation occurred at exon-intron boundaries. Gene-body methylation is highly conserved across eukaryotes; in contrast to promoter methylation, gene-body methylation is positively associated with transcription but not linked to gene silencing [37,38]. DNA methylation in gene bodies may facilitate transcription elongation and/or co-transcriptional splicing and may repress intragenic cryptic promoters [39]. CCCTC binding factor, a DNA-binding protein that is bound to exon 5 of the CD45 gene, may be inhibited by CD45 gene-body DNA methylation, which results in local RNA polymerase II pause, reducing the transcription elongation rate, and favoring the co-transcriptional spliceosome assembly to promote weak upstream exon inclusion and alterative splicing [40]. Intragenic DNA methylation can modulate alternative splicing by recruiting methyl-CpG-binding protein to promote exon recognition [41]. Xu et al. [20] examined whether rs12411216 affects MUC1 promoter methylation and found that none of the CpG sites within the CpG island overlaying the MUC1 promoter region showed an association with the rs12411216 genotype. Our data also revealed no significant association between two functional MUC1 variants and the methylation status of 12 MUC1 promoter CpG sites. By contrast, we observed genome-wide significant associations between both rs4072037 and rs12411216 genotypes and four DNA methylation sites at the gene body of MUC1. Xu et al. [20] showed that the risk allele A at rs4072037 may increase the expression levels of the alternatively spliced MUC1 isoform. In our study, this risk allele was significantly associated with a highly methylation status at the MCU1 gene body. Thus, the possibility that the association between rs4072037 and the total or alternative spliced form of gene expression may be secondary to MUC1 gene-body methylation status should be highly considered. However, epigenetic effects are tissue-specific, and additional studies may be necessary to confirm the conclusion in different tissues and disease statuses.

Role of Nearby Gene Region Variants Determined by Performing Regional Plot Association Studies
By performing regional plot association studies, we tested whether the association between MUC1 functional variants and various phenotypes may be due to LD with adjacent gene variants. When the data of the lead SNPs and two functional MUC1 variants were combined in the regional plot analysis, we observed that a major association with the study phenotypes was focused on variants at the nearby TRIM46-MUC1-THBS3-MTX1 gene region with at least moderate LD between SNPs ( Figure 4A,B). The associations of the A allele of the rs4072037 genotype and the C allele of the rs12411216 genotype with lower systolic blood pressure, uric acid levels, and gout risk and a higher eGFR and hematocrit may be beneficial to participants, whereas the association with higher waist circumference, HbA1C, albuminuria, gastric cancer risk, and lower bone mineral density may be unfavorable to participants.

Limitations
Because only urine albumin but not urine creatinine was measured, the UACR could not be determined in our study. However, the association between MUC1 variants and albuminuria has been confirmed by previous studies [21,25].

TWB Participants
The study cohort for the regional plot association study consisted of participants with samples in the TWB; they were recruited from centers across Taiwan between 2008 and 2020. A total of 107,494 participants with no history of cancer and with genotyping from Axiom Genome-Wide CHB 1 or 2 Array were recruited, and 26,766 participants were excluded from the analysis based on the following criteria: no imputation data (12,289 participants), quality control (QC) for the GWAS (10,956 participants), fasting for <6 h (2862 participants), failure to genotype the rs4072037 or rs12411216 polymorphism (371 participants), and the absence of any study phenotype (288 participants). During QC for the GWAS, all of the participants excluded were due to relative pairs of 2nd degree relatives or closer with identity by descent (IBD) > 0.187. While examining the serum uric acid level, blood pressure status, and lipid and glucose metabolism parameters, we excluded participants with a history of gout, hypertension, hyperlipidemia, and diabetes mellitus, respectively. Figure 5 depicts the flowchart of participant enrollment. The definitions of hypertension, diabetes mellitus, obesity, hyperlipidemia, and current smoking are provided in Supplementary Method S1. Ethical approval was received from the Research Ethics Committee of Taipei Tzu Chi Hospital, Buddhist Tzu Chi Medical Foundation (approval number: 08-XD-005), and Ethics and Governance Council of the Taiwan Biobank (approval number: TWBR10908-01). Each participant signed an informed consent form before participation in the study.

Clinical Phenotypes and Laboratory Examinations
We examined the following clinical phenotypes: body height; body weight; waist circumference; waist-hip ratio; body mass index (BMI); and systolic, mean, and diastolic blood pressure. In addition, we collected the following biochemistry data: lipid profiles such as total, high-density lipoprotein (HDL), and low-density lipoprotein (LDL) cholesterol and triglyceride levels, glucose metabolism parameters such as fasting plasma glucose and hemoglobin A1C (HbA1C) levels, and liver and renal functional test-related parameters such as blood urea nitrogen (BUN), serum creatinine, uric acid, aspartate aminotransferase (AST), alanine aminotransferase (ALT), γ-glutamyl transferase (γ-GT), albumin, and total bilirubin levels. The BMI and eGFR were calculated as reported previously [42]. The hematological parameters analyzed included white and red blood cell counts, platelet counts, and hematocrit and hemoglobin levels. Because of the absence of the urine creatinine level, only the spot urine albumin level was used to evaluate albuminuria. Microalbuminuria was defined as a urine albumin level of ≥30 mg/L.

Selection of Functional MUC1 Variants and Genotyping
For TWB participants, DNA was isolated from blood samples by using a PerkinElmer Chemagic 360 instrument (PerkinElmer, Waltham, MA, USA) following the manufacturer's instructions. SNP genotyping was conducted using custom TWB chips and performed on the Axiom genome-wide array plate system (Affymetrix, Santa Clara, CA, USA). Two previously reported functional MUC1 variants, rs4072037 and rs12411216, were used for the analysis [15,31] (Table S9).

DNA Methylation Analysis
DNA methylation was assessed using sodium-bisulfite-treated DNA from whole blood by using the Infinium MethylationEPIC BeadChipEPIC array (Illumina Inc. San Diego, CA, USA). Four quality control measures were used: correction for dye bias across batches by normalization, removal of background signals, elimination of outliers by the median absolute deviation method, and elimination of probes with poor detection (p > 0.05) and those whose bead counts were <3. A total of 1686 participants were enrolled for the analysis of MUC1 DNA methylation.

Clinical Phenotypes and Laboratory Examinations
We examined the following clinical phenotypes: body height; body weight; waist circumference; waist-hip ratio; body mass index (BMI); and systolic, mean, and diastolic blood pressure. In addition, we collected the following biochemistry data: lipid profiles such as total, high-density lipoprotein (HDL), and low-density lipoprotein (LDL) cholesterol and triglyceride levels, glucose metabolism parameters such as fasting plasma glucose and hemoglobin A1C (HbA1C) levels, and liver and renal functional test-related parameters such as blood urea nitrogen (BUN), serum creatinine, uric acid, aspartate aminotransferase (AST), alanine aminotransferase (ALT), γ-glutamyl transferase (γ-GT), albumin, and total bilirubin levels. The BMI and eGFR were calculated as reported previously [42]. The hematological parameters analyzed included white and red blood cell counts, platelet counts, and hematocrit and hemoglobin levels. Because of the absence of Other phenotypes include: age, body mass index (BMI), waist circumference, waist-hip ratio, aspartate aminotransferase (AST), alanine aminotransferase (ALT), γ-Glutamyl transferase (γ-GT), serum creatinine level, estimated glomerular filtration rate (eGFR), serum albumin, total bilirubin, hemoglobin, hematocrit, red blood cell, leukocyte, and platelet counts, blood urea nitrogen (BUN), albuminuria, microalbuminuria, current smoking status, metabolic syndrome.

Regional Plot Association Analysis
To determine the lead SNP variants around the MUC1 region for various studied phenotypes, we first performed a quality check (QC) for GWAS. The Axiom genome-wide CHB 1 and 2 array plates (Affymetrix, Inc., Santa Clara, CA, USA), each with 27,720 and 67,485 participants, and comprising 611,656 and 640,160 SNPs, respectively, underwent the imputation analysis. By using the East Asian population from the 1000 Genome Project Phase 3 study as the reference panel, genome-wide genotype imputation was performed using SHAPEIT and IMPUTE2. After imputation, QC was performed by filtering SNPs with an IMPUTE2 imputation quality score of >0.3. Indels were removed using VCFtools. All samples enrolled for the analysis had a call rate of ≥97% with IBD ≥ 0.187. For SNP QC, an SNP call rate of <3%, a minor allele frequency of <0.01, and the violation of the Hardy-Weinberg equilibrium (p < 10 −6 ) were the criteria for exclusion from subsequent analyses. Finally, after QC, 84,249 participants and 3639,842 SNPs remained in the whole genome and 236 SNPs at positions between 155.0 and 155.5 Mb on chromosome 1q22 were enrolled for the regional plot association analysis.

Statistical Analysis
Continuous variables are expressed as the mean ± standard deviation. When the distribution was strongly skewed, median and interquartile ranges are presented that were examined using a two-sample t test or analysis of variance. Differences in categorical data distribution were examined using a chi-squared test or chi-squared test for trend. Before analysis, all study parameters were logarithmically transformed to adhere to a normality assumption. We assumed the genetic effect to be additive after adjustment for age, sex, BMI, and current smoking status, and a generalized linear model was used to analyze studied phenotypes in relation to the predictors of investigated genotypes and confounders. Regional plot association studies were performed using the analysis software package PLINK (version 1.07, Shaun Purcell, Cambridge, MA, USA, http://pngu.mgh.harvard.edu/purcell/plink, accessed on 14 August 2021). Genome-wide significance was defined as a p value of <5 × 10 −8 . The LDmatrix software (https://analysistools.nci.nih.gov/LDlink/?tab= ldmatrix, accessed on 19 April 2021) was used for the analysis of linkage disequilibrium (LD). All statistical analyses were performed using SPSS (version 22; SPSS, Chicago, IL, USA). A two-sided p value of <0.05 was considered to be statistically significant.

Conclusions
Our data confirmed the presence of pleiotropy at the TRIM46-MUC-THBS3-MTX1 gene region, especially at MUC1 gene functional variants, with multiple phenotypes being associated with candidate variants. Our data linked rs4072037 to MUC1 gene-body methylation, which may affect alternative splicing and increase the expression of MUC1 alternative splicing isoforms. The bidirectional effect of functional MUC1 gene variants on health should be examined in future studies. These results can provide insights into the critical role of MUC1 in the pathogenesis of cardiometabolic, renal, and hematological disorders.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Acknowledgments:
We greatly appreciate technical support from the Core Laboratory of the Taipei Tzu Chi Hospital, Buddhist Tzu Chi Medical Foundation, and expert statistical analysis assistance from Tsung-Han Hsieh.

Conflicts of Interest:
No potential conflict of interest relevant to this article was reported.