Close association of PFASs exposure with hepatic fibrosis than steatosis: evidences from NHANES 2017–2018

Abstract

Multiple animals and in vitro studies have demonstrated that perfluoroalkyl and polyfluoroalkyl substances (PFASs) exposure causes liver damage associated with fat metabolism. However, it is lack of population evidence for the correlation between PFAS exposure and nonalcoholic fatty liver disease (NAFLD). A cross-sectional analysis was performed of 1150 participants aged over 20 from the US. Liver ultrasound transient elastography was to identify the participants with NAFLD and multiple biomarkers were the indicators for hepatic steatosis and hepatic fibrosis. Logistics regression and restricted cubic splines models were used to estimate the association between PFASs and NAFLD. PFASs had not a significant association with NAFLD after adjustment. The hepatic steatosis indicators including fatty liver index, NAFLD liver fat score, and Framingham steatosis index were almost not significantly correlated with PFASs exposure respectively. But fibrosis indicators including fibrosis-4 index (FIB-4), NAFLD fibrosis score, and Hepamet fibrosis score were positively correlated with each type of PFASs exposure. After adjustment by gender, age, race, education, and poverty income rate, there was also a significant correlation between PFOS and FIB-4 with 0.07 (0.01, 0.13). The mixed PFASs were associated with FIB-4, with PFOS contributing the most (PIP = 1.000) by the Bayesian kernel machine regression model. The results suggested PFASs exposure appeared to be more closely associated with hepatic fibrosis than steatosis, and PFOS might be the main cause of PFASs associated with hepatic fibrosis.

Key messages

1. Current exposure doses of PFAS did not significantly change the risk of developing NAFLD.

2. PFASs exposure appeared to be more closely associated with hepatic fibrosis than steatosis.

3. PFOS might be the main cause of PFASs associated with hepatic fibrosis.

Keywords:

• Nonalcoholic fatty liver disease (NAFLD)
• Perfluoroalkyl and polyfluoroalkyl substances (PFASs)
• Generalized linear model (GLM)
• Restricted cubic spline (RCS)
• Bayesian kernel machine regression (BKMR)

Introduction

Due to the rapid development of the economy and the change of lifestyle, nonalcoholic fatty liver disease (NAFLD) has become the main cause of liver disease in the world, which may develop into steatohepatitis, cirrhosis and even liver cancer, threatening human health [1]. The estimated global prevalence of NAFLD was about 30% in recent years, and it has gradually increased and become a major concern [2,3]. Exposure to environmental pollutants is also an important cause of NAFLD as the liver is the most important metabolic organ.

Perfluoroalkyl and polyfluoroalkyl substances (PFASs) are not only widely found in food packaging [4], nonstick cookware [5], outdoor apparel [6], and other consumer products because of their chemical inertness and high-temperature resistance, but also in surface water [7] and aquatic biota [8] due to their strong resistance to degradation. The Centers for Disease Control and Prevention of the U.S. added them to the routine health testing items decades ago, and 98% of U.S. adults had residual PFASs in their blood in 2014 [9].

PFAS exposure and NAFLD may have a link. A cross-sectional study from National Health and Nutrition Examination Survey (NHANES) (2011 ∼ 2014) on US adults, pointed out that perfluorooctane sulfonic acid (PFOS), perfluorohexane sulfonic acid (PFHxS) and perfluorononanoic acid (PFNA) serum concentrations were positively associated with alanine aminotransferase (ALT) and gamma-glutamyl transferase (GGT) in obese participants aged ≥20 years [10]. Another study with 200 participants aged between 40 to 70 from six water districts were contaminated with perfluorooctanoic acid (PFOA) in the mid-Ohio Valley of Ohio and West Virginia from 2005 to 2006 indicated that cytokeratin 18 M30, a promising serologic NAFLD biomarker, positively associated with PFHxS, PFOA and PFNA in serum [11]. All these indirectly indicate a potential association between PFASs exposure and NAFLD. The trials in vivo and in vitro also provide biological evidence that PFAS exposure led to NAFLD. Dietary treatment of mice with 0.002% (w/w) PFOA or 0.005% (w/w) PFOS for 10 days resulted in significant reductions in serum levels of cholesterol and triglyceride (TG), a moderate increase in the serum activity of alkaline phosphatase and hepatomegaly [12]. PFOA and PFOS have been shown showed to induce hepatocyte peroxisome proliferation, liver hypertrophy, vacuolization and hyperplasia in vitro [13]. In addition, long-term feeding of low concentrations of PFASs to mice, such as 0.0001% (w/w) PFOA for 8 weeks and 0.0001% (w/w) PFOS for 6 weeks might have a protective effect against fatty diet-induced hepatic steatosis [14,15]. But the mechanisms behind this paradoxical phenomenon are not yet clear.

Therefore, whether different exposure levels of PFASs in the US population have a direct impact on the occurrence of NAFLD needs further study. In this cross-sectional study, the generalized linear model (GLM), nonlinear restricted cubic spline (RCS) model, and Bayesian kernel machine regression (BKMR) model were used to assess the single effect, nonlinear dose-response relationship, and overall effects of PFASs on NAFLD-related indicators.

2 Methods

2.1. Data source

