Molecular Classification of Hepatocellular Carcinoma Using Wnt–Hippo Signaling Pathway-Related Genes

Simple Summary The characters of Taiwanese hepatocellular carcinoma (HCC) are different from other parts of the world. We characterized the molecular features of HCC using a newly developed classification system based on the expression of the Wnt–Hippo signaling pathway-related genes. We analyzed the data in terms of prognostic value, transcriptome features, immune infiltration, and clinical characteristics, and compared the resulting subclasses with previously published classifications. A new molecular classification method based on a 272 gene panel of Wnt–Hippo pathways that may provide a new target for the treatment. Abstract In Taiwan, a combination of hepatitis B and C infection, economic boom-related food and alcohol overconsumption, and Chinese medicine prescriptions has led to a high rate of hepatocellular carcinoma (HCC). However, the causative factors and underlying tumor biology for this unique HCC environment have not been identified. Wnt and Hippo signaling pathways play an important regulatory role in HCC development, and their functions are generally considered as positive and negative regulators of cell proliferation, respectively. In this study, we characterized the molecular features of HCC using a newly developed classification system based on the expression of the Wnt–Hippo signaling pathway-related genes. RNA sequencing (RNA-Seq) was performed on liver tumor tissues from 100 patients with liver cancer. RNA-Seq data for 272 previously characterized Wnt–Hippo signaling pathway-related genes were used for hierarchical clustering. We analyzed the data in terms of prognostic value, transcriptome features, immune infiltration, and clinical characteristics, and compared the resulting subclasses with previously published classifications. Four subclasses of HCC (HCCW1–4) were identified. Subclass HCCW1 displayed the highest PCDHB4 expression. Subclass HCCW2 displayed lower Edmondson–Steiner grades (I and II) and CTNNB1 mutation frequencies. Subclass HCCW3 was associated with a good prognosis, the highest PCDHGB7 expression, high CD8+ naïve T cells abundance, and relatively low TP53 mutation rates. Subclass HCCW4 was associated with a poor prognosis, the highest PCDHB2 and PCDHB6 expression, a relatively high abundance of Th1 cells, NKT and class-switched memory B cells, relatively low enrichment of cDC, iDC, and CD4+ memory T cells, and high Edmondson–Steiner grades (III and IV). We also identified Wnt–Hippo signaling pathway-related genes that may influence immune cell infiltration. We developed a panel of 272 Wnt–Hippo signaling pathway-related genes to classify HCC into four groups based on Taiwanese HCC and The Cancer Genome Atlas Liver Hepatocellular Carcinoma datasets. This novel molecular classification system may aid the treatment of HCC.


