The GWAS summary statistics were downloaded from the website of the Japanese ENcyclopedia of GEnetic associations by Riken,¹ including the datasets of chronic hepatitis B and chronic hepatitis C. The GWAS summary statistics of the Pan UK Biobank were downloaded from https://pan.ukbb.broadinstitute.org/, including two datasets: viral hepatitis B and viral hepatitis C (Supplementary Table 1).

Quality control (QC) filtering was performed on the summary statistics data of each cohort. Any variant with low confidence or minor allele frequency (MAF) < 0.01 was removed from further analysis. The number of variants before and after QC is shown in Supplementary Table 1.

After data QC, meta-analyses based on the above GWAS summary statistics of hepatitis infection were performed for each population. The Z-score-based method implemented in the software METAL (Willer et al., 2010) was adopted with sample overlap correction. This approach computed a weighted sum of Z-score for each effect allele of each single-nucleotide polymorphism (SNP) across cohorts, with weights defined by the square root of the effective sample size for each cohort. It also takes the direction of effect into account (Willer et al., 2010).

The associated SNPs were clumped into linkage disequilibrium (LD)-independent genomic regions using the software PLINK (Purcell et al., 2007) with the following parameter settings: –clump-r² = 0.4, –clump-kb = 500, –clump-p1 = 5e-08, and –clump-p² = 5e-02. LD information from the 1,000 Genomes phase 3 data (Auton et al., 2015) of European and East Asian ancestry was used in the analysis. The MHC region (chr6: 25–35 MB) was considered as one locus. The resulting independent genomic loci were compared to those reported in the GWAS catalog for hepatitis C and hepatitis B diseases. Any locus within 1 MB of a reported locus was regarded as a known locus. By excluding known loci, the remaining genome-wide significant (GWS) loci were considered as novel association loci.

Fine mapping via the software FINEMAP (Benner et al., 2016) was conducted to identify variants in the 95% credible set of each novel GWS locus, with LD score calculation based on the EUR and EAS samples of the 1,000 Genome phase3 panel. The maximum number of allowed causal SNPs is set as 10 in the analysis. Causal variants were defined as those with r² > 0.6 with the index SNP and the association P-value < 1 × 10^–4 in the meta-analysis.

Functional annotation of each GWS locus was performed with FUMA GWAS (Watanabe et al., 2017). Annotation of the expression quantitative trait loci (eQTLs) was based on scRNA eQTLs (van der Wijst et al., 2018), DICE (Schmiedel et al., 2018), eQTLGen, blood eQTLs (Westra et al., 2013), and GTEx v8 (Battle et al., 2017). The threshold for significant eQTL SNP–gene pairs was defined as a false discovery rate (FDR) < 0.05. To annotate chromatin interactions, the liver Hi-C datasets were analyzed by using the 3D Genome Browser and 3D-genome Interaction Viewer (Lieberman-Aiden et al., 2009; Wang et al., 2018; Yang et al., 2018).

The microarray dataset GSE49824 was downloaded from the GEO database (Huebener et al., 2014). The dataset contains gene expression profiling data from mouse hepatic stellate cells from four mice with hepatocyte-specific deletion of Hmgb1 and four control mice. The R package GEOquery was used for data download (Davis and Meltzer, 2007). Calibration, quantile normalization, and differential expression analysis were conducted with the R package limma (Ritchie et al., 2015). The significant genes were defined as log₂| fold change (FC)| > 0.5 and P-value < 0.05.

Count data were downloaded from the GEO dataset GSE148301 (Xu et al., 2021). Then, the variance-stabilizing transformation (vst) normalization and differential expression were performed via the R package DESeq2 (Love et al., 2014). The significant genes were defined as log2| fold change (FC)| > 1 and P_adj < 0.05.

The online web portal DAVID (Huang et al., 2007; Huang da et al., 2009b) was used for pathway enrichment analysis with differentially expressed genes as input. The collection of GO and KEGG pathway datasets was used to find significantly enriched pathways (Huang da et al., 2009a).

The protein–protein interaction (PPI) network was predicted using the Search Tool for the Retrieval of Interacting Genes (STRING)² online database (Franceschini et al., 2013), removing interactions with a combined score < 0.4. Then, the PPI network was drawn using Cytoscape (version 3.7.2). The significant interactions were identified using Degree (Smoot et al., 2011). Finally, DAVID was used to conduct pathway enrichment analysis on the 19 genes that were most closely related in the network (Huang da et al., 2009a).

