TBX5 Variants are Associated with Susceptibility to and the Incidence of Liver Cirrhosis and Hepatocellular Carcinoma in the Chinese Population: A Multicenter and Follow-Up Study

Introduction

Chronic hepatitis B virus (HBV) infection occurs worldwide but is endemic in East Asia and Africa, and the progression of chronic HBV infection includes liver cirrhosis (LC) and hepatocellular carcinoma (HCC), which account for 500000 to 1.2 million deaths per year worldwide.¹ Despite the identical cause, the progression of chronic infection and outcome of primary infection are highly variable. More than 350 million people have chronic HBV infection, but only a subset of these individuals, 50% of male carriers and 14% of female carriers,^2,3 will eventually die due to progression to LC and HCC.

Progression is a multifactorial and interindividual-variation disorder involving treatment for HBV infection, coinfection with other hepatotropic viruses⁴ and human immunodeficiency virus (HIV); sex; age;⁵ alcohol; smoking; and human genetic variants. The association of human genetic variants with progression and persistent infection has been recently emphasized in Asia,^6–10 thereby broadening our understanding of genetic predisposition. The familial clustering of hepatocellular carcinoma in HBV endemic areas also consolidated the genetic role in disease progression.¹¹ LC is a precancerous lesion, and some cases of HCC are accompanied by concurrent LC, indicating that there might be some common genetic susceptibility and suggesting that the internal genetic relationship might be an informative and unnoticed part of previous association studies. Therefore, it is of great interest to identify clear susceptibility loci by conducting a parallel comparison among progressive subgroups and primary infection.

A joint analysis strategy was used, including a DNA pooling screen in the genome-wide analysis (GWAS) and subsequent independent replications, to explore genetic susceptibility to LC and HCC, to examine their relation to outcomes of primary infection, explore significance of risk single nucleotide polymorphism (SNP) and to further assess the effect of risk loci on the incidence of progression in follow-up analysis.

Methods

Genetic Association Design and Subject Enrollment

More than 906,600 SNPs in the screen phase stage were initially genotyped in DNA pools, including the pool of HCC (case), the pool of LC (case) and the pool of asymptomatic chronic HBV carriers (AsC) (control), to identify the progressive genetic loci. Associations detected in this pool-based GWAS phase were validated in four independent cohorts using individual genotyping, and the combined samples from four independent cohorts were used to evaluate the genotypic effects on progression.¹² Additionally, the risk loci from progressive cohorts were assessed for susceptibility to consequences of primary HBV infection in three persistent infection cohorts, consisting of subjects with persistent infection (case), subjects who were naivety to HBV infection (control) and subjects who have spontaneously recovered (SR) (control) in each cohort.

The subgroup of persistent HBV infection in Hubei Province was followed-up to weigh the contribution of risk genotypes on the incidence of progressive outcomes from enrollment to December 31, 2017.

All unrelated Han Chinese subjects were recruited at four tertiary hospitals in Hubei Province, central China; three tertiary hospitals in Guangdong; and three tertiary hospitals in Hainan, southern China. Replication 1 progressive cohort, replication 2 progressive cohort and persistent infection cohort 1 were from Hubei Province; replication 3 progressive cohort and persistent infection cohort 2 comprised subjects from Guangdong; and replication 4 progressive cohort and persistent infection cohort 3 comprised subjects from Hainan.

Participants were enrolled according to the uniform predesigned protocol (Supplemental Table S1) in these centers. To maintain homogeneity in progression, the natural history and diagnosis were confirmed in detail by reviewing medical records and obtaining the medical history of all subjects. The patients with persistent infection, LC, HCC and AsC had no serologic evidence for coinfection with hepatitis C virus (HCV), hepatitis E virus (HEV), hepatitis D virus (HDV), cytomegalovirus (CMV), HIV, or medications associated with liver injuries.

The study was approved by the Ethics Committee of Hainan General Hospital and was performed in accordance with the Helsinki Declaration II. Written informed consent was obtained from patients and their next of kin for patients.