Introduction
Liver cancer is the fifth most common cancer and the second common cause of cancerrelated death worldwide, and hepatocellular carcinoma (HCC) is its predominant form [1,2]. In addition to genetic factors, HCC development is associated with chronic viral infection, alcohol abuse, diabetes mellitus (DM), obesity, metabolic diseases, hemochromatosis, and autoimmune hepatitis [3,4]. The incidence of noninfectious HCC in developed countries is increasing due to increased rates of obesity, DM, and metabolic diseases [5,6]. These conditions induce liver injury and progressive inflammation, such that liver cells show a cycle of necrosis, regeneration, somatic mutations, and chromosomal instability [7][8][9]. Recent studies have revealed many genomic alterations of HCCs, and numerous frequently altered genes including TP53, CTNNB1, and TERT [10][11][12][13]. Surgery resection is the major treatment for HCC; however, many HCCs found at the unresected stage require chemotherapy or targeted therapy. Despite many potential therapeutic targets, few drugs have shown clinically promising effects, such as the multikinase inhibitors sorafenib and regorafenib, and immune checkpoint inhibitors; however, most of these drugs increased survival by only a few months, indicating the need for new drugs to treat HCC [4,6].
Nationwide programs have been initiated in Taiwan to reduce hepatitis B virus (HBV)and hepatitis C virus (HCV)-related HCC, beginning in 1984; aggressive countermeasures against HBV and HCV infections have also been implemented by the government. These measures have resulted in a reduction in the number of HBV carriers and a decrease in the rate of HBV-related HCC among younger individuals [14][15][16]. There is no vaccine available for hepatitis C even nowadays. Chinese herbs are also prescribed for the treatment of diseases, including chronic hepatitis, under the national health system; these herbs may contain compounds that modulate aristolochic acid or liver toxins [17]. Moreover, an economic boom was associated with higher rates of DM, metabolic disorders, and obesity [14,18,19]. HCC in Taiwan is unique due to the higher prevalence of chronic hepatitis B than in the Western world. Nonalcoholic fatty liver disease, alcoholic liver disease, and chronic hepatitis C are more common in Western countries.
Aberrant Wnt signaling promotes the development and/or progression of HCC. Mutations of CTNNB1, a gene that codes for β-catenin in the Wnt signaling pathway, are common in HCC [20]. Aberrant activation of Hippo signaling pathway components has also been observed in HCC, especially the Yes-associated protein and transcriptional co-activator with PDZ-binding motif transcription factors [21]. Inactivation of the Hippo pathway is significantly associated with poor prognosis in human HCC [22]. Molecular classification of HCC plays an important role in the selection of therapeutic strategies and evaluation of therapeutic responses and prognoses [10][11][12][13][23][24][25][26]. For example, a recently developed metabolic network classified HCC patients into three subgroups (iHCC1-iHCC3), where the iHCC2 phenotype was associated with aberrant Wnt signaling: 75% of iHCC2 tumors found to carry CTNNB1 mutations and displayed upregulated expression of β-catenin target genes [24]. Another study identified three HCC subtypes (S1-S3), where S1 tumors exhibited Wnt pathway activation through two mechanisms: CTNNB1 mutation and TGF-β activation [25].
In this study, we applied bioinformatic tools for RNA sequencing (RNA-Seq) analysis of 100 Taiwanese HCCs, and then used several approaches to develop a survival-related molecular classification method. This is the first study to analyze the Wnt-Hippo pathways for the molecular classification of HCCs in Taiwan. Our findings improve our understanding of HCC etiology and provide a basis for the identification of new treatment targets. from patients with more than 20 cancer types. This study was approved by the Ethics Committee of the CMUH (CMUH 109-REC3-055), and written informed consent was obtained from all participants according to the standard procedure of the CMUH Tissue Bank.

RNA Extraction and RNA-Seq
RNA was extracted from tissue samples using the NucleoSpin RNA Kit (Macherey-Nagel GmBH, Duren, Germany) following the manufacturer's instructions. The quality, quantity, and integrity of the RNA were evaluated using a NanoDrop1000 spectrophotometer and Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA).
RNA-Seq was performed as described previously [27]. Briefly, samples with an RNA integrity number > 6.0 were used for RNA-Seq. An mRNA-focused, barcoded library was generated using the TruSeq strand mRNA Library Preparation Kit (Illumina, San Diego, CA, USA). The libraries were sequenced on the Nova Seq 6000 instrument (Illumina), using 2 × 151-bp paired-end sequencing flow cells following the manufacturer's instructions.

RNA-Seq Data Analysis
The RNA-Seq data were analyzed as described previously [28]. Briefly, data quality control at the Q20 level was performed using the Trimmomatic tool [29], read alignment to the GRCh38 human genome was conducted using the HISAT2 alignment program [30], expression was quantified with reference to GENCODE v35 (excluding mitochondrial genes), and transcripts were normalized to transcripts per million (TPM) using the StringTie assembler [31].