We aimed to understand the common pathogenesis mechanism of HBV and HCV infection, to aid in protection from harmful coinfection of these viruses. We carried out a meta-analysis between the GWAS results of HBV and HCV within populations of European ancestry and Asian ancestry, respectively (Supplementary Table 1). In the Asian population, the HLA locus on chromosome 6 showed strong association and another previously reported locus on chromosome 19 at the 5’ of IFNL4 also reached GWS (Figure 1A). In addition, we identified a novel GWS locus at 13q12.3, located in 5’ of gene MEDAG (Figures 1A, 2A). From the meta-analysis of the GWAS summary statistics of HBV and HCV in the Pan UK Biobank, we found two novel GWS loci at 1q32.3 and 5p13.3 (P-value < 5 × 10^–8). These loci are in the intron of gene INTS7 and 3’ of gene CDH6, respectively (Table 1 and Figures 1B, 2B,C). INTS7 encodes integrator complex subunit 7, and CDH6 is a type II cadherin. Looking into the association of the top SNP at the novel GWS loci with HBV and HCV infections separately, we found that the top SNP at each novel locus was associated at nominal significance level (P-value < 0.05) in the original GWAS statistics and showed consistent direction of effects (Table 1). They reached GWS level in the meta-analysis due to the synergistic contribution from both disease cohorts instead of being driven by HBV or HCV specifically.

Next, we tried to understand how the novel loci in each population contribute to the development of HBV/HCV infection by performing functional annotation of SNPs at the novel loci and examining potential target genes. It is known that the causal variant at each locus may not be the lead SNP, but a SNP in LD. Therefore, we performed fine mapping to identify likely causal variants satisfying the criteria of LD, association in meta-analysis, and Bayesian statistics posterior probability (Table 2). The results show that lead SNPs rs142568263 at 1q32.3 and rs76561592 at 5p13.3 are in the credible set of potential causal variants. Additionally, four SNPs and two SNPs are in the credible set of these two loci, respectively, each with a posterior probability of 1. For locus 13q12.3, there are only three SNPs with a posterior probability > 0.2. Nine of these SNPs overlap with histone modification marks, such as H3K4me1, H3K4me3, H3K27ac, and H3K9ac in liver tissue, suggesting their potential role as cis-regulatory elements regulating target gene expression (Supplementary Table 3).

The SNPs at each GWAS locus could regulate the expression of distant genes through chromatin interaction, in addition to the regulation of its closest gene (Smemo et al., 2014; Claussnitzer et al., 2015; Xia et al., 2016). To identify target genes of each novel association locus, we examined eQTL genes in liver or immune cell types. Gene LPGAT1 is an eQTL gene for SNP rs150857619 at locus 1q32.3 with FDR < 0.05 (Supplementary Table 4). LPGAT1 encodes lysophosphatidylglycerol acyltransferase 1, which may be involved in triacylglycerol synthesis and secretion in the liver. Mouse Lpgat1 knockout resulted in phosphatidylglycerol remodeling defects and led to oxidative stress, mitochondrial DNA depletion, and mitochondrial dysfunction. The Lpgat1^–/– mice developed diet-induced obesity and hepatopathy phenotypes (Zhang et al., 2019a). These results indicate the role of LPGAT1 in mitophagy, which is an important type of selective autophagy to maintain cellular homeostasis upon virus invasion. The expression level of additional genes was found to have nominal significant association with the top associated SNP or causal variants at each novel GWS locus (P-value < 0.05) (Supplementary Table 4).

We also examined genes having chromatin interactions with these loci in Hi-C datasets (Supplementary Table 5). The Hi-C data suggest that there is a chromatin interaction between the HBV/HCV susceptibility loci and the promoter region of multiple genes, which have been reported to be involved in autophagy and immune signaling pathways. For example, the target genes of locus 5p13.3: CDH6 inhibits autophagy and promotes re-organization of the mitochondrial network (Gugnoni et al., 2017); Drosha is involved in regulating hepatocyte growth (Fabris et al., 2019) and HBV replication control (Yang et al., 2017). GOLPH3 is another targeted gene of locus 5p13.3 functioning in selective autophagy (Lu et al., 2020) and hepatitis C virus secretion (Bishe et al., 2012). Similarly, the target genes at locus 1q32.3: TRAF3IP3 enhances autophagy through interaction with an ATG16L1-binding motif (Peng et al., 2015), and ATF3 is a negative regulator of autophagy (Sood et al., 2017). Another gene of particular interest is HMGB1, which is in the same topologically associated domain (TAD) as SNP rs34780238 at 13q12.3 (Figure 3 and Supplementary Figure 1). Twenty three out of 54 Hi-C genes at these three loci exhibited nominally significant eQTL association with the three GWS loci, including the above autophagy-related genes CDH6, ATF3, and HMGB1 (Supplementary Table 5).