DNA Isolation and Progressive DNA Pooling Construction

Genomic DNA was isolated from whole blood collected in EDTA tubes using a QuickGene DNA whole blood kit S (DB-S) with QuickGene-Mini80 equipment (Fujifilm, Tokyo, Japan). The yield and purity of DNA was qualified by ultraviolet light spectroscopy and divided into three aliquots at −70°C for DNA analysis.

Eighty-eight patients in the LC group, 90 patients in the HCC group, and 66 patients in the AsC group from replication were used to create DNA pools, and the detailed process of DNA pooling has been released in GEO (GSE26034). In accordance with the protocol of DNA pooling,¹³ the DNA concentration of all samples was quantified four times using a Quant-iT PicoGreen dsDNA assay kit (Invitrogen, Carlsbad, CA) and tested by real-time PCR. Each individual contributed a total of 120 ng DNA to create the DNA pooling and de novo a total of three times.^12,13 A total of 12 DNA pools were established, and once created, each pool was diluted to 50 ng/µL for GWAS.¹⁴

DNA-Pool Analysis and Significant SNP Selection

The DNA pools were genotyped according to the Genome-Wide Human SNP 6.0 guide and workflow (Affymetrix, Santa Clara, CA). Sex chromosome markers were excluded in the analysis stage. The silhouette method was employed to find frequency differences between pools by calculating the ratio of allele signal intensities PMA/(PMA+PMB).¹⁵ A silhouette value greater than 0.6 was used as a threshold for further replication, thereby minimizing the possibility of false-positives at the apparent cast of false-negatives. SNPs with MAF below 0.05 were excluded because of the limited power to detect a significant association and the decrease in the accuracy of DNA pooling with rear variants. Two statistical strategies were used to select the top risk variants. For the first strategy, the SNP was significant in single-point (single method) of all three progressive outcomes with a Silhouette value >0.6 as the cutoff. This strategy was proposed by Kaufman and Rousseeuw, and they posited that a Silhouette value greater than 0.5 indicated an acceptable association. For the second strategy, the SNP was significant in window-points (sliding-window method) within the range of 7.5 kb,10 kb,12.5 kb and 15 kb throughout whole genome.

Replicated Individual Genotyping

All subjects were individually genotyped according to the manufacturer’s instructions using the TaqMan 7900HT sequence Detection System (Applied Biosystems, Foster City, Calif.). Each assay was carried out using 10 ng DNA in a 5-µL reaction consisting of TaqMan universal polymerase chain reaction master mix (Applied Biosystems, Foster City, Calif.), forward and reverse primers, and 6-carboxyfluorescein (FAM) and 4,7,2’-trichloro-7’-phenyl-6-carboxyfluorescein (VIC)-labeled probes designed by Applied Biosystems (ABI Assays-on-demand). Allelic discrimination was measured automatically using Sequence Detection Systems 2.1 software (automatic confidence level 95%). Approximately 8% of all genotypes were repeated in independent polymerase chain reactions to check for consistency and to ensure intraplate and interplate genotyping quality control. No genotyping discrepancies were detected between the repeated samples (Figure 1).

Figure 1 The workflow of association analysis.

Genotypic Function Exploration by Web-Based Bioinformatics