Molecular Classification of HCC
Molecular classification of HCC was performed using transcriptomics-based analysis integrating Wnt [32] and Hippo [33] signaling pathways-related gene panels, coupled with survival analysis, and then validated using data from The Cancer Genome Atlas (TCGA). Subjects were also classified according to molecular signatures reported in other HCC studies [24,25]. For the molecular classification, hierarchical clustering was performed using Spearman rank correlation for sample distance calculation and the "average" linkage method using the Morpheus web tool (https://software.broadinstitute.org/morpheus/, accessed on 13 April 2022). Important nodes and subnetworks were predicted and explored using CytoHubba [34], a Cytoscape plugin [35].

Cell Enrichment Analysis
The xCell [36] was used to examine the enrichment of various immune cells in the tumors and to calculate an immune score from the TPM expression matrix.

Statistical Analyses
The Kaplan-Meier method was used to construct overall survival curves, and the p-values were obtained using a log-rank test. The Kruskal-Wallis test and Mann-Whitney U test were used to compare groups in terms of non-normally distributed variables. Contingency table variables were analyzed using Fisher's exact test. All statistical analyses were performed in the Python programming environment. A two-tailed p < 0.05 was considered statistically significant.

Molecular Subtypes of Taiwanese HCC Based on the Expression Profiles of 254 Wnt Pathway Genes
Hierarchical clustering analysis identified four subclasses (A-D) of Taiwanese HCC patients ( Figure 1A). Kaplan-Meier plots of the results of individual and grouped pairwise analysis showed significant differences in survival between subclasses B and D (p = 0.019) and C and D (p = 0.035), but there was no significant difference in a comparison of all four groups (p = 0.077) ( Figure 1B-H). Moreover, we investigated the expression of CTNNB1 up/downstream genes for CTNNB1 overexpression or down expression. FZD1, FZD3, FZD4, FZD5, FZD6, FZD7, GSK3B, TCF7L2, TLE4, LRP5, TLE2, TLE5, and TLE6 were found to be differentially expressed in CTNNB1 overexpression or down expression (Supplementary Figure S1). We found positive correlation for the expression of CTNNB1 with the expression of FZD1, FZD3, FZD4, FZD5, FZD6, FZD7, GSK3B, TCF7L2, and TLE4. However, the expression of CTNNB1 was inversely correlated with the expression of LRP5, TLE2, TLE5, and TLE6.

Molecular Subtypes of Taiwanese HCC Based on the Expression Profiles of 254 Wnt Pathway Genes
Hierarchical clustering analysis identified four subclasses (A-D) of Taiwanese patients ( Figure 1A). Kaplan-Meier plots of the results of individual and grouped wise analysis showed significant differences in survival between subclasses B and 0.019) and C and D (p = 0.035), but there was no significant difference in a compari all four groups (p = 0.077) ( Figure 1B-H). Moreover, we investigated the express CTNNB1 up/downstream genes for CTNNB1 overexpression or down expression. FZD3, FZD4, FZD5, FZD6, FZD7, GSK3B, TCF7L2, TLE4, LRP5, TLE2, TLE5, and were found to be differentially expressed in CTNNB1 overexpression or down expr (Supplementary Figure S1). We found positive correlation for the expression of CT with the expression of FZD1, FZD3, FZD4, FZD5, FZD6, FZD7, GSK3B, TCF7L2, and However, the expression of CTNNB1 was inversely correlated with the express LRP5, TLE2, TLE5, and TLE6.

Molecular Subtypes of Taiwanese HCC Based on the Expression Profiles of 272 Wnt-H Pathways Genes
A panel of Wnt pathway-related genes showed that some were corrected wit vival; therefore, we added several other functional gene panels to improve the pow the analysis. Finally, among 18 Hippo pathwayand 254 Wnt pathway-related gene identified four molecular subclasses (HCCW1-4) ( Figure 3A) that were correlated

Molecular Subtypes of Taiwanese HCC Based on the Expression Profiles of 272 Wnt-Hippo Pathways Genes
A panel of Wnt pathway-related genes showed that some were corrected with survival; therefore, we added several other functional gene panels to improve the power of the analysis. Finally, among 18 Hippo pathway-and 254 Wnt pathway-related genes, we identified four molecular subclasses (HCCW1-4) ( Figure 3A) that were correlated with patient's survival (p = 0.040) ( Figure 3B). The average survival time was significantly shorter for HCCW4 (2575 days; 95% confidence interval (CI): 1953-3197, n = 18) than for HCCW2 (3075 days; 95% CI: 2649-3502, n = 30), HCCW1 (3446 days; 95% CI: 2792-4101, n = 18), and HCCW3 (3556 days; 95% CI: 3143-3968, n = 34). HCCW3 was associated with a statistically significantly better prognosis than HCCW4 based on a two-sided log-rank test (p = 0.005) ( Figure 3C). HCC (LIHC) data from TCGA and obtained similar results (p = 0.008; Figure 4A,B). HCCW4 showed significantly lower survival than HCCW2 (p = 0.003) ( Figure 4C) and HCCW3 (p = 0.007) ( Figure 4D), with a more distinct grouping than that seen for Taiwanese HCC patients, perhaps due to a larger sample size. The Bidkhori classification method [24] applied to our cases based on a panel of 51 metabolic genes revealed no correlation with survival (p = 0.943) ( Figure 3D), whereas the results of the Hoshida classification [25] showed a correlation with survival (p = 0.039) ( Figure 3E).  To validate the findings, we applied a similar molecular classification approach using HCC (LIHC) data from TCGA and obtained similar results (p = 0.008; Figure 4A,B). HCCW4 showed significantly lower survival than HCCW2 (p = 0.003) ( Figure 4C) and HCCW3 (p = 0.007) ( Figure 4D), with a more distinct grouping than that seen for Taiwanese HCC patients, perhaps due to a larger sample size. The Bidkhori classification method [24] applied to our cases based on a panel of 51 metabolic genes revealed no correlation with survival (p = 0.943) ( Figure 3D), whereas the results of the Hoshida classification [25] showed a correlation with survival (p = 0.039) ( Figure 3E).