To further understand the function and role of HMGB1 in the HBV/HCV infection, we analyzed a microarray dataset from the liver of Hmgb1 knockout mice. Our results yielded 194 differentially expressed genes, including 82 upregulated genes and 112 downregulated genes(log2| FC| > 0.5, P-value < 0.05) (Figure 4A and Supplementary Table 6). Consistent with the original report, we did not observe significant changes in expression levels of genes involved in autophagosome formation. However, the key gene Atp6v0d2 functioning at the step of fusion between autophagosome and lysosome showed increased expression in Hmgb1 KO mice (Figure 4B) (FC = 1.58, P-value = 9.11E-04). We further conducted pathway enrichment analysis to examine the signaling pathways affected by Hmgb1 deletion. Five pathways related to virus infection and immune response were significantly enriched among the 194 differentially expressed genes (DEGs), including “innate immune response” and “response to virus” pathways (Figure 4C and Supplementary Table 7). We also constructed a PPI network analysis among the DEGs, revealing 19 significantly associated genes (Figure 4D). Pathway analysis of the significantly associated genes demonstrated significant enrichment of “immune system process”, “defense response to virus”, “negative regulation of viral genome replication” and other immune function-related gene sets (Supplementary Table 8). These findings further suggest that Hmgb1 may regulate hepatitis virus infection.

We also looked into the function of a potential target gene ATF3 at locus 1q32.3, which exhibited chromatin interaction with GWS SNPs at this locus and showed nominally significant association with these SNPs (Supplementary Tables 4, 5). ATF3 is a negative regulator of autophagy and cellular antiviral signaling by interacting with the promoter of Atg5, Stat1, Irf9, and Isg15 (Sood et al., 2017). A study also showed that ATF3 promotes hepatitis B virus x mRNA degradation by stimulating Ski2 expression (Shiromoto et al., 2018). We examined a GEO dataset from C57BL/6J mice injected with adeno-associated virus (AAV) expressing human ATF3 (Xu et al., 2021). We observed the significant upregulation of 389 genes and downregulation of 315 genes in mouse liver tissue (Figure 5A and Supplementary Table 9). Interestingly, the autophagy gene Atg16l2 showed significant upregulation (logFC = 2.89, adj. P-value = 8.07E-75). Atg16l2 is an isoform of Atg16l and can interact with Atg12-conjugated Atg5 to form a complex in vivo. Atg16l2 can function as an E3-like enzyme to catalyze the conjugation of LC3 to phosphatidylethanolamine (Ishibashi et al., 2011). The DEGs were enriched in metabolic signaling pathways including mitophagy-related pathways, such as oxidoreductase activity, mitochondrion morphogenesis, and endoplasmic reticulum, and they interconnect with each other (Supplementary Table 10 and Figure 5B). Consistently, we also found enrichment of the organelle membrane pathway, which is involved in diverse biological functions including autophagy (Lin et al., 2021). It is known that viruses may trigger and modulate mitophagy to affect the innate immune response of the host. ATF3 ChIP-seq data from the ENCyclopedia of DNA Elements (ENCODE) showed a strong signal at the genomic region including gene ATG16L2, suggesting the likely regulation of ATG16L2 by ATF3 (Figure 5C).

To validate the GWS loci identified in our study, we conducted a replication study with data extracted from the biobank of CAG, CHOP. The replication cohort of European ancestry after QC filtering includes 67 cases with the diagnosis of HBV or HCV and 649 controls (Supplementary Tables 2, 11). Both the lead SNP rs142568263 at locus 1q32.3 and the top SNP rs76561592 at locus 5p13.3 are associated with HBV/HCV infection in the replication cohort (rs142568263 P-value = 0.00874, rs76561592 P-value = 0.0139) (Supplementary Table 12); thus, the two loci were replicated with a P-value < 0.025. The number of genotyped HBV/HCV cases of East Asian ancestry that we could find is too small for a replication study; however, the peak SNP rs34780238 at locus 13q12.3 is associated with phenotype HBc antigen for hepatitis B virus (phenotype code 23016, P-value = 7.757E-5) in UKBB reported in Global Biobank Engine (Stanford, CA)⁵ (McInnes et al., 2019). Variants mapped to HMGB1 are associated with multiple blood traits, like mean corpuscular volume and eosinophil counts at GWS level, as reported in literature (Astle et al., 2016; Chen et al., 2020). The results from the replication cohort, the supporting evidence from literature, and the known biological function of the candidate genes at each locus suggest that these loci are likely to have bona fide HBV/HCV associations.