NHANES is a large cross-sectional survey, that combines interviews and physical examinations to assess the health and nutritional status of adults and children in the United States (https://www.cdc.gov/nchs/nhanes/index.htm). The NHANES 2017–2018 survey is the only one that includes the liver ultrasound transient elastography examination which use to identify participants with fatty liver.

2.2. Study participants

There were 9254 participants in NHANES 2017–2018, of whom 5492 (59.3%) had a completed liver ultrasound transient elastography examination. The 983 participants were excluded for age under 20, and 181 were excluded for viral hepatitis (43 for HCV RNA+; 27 for HBsAg+) and 73 for excess alcohol (more than two drinks per day in women and more than three drinks per day in men) [16]. At last, 99 participants were excluded for body mass index (BMI) < 18.5, as the link between undernutrition and fatty liver was more complex and less understood [17]. Figure 1 showed the flowchart of participant selection. The tests of HBV surface antigen and HCV RNA were detailed in the supplementary document.

Figure 1. Flowchart of participants selected from the NHANES 2017–2018.

2.3. Vibration-Controlled transient elastography (VCTE)

VCTE involves an ultrasound machine (FibroScan model 502 V2 Touch) with a handheld transducer that taps, in rapid succession, on the participant’s abdomen. The controlled attenuation parameter (CAP) could be quantified liver fat [18]. CAP values ranged from 100 to 400 dB/m, with higher values indicating higher amounts of liver fat [19]. Participants with NAFLD were diagnosed as having a median CAP score of ≥285 dB/m–optimized sensitivity of 80% and specificity of 77% for detecting hepatic steatosis [20]. Although not as reliable as liver biopsy, VCTE has been used as an identification method for NAFLD in several studies because it is noninvasive and has higher accuracy than common liver color ultrasounds [17,21–23].

2.4. PFAS detection

The blood samples were collected, and serum specimens were stored at −30 °C until being detected. Online Solid Phase Extraction-High Performance Liquid Chromatography- Turbo Ion Spray-Tandem Mass Spectrometry (online SPE-HPLC-TIS- MS/MS) was used for the quantitative detection of PFAS. A total of 5 PFAS with the highest rate of detection (>80%) in the US population were considered, including perfluorodecanoic acid (PFDeA), PFHxS, PFNA, PFOA and PFOS. The concentration of PFOA was the sum of the concentrations of linear PFOA and the branched PFOA isomers. The concentration of PFOS was the sum of the linear PFOS and monomethyl-branched PFOS isomers. The lower limit of detection (LLOD in) for each PFAS was 0.10 ng/mL and replaced with LLOD divided by the square root of 2, which was 0.07 ng/mL [24]. Total relative standard deviations were ≤ 15% (CV ≤ 15%).

2.5. Hepatic steatosis-related biomarkers detection

2.5.1. Fatty liver index (FLI)

FLI was to predict hepatic steatosis, performed by BMI, waist circumference, TG and GGT. The area under the ROC curve (AUROC curve) was 0.84 with ultrasonography as the standard in diagnosing fatty liver [25]. The equation for FLI is: $$FLI=\frac{e^{0.953\times \ln \left[TG\left(mg/dL\right)\right]+0.139\times BMI\left(kg/m^2\right)+0.718\times \ln [\text{GGT}\left(U/L\right)]+0.053\times waistcircumference\left(cm\right)-15.745}}{1+e^{0.953\times \ln \left[TG\left(mg/dL\right)\right]+0.139\times BMI\left(kg/m^2\right)+0.718\times \ln [\text{GGT}\left(U/L\right)+0.053\times waistcircumference\left(cm\right)-15.745}}\times 100$$

2.5.2. NAFLD liver fat score (LFS)

NAFLD LFS was used to estimate liver fat based on metabolic syndrome and type 2 diabetes, serum insulin, aspartate aminotransferase (AST) and the AST/ALT ratio. The AUROC curve was 0.86 with proton magnetic resonance spectroscopy as the standard in diagnosing fatty liver [26]. The equation for LFS is: $$\mathit{NAFLD}\mathrm{ }\mathit{LFS}=-2.89+1.18\times \mathit{metablolic}\mathrm{ }\mathit{syndrome}\mathrm{ }\left(\mathit{yes}=1,no=0\right)+0.45\times \mathit{type}\mathrm{ }2\mathrm{ }\mathit{diabetes}\mathrm{ }\left(\mathit{yes}=2,no=0\right)+0.15\times \mathit{insulin}\mathrm{ }\left(mU/L\right)+0.04\times \mathit{AST}\left(U/L\right)-0.94\times \mathit{AST}/\mathit{ALT}$$

The detailed definition of metabolic syndrome is in the supplementary materials.

2.5.3. Framingham steatosis index (FSI)

FSI was used to identify hepatic steatosis. The AUROC curve of FSI was 0.845 with abdominal Computed Tomography scans as the standard in diagnosing fatty livers [27]. Through constructing a stepwise regression model, FSI was derived from the indexes which were strongly associated with steatosis including age, sex, BMI, TG, hypertension, diabetes and the ALT/AST ratio. The equation for FSI is: $$\mathit{FSI}=\mathrm{ }-7.981+0.011\times \mathit{age}\mathrm{ }\left(\mathit{years}\right)-0.146\times \mathit{sex}\mathrm{ }\left(\mathit{female}=1,\mathrm{ }\mathit{male}=0\right)+0.173\times \mathit{BMI}\mathrm{ }\left(kg/m^2\right)+0.007\times TG\mathrm{ }\left(mg/dL\right)+0.593\times \mathit{hypertension}\mathrm{ }\left(\mathit{yes}=1,no=0\right)+0.789\times \mathit{diabetes}\left(\mathit{yes}=1,\mathrm{ }no=\mathrm{ }0\right)+1.1\times [\mathit{ALT}/\mathit{AST}\ge 1.33\mathrm{ }(\mathit{yes}=1,\mathrm{ }no=0)]$$

2.6. Hepatic fibrosis-related biomarkers

2.6.1. Fibrosis-4 index (FIB-4)

The FIB-4 was originally developed for patients co-infected with HCV/HIV in 2006 [28]. As a widely used non-invasive biomarker, FIB-4 used four indicators associated with liver fibrosis including platelets, ALT, AST and age to predict the level of liver fibrosis in many liver diseases including NAFLD [17,29,30]. The overall statistical consistency with liver biopsy in the classification of hepatic fibrosis was 0.87 [28]. The equation for FIB-4 is [28]: $$\mathit{FIB}-4=\frac{\mathit{Age}\left(\mathit{years}\right)\times \mathit{AST}\left(U/L\right)}{\mathit{Platelet}\mathrm{ }\mathit{count}\mathrm{ }\left(1000\mathit{cells}/\mu L\right)\times \sqrt{ALT}}.$$

2.6.2. NAFLD fibrosis score (NFS)

NFS was constructed based on age, hyperglycemia, BMI, platelet count, albumin, and AST/ALT ratio to separate NAFLD patients with and without advanced fibrosis in 2007. The statistical consistency with liver biopsy in diagnosing advanced fibrosis was 0.82 [31]. The equation for NFS is: $$\begin{array}{c}NFS=-1.675+[0.037\times age\left(years\right)]\\ +[0.094\times BMI\left(kg/m^2\right)]\\ +[1.13\times impairedfastingglucose\left(ordiabetes\right)\\ \left(yes=1,no=2\right)]+\left(0.99\times AST/ALT\right)\\ -[0.013\times \text{platelet}\text{count}\left(\times {10}^9/L\right)]\\ -0.66\times \text{Albumin}\left(\text{g}/dL\right)\end{array}$$

2.6.3. Hepamet fibrosis score (HFS)

HFS is a new fibrosis scoring system, developed in 2020. Variables used to determine HFS were sex, age, homeostatic model assessment score, presence of diabetes, levels of aspartate aminotransferase, and albumin, and platelet counts. HFS discriminated between patients with and without advanced fibrosis with an AUROC curve value of 0.85 [32]. The equation for HFS is: $$\mathit{HFS}=1/\{1+ê[5.390-0.986\times \mathit{Age}(45-64\mathrm{ }\mathit{years}\mathrm{ }of\mathrm{ }\mathit{age})-1.719\times \mathit{Age}(\ge 65\mathrm{ }\mathit{years}\mathrm{ }of\mathrm{ }\mathit{age})+0.875\times \mathit{Male}\mathrm{ }\mathit{sex}\mathrm{ }-0.896\times \mathit{AST}(35-69IU/L)-2.126\times \mathit{AST}(\ge 70IU/L)-0.027\times \mathit{Albumin}\mathrm{ }(4-4.49\mathrm{ }g/dL)-0.897\times \mathit{Albumin}(<4g/dL)-0.899\times \mathit{HOMA}(2-3.99\mathrm{ }\mathit{with}\mathrm{ }no\mathrm{ }\mathit{Diabetes}\mathrm{ }\mathit{Mellitus})-1.497\times \mathit{HOMA}(\ge 4\mathrm{ }\mathit{with}\mathrm{ }no\mathrm{ }\mathit{Diabetes}\mathrm{ }\mathit{Mel}\mathrm{l}\mathit{itus})-2.184\times \mathit{Diabetes}\mathrm{ }\mathit{Mellitus}-0.882\times \mathit{plateles}\times 1000/\mu L\left(155-219\right)-2.233\times \mathit{platelets}\times 1000/\mu L(<155)\left]\right\}$$

Homeostatic model assessment of insulin resistance (HOMA-IR) was calculated based on insulin and glucose [fasting insulin(µU/mL) x fasting glucose(mg/dL)/405].

The detailed methods for the detection of biomarkers related to NAFLD, including BMI, blood pressure (BP), ALT, AST, GGT, albumin, TG, high-density lipoprotein (HDL), fasting plasma glucose (FPG), insulin and platelet count are in the supplementary materials.

2.7. Covariates selection

A directed acyclic graph (DAG) was made based on the results of baseline characteristics comparison and prior knowledge. And the DAGitty v3.0 was used to visualize DAG and calculated the minimal sufficient adjustment sets to recognize the covariates.

The baseline characteristics included socio-economic factors (gender, age, race, education, marital status, poverty income rate (PIR)), pathophysiological factors (cancer, hypertensive, diabetic, liver diseases, creatinine level which might be harmful to the livers and affect chemical concentrations among participants), and behavior factors (drug use, smoking status and physical activity). The methods for information collection of baseline characteristics are shown in the supplementary materials.

2.8. Statistical analysis

2.8.1. Description of distribution

The total population was divided into NAFLD and non-NAFLD. The PFASs value was divided into quartiles. To compare the distribution differences of baseline characters, the Chi-square tests were used for categorical data and the Kruskal-Wallis tests were used for continuous data. Kruskal-Wallis tests were used to compare the distribution of NAFLD-related biochemical values in quartiles of PFASs.

2.8.2. The association between PFASs exposure and NAFLD

Two logistic models were used to study the association between PFASs exposure and NAFLD. Model 1 was a univariate model, and model 2 was a multivariate model adjusted for covariates. The PFASs levels were divided into quartiles and compared three other quartiles with the first to calculate the odds ratio (OR) and 95% confidence interval (95% CI) to estimate the risk of NAFLD in the population exposed to PFASs at each quantile level. Besides, the median exposure level of each group of elements was taken as a continuous variable to conduct a trend test. For models with significant trends, RCS was employed to reveal the dose-response relationship between PFASs and NAFLD. RCS was suitable for nonlinear models and could supplement the shortcomings of GLM [33]. The generation of RCS is detailed in the supplementary material.

2.8.3. Correlations among PFASs and between PFASs and NAFLD-related biomarkers

Logarithmic transformation of PFASs values was performed to approximate the normal distribution, as shown in the quantile-quantile plot (Figure S1). The Spearman’s test was used to analyze the correlation between PFASs and NAFLD-related biomarkers values. And Pearson’s test was used to analyze the correlation among PFASs.

2.8.4. Association of FIB-4 with PFAS

Two GLMs included a univariate model and a multivariate model adjusted for covariates to estimate the increased risk of FIB-4 in the population exposed to single PFASs.

And BKMR was used to explore the mixed exposure effect [34]. The BMKR model is Y_i = ℎ(PFDeA, PFHxS, PFNA, PFOA, PFOS) + βx_i + ε_i, where $Y_i$ denotes the FIB-4 for individual $i\left(i=1,\dots ,n\right)$, $h$ denotes the unknown exposure-response function to be estimated, $\beta$ represents the effect of the covariates ($x_i$ is a vector), and ${\epsilon }_i$ is the residual. The covariates were consistent with those mentioned earlier. The number of iterations of the Markov chain Monte Carlo sampler was 10,000. A BMKR model without variable selection and a BMKR model with gamma prior distribution was fitted to demonstrate the robustness of the results [35].

To account for the complex survey sampling, the strata (sdmvstra), primary sampling units (sdmvpsu) information, and the weight was used to adjust the results of logistics regression. The R ‘Survey’ package was used to adjust the results according to the NHANES’s analytic guideline (https://wwwn.cdc.gov/nchs/nhanes/ analyticguidelines.aspx). The statistically significant difference was defined as a two-sided p < 0.05.

3 Results

3.1. Effect of PFASs exposure on the risk of developing NAFLD

Among those with completed ultrasound transient elastography of the liver, those who were younger than 20 years of age, with excessive alcohol consumption, hepatitis B/C, and a BMI of less than 18.5 were excluded. A total of 4267 people were included and divided into non-alcoholic liver patients and healthy controls. According to NHANES analytic guidelines, the samples were weighted, and then the weighted sample represented over 197.6 million Americans with 36.8% NAFLD, as shown in Table 1.

Table 1. Baseline characteristics among the healthy control and NAFLD.

Variables	Unweighted			Weighted
	Healthy control	NAFLD	P	Healthy control	NAFLD	P
N	2627	1640		125005099	72633232.3
Males (%)	1166 (44.4)	916 (55.9)	<0.001	55436488.3 (44.3)	41007286.0 (56.5)	<0.001
Age (%)			<0.001			<0.001
20-39	952 (36.2)	361 (22.0)		53074047.0 (42.5)	20119812.6 (27.7)
40-59	759 (28.9)	593 (36.2)		39647923.4 (31.7)	28426682.7 (39.1)
≥60	916 (34.9)	686 (41.8)		32283128.5 (25.8)	24086737.0 (33.2)
Race (%)			<0.001			<0.001
Mexican American	278 (10.6)	311 (19.0)		8509323.7 (6.8)	9270185.3 (12.8)
Other Hispanic	256 (9.7)	154 (9.4)		9343918.4 (7.5)	4546665.1 (6.3)
Non-Hispanic White	857 (32.6)	589 (35.9)		77454152.1 (62.0)	45819021.1 (63.1)
Non-Hispanic Black	681 (25.9)	293 (17.9)		15852537.2 (12.7)	5963463.6 (8.2)
Other Race	555 (21.1)	293 (17.9)		13845167.4 (11.1)	7033897.2 (9.7)
College (%)	1555 (59.3)	899 (55.0)	0.007	81591653.2 (65.3)	42616961.3 (58.7)	0.035
Marital status (%)			<0.001			<0.001
Married/Living with partner	1500 (57.2)	1063 (64.9)		73756033.9 (59.0)	50690283.7 (69.8)
Never married	553 (21.1)	220 (13.4)		27969902.5 (22.4)	9639407.8 (13.3)
Widowed/Divorced/Separated	571 (21.8)	356 (21.7)		23213046.3 (18.6)	12293040.3 (16.9)
PIR (%)			0.529			0.201
<2	1020 (44.4)	653 (45.5)		35851661.4 (32.0)	20271865.6 (31.0)
2-4	641 (27.9)	409 (28.5)		30511085.9 (27.3)	20709853.8 (31.7)
>4	634 (27.6)	372 (25.9)		45571881.6 (40.7)	24354910.5 (37.3)
History of cancer (%)	250 (9.5)	183 (11.2)	0.094	12056300.5 (9.6)	8743020.9 (12.0)	0.015
History of hypertension (%)	808 (30.8)	800 (48.8)	<0.001	29570442.6 (23.7)	33164449.0 (45.7)	<0.001
History of diabetes (%)	257 (9.8)	395 (24.1)	<0.001	6468281.9 (5.2)	14718167.8 (20.3)	<0.001
History of liver diseases (%)	93 (3.5)	117 (7.1)	<0.001	3886315.5 (3.1)	4676040.1 (6.4)	0.002
Ever take drug (%)	881 (55.5)	462 (51.7)	0.075	54672680.0 (62.6)	26876405.4 (58.4)	0.219
Smoking Status (%)			<0.001			0.013
None	1615 (61.5)	917 (55.9)		76069809.2 (60.9)	40579922.1 (55.9)
Former	564 (21.5)	462 (28.2)		28878849.4 (23.1)	21042004.4 (29.0)
Current	448 (17.1)	261 (15.9)		20056440.3 (16.0)	11011305.8 (15.2)
Physical Activity (%)			<0.001			<0.001
Inactive	1254 (47.7)	950 (57.9)		49117062.0 (39.3)	38273735.0 (52.7)
Moderate	616 (23.4)	395 (24.1)		31854932.6 (25.5)	19775842.3 (27.2)
Active	757 (28.8)	295 (18.0)		44033104.3 (35.2)	14583655.0 (20.1)
Urine creatinine (mg/dL) (median [IQR])	112 [62, 172]	118 [70, 175]	0.038	112 [59, 168]	115 [72, 173]	0.130

The baseline features were compared between normal and NAFLD group, including socio-economic factors (gender, age, race, education, marital status, PIR), pathophysiological factors (cancer, hypertensive, diabetic, liver diseases), behavior factors (drug use, smoking status and physical activity) and urine creatine level. There were statistical differences in socio-economic characteristics including gender (p < 0.001), age (p < 0.001), race (p < 0.001), education levels (p = 0.035) and marital status (p < 0.001) between normal and NAFLD group while there was no significant difference in PIR. In pathophysiological characteristics, there were significant differences in the history of cancer (p = 0.015), hypertension (p < 0.001), diabetes (p < 0.001), and liver disease (p = 0.002), but not urinary creatinine levels. And in behavior characteristics, there were significant differences in smoking (p = 0.013), and physical activity status (p < 0.001) between NAFLD and normal while there was no significant difference in drug use. Males, aging, Mexicans, married or living with a partner, lower education level, a history of cancer, hypertension, diabetes, liver diseases, smoking and lower physical activity were more common among NAFLD patients.

A total of 1318 samples were successfully merged with PFASs laboratory data. The populations were grouped according to quartile of PFAS exposure concentration and compared baseline characteristics. As shown in Table 2, there were significant differences in the distribution of gender, age, race and PIR in people exposed to different concentrations of the same PFASs. The concentrations of PFNA and PFOA exposure differed between those with and without a college education. And the concentrations of PFOS exposure differed between those with and without a history of hypertension. Marital status, history of cancer, diabetes, liver disease, urinary creatinine levels, drug use, smoking, and physical activity were not associated with PFASs exposure. The details were shown in Table S2.

Table 2. The P value of Chi-square tests conducted by PFASs quartile grouping.

Variables	PFDeA	PFHxS	PFNA	PFOA	PFOS
Gender	0.063	<0.001***	0.062	<0.001***	<0.001***
Age cohort	0.017*	<0.001***	<0.001***	<0.001***	<0.001***
Race	0.035*	0.006***	0.015*	0.005**	0.003**
College	0.094	0.067	0.006**	<0.001***	0.130
Marital status	0.239	0.280	0.456	0.362	0.127
PIR	0.001**	0.005**	0.007**	<0.001***	0.023*
History of cancer	0.692	0.780	0.285	0.433	0.267
History of hypertension	0.659	0.339	0.150	0.737	0.002
History of diabetes	0.487	0.898	0.245	0.976	0.180
History of liver diseases	0.890	0.729	0.977	0.897	0.483
Ever take drug	0.280	0.255	0.367	0.367	0.454
Smoking status	0.182	0.062	0.249	0.227	0.210
Physical activity	0.241	0.223	0.755	0.129	0.149
Urine creatinine	0.150	0.376	0.334	0.286	0.686

***p < 0.001; **p < 0.010; *p < 0.050.

The DAG graphs were formed through baseline characteristic comparison results and prior knowledge, as shown in Figure 2. The covariates were recognized as Gender, Age, Race, Education and PIR.

Figure 2. DAG was made based on the results of baseline characteristics comparison and prior knowledge.

After removing 3 participants with education information missing and 165 participants with PIR information missing, a total of 1150 participants were used to construct two logistic models to study the association between PFASs exposure and NAFLD. Model 1 was a univariate model, and model 2 was a multivariate model adjusted for covariates of gender, age, race, education and PIR. The results were shown in Figure 3(A). Taking the first interval as the reference, all intervals of PFHxS, PFNA, PFOA, and PFOS had no significant association with NAFLD in models 1 and 2. And in PFDeA, the OR of quartile 4 was 0.47 (0.26-0.85) in model 1, but without a significant association in model 2. Besides, the trend tests of all the models were also performed. The trend tests for PFDeA were significant in both model 1 (p = 0.013) and model 2 (p = 0.038). The nonlinear relationship was explored between PFDeA and OR by RCS, as shown in Figure 3(B). No statistically significant non-linear association was found in PFDeA with NAFLD in model 1 (p = 0.981 for the non-linear test) and model 2 (p = 0.597 for the non-linear test).

Figure 3. Effect of PFASs exposure on the risk of developing NAFLD. (A) Logistic regression result of NAFLD and single PFASs. Odds ratios (OR)1 and P1 for model 1 without adjusted; OR2 and P2 for model 2 with adjusted by gender, age, race and education and PIR. ‘#’ was the result of P for trend. (B) RCS results of NAFLD and ln (PFDeA). The median value of ln (PFDeA) was used as a reference point. Blue was the density of the ln (PFDeA) distribution.

3.2. Effects of PFAS exposure on NAFLD-related biomarkers

Liver ultrasound transient elastography is not the gold standard for NAFLD, so the association between PFASs with other NAFLD-related biomarkers such as body measure, biochemical indicators, and FLI, LFS, FSI, FIB-4, NFS, HFS were further analyzed. Table 3 showed the biomarkers have significant differences among different PFASs exposure. The details were shown in Table S3. The concentration (ng/mL) for PFDeA was 0.20 [0.10, 0.30] (median [IQR]), for PFHxS was 1.20 [0.70, 1.90], for PFNA was 0.50 [0.30,0.70], PFOA was 1.47 [0.97, 2.17], for PFOS was 4.60 [2.70, 8.06]. The result suggests BMI changed significantly in all PFASs exposures while waist circumference significantly changed among PFDeA, PFHxS and PFNA exposure. AST was significantly related to PFHxS, PFOA and PFOS and seemed to be most closely related to PFASs exposure among the transaminases. ALT was related to the PFHxS and PFOS while the GGT was related to PFOS. Serum albumin was a critical plasma protein produced by the liver. Advanced cirrhosis was characterized by reduced albumin concentration as well as impaired albumin function [36]. The result showed that differences in albumin might be related to PFHxS, PFOA and PFOS exposure. TG and HDL were related to PFOS and PFDeA respectively. Insulin resistance was the major hallmark of NAFLD, indirectly influencing hepatic glucose and lipid metabolism by increasing the flux of substrates that promoted lipogenesis and gluconeogenesis [37]. PFOS was closely related to fasting plasma glucose (FPG), insulin and HOMA-IR. Liver damage was associated with widespread coagulation dysfunction. The latest research also showed that platelets accelerated liver regeneration and promoted liver fibrosis [38]. In this result, platelet count was significantly related to PFHxS, PFNA, PFOA, and PFOS respectively.

Table 3. The P value of the Kruskal-Wallis tests which were conducted by PFASs quartile grouping.

Biomarkers	PFDeA	PFHxS	PFNA	PFOA	PFOS
BMI (kg/m^2)	0.037*	0.003**	0.024*	0.026*	0.038*
Waist circumference (cm)	0.031*	0.042*	0.028*	0.525	0.313
ALT (U/L)	0.642	0.030*	0.973	0.110	0.014*
AST (U/L)	0.611	0.005**	0.854	0.008**	<0.001***
GGT (U/L)	0.918	0.097	0.158	0.183	0.034*
Albumin (g/dL)	0.224	0.011*	0.409	0.012*	0.028*
TG (mg/dL)	0.944	0.198	0.119	0.116	0.025*
HDL (mg/dL)	0.008**	0.105	0.155	0.311	0.121
FPG (mg/dL)	0.616	0.644	0.350	0.688	0.009**
Insulin (μU/mL)	0.231	0.115	0.060	0.066	0.004**
HOMA-IR	0.310	0.198	0.062	0.159	0.004**
Platelet (x10^3 cells/μL)	0.539	0.010*	0.005**	0.018*	<0.001**
FLI	0.321	0.138	0.031*	0.816	0.127
LFS	0.494	0.724	0.192	0.757	0.006**
FSI	0.683	0.387	0.074	0.866	0.012*
FIB4	0.003**	0.001**	0.001**	<0.001***	<0.001***
NFS	0.426	0.014*	0.118	0.010*	0.014*
HFS	0.310	0.176	0.013*	0.008**	0.005**

***p < 0.001; **p < 0.010; *p < 0.050.

FLI, LFS and FSI were to assess the relationship between PFASs and hepatic steatosis. FLI was significantly related to PFNA, while LFS and FSI were related to PFOS. PFASs exposure appeared to be more closely associated with liver fibrosis than with steatosis. NFS was significantly related to PFHxS, PFOA and PFOS while HFS was significantly related to PFNA, PFOA and PFOS separately. More importantly, FIB4 changes occurred in all the types of PFAS exposed.

As shown in Figure 4(A), through Spearman’s correlation analysis, albumin, FPG, platelet count, FIB4 and HFS had significant correlations with all of PFASs. Platelet count had negative relations with PFASs, while albumin, FPG, FIB4 and HFS had positive relations. Interestingly, three hepatic steatosis indicators FLI, LFS and FSI were almost not significantly correlated with PFASs exposure. But fibrosis indicators FIB-4, NFS and HFS were positively correlated with almost all PFASs exposure. Therefore, PFASs exposure is not closely associated with hepatic steatosis but may aggravate liver fibrosis.

Figure 4. Correlation analysis (A) Spearman correlation coefficients between PFASs and NAFLD biomarkers. (B) Pearman correlation coefficients among PFASs logarithm. The correlation coefficient was shown in color. Blue was a positive correlation and red was a negative correlation. Color darker, circle bigger implied a stronger correlation. Demerit marks showed p > 0.05.

In addition, we also constructed a Pearson’s correlation analysis to explore the correlation between PFASs after logarithmic transformation, as shown in Figure 4(B). All PFASs exposures were positively correlated. The correlation order was PFDeA vs PFNA (r = 0.74), PFNA vs PFOS (r = 0.72), PFHxS vs PFOS (r = 0.71), PFNA vs PFOA (r = 0.67), PFOA vs PFHxS (r = 0.63), PFDeA vs PFOS (r = 0.61), PFOA vs PFOS (r = 0.59), PFOA vs PFDeA (r = 0.51), PFHxS vs PFNA (r = 0.49), PFHxS vs PFDeA (r = 0.28). The P for all correlation tests of significance was less than 0.001.

3.3. Associations of hepatic fibrosis-related biomarkers with PFASs by GLM models

The GLM was used to further verify whether there were significant correlations between PFAS and liver fibrosis indicators. The result was shown in Figure 5. In the univariate model, PFDeA (p = 0.005), PFHxS (p < 0.001), PFNA (p < 0.001), PFOA (p < 0.001) and PFOS (p < 0.001) were all positively related with FIB-4; PFHxS (p = 0.026), PFNA (p = 0.047), PFOA (p = 0.026) and PFOS (p = 0.002) were positively related with NFS; and PFNA (p = 0.009), PFOA (p = 0.007), PFOS (p = 0.010) were positively related to HFS. However, after the adjusted gender, age, race, education, and PIR, only the correlation between PFOS and FIB-4 was significant (p = 0.025).

Figure 5. The association between PFASs and hepatic fibrosis-related biomarkers was analyzed by GLM. Model 1 was unadjusted; model 2 was adjusted by gender, age, race and education and PIR.

The findings of FIB-4 remained robust in our sensitivity analysis using the study cohort without participant selection (N = 1929, PFASs test completed), as shown in Figure S2.

3.4. Multiple exposures of PFAS affected on FIB-4 by BKMR model

The BKMR model was used to explore the overall and individual of the PFASs by gender, age, race, education and PIR adjustment. The posterior inclusion probability (PIP) of each PFASs was calculated. The PIP values for PFASs from high to low were 1.000 for PFOS, 0.441 for PFHxS, 0.025 for PFOA, 0.005 for PFDeA, and 0.004 for PFNA. It suggested that PFOS was significantly correlated to FIB-4 in the BKMR model. Besides, the global effect of the mixed PFASs on FIB-4 was significant. As compared to all of the PFASs were fixed at their median values, when all of the PFASs were at their 75th percentile the FIB-4 increased by 0.143 units, while at their 25th percentile, the FIB-4 decreased by 0.102 units, as shown in Figure 6(A). The higher exposure to the PFASs was accompanied by a greater risk of FIB-4. The interaction of the mixed PFASs was also assessed. For PFOS, it all showed significant and positive effects on the risk of FIB-4 increase whatever all of the other PFASs were fixed at their 25th, 50th, or 75th percentiles, as shown in Figure 6(B). But other PFASs had no significant effect on FIB-4. Therefore, PFOS was the main reason that PFASs increased the risk of FIB-4. Models from sensitivity analyses were very similar to our preliminary results, which proved the robustness of our results, as shown in Figure S3.

Figure 6. The overall and individual exposure effect of PFASs on FIB-4. (A) The overall effect of the PFASs mixture being fixed to different percentiles as compared to their 50 percentiles. (B) The estimated values of the individual effect were calculated by comparing the FIB-4 increase risk when a single PFAS was at its 75th percentile as compared to when that was at its 25th percentile, where all of the remaining PFAS were fixed at their 25th (red line), 50th (green line), or 75th (blue line) percentile.

Discussion

In this study, the associations were investigated between several serum PFASs and NAFLD based on liver ultrasound transient elastography and other non-invasive biomarkers. The GLM and RCS were constructed to reveal there was no significant association between PFASs and CAP. The results of hepatic steatosis-related indicators FLI, LFS and FSI also demonstrated that PFASs had no significant correlation with hepatic steatosis based on ANOVA and Spearman’s tests. But PFASs exposure was indeed closely related to glucose and lipid metabolism. Serum TG has positive correlations with PFHxS, PFNA, PFOA and PFOS while FPG had positive correlations with all of the PFASs in the study. PFASs structurally resembled fatty acids and activated the peroxisome proliferators-activated receptors (PPARs) signaling pathway, responsible for perturbations of glucose and lipid homeostasis partially [39]. A new study showed serum albumin mediated the effects of PFASs on serum lipid levels [40]. PFASs had shown a strong binding affinity to serum proteins and were positively associated with serum albumin [41], which was bound to PFASs and transported to target organs that regulate serum lipid levels [40]. In addition, we also found that PFASs were negatively correlated with BMI and positively correlated with HDL which was also reported in a previous study among European teenagers [42]. But a high level of PFOA (median level = 1635.96 ng/mL) was reported to negatively associate with HDL in workers of a fluorochemical plant [43]. These evidences suggested that the effects of PFASs on hepatic steatosis appeared to be a double-edged sword and associated with the exposure concentration. The activation of PPARs induced genes in both lipogenesis/uptake and fatty acid oxidation pathways. But targeting PPARs was also an emerging strategy for the treatment of steatosis and inflammation. The agonists of PPARα (fenofibrate) and PPARγ (pioglitazone, lobeglitazone) had been shown to reduce hepatic steatosis in both animal studies and clinical trials. Additionally, the lipid metabolism pathways might be influenced by the feedback loop [44]. As previously mentioned, PFOS and PFOA were found to relieve steatosis in mice with diet induce-fatty liver [14, 15]. The molecular mechanisms of the switch between lipogenic and lipid oxidation effects induced by PFASs are not fully understood.

FIB-4, NFS and HFS were used to assess the association between PFASs exposure and hepatic fibrosis. The results of Spearman’s test showed FIB-4, NFS and HFS were positively correlated with PFASs. Elevated transaminases have been found under multiple PFASs exposures. It is consistent with previous studies [45–48]. More importantly, platelet accounts significantly decreased under PFASs exposure. It was also found in the previous study among young and middle-aged Taiwanese populations that platelet count decreased significantly with increasing quartiles of PFOA and PFOS [49]. Platelet was proven to decrease hepatic fibrosis. Platelets contacted hepatic stellate cells (HSCs) and released adenine nucleotides, which subsequently led to the production of adenosine through the degradation by HSCs. Activated HSCs were inactivated by adenosine and had a decreased ability to produce pro-fibrotic growth factor TGF-β and secrete extracellular matrix [38]. The decrease in platelet count associated with PFASs exposure might be the main reason why PFASs showed an association with hepatic fibrosis but not steatosis. So, PFASs exposure appeared to be more closely associated with hepatic fibrosis than steatosis. The reduction of platelets caused by PFASs might be the key to the changes in liver fibrosis. Regular examination for platelet in people with high PFASs exposure might be effective in preventing the development of liver fibrosis.

The GLM was used to verify the association between PFASs and hepatic fibrosis-related indicators. PFOS exposure was linearly significant with FIB-4, both before and after adjustment. Due to the interaction between PFASs, we also adopted the BKMR model to verify the mixed effects and single effects of PFAS. The results showed that the higher the mixed PFASs concentration, the higher the FIB-4 of the participants. PFOS was the main cause of FIB-4 elevation. This result was also verified by the BKMR model under different priors.

Exposure to PFAS has generated publicity in some specific occupational populations. The median concentration of PFOA in the serum of professional ski waxers was 50 ng/ml [50]. In the serum of fluorochemical plant workers, the median concentration of PFOA was 1635.96 ng/mL and PFOS was 33.46 ng/mL [43]. And in the serum of firefighters, the median concentration of PFOA was 31.50 ng/mL and PFOS was 27.85 ng/mL [51]. These were well above the normal levels in the general population in our study. These specific occupational populations are likely to be at high risk for liver fibrosis. Although US environmental protection agency has proposed maximum contaminant level goals in drinking water for PFASs [52] and the US Conference for governmental industrial hygienists has established threshold limit values for three PFAS in the air: perfluoroisobutylene, perfluorobutyl ethylene, and ammonium perfluorooctanoate [53]. The occupational protection administration for PFASs is far from adequate. Further studies of quantitative individual exposure and biomarker data in populations with occupational exposure to PFAS are urgently needed.

This study has several strengths. It is the first study using multiple NAFLD-related biomarkers to analyze the association between PFASs and NAFLD. Then, a variety of models including GLM, RCS, and BKMR were used to assess the single effect, nonlinear dose-response relationship, and overall effects of PFASs on NAFLD. Additionally, the present study used a nationally representative sample of the US population, with a relatively large population, and adjusted for potential confounders, which reduced the potential bias associated with the selection of study subjects and rendered the findings generalizable. Inevitably, the present study also has several deficiencies. The results were not based on the exact liver biopsy. Because of the limitations of a cross-sectional study itself, laboratory research and prospective cohort studies are needed to verify in the future.

Conclusion

Our study indicated that current exposure doses of PFAS did not significantly change the risk of developing NAFLD based on the liver ultrasound transient elastography. From the changes of NAFLD-related biomarkers, including FLI, LFS, FSI, FIB-4, NFS, HFS, PFASs exposure appeared to be more closely associated with hepatic fibrosis than steatosis. PFOS might be the main cause of PFASs associated with hepatic fibrosis.

Ethical approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Author contributions

Conception and writing original draft preparation: Wenli Cheng; Methodology and formal analysis: Min Li, Luyun Zhang and Cheng Zhou; Data acquisition: Xinyu Zhang, Chenyu Zhu and Luyi Tan; Interpretation of the data: Hui Lin; Design of the work and funding acquisition: Wenjuan Zhang and Wenji Zhang. All authors provided critical revisions and approved the final manuscript.

Acknowledgements

Not applicable.

Disclosure statement

The authors declare that they do not have any conflicts of interest related to this study. This manuscript has been read and approved by all the authors and has not been submitted to or is not under consideration for publication elsewhere.

Data availability statement

Data will be made available on request.

Additional information

Funding

This work was funded by the Guangdong Basic and Applied Basic Research Foundation (2021A1515011220, 2023A1515010353); the Foundation of Director of Crops Research Institute, Guangdong Academy of Agricultural Sciences (202205); Special fund for scientific innovation strategy-construction of high-level Academy of Agriculture Science (R2018YJ-YB3002, R2022PY-QY004) and the National Natural Science Foundation of China (81903319).