Transcriptomes of the HCC Subclasses
We further analyzed the expression of genes that may influence the molecular classification of Taiwanese HCC patients, the results showed differential expression of many genes among the four groups; however, we detected no correlations with survival (Supplementary  Table S1). We further analyzed survival-related expression genes and found that PCDHB6, PCDHGB7, PCDHB2, and PCDHB4 were correlated with the molecular classification and patient survival (Figure 5A-D and Supplementary Table S1). Patients with high PCDHB6, PCDHB2, and PCDHB4 mRNA expression levels had lower overall survival. By contrast, high PCDHGB7 mRNA expression was associated with longer survival. We used similar approaches to analyze the TCGA data, and the results showed that 56 genes were correlated with both the molecular classification and survival (Supplementary Table S1).

Transcriptomes of the HCC Subclasses
We further analyzed the expression of genes that may influence the molecular classification of Taiwanese HCC patients, the results showed differential expression of many genes among the four groups; however, we detected no correlations with survival (Supplementary Table S1). We further analyzed survival-related expression genes and found that PCDHB6, PCDHGB7, PCDHB2, and PCDHB4 were correlated with the molecular classification and patient survival ( Figure 5A-D and Supplementary Table S1). Patients with high PCDHB6, PCDHB2, and PCDHB4 mRNA expression levels had lower overall survival. By contrast, high PCDHGB7 mRNA expression was associated with longer survival. We used similar approaches to analyze the TCGA data, and the results showed that 56 genes were correlated with both the molecular classification and survival (Supplementary Table S1).  Table S2). The most important nodes and subnetworks of protein-protein interactions among the 110 upregulated genes in HCCW4 were predicted and explored using CytoHubba; the 10 most significant node genes were WNT8B, SMARCA4, ARID1A, SMARCB1, SMARCC2, SMARCD1, SMARCD2, SMARCD3, FZD2, and WNT9A ( Figure 6).
We also identified 26 and 11 differentially expressed genes that were up-and downregulated, respectively, in HCCW3 compared with those in HCCW1, HCCW2, and HCCW4 (Supplementary Table S3).

Correlation of HCC Subclasses with Immune Infiltration
We investigated the associations between the subclasses and the expression levels of 34 immune cells. Significant differences were observed between HCCW4 and the other three subclasses, with higher abundance seen for three immune cell populations (Th1, natural killer T (NKT), and class-switched memory B-cells) in HCCW4 than compared with HCCW1-3 ( Figure 7A), and lower enrichment for three immune cell populations (cDC, iDC, and CD4+ memory T-cells) ( Figure 7B). In the nonhematopoietic group, hepatocytes were significantly less enriched in HCCW4 ( Figure 7C). In addition, CD8+ naïve T-cells were more abundant in HCCW3 ( Figure 7D). Among the significantly differentially expressed genes related to the Wnt-Hippo pathway, 110 and 14 were up-and downregulated, respectively, in HCCW4 compared with those in HCCW1-3 (Supplementary Table S2). The most important nodes and subnetworks of protein-protein interactions among the 110 upregulated genes in HCCW4 were predicted and explored using CytoHubba; the 10 most significant node genes were WNT8B, SMARCA4, ARID1A, SMARCB1, SMARCC2, SMARCD1, SMARCD2, SMARCD3, FZD2, and WNT9A ( Figure 6).    Table S2). The most important nodes and sub networks of protein-protein interactions among the 110 upregulated genes in HCCW were predicted and explored using CytoHubba; the 10 most significant node genes wer WNT8B, SMARCA4, ARID1A, SMARCB1, SMARCC2, SMARCD1, SMARCD2, SMARCD3 FZD2, and WNT9A ( Figure 6).  34 immune cells. Significant differences were observed between HCCW4 and the other three subclasses, with higher abundance seen for three immune cell populations (Th1, natural killer T (NKT), and class-switched memory B-cells) in HCCW4 than compared with HCCW1-3 ( Figure 7A), and lower enrichment for three immune cell populations (cDC, iDC, and CD4+ memory T-cells) ( Figure 7B). In the nonhematopoietic group, hepatocytes were significantly less enriched in HCCW4 ( Figure 7C). In addition, CD8+ naïve T-cells were more abundant in HCCW3 ( Figure 7D).

