{'fraction_non_alphanumeric': 0.03264751748075859, 'fraction_numerical': 0.026434931334574172, 'mean_word_length': 5.365513928914505, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 8, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The influence of host genetics on liver microbiome composition in patients with NAFLD

Summary Background Human body microbiotas are influenced by several factors, including the interaction of the host with the environment and dietary preferences. The role of host genetics in modulating the liver microbiota in the context of NAFLD remains unknown. To address this gap, we examined the interplay between the liver metataxonomic profile and host genetics. Methods We obtained 16S rRNA gene sequences from liver biopsies and genotypes by Taqman-assays in 116 individuals. We compared taxon abundance at the genus level across host genotypes using dominant models of inheritance. We focused the analysis on variants influencing the risk/ protection against NAFLD-histological severity (PNPLA3-rs738409, TM6SF2-rs58542926, MBOAT7-rs641738, and HSD17B13-rs72613567) and a variant influencing macronutrient intake (FGF21-rs838133). We also explored the variants' combined effect via a polygenic risk score (PRS). Findings We identified at least 18 bacterial taxa associated with variants in the selected loci. Members of the Gammaproteobacteria class were significantly enriched in carriers of the rs738409 and rs58542926 risk-alleles, including Enterobacter (fold change [FC]=6.2) and Pseudoalteromonas (FC=2) genera, respectively. Lawsonella (1.6-FC), Prevotella_9 (FC=1.5), and Staphylococcus (FC=1.3) genera were enriched in rs838133-minor allele carriers, which is linked to sugar consumption and carbohydrate intake. Tyzzerella abundance (FC=2.64) exhibited the strongest association (p = 0.0019) with high PRS values (>4 risk alleles). The percentage of genus-level taxa variation explained by the PRS was ∼7.4%, independently of liver steatosis score and obesity. Interpretation We provided evidence that genetic variation may influence the liver microbial DNA composition. These observations may represent potentially actionable mechanisms of disease.


Introduction
Non-alcoholic fatty liver disease (NAFLD) is the leading global cause of chronic liver disease 1 . The disease pathogenesis involves a complex interrelationship among lifestyle (including excessive intake of energy-dense foods and diet rich in highly refined carbohydrates), systemic metabolic dysregulation associated with related comorbidities such as obesity, type 2 diabetes, and cardiovascular risk, 1,2 genetic predisposition involving nuclear and mitochondrial genomes, 3,4 and epigenetic mechanisms. 5,6 The NAFLD natural history is characterized by not only high inter-individual variability but also phenotypic complexity, with histological stages ranging from nonalcoholic fatty liver (NAFL) to steatohepatitis (NASH), NASH-fibrosis, cirrhosis, and even hepatocellular carcinoma 1 .
Owing to the advances in molecular techniques, including high-throughput OMIC approaches, it has been established that a significant proportion of the variability in NAFLD-histologic traits can be explained by a highly interconnected and dynamic network of factors 7 that includes the gut microbiome. [8][9][10][11][12] Seminal studies examining the biogeography of bacterial communities of the human body 13 revealed that tissues' microbial composition plays a vital role in maintaining not only individual organ physiology but also systemic metabolic homeostasis. We recently found that the liver metataxonomic signature may explain not only differences in the NAFLD pathogenic mechanisms and phenotypic histological complexity but also physiological functions of the host. 14 Several lines of evidence suggest that the microbial composition of the body is shaped by diverse factors, including host genetics. 15 Xie and co-workers demonstrated the heritability of many microbial taxa in the gut microbiome, including those associated with diseases. 16 More importantly, a recent study conducted by Kolde showed that genetic attributes and single genetic variant effects on the microbial features are specific to particular body sites. 17 While it is widely accepted that genetic susceptibility plays a significant role in the biology of NAFLD, 4,7 the role of genetics in modulating the liver microbiota composition associated with the disease remains unknown. To explore whether the liver microbiota is shaped by host genetics, we conducted candidate gene association analyses using as a proxy of NAFLD histological severity common variants in loci that impose either risk or protection against the disease (PNPLA3-rs738409, 18,19 TM6SF2-rs58542926, [20][21][22] MBOAT7-rs641738, 23 and HSD17B13-rs72613567 [24][25][26] ). In addition, as a proxy of major dietary modifiers, we used a variant in a locus that has been reproducibly linked to high carbohydrate intake and macronutrient preferences (FGF21-rs838133). 27,28 Methods

Studied population
This study involves a secondary analysis of microbial 16S rRNA reads from the liver of 116 individuals, categorized as non-NAFLD patients (n = 19) and patients with NAFL (n = 44) and NASH (n = 53), as explained

Research in context
Evidence before this study NAFLD is a multifactorial disease that has become the most prevalent chronic liver disease worldwide. The disease natural history is variable among individuals. NAFLD phenotypic complexity is determined by multiple disease-triggering mechanisms, including environmental factors, genetic makeup, gut microbiota, and lifestyle factors. Pioneering studies on the liver tissue metataxonomic profile of patients with NAFLD revealed that the liver contains a diverse repertoire of bacterial DNA that may be involved in the disease severity. Nevertheless, the potential interaction between the liver microbial DNA composition and the host genetics is not entirely understood.

Added value of this study
We explored the interrelationship between the liver metataxonomic profile in patients with NAFLD across the entire spectrum of the disease severity and variants modifying either risk or protection against the disease and variants involved in carbohydrate intake and macronutrient preferences. Together, we provide evidence that the liver microbiota depends in part on the host genetic background. By combining the effects of traitassociated risk alleles, we demonstrated that »7.4% of genus-level taxa variation is explained by a polygenic risk score. This effect was independent of key covariate parameters, including obesity and the liver amount of fat. Furthermore, we explored the functional consequences of taxa linked to genetic variation, and we detected differentially abundant pathways and predicted bacteria-derived metabolites that may serve to understand the mechanism (s) of disease.

Implications of all the available evidence
We present a proof-of-concept study demonstrating that host genetics plays an essential role in modulating the liver microbial DNA composition. Our results show, for instance, that carriers of the PNPLA3-rs738409 risk allele may present an overabundance of DNA derived from members of the Gammaproteobacteria class, some of which are capable of metabolizing alkane-derived carbons to yield primary alcohols. Collectively, we provide a novel perspective in understanding disease biology.
elsewhere. 14 The cohort of patients with NAFLD included moderately obese or obese subjects (BMI < 40 kg/m 2 , n = 47) and severely/ morbidly obese patients (BMI > 40 kg/m 2 ) elected for bariatric surgery (n = 50). By adopting a matched cohort study design, we ensured that the non-NAFLD liver samples in each group matched patients' features. 14 As non-NAFLD liver tissue we included samples from patients without either NAFLD or features of the metabolic syndrome who were selected from patients attending the Liver Unit, whose age and sex matched the NAFLD patients. These patients presented near-normal liver histology in the liver specimens obtained by percutaneous liver biopsy, and the reason for performing a liver biopsy in these subjects was based on the presence of persistently mildly elevated serum liver enzymes activity. In all the non-NAFLD subjects, all causes of common liver disease were ruled out, and they were included in the study if they did not have histological evidence of fatty change; the histological diagnosis of non-NAFLD livers was minimal changes. In the population of morbid obese patients, non-NAFLD subjects were obese patients who also underwent bariatric surgery and had not features of NAFLD demonstrated in the liver biopsy.
Patients were recruited prospectively and were considered for inclusion if they had histopathologic evidence of fatty liver disease, either NAFL or NASH, on liver biopsy performed within the study period. Selection of biological samples included in the present study was based on sufficient amounts of high-quality nucleic acids (DNA) for performing metagenomic explorations. Exclusion criteria: Secondary causes of steatosis, including alcohol abuse ( 30 g for men and 20 g for women, of alcohol daily), total parenteral nutrition, hepatitis B and hepatitis C virus infection, and the use of drugs known to precipitate steatosis were excluded. In addition, patients with any of the following diseases were excluded: autoimmune liver disease, metabolic liver disease, Wilson's disease, and a-1-antitrypsin deficiency. Patients under treatment with antibiotics, immunosuppressive medication or proton-pump bomb inhibitors were also excluded.
Biological specimens from NAFLD-participants and non-NAFLD subjects were consecutively selected during the same study period from the same population of patients attending the participating institutions located in Argentina, ensuring that all shared the same demographic characteristics (occupation, educational level, place of residence, and ethnicity) as the matched patients.
Arterial hypertension was defined as systolic blood pressure 140 mm Hg and/or diastolic blood pressure 90 mm Hg following repeated examination and considering the average of 2 or more blood pressure measurements on at least 2 subsequent visits. 29 All liver specimens were obtained by liver biopsy that was performed before any intervention (if required) with ultrasound guidance using a modified 1.4-mmdiameter Menghini needle (Hepafix, Braun, Germany) under local anaesthesia on an outpatient basis, or during bariatric surgery (surgically excised samples from the left lobe were immediately collected after the abdomen was opened and before organs were manipulated). A portion of each liver biopsy specimen was routinely fixed in 40 g/l formaldehyde (pH 7.4), embedded in paraffin, and stained with haematoxylin and eosin, Masson trichrome, and silver impregnation for reticular fibres. All biopsies were at least 3 cm in length and contained a minimum of eight portal tracts. 14

Ethics
Biological specimens, including blood samples and liver biopsies from all subjects included in this study, were obtained with written, informed consent under the Institutional Review Board-approved protocols (protocol numbers: 104/HGAZ/09, 89/100, 1204/2012, and updated DI-2019-376-GCABA-HGAZ). Protocol # DI-2019-376-GCABA-HGAZ ("Genetics of steatohepatitis and its association with the liver tissue microbiome") was approved under the Argentinean Law No. 3301 on Protection of the Rights of Subjects in Health Research, regulated by Decree No. 58/2011, which establishes the regime for health research activity with human beings. All data were de-identified prior to use. All investigations performed as a part of the present study were conducted under the guidelines of the 1975 Declaration of Helsinki, as revised in 1993.

Host genetics: selection of variants and genotyping strategy
We selected variants in genes that serve as valid instruments for assessing the interrelationship with the liver microbial DNA composition. We focused the variant selection on two axes, including variants influencing the risk and /or protection against severe NAFLD-histological phenotypes and key external factors involved in the disease pathogenesis, such as dietary preferences. Hence, variant selection involved: (i) major genetic modifiers of NAFLD and NASH natural history that have been consistently replicated across ethnicities 4 , including coding single nucleotide polymorphisms in PNPLA3 (rs738409), TM6SF2 (rs58542926), MBOAT7 (rs641738), and the splice donor variant in HSD17B13 (rs72613567), and (ii) the synonymous rs838133 variant in the first exon of FGF21, which is linked to food intake regulation, macronutrient preference (high carbohydrate intake), and central reward pathways. 27,28 The genetic analyses were performed on genomic DNA extracted from white blood cells by a standard method. 30 For rs738409 genotyping, a high-throughput genotyping method involving PCR amplification of genomic DNA with two-tailed allele-specific primers was implemented as previously described. 31 Genotyping of rs58542926, rs641738, rs72613567, and rs838133 was performed using TaqMan genotyping assays (dbSNP rs58542926 assay C__89463510_10, #4351379, dbSNP rs641738 assay C___8716820_10, # 4351379, rs72613567 custom assay #4332072, and dbSNP rs838133 assay C___8832415_10 Cat. # 4351379, respectively, Applied Biosystems, Foster City, CA) according to the manufacturer's instructions, as explained in our previous publications. 22,26,32 To ensure genotyping quality, we included DNA samples as internal controls, hidden samples of known genotype, and negative controls (water). The overall genotype completion rate was 100%. To account for possible population stratification, we used a collection of 13 SNPs at different loci (located in chromosomes 4, 15, 17, 13, 1, and 3) and then analysed the data with the Structure program Version 2 33 as we explained elsewhere. 22,26,32 We found no evidence of stratification in our sample because the cases and the non-NAFLD subjects showed similar Q values and the Structure program assigned a similar distance to clusters, with no further improvement in the fitting model by adding up to four clusters (the greatest ln of likelihood was obtained for K = 1).
In addition, we explored the combined effect/s of variants by examining a polygenic risk score (PRS) calculated as the sum of trait-associated risk alleles across the selected loci. In the case of the HSD17B13 (rs72613567) variant, we considered the deletion allele (-) as the risk allele.

Microbial DNA data collection
Using high-throughput 16S rRNA gene sequencing, we obtained pertinent liver metataxonomic information from the tissue samples of the 116 individuals as explained earlier. 14 In brief, microbial DNA was isolated from fresh liver specimens by a manual protocol, after which extracted bulk DNA samples were amplified with barcoded primers, DNA libraries were constructed and high-throughput sequencing was performed on an Illu-minaMiSeq platform by Macrogen Inc. (NGS Division) Seoul, South Korea. A negative (blank) control was used to test the potential presence of contaminant DNA and/ or cross-contamination; the sample showed no product and did not pass the quality control analysis. 14 The 16S rRNA gene hypervariable region 3 (V3) and V4 amplicons were generated via PCR amplification using primers as reported elsewhere. 14 The 16S rRNA gene sequencing data were filtered, denoised, and processed on the QIIME2 (version 2018.11) platform (http:// qiime2.org/index.html). To account for the generally non-normal distribution of microbial community composition data, 16S rRNA abundances underwent normalization by total sum combined with squared root transformation. High-quality sequences were assigned to operational taxonomic units (OTUs) using the QIIME pipeline 34 whereby default parameters were used in the selection of OTUs for constructing the OTU table. High-quality amplicon sequence variants were classified using the vsearch algorithm with default parameters and SILVA 16S-only 99% identity database (release 137) to build our BIOM feature table of OTUs. Assessment of differences in the liver bacterial DNA composition across the entire disease severity spectrum was further confirmed by an independent molecular approach as explained earlier. 14

Functional profiling of liver microbial DNA
To predict metagenome functional content from the OTU table, we utilized the bioinformatics software package PICRUSt2: Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (https://github.com/picrust/picrust2/). To additionally detect differentially abundant KEGG ortholog pathways, we used the Piphillin metagenomics inference tool. 35 The normalized sequence abundance table and the weighted nearest-sequenced taxon index values persample were used to predict pathways. The resulting output was a list of MetaCyc pathway abundances according to the MetaCyc database (https://metacyc. org/). The STAMP software package (version 2.1.3), which relies on the concept of biological relevance in the form of confidence intervals, was employed to determine differentially enriched metabolic pathways (p < 0.05) and their effect sizes (h 2 ) (http://kiwi.cs.dal.ca/ Software/STAMP). 36 In addition, we reported p-values corrected through the Benjamini-Hochberg procedure for multiple testing corrections.

Statistical analysis
We compared taxon abundances at the genus level across genotypes assuming a dominant model of inheritance and comparing YY with XY + XX, in which X is the risk allele. The use of a dominant model for all tested variants is justified as we performed pairwise comparisons of the relative abundance of each taxon. Differences between genotype groups were assessed using a pairwise Wilcoxon rank-sum test implemented in the R-base function and the default method for pvalue correction using the web-application Calypso (FDR < 0.05). 37 To identify the potential discriminatory power of individual bacterial DNA taxa over PRS, we used the area under the receiver operating characteristic (ROC) curve (AUC) and fold change. To adjust for relevant covariates (sex, obesity degree, and liver disease severity), we used non-parametric ROC estimation using bootstrap to obtain standard errors and confidence intervals along with linear regression to model the covariates.

Articles
We calculated the fold change in the abundance of taxon of interest as the ratio of the mean of sequencing reads assigned to the taxon normalized by the total number of sequencing reads combining with square root transformation in each group. Measurements across genotype groups were compared by non-parametric statistical tests (Wilcoxon rank test), with p < 0.05 indicating statistical significance. As indicated in some cases, pair-wise comparisons between genotype groups were performed through Student's t-test, and relative abundance was calculated as the square root of each feature read count divided by the total number of features reads in each sample, and the standard error (SE) is indicated by error bars.
The R-squared in ANOVA estimated the genetic contribution to the variance in the abundance of different genera.
We used ordinary linear regression to study the correlation between bacterial DNA abundance as dependent variable and independent variables, in particular, PRS as a continuous or dichotomic variable and relevant covariates such as sex (0=female, 1=male), steatosis score (0, 1, 2, and 3 grades), obesity degree (obese/ overweight vs. severe obesity) and liver disease severity (non-NAFLD, NAFL, and NASH). After fit the model, margins subroutine was used to estimate predicted media at fixed values of some covariates as implemented in the STATA v16.0 package (StataCorp LLC, College Station, Texas, USA).
For the assessment of the association between blood octane levels and liver fibrosis in a population-based study, we used the svyset (survey design for dataset) command (STATA 16.0). To avoid skewed distributions, continuous variables were log-transformed (Log 10 ). Statistical analyses were performed by linear regression among blood octane concentration as the dependent variable and quantiles of the median stiffness as the independent variable. Age (years), gender (women= 0, men= 1), body mass index (kg/m 2 ), diabetes (no = 0, yes = 1), and alcohol consumption (g per day) were also included as covariates.
For the comparison of clinical, biochemical, and histological characteristics we used Mann-Whitney U test, except for female/male proportion between studied groups that was assessed by a Chi-square test.
Sample size and power calculations indicate that a sample with »50 participants per genotype group has a power of »0.80 to detect media differences representing »2.0 fold-changes with a probability of <0.05 for type I error.

Role of funders
The funders had no role in the conceptualization, study design, data collection, analysis, interpretation of data, in writing the paper, or in the decision to submit the paper for publication.

Results
Host genetics and hepatic microbial DNA associations Detailed clinical, biochemical, and histological characteristics of the study participants is shown in Table 1. Briefly, the sample comprised of 55.17% women, with a mean age of 46 § 2.3 years. Genotype frequencies were in Hardy-Weinberg equilibrium ( Table 2). Visualization of the hierarchical relationship of liver microbial DNA profile in the whole population for the 950 most abundant taxa is depicted in Figure 1.
We next tested for the association between the relative abundance of bacterial DNA taxa and minor alleles of the selected variants. We found that two genera were significantly enriched in carriers of the PNPLA3-rs738409 G-risk allele, including Enterobacter and Marivita.
Genera Butyricicoccus and Streptococcus were significantly depleted, and Tyzzerella was increased in carriers of the MBOAT7-rs641738 T-risk allele. Among carriers of the protective A-INS allele of the HSD17B13-rs72613567 variant, the abundance of genera Fusobacterium and Parasutterella was significantly decreased while that of the Methylotenera genus was enriched.
Regarding the variant influencing dietary preferences, we found that the abundance of genera Lawsonella, Prevotella_9, and Staphylococcus was higher, and Lactobacillus quantity was diminished in carriers of the FGF21-rs838133 minor A-allele.
The direction of effects, including fold changes across all variants, is shown in Figure 2.
Combined effect of risk alleles on the liver microbial DNA composition: Taxon abundance and variance explained by the host genetics We further explored the effect of the combined effect of NAFLD-risk/protection alleles on the composition of the liver microbiota by a PRS. We first investigated whether the genetic risk of NAFLD was associated with alpha (within sample) diversity, and we found no detectable association between PRS for NAFLD and Shannon index (p = 0. 93) or Simpson index (p = 0.8).
Then, using univariate analysis, we compared the taxon abundances across two groups of individuals, including subjects carrying 4 risk alleles and subjects carrying >4 risk alleles, respectively. As shown in Figure 3a, among all the microbial DNA taxa present in the liver, certain taxa showed significant associations with the PRS. Specifically, the abundance of the Tyzzerella genus-a member of the Firmicutes phylum and the Clostridia class-showed the strongest association with high PRS values (>4 risk alleles), as indicated by 2.64-fold change differential abundance (p = 0.0019, FDR < 0.05) and shown in Figure 3b. Conversely, the Lactobacillus genus-a member of the Firmicutes phylum and Bacilli class-exhibited the strongest association with low PRS values (4 risk alleles) (0.89-fold change, p = 0.033, FDR < 0.05, Figure 3b).
Next, we calculated the percentage of the liver microbial taxa variance that could be explained by the PRS, which was built with proxy-based variants of NAFLD histological severity and macronutrient intake. The proportion of genus-level taxa variance explained by the host genetics (PRS) was »7.4% (adjusted R 2 = 0.0738, p < 0.005). This effect was independent of the presence of fatty liver and the obesity status of the studied population. Figure 3c shows the predicted marginal means of Tyzzerella abundance by PRS according to subjects' liver steatosis scores in obese/overweight and severely obese groups and across the two main stages of the disease severity, namely NAFL and NASH. The Margins algorithm performed the estimation after linear regression with Tyzzerella abundance and PRS as dependent and independent continuous variables, respectively, and sex, steatosis score, and BMI-classified groups as covariates. We found that correlation between Tyzzerella abundance and the number of risk alleles (PRS) was only significant in NASH (Coefficient: 5.42, 95% CI: 0.92 -9.93, p value = 0.02) as opposed to in NAFL (Coefficient: 2.89, 95% CI: -2.97 -8.76, p value = 0.32) patients.
Because of previous observations in well-characterized cohorts addressing the relationship between the gut microbiota composition and cardiovascular disease (CVD) risk demonstrated a potential link between Tyzzerella and Tyzzerella 4 abundance and high CVD risk profile, 38 we explored the possible association between Tyzzerella enrichment and the presence of arterial hypertension. We found a trend toward a relative increased abundance of Tyzzerella in the liver of subjects with arterial hypertension (mean: 13.85 § SE: 2.85, n = 53) compared to non-hypertensive individuals (mean: 10.03 § SE: 2.80, n = 60) (p value = 0.056 [Mann-Whitney test]). However, the association was not significant after adjustment for the relevant covariates.

Functional annotation enrichment shows distinctive predicted functional signatures
After demonstrating the host genetic impact on liver microbial taxa, we examined genetic contributions to liver microbial functions. We found significantly enriched specialized pathways in carriers of the Bars depict fold changes of liver bacterial DNA taxa abundance (whereby >1 indicates enrichment and <1 depletion) in carriers of variants in target loci, including NAFLD/NASH-risk loci (PNPLA3 rs738409 CC n=31, CG+GG n=85), TM6SF2 rs58542926 CC n=89, CT+TT n=27, MBOAT7 rs641738 CC n=44, CT+TT n=72, and HSD17B13 rs72613567 -/-n=83, -A+AA n=31), and a locus influencing dietary macronutrient intake (FGF21 rs838133 GG n=61, GA+AA n=55). Measurements across sample groups are compared by nonparametric statistical tests [Wilcoxon rank test], with p < 0.05 indicating statistical significance.
According to the NCBI taxonomy database, taxa are shown at the genus level (p_ stands for phyla and c_ for class). We calculated the relative abundance (fold change) of the taxon as the ratio of the media per group of sequencing reads assigned to the taxon normalized by the total number of sequencing reads combining with square root transformation. 37 PNPLA3-rs738409 G-risk allele. The most remarkable ones include octane oxidation and fatty acid salvage, also known as b-oxidation degradation de novo fatty acid biosynthesis shunt (Figure 4a). As shown in Figure 4b, among the pathways associated with the TM6SF2-rs58542926 T-risk allele, we found L-tryptophan degradation to 2-amino-3-carboxymuconate semialdehyde and NAD biosynthesis II from tryptophan (one of the essential coenzymes in redox reactions in the cell). In carriers of the HSD17B13-rs72613567 A-INS allele, we found that phosphatidylglycerol biosynthesis I and gondoate biosynthesis microbial pathways were under-represented and the super pathway of adenosylcobalamin salvage from cobinamide I was over-represented ( Figure 4c). Specialized pathways observed in carriers of the MBOAT7-rs641738 T-risk allele is shown in Figure 4d, indicating that the superclass of nucleoside and nucleotide and purine biosynthesis serving as generating cellular energy are enriched. Finally, predicted functionality associated with the PRS showed a trend toward the association with lipopolysaccharide biosynthesis (Figure 3d).

Discussion
Our results indicate that host genetic determinants of the risk/protection for NAFLD/NASH and dietary macronutrient intake are related to the liver microbial DNA b. Bar chart shows significantly different taxa according to the polygenic risk score (PRS) (p < 0.05, [ANOVA]). Relative abundance is calculated as the square root of each feature read count divided by the total number of features reads in each sample. The standard error is indicated by error bars. Pair-wise comparisons are performed through t-test and are annotated as *: p < 0.05, **: p < 0.01, and ***: p < 0.001.
c. Plots show the effects of host genetics (PRS, the number of risk allele carried) on Tyzzerella DNA abundance [predicted by linear regression and Margins subroutine] in NAFLD patients with steatosis scores in the 1À3 range, comparing those who are overweight/moderately obese vs. severely obese patients and after adjusting for confounding factors (age and gender, in addition to steatosis score and obesity degree); the analysis was performed in patients stratified by liver disease severity, separately (NAFL in the left panels vs. NASH in the right panels).
d. Predicted functionality associated with the PRS. The functional inference was explored based on KEGG KO pathways by applying the Piphillin metagenomics inference tool 35 . Bars indicate mean proportions (%), confidence intervals, and their associated pvalues for the top-ranked metabolic pathways corresponding to both groups with 4 PRS > 4. STAMP was used to determine differentially enriched metabolic pathways (corrected p < 0.05) and their effect sizes (h 2 ). The Welch's test was performed to assess statistical significance and corrected p values were calculated using Storey's FDR approach.
composition. The strongest associations were found between members of the Gammaproteobacteria class and NAFLD/NASH risk alleles of the rs738409 and rs58542926 variants, including Enterobacter and Pseudoalteromonas genera, respectively. These findings are consistent with our previous observation indicating that DNA derived from Gammaproteobacteria was associated with more severe forms of liver disease, including high scores of lobular and portal inflammations. 14 It is also known that the body microbiome is sensitive to the host's diet and lifestyle. In this study, we found that carriers of the FGF21-rs838133-minor allele presented higher abundances of Lawsonella (a member of the Actinobacteria class), Prevotella_9, and Staphylococcus genera, and a decreased abundance of Lactobacillus. This finding suggests that dietary sugar consumption and carbohydrate intake favor certain bacterial taxa, such as Prevotella_9 that has been associated with the risk of major common diseases, including cardiovascular disease, 38 colon cancer, 39 and even immunologic diseases such as rheumatoid arthritis. 40 Most importantly, a PRS generated to assess the combined effect of NAFLD/NASH associated variants' risk alleles and the FGF21-rs838133 minor allele revealed that the abundance of the Tyzzerella taxa showed the strongest association with high PRS values (>4 risk alleles). Findings yielded by the analyses of the gut microbiome of the participants of the Bogalusa Heart Study-a long-term epidemiologic study investigating the natural history of atherosclerosis-indicate that the genus Tyzzerella is associated with lifetime cardiovascular risk profile. 38 Collectively, our findings might serve to explain, at least in part, the link between NASH/ NAFLD and CVD. 41,42 Although several factors influence the human (and mammalian) body microbiota composition, it appears that host genetic factors account for a considerable proportion of variance. For example, seminal studies in twins showed that microbial gut composition is »40% heritable. 15,16 Recent evidence that leveraged 16,234 gut microbiome profiles collected over 14 years from 585 wild baboons suggests that, after controlling for diet, age, and socioecological variation, 97% of the gut microbiome phenotypes are significantly heritable. 43 Of note, the percentage of genus-level taxa variation explained by our NAFLD/NASH PRS was »7.4%, which was independent of key covariate parameters that may affect this effect, including the steatosis score and obesity degree. Even age and gender were not contributing factors. Considering that the PRS based on selected variants in few genes was significant, it is justified to consider that the observed variance is higher than expected by chance. It should be highlighted that host genetics is also relevant in modelling hostÀgut microbe interactions in health and disease 44 , including NAFLD in obese youth as reported recently. 45 Finally, we aimed to detect differentially abundant pathways and predicted bacteria-derived metabolites from the taxa associated with the tested variants. We reasoned that exploring putative functional consequences of taxa linked to genetic variation may aid in better understanding of the disease mechanism(s). We found multiple altered microbial metabolic pathways associated with host genetics. For example, in carriers of the rs738409 G-risk allele, octane oxidation and the novo fatty acid biosynthesis shunt were over-represented. The key process in octane oxidation is the alkane hydroxylase system that introduces molecular oxygen in the C1 atom of the hydrocarbons at the expense of NADH to yield primary alcohols. 46 Previous studies demonstrated that several alkane-degrading bacteria can use diverse compounds as a carbon source in addition to alkanes, 47 which are further oxidated to fatty acids via bacterial b-oxidation pathway (BioCyc ID: P221-PWY). While the alkB domain has been found in many bacterial genomes, including Proteobacteria, Actinobacteria, and Bacteroidetes, it appears that the alkane hydroxylase system is particularly relevant in the Gammaproteobacteria class, including those from Pseudomonas aeruginosa. 48 Volatile organic compounds, including octane, were found either in the exhalation or faeces of human subjects with diverse medical conditions associated with oxidative stress and chronic inflammation, including lung cancer, 49 obstructive sleep apnoea, 50 gastrointestinal diseases, 51 and NAFLD. 52 Notably, data analysis based on NHANES 2017À2018 population-based survey of non-occupationally exposed individuals that included liver fibrosis assessment via transient elastography (FibroScan) examination shows a significant association between blood octane concentrations and liver stiffness (beta 0.0016 § 0.0007, p = 0.031 adjusted by age, gender, diabetes, obesity [BMI], and alcohol consumption) (Supplementary Figure 1).
Collectively, primary alcohols can be synthesized either via fatty acid or amino acid pathways in diverse bacteria. The link between the gut microbiota and alcohol (ethanol) production has been described by several authors. 53,54 However, our study provides novel evidence regarding the association between PNPLA3-rs738409 and taxa capable of metabolizing alkanederived carbons to yield primary alcohols.
Likewise, the rs738409 G-risk allele was found to be associated with a metabolic pathway involved in b-oxidation degradation, also known as de novo fatty acid biosynthesis shunt. Fatty acid biosynthesis shunt is initiated by the b-ketoacyl-[acyl-carrier-protein] synthase that catalyses the condensation of acetyl-CoA with malonyl-[acp], producing acetoacetyl-[acp] and that can be transformed into longer fatty acids via the fatty acid elongation cycle. This might be a new mechanism to explain why carriers of the PNPLA3 variant are exposed to 3.28-fold higher risk of liver fat accumulation. 18,19 The comparison of functional prediction between carriers and non-carriers of the HSD17B13-72613567 protective A-INS allele yielded results that might expand our knowledge of the protective function of this variant against the most detrimental NAFLD-histological outcomes. 24,26 Specifically, we found over-representation of the super pathway of adenosylcobalamin salvage from cobinamide I-a pathway involved in vitamin and cofactor biosynthesis-in A-INS allele carriers, along with a lower proportion of taxa linked to fatty acid metabolism, such as phosphatidylglycerol biosynthesis and components of the gondoate pathway, which is essential for the production of unsaturated fatty acids, including arachidonate. Arachidonate production is the rate-limiting step in the synthesis of prostaglandins and is associated with NAFLD 55 and also participates in the Land's cycle.
In addition, the predicted functional consequence/s of associated taxa with the PRS suggested that carriers of >4 risk alleles might be exposed to some genera with high potential for lipopolysaccharide biosynthesis. This finding is consistent with our earlier observation that lipopolysaccharides derived from Gram-negative bacteria localize in the portal tract of patients with more severe disease. 14 However, when interpreting these findings, some study limitations should be considered. In particular, our results are based on the examination of the bacterial DNA profile of liver specimens in association with a few genetic variants. Therefore, much greater (by many orders of magnitude) sample sizes will be necessary to fully explore the influence of the whole genome on the liver microbial DNA composition. Nonetheless, we were able to overcome several challenges, including explorations of microbiome profiles in liver tissue, which is very difficult to obtain and analyze for relevant data.
It is also plausible that different results may arise in other populations, as the inter-relationship between genetics and the liver microbiome can be modified by environmental factors not explored in the present study. These potential differences can be explained by the presence of population structure across individuals from different parts of the world, as well as by different biogeography of bacterial taxa in different parts of the body that are determined not only by host genetics, but are also modified by factors that impose selective pressure on the body microbiotas, such as antibiotic usage, travel, diet, microbes in the soil, vertical transmission from mothers to children, and diverse causes of horizontal transmission between humans and/or animals. 56 In conclusion, we studied the intricate link between host genetics and the liver microbiome composition in patients with NAFLD. We were able to show the most significant liver bacterial DNA that may be affected by host variation in candidate genes relevant to the disease biology. The most striking result to emerge from our study is that host genetics may shape the metabolic and biological environment in which the NAFLD-associated liver tissue-microbial DNA resides.

Declaration of interests
Nothing to declare.

Contributors
All authors (CJP; AS; MFQ; GOC; MG; and SS) read and approved the final version of the manuscript. CJP: study concept and design; biological material collection, data acquisition; data analysis and interpretation; genetic studies; general study supervision; statistical analysis; manuscript drafting; and securing funding. AS: bioinformatics analysis of sequencing data and data processing and presentation. MFQ: sample processing. GOC and MG: performed liver biopsies and collected biological samples. SS: study concept and design; data acquisition; performed liver biopsies and collected biological material; histological evaluation; genetic studies; data analysis and interpretation; general study supervision; manuscript drafting; securing funding. All authors have read and approved the final version of the manuscript. In addition, CJP, GOC, and SS have verified the underlying data.