We used RNAfold (http://rna.tbi.univie.ac.at//cgi-bin/RNAWebSuite/RNAfold.cgi) to predict the biological effect of the significant SNPs on RNA structure. RNAfold is a classic database to predict RNA structure and energy change in RNA formation according to RNA sequence. Free energy represents the energy required to change the secondary structure from the current RNA structure. Additionally, we further explored the effects of SNPs on gene expression by investigating a public database of the GTEx portal (http://www.gtexportal.org/home/).

Follow-Up Outcomes

The follow-up time was calculated from the date of enrollment to either the date of patients with LC or HCC occurrence or December 30, 2017, for those without LC or HCC occurrence. Progressive events were ascertained at an average of three months by phone, out-of-hospital and hospitalized medical records. Patients will be diagnosed with HCC if they have LC and HCC at the same time. Family history defined the direct relatives who had HBV-related complications. LC was diagnosed by liver ultrasonography and abdominal computerized tomography features suggestive of cirrhosis based on a quantitative scoring system derived from the appearance of the liver margins, portal vein calibre and spleen diameter, supplemented with the cirrhotic complications (portal hypertension or bleeding in the esophageal varices or gastric varices, swelling in the legs and abdomen and splenomegaly). HCC was diagnosed based on pathologic confirmation or a positive lesion detected by imaging techniques (liver ultrasonography and abdominal computerized tomography) and a-fetoprotein levels exceeding 400 mg/l.

Statistical Analysis

The silhouette value was used to estimate genotypic cluster differences between cases and controls by the R package (www.r-project.org). This process was performed in a blinded manner, so the analyzer did not know the case/control status of each subgroup. The replication analysis was calculated using STATA software (version 10.0/SE, Stata, College Station, Tex.). Hardy-Weinberg equilibrium (HWE) assumptions were independently tested in each persistent infection cohort, chi-square analysis was employed to examine the differences in allele frequencies and genotype, the Cochran-Armitage test was used for genotypic trend analysis, and three genetic models fit the genetically susceptible type.

Unconditional logistic analysis was employed to estimate the odds ratio (OR) and 95% CI. In the combined analysis, the ORs and 95% CIs were adjusted by sex, age and cohort.

In follow-up analysis, the effects of risk genotype(s) were assessed by Kaplan‒Meier analyses and Cox regression analysis. In multivariate analysis, all available baseline covariates with P < 0.20 in unadjusted analysis were considered potential independent predictors. By using a backward stepwise approach with factors sequentially eliminated according to the results of the likelihood ratio test, we constructed Cox proportional hazards models to examine the association between the risk SNPs and progressive events.

Results

Patient Characteristics

All subjects were self-reported Han Chinese ancestry. A total of 8906 participants were recruited from June 2007 to June 2020, comprising 1610 HCC patients, 1430 LC patients, 1540 AsC patients, 1406 persistent infection patients, 1406 naïve to HBV infection participants and 1514 natural clearance participants. The subjects naïve to HBV infection and natural clearance in each replication showed no significant deviation from HWE. The characteristics of the participants in each cohort are shown in Table 1.

Table 1 Characteristics of Participants in Each Cohort

Pooled GWAS Data for Screening Risk SNPs for Replication

After quality control and silhouette score ranking, 1848 SNPs with average silhouette widths were more than 0.7 in the single-point method (data not shown). According to predesigned SNP selection strategies, 7 SNPs (Table 2) met the selection criteria and were prioritized for further replications: rs11119081, rs2638121, rs7021235, rs3825214, rs6457617 and rs71699 from the single method, and rs2996466 from the sliding-window method.

Table 2 The Selected SNPs for Replication Genotyping

Replication 1 Panel for Progression

To test the potential association of 7 SNPs, 1009 subjects consisting of 234 AsC patients, 300 LC patients and 300 HCC patients were individually genotyped in the screening stage. We compared the frequency of both allele and genotype between AsC versus LC and AsC versus HCC. The P_allele and P_genotype values for each SNP are shown in Table 3. Of these SNPs, rs3825214 within TBX5 (T-box transcription factor 5) was strongly associated with LC (P_dominant=0.004) and HCC (P_dominant= 0.0009). The remaining SNPs showed no significant associations with either LC or HCC (Table 3).

Table 3 Association of rs3825214 in Each Progressive Cohort

Replication 2 Panel for Progression

To confirm the significance in additional Hubei cohorts, 1383 subjects consisting of 435 patients with AsC, 394 patients with LC and 454 patients with HCC were further assessed. The P_allele and P_genotype values for rs3825214 are shown in Table 3. Again, the association of rs3825214 with LC was observed (P_dominant=0.0027), and the association of rs3825214 with HCC was discovered. (P_dominant=0.0089). However, the remaining four SNPs were not nominally significant.

Replication 3 Panel for Progression

To validate the results, 1381 subjects from Guangdong, including 450 patients with AsC, 418 patients with LC and 392 patients with HCC, were evaluated. The P_allele and P_genotype values for rs3825214 are presented in Table 3. Consistently, rs3825214 of TBX5 showed an association in both LC (P_dominant<0.001) and HCC (P_dominant=0.0014).

Replication 4 Panel for Progression

To further validate the results, 1151 subjects from Hainan, including 421 patients with AsC, 318 patients with LC and 412 patients with HCC, were genotyped. The P_allele and P_genotype values for rs3825214 are shown in Table 3. Consistently, rs3825214 of TBX5 showed evidence of association in both LC (P_dominant<0.001) and HCC (P_dominant=0.0014).

Combined Data Analysis for Progression

For all the subjects of Han Chinese ancestry, rs3825214 within TBX5 was in the same direction among the cohorts, and the genotype frequencies of the risk loci in persistent infection cohorts conformed to HWE expectations. Therefore, we combined the samples across all progressive cohorts. In the combined analysis, rs3825214GG was associated with an increased risk of LC (versus AA/AG P < 0.001; OR = 1.98; 95% CI: 1.65–2.38) and HCC (P < 0.001; OR, 1.68; 95% CI: 1.4–2.02).

Association of rs3825214 with Persistent HBV Infection

The susceptible relationship of rs3825214 with persistent HBV infection, natural clearance and naivety to HBV infection was assessed in three independent panels: persistent infection cohort 1 comprised 666 individuals with persistent HBV infection, 757 individuals who were naïve to HBV infection, and 881 individuals with natural clearance; persistent infection cohort 2 comprised individuals 340 with persistent HBV infection, 320 individuals who were naivety to HBV infection, and 321 individuals with natural clearance; and persistent infection cohort 3 included 409 individuals with persistent HBV infection, 332 individuals who were naivety to HBV infection, and 315 individuals with natural clearance.

The genotype and allele of rs3825214 presented no significant association in the two-two comparison among persistent HBV infection, natural clearance and naïve to HBV infection (Table 4).

Table 4 Association of rs3825214 with Persistent HBV Infection, Natural Clearance and Naïve to HBV Infection

Function Exploration and Expression Alternation

We used RNAfold and the GTEx portal to predict the function of rs3825214 (G>A). In the RNAfold analysis, the centroid secondary and minimum free energy (minimum free energy, MFE) structures of rs3825214 are shown in Figure 2. The genotype of rs3825214 alters the minimum free energy, free energy of the thermodynamic ensemble, frequency of the MFE structure in the ensemble and ensemble diversity, minimum free energy (G: −41.60 kcal/mol vs A: −39.50 kcal/mol), free energy of the thermodynamic ensemble (G: −46.13 kcal/mol vs A: −44.78 kcal/mol), frequency of the MFE structure in the ensemble (G: 0.06% vs A: 0.02%), and ensemble diversity (G: 44.91 vs A: 55.83) (Figure 2).

Figure 2 The RNAfold algorithm predicts the genotypic impact of rs3825214 on centroid secondary and minimum free energy. (A) MFE and centroid secondary structure of rs3825214-A. (B) MFE and centroid secondary structure of rs3825214-G.

In expression analysis, we found that rs3825214 contributed to expression quantitative trait loci (eQTL) based on the public database GTEx Portal. The rs3825214 genotype significantly altered the intron excision ratio of TBX5 in the heart (P=1.924e-8) (Figure 3).

Figure 3 Violin plots showing the association between normalized intron-excision ratio and rs3825214 genotype. The GG allele was significantly correlated with increased intron excision in the heart. Data were generated using the Single‐Tissue eQTL function in the GTEx Data Portal.

Multivariate Analysis of Follow-Up Outcomes

A total of 652 persistent infection patients at baseline were included for follow-up in persistent infection cohort 1. Twenty-one (3.2%) patients lost their phone connection and were without hospitalized register records, 31 (4.7%) declined to offer disease information, and 29 (4.4%) patients refused follow-up. To preclude the potential bias created by loss to follow-up, the 81 subjects were compared, and there was no significant difference in genotype distribution of rs3825214 (MM:Mm:mm 25:44:12/190:279:106, P=0.58), sex (male 62 versus 496, P=0.51), or age (40.2 versus 38.6, P=0.21) compared with subjects who did complete follow-up examinations.

A total of 571 patients were included in the follow-up analysis. During a median of 5.1 years (range 0.098–7.31 years) of follow-up with a total of 2876.2 person-years, 93 (16.29%) patients had LC, and 74 (12.96%) had HCC.

In univariate Cox analysis, the risk factors independently predicting LC included age (P<0.001, OR = 1.04; 95% CI: 1.02–1.06), alcohol never (versus former, P=0.002, OR = 2.3; 95% CI: 1.34–3.92; versus current, P=0.009, OR = 1.97; 95% CI: 1.19–3.27), INR (P<0.001, OR = 3.13; 95% CI: 1.92, 5.11), HBeAg persistence (P<0.001, OR = 2.98; 95% CI: 1.94–4.59), HBV positive (P<0.001, OR = 2.21; 95% CI: 1.44–3.37), and rs3825214GG (versus AG+AA, P<0.001, OR = 3.03; 95% CI: 1.95–4.70) (Table 5).

Table 5 Univariate Cox Analysis of Baseline Characteristics of Follow-Up Subjects

The predictive risk factors for HCC included age (P<0.001, OR = 1.07; 95% CI: 1.05, 1.09), HBeAg persistence (P<0.001, OR = 3.01; 95% CI: 1.84, 4.91), HBV positivity (P<0.001, OR = 2.14; 95% CI: 1.48, 3.93), and rs3825214GG (versus AG+AA, P<0.001, OR = 4.32; 95% CI: 2.68, 6.89) (Table 6).

Table 6 HBV-Related LC and HCC in Cox Proportional Hazard Regression Models

After accounting for covariates, we also found that rs3825214GG remained associated with the incidence of LC and HCC in Kaplan‒Meier analyses (Figure 4).

Figure 4 Kaplan‒Meier graphs of time to progressive events, (A) rs3825214 GG (21/85, 24.7%) vs AG/AA (25/444, 5.6%) with log-rank P value <0.0001 in HCC incidence; (B) rs3825214 GG (24/82, 29.3%) vs AG/AA (46/423, 10.9%) with log-rank P value <0.0001 in LC incidence.

Discussion

By investigating a large number of HBV-related subjects in the Chinese Han population from multiple centers, rs3825214 within TBX5 was specifically associated with LC and HCC. With at least 5 years of follow-up of persistent HBV infection patients, we observed that the susceptibility variants increased the incidence of LC and HCC, which highlights the genotype-phenotype link between risk genotypes and progressive events.

These results underscore several important principles. First, consistent with previous analysis, susceptibility exploration on a genome-wide scale has the potential to detect novel risk loci and provide independent information beyond candidate genes for progression, shedding light on the pathophysiology/function of HBV-related progression. Second, the associations are detectable in both LC and HCC, considering etiologic homogeneity, suggesting that the progressive phenotypes might share the pathological pathway, which indicates the transferability of preventive and therapeutic options for progression. Third, the genotype(s) modified the risk of LC and HCC provide incremental clues about genetic risk discrimination for chronic HBV-infected individuals. Fourth, the risk polymorphism was related to the phenotype and outcomes in long-term observation.

The rs3825214 is located in the intron region of TBX5 at 12q24.1, harboring a 5-kb linkage disequilibrium region at the 3’-end region, indicating that the 5-kb haplotype block including the promoter region as a whole passes genetic information into the next generation (Supplementary Figure 1). The LD structure and biologic rationale are supportive clues for TBX5 genetic variant involvement in disease progression. TBX5 encodes a transcription factor involved in multiple organ differentiation at the embryogenesis stage, and mutations in TBX5 cause Holt-Oram syndrome in humans, characterized by congenital cardiac defects and forelimb malformations. The genetic variants in TBX5 have also been correlated with blood pressure¹⁶ and cardiac arrhythmias¹⁷ in GWAS for adults. Additionally, genotype rs3825214 was related to lone atrial fibrillation in the Chinese Han population,¹⁸ indicating that genetic variations in TBX5 are likely related to organ function defects. The TBX5 rs3825214 genotype from GTEx data regulated alternative splicing to affect heart development (Figure 3 and Supplementary Figure 1). Additionally, RNAfold analysis showed that the TBX5 rs3825215 genotype affected RNA formation energy and stability, so functionally, the TBX5 rs3825215 genotype influenced epigenetics. Our results in combination with existing evidence suggest that genetic variants in TBX5 could result in deficient liver tissue development and enhance vulnerability to HBV-related liver immune injuries and/or the inability to repair damaged cells; thus, enduring injuries are likely attributed to LC, which also increases the risk of HCC. Interestingly, the current study reported that hypermethylation of the TBX5 promoter region was associated with HCC in the Japanese population,¹⁹ posing a potential link between genetic variations and methylation status, an alternative explanation for the influence of genetic variations on HCC development.

Our previous studies demonstrated the two susceptible directions to either persistent infection²⁰ or progression.²¹ In interest of specificity, the association of the risk loci with persistent infection, natural clearance and naïve to HBV infection was extensively evaluated in a well-defined population collected from the same region of a progressive cohort. Given that HBV vaccination was available in mainland China in 1980, an age above 40 years²² and no vaccination history are preconditions for natural exposure to HBV. Based on preconditions, community-based subject enrollment, including vertical and horizontal exposure, is representative of the general population in HBV endemic regions. This exploration confirmed the specific association; rs3825214-GG demonstrated no progressive characteristics in the general population once the rs3825214-GG individual had chronic HBV infection that increased the risk of developing LC or/and HCC.

The risk loci have been related to the clinical severity and HCC development in Asian chronic HBV infection patients,²³ while the relationship of the risk genotypes to HBV-related progressive events has not previously been defined. In our hospital-based follow-up, 167 of 575 patients developed progressive outcomes overall. The exploration of this is significant, as the findings not only confirm that rs3825214 is associated with LC and HCC but also suggests that risk genotypes from other association studies might have an effect on prognosis when no evidence supports biological differences. Accumulating genetic data improve the prediction accuracy of progressive outcomes.

Several limitations of our study need to be mentioned, since GWAS based on the DNA pool was applied as the first step to obtain major genetic predisposition in overview and lost its susceptibility in detail. There are potential biases in statistical methods, and different references for imputation and prioritizing strategies of SNP selection might yield a different association. The mechanisms underlying the association of TBX5 with LC and HCC need to be systematically elucidated. Moreover, in the functional annotation, the condition should keep in mind that the risk genotypes that we tested may not actually be causal but rather correlate the causal genotypes by LD linkage relations; thus, further studies are first needed to fine-map the risk region of TBX5 and clarify the influence of risk genotype(s) or haplotype(s) on the gene regulation of expression and function processes, which in turn provide mechanistic plausibility. The follow-up study is limited to hospital persistent infection patients, which might generate bias.

To the best of our knowledge, our report has extended our understanding of the genetic architecture underlying LC and HCC. Accumulating genetic information in an unbiased manner indicated that susceptibility might be higher; thus, it remains to be determined whether the new high-resolution genotype technology can systematically and specifically reveal the progressive genetic profile. In aggregate, it is important to gain a better understanding of the genetic profile that mediates prognostic differences in addition to clinical risk factors.