Correlation of the HCC Subclasses with Clinical Characteristics
Next, we explored the tumor-related clinicopathological variables associated with our classification results. Fisher's exact test revealed significant correlations between clinicopathological features and HCC subclasses in our cohort. A lack of advanced Edmondson-Steiner grade (p = 0.0007) and lower CTNNB1 mutation frequency (p = 0.0586) were associated with the HCCW2 subclass, and Edmondson-Steiner grade III/IV (p = 0.0203)

Correlation of the HCC Subclasses with Clinical Characteristics
Next, we explored the tumor-related clinicopathological variables associated with our classification results. Fisher's exact test revealed significant correlations between clinicopathological features and HCC subclasses in our cohort. A lack of advanced Edmondson-Steiner grade (p = 0.0007) and lower CTNNB1 mutation frequency (p = 0.0586) were associated with the HCCW2 subclass, and Edmondson-Steiner grade III/IV (p = 0.0203) was associated with HCCW4 subclass. HCCW3 was negatively correlated with the TP53 mutation rate (p = 0.0049) ( Table 1). We also compared our classification results with those of the Bidkhor (iHCC1-3) and Hoshida (S1-S3) classification systems. HCCW3 was marginally significantly associated with the Hoshida S3 class (p = 0.0526), whereas HCCW4 was significantly associated with the Bidkhor iHCC3 class (p < 0.0001) and Hoshida S1 and S2 classes (p = 0.0039) ( Table 1).

Discussion
Comprehensive integrated analysis of HCC is a powerful tool for understanding molecular events relevant to HCC [10][11][12][13]24,25]. In this study, we used a functional gene panel approach to explore correlations of molecular events with molecular changes and survival and devised a new molecular panel including 254 and 18 Wnt and Hippo signaling pathway-related genes, respectively, to classify Taiwanese HCC patients into four molecular subgroups (HCCW1-4). Then, the prognostic value, transcriptome features, immune infiltration, and clinical characteristics of the subgroups were explored.
The expression of PCDHB4 was highest for HCCW1. HCCW2 was associated with relatively low Edmondson-Steiner grade (I-II) and a relatively low CTNNB1 mutation frequency rate. HCCW3 was associated with a good prognosis, the highest expression of PCDHGB7, relatively high CD8+ naïve T-cell abundance, and relatively lower TP53 mutation rates. HCCW4 was associated with a poor prognosis, the highest PCDHB2 and PCDHB6 expression levels, relatively high abundance of four immune cell populations (Th1, NKT, and class-switched memory B-cells), relatively low enrichment of seven immune cell populations (cDC, iDC, and CD4+ memory T-cells), and relatively high Edmondson-Steiner grade (III-IV). Our method was also applied to TCGA HCC data; the results showed that the subgroup with the poorest survival (HCCW4) was well correlated with Hoshida classification groups S1 and S2, and Bidkhori metabolic classification group iHCC3 [24,25].
PCDHs are a group of cell adhesion molecules in the cadherin superfamily that can be divided into clustered (cPCDHs) and nonclustered molecules (ncPCDHs). PCDHs can interact with intracellular molecules including components of the WAVE complex, Wnt pathway, and apoptotic cascades. PCDHs are expressed prominently in the central nervous system, where they play important neurobiological roles. They are also involved in cancer development; loss or dysregulation of PCDHs is associated with multiple types of cancer [37]. Most PCDH genes are considered to be tumor suppressors, and epigenetic mechanisms govern their expression through their unique genomic organization that allows for long-range epigenetic silencing [37,38]. Several studies have shown that PCDH genes are also involved in liver cancer development, including PCDH10, PCDH17, PCDH19, PCDH20, PCDHGC4, and PCDHGC5 [39][40][41][42][43]. Alterations of these PCDHs can occur through several oncogenic pathways; for example, PCDH20 functions as a tumor-suppressor gene by antagonizing the Wnt/β-catenin signaling pathway in HCC [43]. In this study, PCDHB2, PCDHB4, and PCDHB6 were overexpressed in the HCCW4 subgroup, in association with worse survival; this finding contradicts a tumor suppressor role. Our PCDHB6 results were confirmed by the TCGA-LIHC data. Aberrant expression of these genes may have resulted from context-dependent differences among PCDHs. We also found that PCDHGB7 upregulation was correlated with better HCC survival, which suggests a tumor suppressor role for PCDHGB7 in HCC. Moreover, we examined the PCDHB2, PCDHB4, PCDHB6, and PCDHGB7 mRNA and protein expression level by TCGA and Human Protein Atlas [44] databases, respectively. Among them, the RNA expression of PCDHB4 and PCDHGB7 were correlated with the protein level in LIHC. However, PCDHB2 and PCDHB6 expression were detected only at the mRNA level. PCDHB2 and PCDHB6 protein levels were not observed in >50% of liver cancer specimens by immunohistochemistry.
Wnt signaling pathway components play important regulatory roles in immune suppression and immune cell exclusion within the tumor microenvironment [45,46], and include macrophages and DC, T, NK, and myeloid-derived suppressor cells [47]. In a mouse model of HCC, β-catenin pathway activation promotes immune escape and resistance to anti-PD-1 treatment through Wnt/β-catenin signaling; this resulted in defective recruitment of DCs and impaired T cell activity, in turn leading to an impaired anti-tumor immune response [48]. The Wnt and Hippo pathways play important role in HCC development [20,[49][50][51]; our results confirm these findings and suggest that more aggressive and specific treatment options are needed for the poor survival group.
In HCC, TP53 loss of function is associated with high levels of centrosome amplification, aneuploidy cell proliferation, and chromosome instability [10,52], as well as poor prognosis [53]. In the current study, the HCCW3 subclass had lower TP53 mutation frequencies; patients in this group had a better prognosis than those in the other groups. Studies assessing the prognostic impact of TP53 mutations on the outcomes of gastroesophageal adenocarcinomas have demonstrated conflicting results [54][55][56]. From these results, we suggest that TP53 mutations may not completely correlate with poor survival, and other genetic alterations may also play a role for the survival of patient as our two patients with TP53 mutations in the HCCW3 group with a good prognosis.
In this study, we obtained a list of Wnt-Hippo pathways-related genes which may influence the infiltrations of immune cells between the good and poor prognosis groups. We also found higher proportion of Th1, NKT, and class-switched memory B-cells in the poor prognosis group. NKT cells contribute to fibrosis progression in nonalcoholic fatty liver disease [57]. A total of 110 upregulated Wnt-Hippo signaling pathway-related genes were identified in our poor prognosis group. Among these genes, the 10 most significant node genes were WNT8B, SMARCA4, ARID1A, SMARCB1, SMARCC2, SMARCD1, SMARCD2, SMARCD3, FZD2, and WNT9A. These genes may play important regulatory roles in immune cell migration, differentiation, and activation in HCC. Furthermore, we also conducted RNA-Seq in 12 normal samples. The expression of 88 genes were significantly higher in HCC tissues than that in adjacent nontumor tissues (82 genes in Wnt signaling pathway and 6 genes in Hippo signaling pathway), while the expression of 40 genes were significantly lower in HCC tissues than that in adjacent nontumor tissues (37 genes in Wnt signaling pathway and 3 genes in Hippo signaling pathway) (Supplementary Table S4). These alterations may be used as markers for the diagnosis of HCC and need more cases to confirm it.

Conclusions
In summary, comprehensive and integrated approaches, such as that developed in this study, will lead to a better understanding of the molecular mechanisms of HCC, and could in turn lead to the discovery of new therapeutic strategies. The results of this study highlight the uniqueness of Taiwanese HCC patients; this uniqueness arises from complex interactions among genetic, lifestyle, and environmental conditions that have arisen in the past 40 years.

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki, and approved by the Ethics Committee of the China Medical University Hospital (CMUH 109-REC3-055).
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The RNA-Seq data from this study was submitted to the NCBI Sequence Read Archive (SRA) under BioProject accession nos. PRJNA866195.