Genetic variations interact with polybrominated diphenyl ether exposure to alter lipid homeostasis

Domain: chemistry
Published: 2026-03-05
Editorial outcome: minor_revision

Abstract
Polybrominated diphenyl ethers (PBDEs) are implicated in dyslipidemia, but the molecular basis of individual susceptibility remains elusive. Here we report an analysis based on the China National Human Biomonitoring cohort, where we integrate exposome, genomic, and metabolomic data to identify 3,571 genetic variants that interact with PBDE exposure to influence dyslipidemia risk. Metabolomic analysis highlights glycine and glycerophosphate as key mediators. A polygenic risk score derived from these PBDE-interactive variants significantly enhances dyslipidemia prediction in highly exposed individuals. Among these, rs9869609 emerges as a candidate causal variant, showing the strongest association with hypercholesterolemia risk (β=1.18, FDR=0.0078). Further functional validation using single-base CRISPR/Cas9 editing reveals that the rs9869609-G allele downregulates SLC6A20 expression by strengthening BHLHE40 binding, which further impairs glycine transport and promotes cholesterol accumulation, particularly under 2,2′,4,4′-Tetrabromodiphenyl ether exposure. Collectively, our study elucidates a gene-environment interaction mechanism through which genetic variants modulate lipid metabolism in response to PBDE exposure.

Introduction
Metabolic disorders such as obesity, type 2 diabetes (T2D), and dyslipidemia, account for approximately 60% of global mortality, posing a significant public health challenge1. Although genetic susceptibility plays a fundamental role, environmental exposures, particularly widespread pollutants like polybrominated diphenyl ethers (PBDEs), are also considered important contributors to metabolic dysregulation2. PBDEs are widely used as brominated flame retardants and have been implicated in metabolic toxicity due to their potential for bioaccumulation and capacity to disrupt lipid metabolism3. However, individual responses to PBDEs vary considerably, which cannot be explained solely by exposure levels4. This heterogeneity, along with inconsistent epidemiological findings linking PBDEs to metabolic diseases5-8, underscores the pressing need to identify the determinants of individual susceptibility.

While genome-wide association studies (GWPAS) have identified numerous genetic variants associated with metabolic traits, these findings explain less than 30% of disease heritability9,10. This points to the potential role of gene-environment (G × E) interactions, where genetic background modulates responses to environmental exposures11,12. Initial G × E interaction research has focused on established disease-related loci, such as the obesity-related gene FTO, and their interactions with lifestyle factors like diet and physical activity13,14. Emerging evidence also indicates that environmental pollutants, including dioxins and other endocrine disruptors, may mimic endogenous ligands and directly alter proteins encoded by risk variants, thereby disrupting metabolic pathways15,16. Recent genome-wide G × E interaction studies in neurological phenotypes have further highlighted how environmental factors influence epigenetic regulation of gene expression17. However, systematic investigations of G × E interactions involving metabolic traits, particularly in relation to PBDEs, remain limited, impeding a holistic grasp of the dynamics and underlying mechanisms of G × E interactions in metabolic health.

To address this gap, high-throughput metabolomics offers a powerful approach in mechanism insight, as perturbed metabolite profiles can reveal pathway-specific disruptions and serve as early dysfunction biomarkers18. Notably, alterations in amino acid metabolism, such as elevated branched-chain amino acids (BCAAs) and reduced glycine have been implicated in pollutant-induced metabolic dysregulation19,20. BCAAs may exacerbate insulin resistance through mTOR pathway activation, while glycine deficiency can impair glutathione synthesis and cholesterol metabolism, exacerbating oxidative stress and dyslipidemia21-24. Nevertheless, the specific metabolic pathways mediating the interaction between PBDEs and specific genetic variants leading to lipid dysregulation remain unclear.

Against this backdrop, our study is driven by the hypothesis that PBDE exposure interacts with genetic variants to disrupt lipid homeostasis via metabolic pathway alterations. Using the China National Human Biomonitoring (CNHBM) cohort, which integrates exposome, multi-omics, and phenomics profiles in China25, we aimed to systematically identify genetic variants that modulate PBDE-induced dyslipidemia, delineate the metabolic pathways mediating these gene-PBDE interactions, and provide mechanistic insights of individual susceptibility using single-base-edited cellular models (Fig. 1). By bridging epidemiological findings with experimental validation, this study seeks to clarify G × E mechanisms in metabolic disorders, uncover biomarkers for risk stratification, and inform precision intervention strategies to mitigate pollutant-related metabolic risks.

Results
Participant enrollment and basic characteristics
The 2017-2018 CNHBM baseline survey included 21,888 participants aged 3-79 years, with 10% randomly selected for multi-omics profiling (Supplementary Fig. 1a). From this sub-cohort, we restricted our analysis to adults with available lipid profiles and excluded individuals with missing data on PBDE exposure, genomic, metabolomic, or key covariate measurements, resulting in a final analytic dataset comprised 871 individuals (Supplementary Fig. 1b). Concentrations of 14 PBDE congeners were measured, with only 2,2′,4,4′-Tetrabromodiphenyl ether (BDE-47) and 2,2',4,4',5,6'-Hexabromodiphenyl ether (BDE-154) detected in > 50% of samples (Supplementary Fig. 1c). We calculated the sum of all 14 PBDEs (ΣPBDE) to represent overall exposure and stratified participants into high- and low-exposure groups based on the median ΣPBDE concentration (0.816 μg/L). The two exposure groups showed comparable distributions of age, sex, body mass index (BMI), and lifestyle factors (smoking, alcohol consumption). However, the high-exposure group demonstrated significantly increased levels of total cholesterol (TC; P = 0.015) and triglycerides (TG; P = 0.002) and decreased levels of high-density lipoprotein cholesterol (HDL-C; P = 0.004), along with marginally increased low-density lipoprotein cholesterol (LDL-C; P = 0.139), as depicted in Table 1. Association analyses revealed statistically significant relationships between ΣPBDE and the measured lipids except for LDL-C, with similar patterns observed for BDE-47 (Supplementary Fig. 1d).

Interactions of genetic variants and PBDE exposure on lipid profiles
To investigate genetic variants modulating lipid profiles in response to PBDE exposure, we performed genome-wide G × E interaction analysis using adjusted general linear regression models. Manhattan plots depicted the genomic distribution of significant PBDE-interacting single nucleotide polymorphisms (SNPs) at a false discovery rate (FDR) < 0.05 (Fig. 2a). We identified 1,108 SNPs across 178 loci exhibiting significant interactions on TC levels, with rs1392296, annotated to SLC6A20 (Solute Carrier Family 6 Member 20), showing the strongest signal (FDR = 0.0078). SLC6A20 encodes a membrane transporter in the SLC6 family, facilitating molecular transport across cellular boundaries. Compared to individuals with the AA genotype at rs1392296 under low PBDE exposure, those with the GG genotype displayed TC increases of 0.45 unit at low exposure and 1.18 unit at high exposure (Fig. 2b), corresponding to a standardized effect size of 1.13-fold standardized deviation (SD) of TC. Similarly, PBDE exposure interacted with 955 SNPs (177 loci), 769 SNPs (140 loci), and 739 SNPs (127 loci) to influence LDL-C, HDL-C, and TG levels, respectively. The top SNPs from these analyses were annotated to biologically relevant genes: CETP (FDR = 0.0089), which encodes cholesteryl ester transfer protein critical for lipoprotein remodeling26; CPNE4 (FDR = 0.0104), a calcium-dependent phospholipid-binding protein potentially involved in lipid metabolism27; and CAMK1D (FDR = 0.0092), a calcium-dependent protein kinase implicated in hormonal regulation of metabolic processes28. These results suggest that PBDE exposure alters lipid profiles through specific genetic interactions, highlighting the critical role of G × E interactions in lipid metabolism regulation.

Building on these findings, we constructed lipid-specific polygenic risk scores (PRS) incorporating independent PBDE-interacting SNPs to evaluate their predictive utility for dyslipidemia risk. In populations with high PBDE exposure, the PRS exhibited significant positive associations with elevated TC, LDL-C, and TG (P = 0.001, 0.026, and 0.003, respectively) and a negative association with HDL-C (P = 0.001); whereas these associations were absent under low exposure (Fig. 2c). These data suggest that higher PRS values amplify dyslipidemia susceptibility in the context of elevated PBDE exposure, affirming the role of genetic predisposition in modulating pollutant-induced metabolic dysregulation.

Functional characterization of PBDE-interacting SNPs
To elucidate the gene-level mechanisms by which genetic variants affects lipids under PBDE exposure, we characterized significantly interactive SNPs basing on their genomic locations and molecular functions. For comparative analysis, a control set of non-interactive SNPs was generated, matched to the interactive SNPs for linkage disequilibrium, minor allele frequency (MAF), and variant type. Genomic location annotation revealed that approximately 90% of the interactive SNPs resided in non-coding regions, predominantly intronic and intergenic regions (Fig. 3ab). These SNPs exhibited significant enrichment in histone marks linked to epigenetic regulatory elements, including enhancers (H3K27ac, H3K4me1), and other active marks (DNase I, H3K9ac, H3K4me2), while showing depletion in repressive marks (H3K27me3) (Fig. 3c). Analysis of transcription factor (TF) binding sites further indicated significant enrichment of motifs for key TFs, such as BHLHE40, CTCF, FOXA1, and CEBPB (Fig. 3d). These results suggest that interactive SNPs predominantly modulate lipid phenotypes by regulating epigenetic elements and gene transcription.

We subsequently employed an integrative framework, combining nearest-gene, expression quantitative trait loci (eQTLs), and MPAGMA annotations, to identify potential target genes regulated by these SNPs (Fig. 3e). For the four lipid phenotypes, we identified 1,241, 1,261, 957, and 921 interacting genes, respectively (Supplementary Data 1-4). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis delineated distinct metabolic roles: TC-related genes were enriched in cholesterol biosynthesis, LDL-C genes in lipid metabolism regulation via APOB-mediated lipoprotein assembly, HDL-C genes in reverse cholesterol transport, and TG genes in fatty acid β-oxidation (Fig. 3f). Gene Ontology (GO) term enrichment further clarified the functional dimensions of these pathways (Fig. 3g). In biological processes, steroid metabolism and lipid metabolic regulation are significantly enriched. In molecular functions, prominent LDL particle binding and cholesterol transport activity provide a functional context for CPNE4 and CAMK1D in intracellular lipid signaling and metabolism. In cellular components, significant enrichment of lipoprotein particles and cholesterol synthesis complexes supports the structural basis for these processes. Together, these results indicate that PBDE-interacting SNPs exert regulatory influence over gene expression within key metabolic pathways, providing critical mechanistic insights into their contributions to lipid metabolism under PBDE exposure.

Metabolic pathways in dyslipidemia by Gene × PBDE interactions
Metabolomics offers a powerful approach to investigate how metabolic processes in organisms respond to varying environmental and genetic conditions by thoroughly examining shifts in small-molecule metabolites29. In order to identify the metabolic pathways and key metabolites involved in the interplay of PBDE exposure and susceptible variants, we first screened the differential metabolites associated with an interaction index. The index was derived by multiplying PBDE exposure levels by the PRS constructed above for each lipid phenotype. Our analysis revealed 96, 89, 66, and 39 metabolites associated with the G × E interactions for the four lipid traits, respectively, with approximately one- to two- thirds of these metabolites also directly linked to the corresponding lipid phenotypes (Fig. 4a). Interestingly, the number of metabolites implicated to these interactions exceeded those associated solely with PBDE exposure, indicating a broader capture of responding molecules of PBDE exposure (Supplementary Fig. 2ab).

Among the 63 overlapping metabolites identified for TC, citrate, glycine, and bile acid were particularly prominent, pointing to disruptions in energy metabolism and the bile acid cycle that may impair cholesterol synthesis and catabolism. These metabolites also exerted significant effects across the other three lipid types, underscoring their central involvement in G × E-driven lipid dysregulation. For LDL-C, metabolites like arginine/proline, behenoyl dihydrosphingomyelin, and methionine sulfoxide stood out, likely contributing to LDL-C perturbations by influencing oxidative stress and vascular function. In the case of HDL-C, arginine and glycine levels were notably diminished in response to the G × E interaction index, yet they exhibited a positive correlation with HDL-C concentrations. This suggests that G × E interactions may compromise vascular endothelial function, exacerbate oxidative damage, impair reverse cholesterol transport, and reduce HDL functionality. For TG, key metabolites including glycerol 3-phosphate, glycocholate, and purine were identified, potentially worsening hypertriglyceridemia by promoting glycerophospholipid synthesis and purine-related lipid accumulation under PBDE exposure. These findings collectively underscore the complex metabolic alterations driven by the interplay of environmental exposure and genetic predisposition across lipid profiles.

Subsequently, we conducted KEGG enrichment analysis on the overlapping metabolites to pinpoint key metabolic pathways (Fig. 4b). For both TC and LDL-C, the enriched pathways included glycine metabolism, alanine metabolism, and glutathione metabolism, indicating that glycine’s involvement in antioxidant defense and amino acid metabolism may exacerbate cholesterol dysregulation. In the case of TG, enriched pathways such as glycerophospholipid metabolism and the citrate cycle emerged, suggesting that heightened lipid synthesis and energy dysregulation contribute to TG accumulation. We also performed pathway enrichment analysis on all differential metabolites as a supplement. Beyond validating the primary pathways from the overlapping metabolite analysis, the results further identified key pathways enriched across all differential metabolites, notably arachidonic acid metabolism and unsaturated fatty acids biosynthesis (Supplementary Fig. 2c). This further broadens our understanding of the G × E interaction metabolic network.

To further elucidate how specific genetic variants disrupt lipid metabolism through these metabolites, we performed mediation analyses to assess their roles in the associations between PBDE-interacting variants and elevated dyslipidemia risk. This analysis identified 235 significant genetic variants influencing lipid phenotypes via 13 critical metabolites (Fig. 4c). Notably, glycine showed a substantial mediation effect across all four lipid traits, mediating the impact of 11, 12, 8, and 9 SNPs, respectively, including the most significant SNP for TC (rs1392296). Behenoyl dihydrosphingomyelin demonstrated the strongest mediation effect (18.62%) on LPDL-C disruption linked to rs12487683, possibly by altering sphingolipid metabolism and membrane dynamics. Additionally, arginine metabolism, encompassing arginine, proline, leucylhydroxyproline, and related compounds, mediated impairments in TC and HDL-C, potentially through alterations in nitric oxide production and protein metabolism. These key metabolites and pathways also displayed significance in the direct associations between PBDE exposure and dyslipidemia (Supplementary Fig. 2de). Consistent with leading candidates from differential analyses and pathway enrichment, these results collectively emphasized the critical role of metabolic pathways, particularly those involving glycine, citrate, and lipid-related intermediates, in driving lipid disturbances driven by G × E interactions.

Furthermore, we expanded the scope of our mediator analysis to include all G × E-associated differential metabolites as potential mediators. This supplementary analysis identified 14 mediating metabolites (Supplementary Fig. 2f). These metabolites primarily involve previously understudied pathways of inflammation and amino acid metabolism disorders, providing more comprehensive mechanistic insights into G × E interactions affecting lipid metabolism.

Fine mapping of the top variant in hypercholesterolemia
To investigate the molecular mechanisms by which PBDE-interacting variants disrupt lipid metabolism, we first screened 20 top loci (the top five for each lipid phenotype) from Fig.2a and prioritized those with the most significant effects based on two core indicators: FDR and absolute standardized β (interaction β per SD of outcome) (Supplementary Data 5). As a result, the locus with tag SNP rs1392296 outperformed others, exhibiting the lowest FDR (0.0078) and highest |β| (0.6517). The latter indicator suggests a relatively high impact of PBDE-SNP interaction on TC. Within this locus, seven another significant SNPs displayed strong linkage disequilibrium (LD) (r² ≥ 0.6) with the tag SNP (Fig. 5a). We calculated an integrated functional prediction score for these SNPs by combining outputs from CADD, RegulomeDB, and 3DSNP, which identified rs9869609 (r² = 0.991 with the tag SNP, MAF = 6%) as the candidate causal SNP due to its greatest potential for functional regulation (Fig. 5b). Genomic visualization further confirmed that rs9869609 is significantly enriched in DNase hypersensitive sites and active histone modification peaks (H3K4me1, H3K4me3, H3K9ac, and H3K27ac), underscoring its robust regulatory potential (Fig. 5c). Therefore, we prioritized functional validation of this variant to explore how SNPs interact with PBDE to drive lipid dysregulation.

The eQTL analysis based on the GTEx database (v8) indicated that rs9869609 could significantly downregulate gene expression of SLC6A20 in lipid-metabolic tissues, such as skeletal muscle (β = -0.09, P = 0.0098) (Fig. 5d). This variant is located in an intergenic region, approximately 347 kb downstream of the SLC6A20 gene, suggesting its distal cis-enhancer role to regulate the gene expression. We further confirmed the population-level effect of this variant, in which significant interaction effects were observed for TC (FDR = 0.0082) and LDL-C (FDR = 0.0104) levels. Especially, under high PBDE exposure, individuals carrying the GG genotype exhibited substantially altered lipid profiles compared to those with AA genotype, with hypercholesterolemia prevalence increasing from 16.8% to 26.1%, while the change was neglectable under low exposure (Supplementary Fig. 3ab).

To validate this regulatory effect, particularly under PBDE exposure, we conducted dual-luciferase reporter assays in HepG2 cells, comparing the A and G alleles of rs9869609. We used BDE-47, a prevalent PBDE homolog with established epidemiological ties to dyslipidemia, as the exposure agent. Constructs with the rs9869609-G allele exhibited lower transcriptional activity than those with the A allele. Exposure to BDE-47 further reduced activity for both alleles, with a more pronounced decrease observed for the G allele (Fig. 5e), suggesting that BDE-47 amplifies the suppressive effect of the G allele on SLC6A20 transcription.

To delve deeper into the regulatory mechanisms, we considered that allele-specific activity at rs9869609 might arise from differential TF binding affinities30. TF motif analysis using JASPAR identified BHLHE40 as a candidate TF with a preference for binding the G allele (Fig. 5f). BHLHE40, a member of the Basic Helix-Loop-Helix family, is a recognized transcriptional repressor involved in metabolic dysfunction and alcohol-related liver disease31.To test whether BHLHE40 mediates the allele-specific effects of rs9869609 (A > G), we performed reporter gene assays in HepG2 cells following siRNA-mediated knockdown of BHLHE40. The knockdown resulted in increased luciferase activity for constructs carrying either the A or G allele, with minimal differences between them. However, when exposed to BDE-47, the activity of both alleles was suppressed again, with the G allele showing a more significant reduction compared to the A allele (Fig. 5g). Collectively, these results indicate that rs9869609-G functions as an allele-specific silencer, enhancing the suppression of SLC6A20 transcription in the presence of BDE-47 through stronger binding of BHLHE40 to the G allele.

Mechanism validation in BDE-47-exposed rs9869609-edited cells
To investigate the downstream responsive genes of BDE-47 exposure, we performed RNA-seq in HepG2 cells following BDE-47 exposure (12.5 μM, 24 h). As a result, we identified 1,221 significantly differentially expressed genes (DEG), including SLC6A20, CETP and CAMK1D (Supplementary Fig. 4a). KEGG analysis revealed that the DEG were primarily enriched in cholesterol metabolism and lipoprotein biosynthesis pathways, indicating that PBDE exposure disrupts lipid homeostasis through these key metabolic pathways (Supplementary Fig. 4b). Notably, SLC6A20 was among the downregulated genes (log (FC) = -2.21, FDR = 0.0032), pointing to its role in the early cellular response to PBDE exposure.

We then utilized CRISPR/Cas9 genome editing to create HepG2 cell lines with AA, AG, and GG genotypes at rs9869609 (Fig. 6a) and evaluated the interaction effects on SLC6A20 expression. Cells with the GG genotype showed reduced SLC6A20 mRNA and protein levels compared to those with the AA genotype. Upon BDE-47 exposure, SLC6A20 expression reduced in both AA and GG cells, with a more pronounced reduction observed in the GG genotype relative to AA (Fig. 6b). These results suggest that BDE-47 amplifies the suppressive impact of the rs9869609-G allele on SLC6A20 expression.

Drawing from our earlier analyses, rs9869609, along with the tag SNP, interacts with PBDE exposure to disrupt cholesterol levels, probably through metabolites including glycine, which is transported by SLC6A20. To test this effect in gene-edited cells, we measured cholesterol and glycine levels in cells with different rs9869609 genotypes under BDE-47 exposure. GG Cells showed elevated cholesterol levels compared to those with the AA genotype, with a more significant increase following BDE-47 exposure (Fig. 6c). Similarly, intracellular glycine levels were lower in GG cells than in AA cells, with a more substantial decrease observed upon BDE-47 exposure (Fig. 6d). These results confirm the G × E interaction that impairs glycine availability and exacerbates cholesterol dysregulation.

SLC6A20, a key membrane protein for glycine transport, may hinder cholesterol degradation via the bile acid pathway by reducing intracellular glycine, an essential co-factor in forming glycine-conjugated bile acids for cholesterol excretion23,32,33. To clarify SLC6A20’s role in this G × E interaction, we transfected HepG2 cells with siRNAs targeting SLC6A20 (Fig. 6e) and examined the downstream consequences. As expected, SLC6A20 knockdown markedly elevated cholesterol levels compared to control cells. Upon BDE-47 exposure, cholesterol levels increased in both knockdown and control groups, with a more pronounced rise in the knockdown cells relative to their non-exposed counterparts (Fig. 6f). Likewise, intracellular glycine levels dropped significantly in the SLC6A20 knockdown group, and BDE-47 exposure further reduced glycine levels in both knockdown and control cells, with the knockdown group exhibiting the lowest levels compared to non-exposed knockdown cells (Fig. 6g). To evaluate the impact of SLC6A20 on glycine-conjugated bile acids, which are essential for cholesterol excretion, we measured glycocholic acid (GICA) and glycodeoxycholic acid (GDCA) levels in SLC6A20 knockdown HepG2 cells. Ultra-High-Performance Liquid Chromatography-High Resolution Mass Spectrometry (UPLC-HRMS) analysis revealed significantly reduced levels of GCA and GDCA concentrations in siSLC6A20 cells compared to siControl cells. Following BDE-47 exposure, GCA and GDCA levels exhibited a further decrease in siSLC6A20 cells, while siControl cells showed a more moderate reduction (Fig. 6hi). Collectively, these findings suggest that BDE-47 enhances the rs9869609-G allele’s suppression of SLC6A20, impairing glycine transport and driving cholesterol accumulation. This provides a mechanistic framework for how G × E interactions disrupt lipid homeostasis.

Methods
Ethics
The study complies with all relevant regulations and the Declaration of Helsinki and was approved by the Ethical Review Committee of the National Institute of Environmental Health (NIEH), Chinese Center for Disease Control and Prevention (China CDC) as part of the CNHBM cohort (no. 201701). Moreover, written informed consent was obtained from each participant before the study began. For all participants who were minors, written informed consent was obtained from their parent or legal guardian, in accordance with CNHBM ethical protocols and regulations.

Data sources and study participants
This study leverages data from a subset of the China National Human Biomonitoring (CNHBM) baseline cohort, a cross-regional survey that enrolled 21,888 participants aged 3-79 years, representing multiple ethnicities including Han (86.4%), Tibetan, Zhuang, Mongolian, and others. For the present multi-omics analysis, a secondary sample of 1,824 participants was selected from the baseline cohort using stratified sampling across all age groups to maintain demographic representativeness. Peripheral venous blood samples were collected in a standardized manner from all enrolled individuals, immediately processed via centrifugation to separate serum and blood cells. Serum samples were utilized for metabolomics analysis as well as exposure metabolomics, while blood cells were reserved for genomics analysis. Detailed protocols for participant recruitment, sample collection and processing have been described previously25. Following the exclusion of participants under 18 years and those with incomplete data, the final analytical sample consisted of 871 individuals (Fig 1, Supplementary Fig. 1ab).

Multi-omics Data Generation
Exposure Data: Fourteen serum PBDE congeners (BDE-17, 28, 47, 66, 71, 85, 99, 100, 138, 153, 154, 183, 190, 209) were analyzed via the CDC-established method on an Agilent 7890B gas chromatograph coupled to a Waters Autospec Premier dual-focus magnetic mass spectrometer (VSIR mode, resolution 10,000). Serum samples were pretreated with acid-methanol, extracted with organic solvents, and extracts were concentrated, SPE-purified and reconstituted prior to analysis. Instrumental separation was performed on a capillary column (helium carrier gas) with pulsed splitless injection and positive EI (35 eV) ionization. The limit of detection (LOD) of these pollutants ranged from 0.01 to 1.20 μg/L, with non-detectable measurements imputed as LOD/√2. Detailed LODs and median concentrations are presented in Supplementary Data 7. Participants were stratified into high- and low-exposure groups by the median of total PBDE levels (ΣPBDE; 0.8157 μg/L).

Genomics Data: A 1.5-mL peripheral blood sample was collected from each subject at baseline for DNA extraction. Samples were genotyped using the Asian Screening Array-China Health Industry Alliance (ASA-CHIA) BeadChip (Illumina, Inc., San Diego, California, USA), following the Illumina Infinium HTS chip standard operating procedures63. After applying quality control measures, raw genotyping data for ~700,000 SNPs were retained for each sample, with alignment to the human genome Hg19 assembly. Genotype imputation was performed using the Michigan Imputation Server (https://imputationserver.sph.umich.edu). SNPs were excluded based on the following criteria: (1) imputation quality score < 0.2; (2) MAF < 0.05; (3) SNP call rate < 95%; (4) Hardy-Weinberg equilibrium P-value < 1 × 10-6.

Untargeted Metabolomics Data: Metabolomic profiling of serum was conducted on the CalOmics platform (Calibra Scientific Inc.). A total of 1824 human serum samples were analyzed as single injections without technical replicates. Sample preparation involved protein precipitation and metabolite extraction using methanol, and quality control (QC) measures included pooled QC samples, process blanks, and deuterated internal standards. The analysis was conducted using an ACQUITY 2D UPLC system (Waters) coupled to a Q Exactive Quadrupole-Orbitrap mass spectrometer (Thermo Fisher Scientific). Four complementary liquid chromatography- mass spectrometry (LC-MS/MS) methods were employed: two using C18 columns under different mobile phase conditions for electrospray ionization positive mode (ESI⁺, denoted as POSa and POSb) and negative mode (ESI⁻, denoted as NEGa), and one using a HILIC column (BEH Amide) for ESI⁻ mode (denoted as NEGb). Full scans were acquired at 35,000 resolutions (m/z 70–1000). Raw data were processed using the in-house software MetaPro64 for peak picking, alignment, and integration. Metabolites were identified with Metabolomics Standards Initiative (MSI) Level 1 confidence65,66 by matching retention time (±0.1 min), exact mass (error <10 ppm), and MS/MS spectra (match score >75%) against an in-house library of authentic standards.

Phenotypes and Covariates
The primary outcome variables in this study were four blood lipid components: TC, LDL-C, HDL-C, and TG. Blood samples were obtained through venipuncture, immediately centrifuged, and the resulting serum was analyzed using a standardized biochemical analyzer to determine the concentrations of these four lipid markers. To address right-skewed distributions and meet the assumptions of parametric statistical tests, the four lipid phenotypes underwent logarithmic transformation (log10).

To adjust for potential confounding factors, the following covariates were included in all statistical models: sex (male and female), age group (19- 39, 40- 59, and ≥ 60 years), BMI (< 18.5, 18.5- 25, 25- 30, ≥ 30 kg/m2), smoking status (smoked, defined as at least once a week, and never smoked), and alcohol consumption (drinking, defined at least once a week, and never drinking). BMI was derived from measured height and weight (kg/m²). These covariates were included as adjustment factors in the statistical models to minimize potential biases in assessing the relationship between PBDE exposure and lipid outcomes.

Genome-wide G × E interaction analysis
To examine the associations between PBDE exposure, SNPs, and the four lipid profiles (TC, LDL-C, HDL-C, and TG), we utilized a generalized linear regression model, implemented in PLINK (v.1.90), to evaluate the interaction effects of PBDE exposure and SNPs on lipid levels. As demonstrated by recent studies67,68, such G × E interactions allow for the identification of interaction loci that respond to environmental exposure and modulate blood lipid levels, thereby facilitating the characterization of individuals at high risk of dyslipidemia associated with high environmental exposure. We then generated Manhattan plots to illustrate the genomic distribution of significant loci (FDR < 0.05), with the top five SNPs annotated to their corresponding genes; for the most significant SNP in each plot, a stratified analysis of dyslipidemia risk was performed to highlight the interaction effect.

PRS construction for lipid phenotype
To assess the cumulative interaction between lipid-specific genetic predisposition and PBDE exposure on elevated levels of each lipid phenotype, we constructed a PRS using significant SNPs identified from the prior interaction analysis. Independent SNPs were selected for the models through a linkage disequilibrium-driven clumping approach in PLINK (v.1.90), with parameters set as follows: --clump-p 0.0001, --clump-r2 0.20 --clump-kb 250. A total of 172, 168, 130, and 124 independent SNPs were included in the PRS model for the four lipid phenotypes (TC, LDL-C, HDL-C, and TG), respectively. Participants were then scored additively in a weighted manner based on the number of risk alleles they carried, utilizing the clumping-plus-thresholding (C+T) method according to the formula:
PRS = ∑n i=1 βi SNP i (1)
where n represents the total number of variants, SNP i (0, 1, or 2) denotes the count of risk alleles (corresponding to homozygous non-risk allele, heterozygous allele, or homozygous risk allele, respectively), and β i is the effect size of the risk allele69. Using the lowest quintile (bottom 20%) of the PRS as the reference group, increased levels of corresponding lipid were estimated across the remaining four PRS quintiles in linear regression models, stratified by PBDE exposure levels. The dose-response relationship was assessed via trend analysis, with statistical significance defined as P trend < 0.05.

Annotation and functional enrichment of SNPs
To investigate the enriched characteristics of PBDE-interacting SNPs, a control variant set was generated with a 1:1 match to the interactive SNPs for each lipid phenotype, ensuring equivalence in MAF, number of LD SNPs, and variant type. This was achieved using the web-based tool vSampler (version 8.0, http://mulinlab.org/vsampler/). Both interactive SNPs and their control variants were annotated using SnpEff software (version 8) with the GRCh37 reference genome from ENSEMBL. Variants were categorized into 5′ untranslated region (5’UTR), 3′ untranslated region (3′UTR), upstream, downstream, intron, and intergenic region (IGR). Functional annotation data, including histone modification ChIP-seq peaks (H3K27ac, H3K27me3, H3K4me1, H3K4me2 and H3K9ac), DNase I hypersensitive sites, and binding site of 59 transcription factors, were retrived from the ENCODE portal (https://www.encodeproject.org/data/annotations/)70.

To identify the putative interactive genes, we employed an integrated framework for assigning target genes to interactive SNPs (Fig. 3e): (1) Physical distance strategy -- the nearest genes within a 10kb region surrounding the interactive SNPs were included as candidates; (2) eQTL strategy -- eQTL data from the GTEx database were leveraged to pinpoint target genes whose expression levels were affected by interactive SNPs; (3) MAGMA strategy -- target genes were determined by aggregating variant signals to the gene level using Multi-marker Analysis of GenoMic Annotation (MAGMA)71. Genes identified through this approach were designated as PBDE-interacting genes. To explore their biological significance, KEGG and GO term enrichment analyses72 were performed, with the top enriched pathways visualized.

Fine mapping analysis
We next performed a functional annotation for the genetic variants using multiple bioinformatic tools and databases, including the HaploReg database (https://pubs.broadinstitute.org/mammals/haploreg/haploreg.php), Roadmap Epigenomics Consortium (https://www.ncbi.nlm.nih.gov/geo/roadmap/epigenomics/), RegulomeDB database (https://www.regulomedb.org/regulome-search), CADD database (https://cadd.gs.washington.edu/score), 3DSNP (https://omic.tech/3dsnpv2/), which integrated multiple histone modification chromatin immunoprecipitation (ChIP)-seq peaks, transcription factor (TF) ChIP-seq peaks. To be specific, (1) regulatory chromatin histone ChIP-Seq were characterized by using the data from the Roadmap Epigenomics Consortium; (2) prioritizing these variants with the scores of regulatory elements by using RegulomeDB, CADD, and 3DSNP prediction databases73. CADD first filters out potentially harmless variants, focusing on functionally significant sites74. RegulomeDB then determines whether these sites reside within regulatory regions and assesses their regulatory potential75. Finally, 3DSNP connects variants within regulatory regions to their actual distal target genes, completing the full biological narrative from variant to target gene76. Specifically, we apply minimum-maximum normalization to the raw outputs of these databases, uniformly converting their scores to the 0–1 range.

Cell culture and cell viability detection
Biological experiments were conducted using the HepG2 cell line, a model widely employed in drug metabolism and hepatotoxicity studies due to its functional resemblance to human liver tissue77. The HepG2 cells were obtained from the China Center for Type Culture Collection (Wuhan, China), and cultured in Dulbecco’s Modified Eagle’s Medium (DMEM; Gibco, USA) supplemented with 10% fetal bovine serum (FBS; Gibco, USA) and 1% antibiotics (100 U/mL penicillin and 0.1 mg/mL streptomycin). Cells were maintained in a humidified incubator at 37 °C with 5% CO2. Cell line authenticity was verified through short tandem repeat profiling (Applied Biosystems), and the absence of Mycoplasma contamination was confirmed using the MycoAlert kit. The most recent testing data was in May 10, 2024.

In the population effects analysis, to better reflect real-world exposure conditions, we used total PBDEs as a quantitative indicator of mixture exposure to assess its impact on blood lipids. Analysis results indicate that the primary effects of BDE-47 in the population align with those of the total PBDEs mixture. In the experimental section, to streamline procedures and ensure representativeness of exposure levels, effect intensity, and mechanisms, we selected BDE-47 as the representative compound for this class of pollutants to infer potential mechanisms of action for similar pollutants78.

Cell viability was assessed using the CCK-8 assay following the manufacturer’s protocol. HepG2 cells were cultured in 96-well plates with exposure to gradient concentrations of BDE-47 (0, 6.25, 12.5, 25, 50, 100 μM) for durations of 24, 48, and 72 hours. Following exposure, the culture medium was replaced with medium containing 10% CCK-8 solution, and cells were incubated for an additional 50 minutes. Optical density (OD) at 450 nm was measured using a microplate reader (BioTek Epoch, USA), and cell viability (%) was calculated relative to the control group. Growth curves were plotted to assess the dose- and time-dependent effects of BDE-47 on HepG2 cell proliferation. Results indicated significant proliferation inhibition at 50 μM after 48-hour treatment (Supplementary Fig. 4c). Therefore, the 48-hour exposure to 25 μM BDE-47 was selected for subsequent toxicological assessments to balance observable effects with maintained cell viability.

Detection of lipids and metabolites
TC was quantified in HepG2 cells from BDE‑47 exposure and control groups (3 biological replicates per group, each with 3 technical replicates) using a commercial total cholesterol assay kit (Solarbio, Beijing, China). Intracellular glycine levels were determined in parallel with a glycine assay kit (Jiancheng Bioengineering Institute, Nanjing, China). Cell processing and analysis were performed according to the respective manufacturer’s protocols. All measurements were normalized to the total protein concentration of each sample to ensure accurate quantification.

Targeted UPLC-HRMS quantified intracellular GCA and GDCA. Methanol-extracted samples (n = 18; BDE-47-exposed vs. control HepG2 cells, 3 biologicals × 3 technical replicates) were separated on a Hypersil GOLD C18 column (0.1% formic acid/acetonitrile, 15-min gradient) and analyzed in negative-ion full-scan mode (m/z 70–1050, resolution 70,000) using a Dionex UPLC/Q Exactive Orbitrap system. External calibration curves (R² ≥ 0.99) were used for quantification, with pooled QC samples injected every 10 runs. Data were acquired with Xcalibur 4.0 and processed using TraceFinder 5.1 for peak integration and concentration calculation.

Dual-luciferase reporter assays
DNA fragments spanning 1,000 base pairs around the SNP rs9869609 (A or G allele) were subcloned into the pGL3-promoter vector (Promega, USA) and synthesized by Tsingke Biological Technology (Wuhan, China). For transient transfection assays, HepG2 cells were cultured in 96-well plates and co-transfected with constructed luciferase vector (containing either the rs9869609-A or G allele) and the pRL-SV40 Renilla luciferase plasmid (Promega, USA) using Lipofectamine 3000 (Invitrogen, USA). Twenty-four hours post-transfection, the control group was maintained in standard medium, while the experimental group was treated with medium containing 12.5 μM BDE-47 for an additional 24 hours. Luciferase activity was measured using the Dual-Luciferase Reporter Assay Kit (Promega, USA) according to the manufacturer’s protocol, with Renilla luciferase luminescence normalized to that of firefly luciferase for each sample.

RNA interference
For RNA interference, small interfering RNA (siRNA) oligonucleotides targeting BHLHE40 and SLC6A20, along with a non-targeting siRNA control, were obtained from Tsingke Biological Technology (Wuhan, China), with the sequences provided in Supplementary Data 8. Transfection was performed using Lipofectamine RNAiMAX (Invitrogen, USA), and the knockdown efficiency was validated by quantitative reverse transcription PCR (qRT-PCR) (Supplementary Fig. 4de).

RNA-seq
HepG2 cells were treated with 12.5 μM BDE‑47 or vehicle control for 24 hours (n = 3 biological replicates per group). Total RNA was extracted using TRIzol reagent (Invitrogen), and RNA integrity was verified (RIN > 8.0) on an Agilent Bioanalyzer. Ribosomal RNA was depleted, and stranded cDNA libraries were constructed following the manufacturer’s protocol (Vazyme). Libraries were sequenced on a DNBSEQ‑T7 platform (MGI Tech) in 150‑bp paired‑end mode. Raw reads were quality‑trimmed with fastp (v0.23.2) and aligned to the human reference genome (GRCh38) using HISAT2 (v2.2.1). Gene‑level counts were obtained with featureCounts (v2.0.1). Differential expression analysis was performed with edgeR (v3.40.2), applying a threshold of |log₂FC| > 1 and FDR < 0.05. The raw sequencing data have been deposited in the Genome Sequence Archive for Human (GSA-Human) under accession code HRA014233.

CRISPR/Cas9 mediated single-base mutation
CRISPR/Cas9-mediated editing of the rs9869609 variant was conducted in the HepG2 cell line. The single-guide RNA (sgRNA) was designed using the established CRISPR design tool (https://cctop.cos.uni-heidelberg.de/).The sgRNA was subcloned downstream of the human U6 promoter into the plasmid pSpCas9(BB)-2A-Puro (PX459) V2.0 (Cat# 48138, GeneAdv, Suzhou) via BbsI restriction sites. Additionally, a double-stranded donor (dsDonor) sequence (1.6 kb) containing the rs9869609 variant (A > G) was designed based on the highly efficient sgRNA and inserted into the pUC57 vector through the EcoRV restriction site. All oligonucleotides used are listed in Supplementary Data 8. Transfection was performed in HepG2 cells with 80% confluency. A total of 2.5 μg of the Cas9 plasmid (PX459-sgRNA) and 1 μg of dsDonor were co-transfected into cells using Lipofectamine 3000. After 48 hours, 0.8 mg/mL puromycin was added to select for transfected cells. Surviving cells were then trypsinized and seeded into 96-well plates via fluorescence-activated cell sorting (FACS) to isolate single-cell clones. These clones were subsequently subcultured and verified by Sanger sequencing over a period of 14-21 days.

QRT-PCR and western blotting (WB)
Total RNA was extracted from cells with Trizol Reagent (Thermo Fisher Scientific, USA). Reverse transcription was carried out with PrimeScriptTM RT-PCR Kit (Takara, Japan), followed by qPCR on a QuantStudio 5 qPCR system (Applied Biosystems, Thermo Fisher) using SYBRTM Green PCR Master Mix (Takara, Japan). Gene expression levels were normalized to GAPDH, which served as an internal reference. The corresponding primers were listed in Supplementary Data 8.

Total protein was harvested from HepG2 cells of BDE47 exposure and control groups (3 biological replicates per group) using RIPA lysis buffer (Beyotime, China) supplemented with protease inhibitors PMSF (Beyotime, China), cOmplete cocktail (Sigma, USA) and PhosSTOP (Sigma, USA). Protein was quantified using BCA reaction (Beyotime, China) and denaturated at 99 °C for 5 min. Equal amount of protein was separated by electrophoresis in 10% SDS-PAGE gels and transferred onto 0.45 mm PVDF membrane. Protein was incubated with antibodies against SLC6A20 (1:1000, Proteintech, Cat# 30843-1-AP) and β-actin (1:1000, Proteintech, Cat# 20536-1-AP) at 4 °C overnight. HRP-conjugated anti-rabbit IgG (1:5,000, Proteintech, Cat# SA00001-2) was used as secondary antibody. Chemiluminescence signal was developed with enhanced chemiluminescence substrate (Beyotime, China) by Image Lab software. Unprocessed scans of the key blots (for Fig. 6B) are supplied in the Source Data file, and the detail information of antibodies was provided in the Reporting summary.

Statistics & Reproducibility
This study employed R (version 4.4.1) to perform the statistical analysis of epidemiological and experimental data. The generalized linear regression model (R package ‘stats’), adjusted by covariates, was employed to analyze the metabolomic data, identifying metabolites associated with the G × E interaction index and corresponding lipid levels. The interaction index, calculated as the product of PBDE exposure and the PRS for each lipid phenotype, represented the overall interaction effect. Overlapping metabolite sets were further enriched using the KEGG database within MetaboAnalyst 3.079.To examine the mediating role of serum differential metabolites in G × E interactions with lipid phenotypes, causal mediation analyses were conducted using the R package ‘mediation’. After covariates adjustment, we calculated the total effect (TE), natural direct effect (NDE), and natural indirect effect (NIE), with the mediation proportion of significant metabolites determined as the ratio of NIE to TE, reflecting the mediation extent. For experimental data, two-group comparisons were analyzed using unpaired t-tests in Prism software (GraphPad). For multiple-group comparisons involving normally distributed data, one-way analysis of variance (ANOVA) was applied, followed by Fisher’s least significant difference (LSD) tests for post-hoc analysis. A P-value < 0.05 was considered statistically significant, with significance levels denoted in the corresponding figures.

Sample sizes were determined based on the predefined sample screening protocol for epidemiological and experimental cohorts. No statistical method was used to predetermine sample size. No data were excluded from the analyses. The experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessment. All key experimental assays were independently replicated at least three times, with consistent results confirming the reproducibility of findings.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Hu N. et al. conducted an interesting study on individual heterogeneity in response to environmental exposure, with a\nspecific focus on PBDE-induced dyslipidemia. By integrating exposome, genomics, and metabolomic data from the China\nNational Human Biomonitoring cohort, the authors identified 3,571 genetic variants that interact with PBDE exposure to\ninfluence dyslipidemia risk with key metabolites serving as mediators. Among these, rs9869609 was pinpointed as the top\ncausal variant, exhibiting the strongest association with hypercholesterolemia risk. Further functional validation through\nsingle-base CRISPR/Cas9 editing revealed that the rs9869609-G allele downregulates SLC6A20 expression by enhancing\nBHLHE40 binding. This, in turn, impaired glycine transport and promoted cholesterol accumulation, particularly under BDE-\n47 exposure.\nThis study provides a comprehensive genome-wide analysis of gene-environment interactions in metabolic diseases. It\nmakes a significant contribution by identifying vulnerable populations and clarifying the mechanistic pathways underlying\nenvironmentally induced metabolic disorders. Additionally, it opens new avenues for research on the health effects of\nenvironmental pollutants. However, several key limitations need to be addressed to improve the quality and impact of the\nmanuscript.\nMajor comments:\n1.The authors constructed a polygenic risk score (PRS) using independent variants to predict lipid phenotypes, finding\nincreased risk only in populations with high PBDE exposure. Why does the PRS not confer risk in low-exposure groups?\nThis observation warrants further discussion\n2.While the study identifies loci interacting with PBDE exposure, it does not clearly differentiate this approach from\ntraditional GWAS for disease-associated loci. Please elaborate on why gene-environment interaction (G×E) loci may offer\ngreater biological or clinical relevance than main-effect GWAS hits, particularly for environmental metabolic disorders.\n3.In metabolomic analysis, the authors focus on overlapping metabolites for pathway enrichment and mediation analysis,\npotentially missing broader mechanistic insights. A supplementary correlation analysis of all significant metabolites could\nreveal additional candidate molecules and pathways.\n4.The choice of rs9869609 (over the most statistically significant SNPs) lacks clear justification. Summing functional scores\nacross different ranges may introduce bias; a weighted or standardized approach would be more rigorous. Please clarify the\nprioritization criteria.\n5.The functional validation uses the maximum non-toxic dose of BDE-47. How does this dose compare to real-world\nexposure levels in the cohort? Are the observed effects plausible at environmentally relevant concentrations? Address\npotential limitations in extrapolating experimental findings to population exposures.\n6.The enrichment of histone marks suggests epigenetic regulation, yet the study focuses solely on TF-mediated gene\nexpression. Could other mechanisms (e.g., chromatin remodeling, splicing modulation) contribute? Please expand the\ndiscussion on alternative mechanisms.\nMinor comments:\n1.SLC6A20 is implicated in glycine transport and cholesterol degradation via bile acid pathways. Were bile acid-related\nmetabolites altered by PBDE exposure? Include relevant data or discuss their potential role.\n2.HepG2 cells were used for validation—are these findings generalizable to primary hepatocytes or in vivo systems? Briefly\njustify the model choice and discuss translational limitations.\n3.The term \"dyslipidemia\" is not defined, and the analysis uses quantitative lipid markers. How do these findings translate to\nclinical dyslipidemia risk? Clarify the connection.\n4.Why total cholesterol was chosen over other lipid traits for experimental validation? Were the loci for different lipids\ncompared? Discuss their interrelationships.\n5.Provide basic genomic details of rs9869609 (e.g., distance to SLC6A20, intronic/intergenic status, MAF) to contextualize\nthe variant’s potential functional impact.\n6.The PBDE-specific PRS lacks internal/external validation. Discuss its utility for future research (e.g., risk prediction,\npersonalized interventions).\n7.The experiments used BDE-47 as exposure, while analyses modeled total PBDEs. Explain this discrepancy and its\nimplications for interpretation.\n8.Other typographical and formatting Issues:\n·In Fig. 2c, specify units for “Beta of blood lipids” in PRS models.\n·In Fig. 6, correct \"gpro\" to \"g pro.\" in metabolite/lipid labels and describe normalization methods.\n·In Fig. 4b, indicate the significance threshold for pathway enrichment.", "reviewer_2": "This study investigates the interaction between genetic variants and polybrominated diphenyl ether (PBDE) exposure in\nmodulating lipid homeostasis, combining epidemiological data from the China National Human Biomonitoring (CNHBM)\ncohort with multi-omics profiling and functional validation. The authors performed gene-environment (G×E) interaction\nanalysis, identifying in SNPs associated with PBDE-induced dyslipidemia, with rs9869609 (near SLC6A20) emerging as\na top candidate (among others). Metabolomic analysis highlighted glycine and glycerophosphate as key mediators, while\nCRISPR/Cas9-edited HepG2 cells demonstrated that the rs9869609-G allele suppresses SLC6A20 expression via\nenhanced BHLHE40 binding, impairing glycine transport and exacerbating cholesterol accumulation under BDE-47\nexposure.\nThe article is highly compelling due to its translational approach, bridging population-level findings with mechanistic\ninsights. The identification of SLC6A20 as a susceptibility marker for PBDE-induced dyslipidemia and its validation through\nknockout experiments in HepG2 cells strengthens the biological plausibility of the results. The study is well-structured,\nclearly written, and addresses an important gap in understanding how environmental pollutants interact with genetic\nbackground to influence metabolic health.\nWhile the findings are robust, some limitations should be noted. The functional validation was conducted only in vitro, and\nfuture studies using in vivo models (e.g., genetically modified mice or humanized animal systems) would further solidify the\nconclusions. Additionally, the generalizability of the findings need to confirm in other cohorts, and the clinical relevance of\nthe polygenic risk score (PRS) warrants further investigation.\nDespite these limitations, the study represents a significant advance in environmental health research, providing novel\ninsights into gene-environment interactions in dyslipidemia. The work is of broad interest to the fields of toxicology,\nmetabolism, and precision medicine. With minor revisions to address, this article is suitable for publication.\nGeneral remarks\n1.Clarification of focus on SLC6A20\n1.1)The authors focus on SLC6A20 due to its strong statistical association with TC (rs1392296, FDR=0.0078) and\ncompelling functional evidence linking its suppression to glycine-mediated cholesterol metabolism. However, other\ncandidate genes (e.g., CTCF, enriched across all lipid traits) are not explored experimentally. While prioritization is\ninevitable, a brief discussion justifying why SLC6A20 was selected over equally plausible candidates—such as CTCF, a\nglobal epigenetic regulator including of lipid-related genes, or CETP, directly involved in lipoprotein remodeling—would\nstrengthen the manuscript. Future studies could address whether these genes act independently or synergistically with\nSLC6A20 in PBDE-induced dyslipidemia.\n1.2)The polygenic risk score incorporates multiple loci, yet only SLC6A20 is functionally tested. A brief explanation of why\nother top genes were excluded (e.g., lack of eQTL support, lower effect sizes) would bolster transparency.\n1.3) While the role of SLC6A20 in glycine transport is compelling, dyslipidemia is multifactorial. The study does not rule out\nadditive effects of other candidate genes. Please discuss\n2. Association between PBDE, SLC6A20 and dyslipidemia\n2.1) The study lacks RNA-seq or single-cell transcriptomics to validate SLC6A20 expression changes in human tissues or\nadditional PBDE-responsive genes.\n2.2) The other top genes identified in genetic variants could had been investigated in HeparG cells.\n2.3) The claim that rs9869609-G is the \"top causal variant\" is premature without: In vivo validation (e.g., knock-in mice).\nPlease temper this in the article.\nEvidence that this variant explains a significant proportion of dyslipidemia risk in the population. The PRS improves\nprediction but has modest effect sizes (β=1.18 for TC).\nHere is some piece of advice that the authors could take into consideration and disuss : Clarify Exposure Relevance: It could\nbe interesting to compare in vitro BDE-47 doses to human urinary levels and assess dose-response relationships.\nAddress Population Bias: Replicate findings in diverse cohorts or perform trans-ethnic meta-analyses.\nStrengthen Mechanistic Claims:\nProvide direct evidence of bile acid conjugation defects in SLC6A20-deficient models.\nIntegrate transcriptomics to identify downstream effects of PBDE-exposed genotypes.\nMinor points\nFigure 2 Color Consistency Issue:\nIn panel (a), colors distinguish chromosomes, but in panels (b) (stratified risk) and (c) (PRS trends), the same/similar colors\nrepresent PBDE exposure groups (low/high). This is confusing.\nUse distinct color schemes (e.g., grayscale for chromosomes, blue/orange for exposure groups) or add explicit labels to\nlegends.\nFigure 4 (Metabolome Profiling)\nThere is a legend Mismatch in Panel (b3) (HDL-C):\nThe legend describes two circle sizes (likely representing effect size or significance), but three sizes appear in the plot\n(small/medium/large).\nPlease align the legend with the visual output\nOn panel c line widths effectively represent mediation strength, but key metabolites (e.g., glycine, citrate) could be\nhighlighted with bold labels or arrows to emphasize central pathways.\nTypos and suggestions\nLine 164: unit intead of units\nLine 168: Suggest rather than demonstrate\nLine 1026 Characterized\nLine 1085: H3K4me1 repeated twice"}, "author_response": "Response to Reviewers\nReviewer #1 (Comments to the Author):\nHu N. et al. conducted an interesting study on individual heterogeneity in response to environmental\nexposure, with a specific focus on PBDE-induced dyslipidemia. By integrating exposome, genomics,\nand metabolomic data from the China National Human Biomonitoring cohort, the authors identified\n3,571 genetic variants that interact with PBDE exposure to influence dyslipidemia risk with key\nmetabolites serving as mediators. Among these, rs9869609 was pinpointed as the top causal variant,\nexhibiting the strongest association with hypercholesterolemia risk. Further functional validation\nthrough single-base CRISPR/Cas9 editing revealed that the rs9869609-G allele downregulates\nSLC6A20 expression by enhancing BHLHE40 binding. This, in turn, impaired glycine transport\nand promoted cholesterol accumulation, particularly under BDE-47 exposure.\nThis study provides a comprehensive genome-wide analysis of gene-environment interactions in\nmetabolic diseases. It makes a significant contribution by identifying vulnerable populations and\nclarifying the mechanistic pathways underlying environmentally induced metabolic disorders.\nAdditionally, it opens new avenues for research on the health effects of environmental pollutants.\nHowever, several key limitations need to be addressed to improve the quality and impact of the\nmanuscript.\nResponse: Thanks for the positive comments. We have addressed the concerns in the following\nresponses.\nMajor comments:\n1. The authors constructed a polygenic risk score (PRS) using independent variants to predict lipid\nphenotypes, finding increased risk only in populations with high PBDE exposure. Why does the\nPRS not confer risk in low-exposure groups? This observation warrants further discussion\nResponse: Thanks for your comments. The PRS was associated with increased dyslipidemia risk\nexclusively in the high-PBDE exposure group, with no comparable effect under low exposure. This\ndemonstrates a clear G×E interaction, where genetic susceptibility becomes phenotypically\nsignificant only above a certain exposure threshold. Elevated PBDE levels appear to amplify the\nfunctional impact of risk variants, whereas low-level exposure remains insufficient to disrupt lipid\nmetabolism. These findings underscore those genetic effects are not uniform but depend on specific\nenvironmental exposure. We have elaborated the discussions on these mechanisms in the revised\nmanuscript.\nIn the Discussion section, page 14, line 453-460, we now state:\n“Notably, this risk was negligible in low-exposure groups, suggesting a threshold effect wherein\nlow PBDE levels are insufficient to disrupt metabolic pathways regulated by the interactive variants\nin the PRS models. In contrast, high PBDE exposure amplifies these disruptions through a\nsynergistic interaction with these variants, significantly increasing the risk of lipid metabolic\ndisorders. These findings demonstrate that genetic susceptibility to dyslipidemia is dynamically\nmodulated by PBDE exposure.”\n2. While the study identifies loci interacting with PBDE exposure, it does not clearly differentiate\nthis approach from traditional GWAS for disease-associated loci. Please elaborate on why gene-\nenvironment interaction (G×E) loci may offer greater biological or clinical relevance than main-\neffect GWAS hits, particularly for environmental metabolic disorders.\nResponse: Thanks for your insightful comment regarding the distinction between our G×E\ninteraction approach and traditional GWAS. G×E interaction loci offer greater biological and\nclinical relevance than main-effect GWAS loci. Specifically, they reveal how PBDEs and other\nenvironmental exposures modulate genetic effects, elucidate specific pathogenic pathways that are\nelusive to GWAS, identify genetically susceptible subgroups for targeted exposure-reduction\ninterventions, and link genetic susceptibility to modifiable environmental factors to provide\nactionable evidence for precision prevention.\nTo further demonstrate this advantage, we performed an overlap analysis between G×E interaction\ngenes and known GWAS loci in the GWAS Catalog (https://www.ebi.ac.uk/gwas/; EFO_0004530;\nEFO_0004612; EFO_0004611; EFO_0004574). Among 1,037 genes associated with PBDE-lipid\ninteractions phenotypes, 210 overlapped with known GWAS loci, while the remainder were novel.\nThese findings highlight exposure-specific biological pathways that traditional GWAS may\noverlook. We have added this discussion in the revised manuscript.\nIn the Discussion section, page 13, line 444-449, we now state:\n“G×E interaction analysis overcomes these shortcomings by revealing conditional genetic effects\ntied to specific environments and capturing more phenotypic variance, which our study also\nconfirms, of the 1037 PBDE-interactive genes, only 210 overlapped with known GWAS genes,\nshowing traditional GWAS overlook exposure-specific pathways.”\n3. In metabolomic analysis, the authors focus on overlapping metabolites for pathway enrichment\nand mediation analysis, potentially missing broader mechanistic insights. A supplementary\ncorrelation analysis of all significant metabolites could reveal additional candidate molecules and\npathways.\nResponse: Thanks for your constructive suggestions. Recognizing the limitation of focusing solely\non overlapping metabolites, we expanded our analysis to include all significant metabolites through\nKEGG pathway and mediation analyses. This approach revealed additional lipid homeostasis\npathways, including arachidonic acid metabolism and unsaturated fatty acid biosynthesis, and\nidentified 14 novel mediating metabolites. These results have been added to the revised manuscript.\nIn the Result section, page 9, line 270-276, we added:\n“We also performed pathway enrichment analysis on all differential metabolites. The results not\nonly validated the primary pathways identified through overlapping metabolite analysis but also\nrevealed novel pathways significantly enriched across all differential metabolites, including\narachidonic acid metabolism and biosynthesis of unsaturated fatty acids (Fig. S2c). This further\nbroadens our understanding of the G×E interaction metabolic network.”\nIn the Result section, page 10, line 297-302, we added:\n“Furthermore, we expanded the scope of our mediator analysis to include all G×E-associated\ndifferential metabolites as potential mediators. This supplementary analysis identified 14 mediating\nmetabolites (Fig. S2f). These metabolites primarily involve previously understudied pathways of\ninflammation and amino acid metabolism disorders, providing more comprehensive mechanistic\ninsights into G×E interactions affecting lipid metabolism.”\nIn the Discussion section, page 16, line 559-563, we added:\n“Additionally, we highlighted the involvement of supplementary pathways such as arachidonic acid\nmetabolism and biosynthesis of unsaturated fatty acid, thereby providing a more comprehensive\nperspective on the metabolic dysregulation induced by PBDEs.”\nFig. S2c, KEGG pathway enrichment analysis for PBDE, TC, LDL-C, HDL-C, and TG; bubble size\ncorresponds to the number of metabolites in the pathway. Dot color and size represent P-values and\ngene counts, respectively.\nFig. S2f, Sankey diagram of the mediating effect of PBDE on lipid characteristics, with lines\nrepresenting mediating effects.\n4. The choice of rs9869609 (over the most statistically significant SNPs) lacks clear justification.\nSumming functional scores across different ranges may introduce bias; a weighted or standardized\napproach would be more rigorous. Please clarify the prioritization criteria.\nResponse: Thank you for your insightful comments. The prioritization of rs9869609 was based on\nthree criteria: (1) its significant association with hypercholesterolemia in G×E analyses, showing\noptimal FDR values; (2) functional annotations indicating the highest regulatory potential across\nthree databases; and (3) its location in a regulatory region of SLC6A20 with significant eQTL effects.\nWe have added corresponding descriptions in the Methods and Results sections of the manuscript.\nIn the Method section, page 21, line 746-766, we now state:\n“We next performed a functional annotation for the genetic variants using multiple bioinformatic\ntools and databases, including the HaploReg database\n(https://pubs.broadinstitute.org/mammals/haploreg/haploreg.php), Roadmap Epigenomics\nConsortium (https://www.ncbi.nlm.nih.gov/geo/roadmap/epigenomics/), RegulomeDB database\n(https://www.regulomedb.org/regulome-search), CADD database\n(https://cadd.gs.washington.edu/score), 3DSNP (https://omic.tech/3dsnpv2/), which integrated\nmultiple histone modification chromatin immunoprecipitation (ChIP)-seq peaks, transcription\nfactor (TF) ChIP-seq peaks. To be specific, (1) regulatory chromatin histone ChIP-Seq were\ncharacterized by using the data from the Roadmap Epigenomics Consortium; (2) prioritizing these\nvariants with the scores of regulatory elements by using RegulomeDB, CADD, and 3DSNP\nprediction databases65. CADD first filters out potentially harmless variants, focusing on functionally\nsignificant sites66. RegulomeDB then determines whether these sites reside within regulatory\nregions and assesses their regulatory potential67. Finally, 3DSNP connects variants within regulatory\nregions to their actual distal target genes, completing the full biological narrative from variant to\ntarget gene68. Specifically, we apply minimum-maximum normalization to the raw outputs of these\ndatabases, uniformly converting their scores to the 0–1 range.”\nIn the Result section, page 10, line 313-317, we now state:\n“We calculated an integrated functional prediction score for these SNPs by combining outputs from\nCADD, RegulomeDB, and 3DSNP, which identified rs9869609 (r² = 0.991 with the tag SNP) as the\ncandidate causal variant due to its greatest potential for functional regulation (Fig. 5b, Table S6).”\nFig. 5b, Bar Plot of SNP Annotation Scores based on different genomic databases (Min-max\nstandardization of RegulomeDB, 3DSNP, and CADD scores).\n5. The functional validation uses the maximum non-toxic dose of BDE-47. How does this dose\ncompare to real-world exposure levels in the cohort? Are the observed effects plausible at\nenvironmentally relevant concentrations? Address potential limitations in extrapolating\nexperimental findings to population exposures.\nResponse: Thank you for the excellent suggestion. To address this, we established dose-response\nrelationships of BDE-47-induced lipid accumulation in both cellular models and human populations\n(Fig. R2A). The benchmark dose (BMD), a widely-used indicator to quantify the minimum\nconcentration for chemical toxicityR1, was 4.12 μM in HepG2 cells. Our in vitro concentration of\n12.5 μM (approximately three times the cellular BMD) based on its consistently induced phenotypic\nalterations and alignment with established toxicological paradigms for acute BDE-47 exposureR2.\nIn the human population (CNHBM cohort), the BMD for the association between BDE-47 exposure\nand dyslipidemia was 0.0711 μg/L (Fig. R2B). Notably, three times the population BMD (0.21 μg/L)\ncorresponds to the third quartile (Q3) of human BDE-47 exposure (Fig. R2C). This consistent dose-\nresponse pattern across experimental and human systems supports the translational relevance of our\nfindings.\nHowever, we acknowledge the limitations in extrapolating from acute in vitro exposure to chronic\nlow-dose human exposure, particularly regarding toxicokinetic differences, and have addressed this\nin the Discussion section.\nIn the Discussion section, page 17, line 575–577, we now state:\n“Third, the extrapolation from acute in vitro exposure to chronic low-dose human exposure remains\nchallenging due to differences in toxicokinetic and the complex interplay of multiple organ systems\nin vivo.”\nFig. R1 Dose-response curves for BDE-47 concentration versus the total cholesterols (TC) in vitro\n(A) and in CNHBM (B), with BDE-47 distribution in CNHBM (C). Note: BMDL: Lower\nbenchmark dose limit; BMD: Benchmark dose; Q1: First quartile; Q3: Third Quartile; M: Median.\n6. The enrichment of histone marks suggests epigenetic regulation, yet the study focuses solely on\nTF-mediated gene expression. Could other mechanisms (e.g., chromatin remodeling, splicing\nmodulation) contribute? Please expand the discussion on alternative mechanisms.\nResponse: Thank you for your insightful suggestions regarding alternative regulatory mechanisms\nbeyond TF-mediated gene expression. We agree that the observed molecular features indicate a\nbroader landscape of gene regulation. In response to your comment, we have expanded the\ndiscussions to include additional regulatory mechanisms.\nIn the Discussion section, page 14, line 466-497, we now state:\n“The intricate interplay between genetic susceptibility and environmental exposures collectively\nconstructs a sophisticated regulatory network that governs metabolic pathways11. Our investigation\nrevealed that the majority of interactive variants are located in non-coding regions enriched with\nregulatory elements such as promoters and enhancers. This indicates that these variants likely\nfunction by interfering with TF binding at promoter or enhancer sites, thereby modulating\ntranscriptional activity. … Additionally, we observed that these interactive variants are significantly\nenriched in epigenetic features. Emerging evidence suggests that environmental pollutants or their\nmetabolites may act as substrate mimics or inhibitors that directly bind to histone-modifying\nenzymes, thereby altering catalytic activity and leading changes in histone modification\npatterns11,17,42. Such epigenetic alterations can subsequently influence downstream RNA processing\nevents. For instance, histone modifications such as H3K36me3 can recruit splicing factors and\ndirectly influence the elongation rate of RNA polymerase II, thereby modulating alternative splicing\nevents43. These potential mechanisms provide clear directions for our future research.”\nMinor comments:\n1. SLC6A20 is implicated in glycine transport and cholesterol degradation via bile acid pathways.\nWere bile acid-related metabolites altered by PBDE exposure? Include relevant data or discuss their\npotential role.\nResponse: Thank you for this suggestion. To strengthen the mechanistic evidence for bile acid\nconjugation defects in SLC6A20 deficiency, we performed LC-MS analysis of bile acid metabolites\nin SLC6A20 knockdown HepG2 cells. The results demonstrated significantly reduced levels of\nglycocholic acid (GCA) and glycochenodeoxycholic acid (GDCA) in SLC6A20 knockdown HepG2\ncells, with a more pronounced decrease under BDE-47 exposure. These findings indicate that both\nSLC6A20 deficiency and BDE-47 exposure impair glycine transport, thereby reducing the synthesis\nof glycine-conjugated bile acids that are critical for cholesterol excretion. We have added a\ncorresponding description in the Result section of the manuscript.\nIn the Result section, page 12, line 408-414, we now state:\n“To evaluate the impact of SLC6A20 on glycine-conjugated bile acids, which are essential for\ncholesterol excretion, we measured GCA and GDCA levels in SLC6A20 knockdown HepG2 cells.\nLC-MS analysis revealed significantly reduced levels of GCA and GDCA concentrations in\nsiSLC6A20 cells compared to siControl cells. Following BDE-47 exposure, GCA and GDCA levels\nexhibited a further decrease in siSLC6A20 cells, while siControl cells showed a more moderate\nreduction (Fig. 6hi).”\nFig.6hi, GCA levels and GDCA levels in HepG2 cells with SLC6A20 knockdown, evaluated under\ncontrol conditions and BDE-47 exposure.\n2. HepG2 cells were used for validation—are these findings generalizable to primary hepatocytes\nor in vivo systems? Briefly justify the model choice and discuss translational limitations.\nResponse: We appreciate this important note. HepG2 cells were selected for initial validation due\nto their retained hepatocyte functions, though we acknowledge their limitations in metabolic\ncompetency compared to primary hepatocytes. Our findings warrant further validation in more\nphysiologically relevant models.\nIn the Discussion section, page 17, line 571-575, we now state:\n“Second, though used as a typical model for lipid metabolism, HepG2 cells exhibit reduced\nmetabolic competency, altered signaling pathways, and lack the full phenotypic complexity of\nprimary human hepatocytes or in vivo systems, which alert for future validations in other cell lines\nand in vivo models.”\n3. The term \"dyslipidemia\" is not defined, and the analysis uses quantitative lipid markers. How do\nthese findings translate to clinical dyslipidemia risk? Clarify the connection.\nResponse: Thank you for this insightful comment. To more precisely analyze the quantitative\nassociation between G×E interactions and lipid phenotypes, our models utilize quantitative lipid\nindicators. The term \"dyslipidemia\" is used in this study as a broad clinical concept referring to\nabnormal lipid levels, including elevated total cholesterol, elevated LDL-C, low HDL-C, or high\ntriglycerides.\nTo illustrate the association between these quantitative indices and disease dyslipidemia, we\nanalyzed the effect of variants rs98696609 on disease risk and found significant G×E interactions\nwith PBDE exposure on TC (FDR = 0.0082) and LDL-C (FDR = 0.0104). Under high exposure,\nGG genotype carriers showed substantially increased hypercholesterolemia prevalence (26.1%)\ncompared to AA carriers (16.8%). We have added a corresponding description in the Result section\nof the manuscript.\nIn the Result section, page 10, line 328-334, we added:\nWe further confirmed the population-level effect of this variant, in which significant interaction\neffects were observed for TC (FDR = 0.0082) and LDL-C (FDR = 0.0104) levels. Especially, under\nhigh PBDE exposure, individuals carrying the GG genotype exhibited substantially altered lipid\nprofiles compared to those with AA genotype, with hypercholesterolemia prevalence increasing\nfrom 16.8% to 26.1%, while the change was neglectable under low exposure (Fig.S3ab).\na, Cumulative distributions of lipid phenotypes (TC, LDL-C, HDL-C, TG) by rs9869609 genotype\nand exposure level. b, Stratified risk analysis of the rs9869609 in lipid phenotypes (TC, LDL-C,\nHDL-C, TG). Abbreviations: |βstd|, absolute standardized β.\n4. Why total cholesterol was chosen over other lipid traits for experimental validation? Were the\nloci for different lipids compared? Discuss their interrelationships.\nResponse: Thank you for your thoughtful comments. Our priority of TC for experimental validation\nis supported by phenotypic associations, G×E interaction analyses and biological relevance of TC.\nPhenotypically, TC showed the strongest positive associations with both ΣPBDEs and BDE-47\nexposure (Fig. S1d). Genetically, TC exhibited the broadest spectrum of G×E interaction variants\namong all lipid traits (Fig.2a). Through systematic screening of the top 20 loci prioritized by\nstatistical significance and interaction effect size (Table S5), the rs1392296 locus emerged as the\nmost compelling candidate, exhibiting the lowest FDR (0.0078) and highest |β | (0.6517). This\nstd\nrobust statistical evidence, combined with the biological relevance of its annotated gene SLC6A20\nin cholesterol metabolism, supported our focus on TC for subsequent mechanistic investigation.\nThese results have been incorporated into the revised manuscript.\nIn the Result section, page 10, line 305-313, we now state:\n“To investigate the molecular mechanisms by which PBDE-interacting variants disrupt lipid\nmetabolism, we first screened 20 top loci (the top five for each lipid phenotype) from Fig.2a and\nprioritized those with the most significant effects based on two core indicators: FDR and absolute\nstandardized β (interaction β per SD of outcome) (Table S5). As a result, the locus with tag SNP\nrs1392296 outperformed others, exhibiting the lowest FDR (0.0078) and highest |β | (0.6517). The\nstd\nlatter indicator suggests a relatively high impact of PBDE-SNP interaction on TC.”\nRegarding whether genetic variants for different lipids should be compared, we analyzed the overlap\nof G×E interaction variants across the four lipid phenotypes. Results showed TC and LDL-C shared\n193 interactive variants (Fig. R2A) and 80 overlapping annotated genes, mostly concentrated in\n8\ncholesterol metabolism pathways (Fig. R2B), while other lipid phenotypes had few common\nvariants. This confirms that validating TC’s mechanism not only clarifies its own response to PBDEs\nbut also provides a shared mechanistic reference for understanding LDL-C and other lipid\ndysregulation, balancing the study’s focus and efficiency. We have added a corresponding\ndiscussion in the revised manuscript.\nIn the Discussion section, page 16, line 534-537, we now state:\n“Notably, the interactive effect between rs9869609 and PBDEs on LDL-C, together with other\nshared variants influencing both TC and LDL-C, providing a mechanistic basis for the coordinated\ndysregulation of LDL-C and TC.”\nFig. R2 Venn diagram of genetic variants (A) and genes (B) associated with TC and LDL-C\n5. Provide basic genomic details of rs9869609 (e.g., distance to SLC6A20, intronic/intergenic status,\nMAF) to contextualize the variant’s potential functional impact.\nResponse: Thank you for this suggestion. As detailed in Table S6, rs9869609 is an intergenic\nvariant located approximately 347 kb downstream of SLC6A20. With a MAF of 6% and high\nintegrated functional scores from CADD, RegulomeDB, and 3DSNP, this variant was prioritized as\na candidate causal variant potentially functioning through distal cis-regulatory mechanisms.\nIn the Result section, page 10, line 313-328, we now state:\n“We calculated an integrated functional prediction score for these SNPs by combining outputs from\nCADD, RegulomeDB, and 3DSNP, which identified rs9869609 (r² = 0.991 with the tag SNP, MAF\n= 6%) as the candidate causal variant due to its greatest potential for functional regulation (Fig. 5b,\nTable S6). …This variant is located in an intergenic region, approximately 347 kb downstream of\nthe SLC6A20 gene, suggesting its distal cis-enhancer role to regulate the gene expression.”\n6. The PBDE-specific PRS lacks internal/external validation. Discuss its utility for future research\n(e.g., risk prediction, personalized interventions).\nResponse: Thank you for this suggestion. We acknowledge the challenges in accessing databases\nthat combine genome-wide genotype data with PBDE exposure measurements, which limited our\nability to validate findings in an independent cohort. This limitation has been explicitly noted in the\nDiscussion section. As suggested, we have also discussed the potential utility of the PRS for future\nresearch directions.\nIn the Discussion section, page 14, line 460-464, we now state:\n9\n“The PRS developed in this study validates polygenic G×E interactions, showing significant\nassociation with total cholesterol and confirming a concerted genetic mechanism. Although not yet\nclinically applicable, this PRS serves as a foundational tool for identifying high-risk individuals\nwho may benefit from targeted exposure mitigation.”\nIn the Discussion section, page 17, line 577-581, we now state:\n“Lastly, the findings’ generalizability may be constrained by the specific demographic and\nenvironmental context of the Chinese cohort studied. Validation across more diverse populations is\nessential to establish the broader applicability of these G×E interactions.”\n7. The experiments used BDE-47 as exposure, while analyses modeled total PBDEs. Explain this\ndiscrepancy and its implications for interpretation.\nResponse: Thank you for this comment. BDE-47 was selected for experimental validation based on\ntwo considerations: it showed the highest detection rate in our cohort and correlated strongly with\ntotal PBDEs levels, representing the primary exposure effect; and as a highly bioactive congener, it\nenables clearer mechanistic investigation in experimental systemsR3. The use of BDE-47 in\nexperiments provides a clear and representative entry point for mechanism studies, while employing\ntotal PBDEs in population analyses better reflects the cumulative effects of actual exposure.\nIn the Method section, page 22, line 780-787, we now state:\n“In the population effects analysis, to better reflect real-world exposure conditions, we used total\nPBDEs as a quantitative indicator of mixture exposure to assess its impact on blood lipids. Analysis\nresults indicate that the primary effects of BDE-47 in the population align with those of the total\nPBDEs mixture. In the experimental section, to streamline procedures and ensure representativeness\nof exposure levels, effect intensity, and mechanisms, we selected BDE-47 as the representative\ncompound for this class of pollutants to infer potential mechanisms of action for similar pollutants75.”\n8. Other typographical and formatting Issues:\nIn Fig. 2c, specify units for “Beta of blood lipids” in PRS models.\nIn Fig. 6, correct \"gpro\" to \"g pro.\" in metabolite/lipid labels and describe normalization methods.\nIn Fig. 4b, indicate the significance threshold for pathway enrichment.\nResponse: We appreciate this suggestion. We have carefully addressed all the typographical and\nformatting issues raised, including the addition of units in Fig. 2c and correction of metabolite labels\nin Fig. 6. Metabolite levels were normalized to total protein (mmol/g pro) using BCA assay. The\nsignificance threshold in Fig. 4b is specified as P < 0.05. All modifications have been incorporated\nadj\ninto the revised manuscript.\nThis study investigates the interaction between genetic variants and polybrominated diphenyl ether\n(PBDE) exposure in modulating lipid homeostasis, combining epidemiological data from the China\nNational Human Biomonitoring (CNHBM) cohort with multi-omics profiling and functional\nvalidation. The authors performed gene-environment (G×E) interaction analysis, identifying in\nSNPs associated with PBDE-induced dyslipidemia, with rs9869609 (near SLC6A20) emerging as\na top candidate (among others). Metabolomic analysis highlighted glycine and glycerophosphate as\nkey mediators, while CRISPR/Cas9-edited HepG2 cells demonstrated that the rs9869609-G allele\nsuppresses SLC6A20 expression via enhanced BHLHE40 binding, impairing glycine transport and\nexacerbating cholesterol accumulation under BDE-47 exposure.\nThe article is highly compelling due to its translational approach, bridging population-level findings\nwith mechanistic insights. The identification of SLC6A20 as a susceptibility marker for PBDE-\ninduced dyslipidemia and its validation through knockout experiments in HepG2 cells strengthens\nthe biological plausibility of the results. The study is well-structured, clearly written, and addresses\nan important gap in understanding how environmental pollutants interact with genetic background\nto influence metabolic health.\nWhile the findings are robust, some limitations should be noted. The functional validation was\nconducted only in vitro, and future studies using in vivo models (e.g., genetically modified mice or\nhumanized animal systems) would further solidify the conclusions. Additionally, the\ngeneralizability of the findings needs to confirm in other cohorts, and the clinical relevance of the\npolygenic risk score (PRS) warrants further investigation.\nDespite these limitations, the study represents a significant advance in environmental health\nresearch, providing novel insights into gene-environment interactions in dyslipidemia. The work is\nof broad interest to the fields of toxicology, metabolism, and precision medicine. With minor\nrevisions to address, this article is suitable for publication.\nResponse: Thanks for the positive comments. We have addressed the concerns in the following\nresponses.\nGeneral remarks\n1. Clarification of focus on SLC6A20\n1.1) The authors focus on SLC6A20 due to its strong statistical association with TC (rs1392296,\nFDR=0.0078) and compelling functional evidence linking its suppression to glycine-mediated\ncholesterol metabolism. However, other candidate genes (e.g., CTCF, enriched across all lipid traits)\nare not explored experimentally. While prioritization is inevitable, a brief discussion justifying why\nSLC6A20 was selected over equally plausible candidates—such as CTCF, a global epigenetic\nregulator including of lipid-related genes, or CETP, directly involved in lipoprotein remodeling—\nwould strengthen the manuscript. Future studies could address whether these genes act\nindependently or synergistically with SLC6A20 in PBDE-induced dyslipidemia.\nResponse: We appreciate the insightful comments and have supplemented additional information\nin response. First, our locus prioritization was based on G×E interaction significance (FDR) and\neffect size (Table S5). The rs1392296 locus, annotated to SLC6A20, exhibited the strongest\nassociation with TC (FDR = 0.0078, |β | = 0.6517). Subsequent fine-mapping revealed rs9869609\nstd\nas the top functional candidate based on CADD, RegulomeDB, and 3DSNP annotations among LD\nSNPs (Table S6). Additionally, we have supplemented functional annotations for other top loci in\nCAMK1D, with GTEx eQTL data supporting its regulatory role (Fig. R1). CAMK1D encodes a\nkinase involved in metabolic regulation R4. These additions have been incorporated into the results\nand discussion sections.\nIn the Result section, page 17, line 305-322, we now state:\n“To investigate the molecular mechanisms by which PBDE-interacting variants disrupt lipid\nmetabolism, we first screened 20 top loci (the top five for each lipid phenotype) from Fig.2a and\nprioritized those with the most significant effects based on two core indicators: FDR and absolute\nstandardized β (interaction β per SD of outcome) (Table S5). As a result, the locus with tag SNP\nrs1392296 outperformed others, exhibiting the lowest FDR (0.0078) and highest |β | (0.6517). The\nstd\nlatter indicator suggests a relatively high impact of PBDE-SNP interaction on TC. Within this locus,\nseven another significant SNPs displayed strong LD (r² ≥ 0.8) with the tag SNP (Fig. 5a). We\ncalculated an integrated functional prediction score for these SNPs by combining outputs from\nCADD, RegulomeDB, and 3DSNP, which identified rs9869609 (r² = 0.991 with the tag SNP, MAF\n= 6%) as the candidate causal SNP due to its greatest potential for functional regulation (Fig.\n5b).Genomic visualization further confirmed that rs9869609 is significantly enriched in DNase\nhypersensitive sites and active histone modification peaks (H3K4me1, H3K4me3, H3K9ac, and\nH3K27ac), underscoring its robust regulatory potential (Fig. 5c). Therefore, we prioritized\nfunctional validation of this variant to explore how SNPs interact with PBDE to drive lipid\ndysregulation.”\nIn the Result section, page 6, line 165-173, we now state:\n“Similarly, PBDE exposure interacted with 955 SNPs (177 loci), 769 SNPs (140 loci), and 739\nSNPs (127 loci) to influence LDL-C, HDL-C, and TG levels, respectively. The top SNPs from these\nanalyses were annotated to biologically relevant genes: CETP (FDR = 0.0089), which encodes\ncholesteryl ester transfer protein critical for lipoprotein remodeling; CPNE4 (FDR = 0.0104), a\ncalcium-dependent phospholipid-binding protein potentially involved in lipid metabolism; and\nCAMK1D (FDR = 0.0092), a calcium -dependent protein kinase implicated in hormonal regulation\nof metabolic processes.”\nIn the Discussion section, page 16, line 537-544, we now state:\n“Beyond SLC6A20, other genes contribute to PBDE-induced dyslipidemia, underscoring the\npolygenic nature of this process CETP mediates cholesterol ester transfer between lipoproteins,\nCTCF acts as an epigenetic regulator of lipid-related genes55, while CPNE4 and CAMK1D may\nparticipate in intracellular lipid signaling and metabolic regulation. These genes represent promising\ntargets for future investigation into the complex interplay between genetic susceptibility and\nenvironmental triggers in dyslipidemia.”\n1.2) The polygenic risk score incorporates multiple loci, yet only SLC6A20 is functionally tested.\nA brief explanation of why other top genes was excluded (e.g., lack of eQTL support, lower effect\nsizes) would bolster transparency.\nResponse: We appreciate your insightful comment. The PRS was constructed using established\nmethods and incorporated all genome-wide significant variants from our G×E interaction analysis.\nAmong these, SLC6A20 was prioritized for experimental validation based on statistical and\nfunctional evidence: it exhibited both the lowest FDR (0.0078) and the largest interaction effect size\n(|β | = 0.6517) for TC, and the regulatory variant rs9869609 within this locus demonstrated the\nstd\nhighest integrated functional score and strong eQTL.\nDue to practical experimental constraints, we focused our validation efforts on the most promising\ncandidate. To ensure transparency, we have provided complete functional annotations for all\ncandidate genes in Tables S5 and S6, including functional prediction scores, and literature-based\nbiological functions for each gene. We have also expanded discussion of other key genes such as\nCETP, CPNE4 and CAMK1D in the manuscript.\nIn the Discussion section, page 16, line 537-544, we now state:\n“Beyond SLC6A20, other genes contribute to PBDE-induced dyslipidemia, underscoring the\npolygenic nature of this process CETP mediates cholesterol ester transfer between lipoproteins54,\nCTCF acts as an epigenetic regulator of lipid-related genes55, while CPNE4 and CAMK1D may\nparticipate in intracellular lipid signaling and metabolic regulation. These genes represent promising\ntargets for future investigation into the complex interplay between genetic susceptibility and\nenvironmental triggers in dyslipidemia.”\n1.3) While the role of SLC6A20 in glycine transport is compelling, dyslipidemia is multifactorial.\nThe study does not rule out additive effects of other candidate genes. Please discuss\nResponse: Thanks for your insightful comment. We acknowledge the polygenic nature of\ndyslipidemia, which motivated the development of a PRS in this study. Our transcriptomic analysis\nin HepG2 cells exposed to BDE-47 identified 1,221 differentially expressed genes, including key\nlipid regulators such as CETP and CAMK1D alongside SLC6A20. These results indicate that PBDE-\ninduced dyslipidemia likely results from coordinated disruption across multiple genetic pathways,\nrather than the effect of a single gene.\nWe agree that investigating additive or interactive effects among these genes would be valuable.\ncompelling candidate. We have acknowledged this limitation and plan to pursue functional studies\nof additional genes in future work.\nIn the Discussion section, page 17, line 568-571, we now state:\n“First, our mechanistic experiments focused solely on the rs9869609 variant’s role in PBDE-induced\nhypercholesterolemia, future studies on other key genes, especially their independently or\nsynergistically act, are needed for comprehensive validation.”\n2. Association between PBDE, SLC6A20 and dyslipidemia\n2.1) The study lacks RNA-seq or single-cell transcriptomics to validate SLC6A20 expression\nchanges in human tissues or additional PBDE-responsive genes.\nResponse: Thank you for your insightful suggestion. In response, we performed RNA-seq in HepG2\ncells exposed to BDE-47 (12.5 μM, 24 h). The analysis identified 1221 significantly differentially\nexpressed genes, including SLC6A20, CETP and CAMK1D. KEGG pathway analysis revealed\nenrichment in cholesterol metabolism and lipoprotein biosynthesis pathways, indicating PBDE\nexposure disrupts lipid homeostasis through these metabolic pathways. Specifically, SLC6A20\nshowed significant downregulation (log₂FC = -2.21, FDR = 0.0032), supporting its involvement in\nresponse to PBDE exposure. These results have been included in the revised manuscript.\nIn the Result section, page 11, line 364-372, we now state:\n“To investigate the downstream responsive genes of BDE-47 exposure, we performed RNA-seq in\nHepG2 cells following BDE-47 exposure (12.5 μM, 24 h). As a result (Fig S4), we identified 1221\nsignificantly differentially expressed genes (DEG), including SLC6A20, CETP and CAMK1D\n(|log2(FC)| > 1, FDR < 0.05). KEGG analysis revealed that the DEG were primarily enriched in\ncholesterol metabolism and lipoprotein biosynthesis pathways, indicating that PBDE exposure\ndisrupts lipid homeostasis through these key metabolic pathways. Notably, SLC6A20 was among\nthe downregulated genes (log2(FC) = -2.21, FDR = 0.0032), pointing to its role in response to PBDE\nexposure.”\nFig. S4. Transcriptome Analysis following BDE-47 exposure in HepG2 cells.\na, Volcano plot of log2(FC) for Differentially expressed genes (DEG) versus false discovery rate\n(−log10(FDR)). b, KEGG pathway enrichment analysis for DEG; Dot color and size represent P-\nvalues and gene counts, respectively.\n2.2) The other top genes identified in genetic variants could had been investigated in HepaRG cells.\nwith comprehensive functional annotations, including prediction scores, and literature-supported\nbiological roles for all candidate genes. We have also expanded the discussion of CETP, CPNE4,\nand CAMK1D in the manuscript. However, due to technical and resource limitations, we focused\nexperimental validation on the most statistically compelling candidate, SLC6A20. We have\nexplicitly acknowledged this constraint in the revised manuscript.\nIn the Discussion section, page 16, line 537-544, we now state:\n“Beyond SLC6A20, other genes contribute to PBDE-induced dyslipidemia, underscoring the\npolygenic nature of this process CETP mediates cholesterol ester transfer between lipoproteins54,\nCTCF acts as an epigenetic regulator of lipid-related genes55, while CPNE4 and CAMK1D may\nparticipate in intracellular lipid signaling and metabolic regulation. These genes represent promising\ntargets for future investigation into the complex interplay between genetic susceptibility and\nenvironmental triggers in dyslipidemia.”\nIn the Discussion section, page 17, line 568-571, we now state:\n“First, our mechanistic experiments focused solely on the rs9869609 variant’s role in PBDE-induced\nhypercholesterolemia, future studies on other key genes, especially their independently or\nsynergistically act, are needed for comprehensive validation.”\n2.3) The claim that rs9869609-G is the \"top causal variant\" is premature without: In vivo validation\n(e.g., knock-in mice). Please temper this in the article.\nEvidence that this variant explains a significant proportion of dyslipidemia risk in the population.\nThe PRS improves prediction but has modest effect sizes (β=1.18 for TC).\nResponse: Thank you for your critical comment. We agree that claiming causality requires in vivo\nvalidationand have revised the terminology to \"candidate causal variant\" throughout the manuscript.\nRegarding the PRS effect size, the observed effect (β = 0.62 mmol/L, standardized β ≈ 0.59) meets\nCohen's criterionR5 for a large effect (standardized β > 0.5), confirming the substantial impact of\nG×E interactions on lipid levels. It should be noted that the single-variant effect (β = 1.18 mmol/L)\nand PRS represent distinct aspects of genetic susceptibility. The former reflects a strong interaction\nunder specific environmental conditions, while the latter captures cumulative risk across multiple\nloci. Both approaches provide complementary insights into the genetic architecture of PBDE-\ninduced dyslipidemia.\n3. Here is some piece of advice that the authors could take into consideration and discuss: Clarify\nExposure Relevance: It could be interesting to compare in vitro BDE-47 doses to human urinary\nlevels and assess dose-response relationships.\nResponse: Thank you for the excellent suggestion. To address this, we established dose-response\nrelationships of BDE-47-induced lipid accumulation in both cellular models and human populations\n(Fig. R2A). The benchmark dose (BMD), a widely-used indicator to quantify the minimum\nconcentration for chemical toxicityR1, was 4.12 μM in HepG2 cells. Our in vitro concentration of\n12.5 μM (approximately three times the cellular BMD) based on its consistently induced phenotypic\nalterations and alignment with established toxicological paradigms for acute BDE-47 exposureR2.\nand dyslipidemia was 0.0711 μg/L (Fig. R2B). Notably, three times the population BMD (0.21 μg/L)\ncorresponds to the third quartile (Q3) of human BDE-47 exposure (Fig. R2C). This consistent dose-\nresponse pattern across experimental and human systems supports the translational relevance of our\nfindings.\nHowever, we acknowledge the limitations in extrapolating from acute in vitro exposure to chronic\nlow-dose human exposure, particularly regarding toxicokinetic differences, and have addressed this\nin the Discussion section.\nIn the Discussion section, page 17, line 575-577, we now state:\n“Third, the extrapolation from acute in vitro exposure to chronic low-dose human exposure remains\nchallenging due to differences in toxicokinetic and the complex interplay of multiple organ systems\nin vivo.”\nFig. R2 Dose-response curves for BDE-47 concentration versus the total cholesterols (TC) in vitro\n(A) and in CNHBM (B), with BDE-47 distribution in CNHBM (C). Note: BMDL: Lower\nbenchmark dose limit; BMD: Benchmark dose; Q1: First quartile; Q3: Third Quartile; M: Median.\n4. Address Population Bias: Replicate findings in diverse cohorts or perform trans-ethnic meta-\nanalyses.\nResponse: We appreciate the suggestion regarding trans-ethnic validation. To clarify, the CNHBM\ncohort provides substantial diversity, encompassing participants from 31 provinces across China\nand representing multiple ethnic groups, as detailed in the Methods section. This geographical and\nethnic heterogeneity strengthens the reliability of our findings within the Chinese population.\nWe acknowledge that external validation in non-East Asian cohorts would further reinforce the\ntrans-ethnic relevance of our results. However, the limited availability of large-scale cohorts with\nintegrated genotype data and detailed PBDE exposure measurements has constrained our ability to\nperform trans-ethnic meta-analyses. We have explicitly recognized this as a limitation in the revised\nmanuscript.\nIn the Discussion section, page 18, line 607-611, we now state:\n“Notably, the CNHBM cohort comprises a multinational population, including Han (18,794, 86.4%),\nTibet (551, 2.5%), Zhuang (392, 1.8%), Mongolia (340, 1.6%), Hui (279, 1.3%), Man (249, 1.1%),\nLi (180, 0.8%), Uyghur (138, 0.6%), Dong (138, 0.6%), and other ethnicities (685, 3.2%), reflecting\n16\nits broad ethnic and geographic representation.”\nIn the Discussion section, page 17, line 577-581, we now state:\n“Lastly, the findings’ generalizability may be constrained by the specific demographic and\nenvironmental context of the Chinese cohort studied. Validation across more diverse populations is\nessential to establish the broader applicability of these G×E interactions.”\n5. Strengthen Mechanistic Claims:\nProvide direct evidence of bile acid conjugation defects in SLC6A20-deficient models.\nResponse: Thank you for this suggestion. To strengthen the mechanistic evidence for bile acid\nconjugation defects in SLC6A20 deficiency, we performed LC-MS analysis of bile acid metabolites\nin SLC6A20 knockdown HepG2 cells. The results demonstrated significantly reduced levels of\nglycocholic acid (GCA) and glycochenodeoxycholic acid (GDCA) in SLC6A20 knockdown HepG2\ncells, with a more pronounced decrease under BDE-47 exposure. These findings indicate that both\nSLC6A20 deficiency and BDE-47 exposure impair glycine transport, thereby reducing the synthesis\nof glycine-conjugated bile acids that are critical for cholesterol excretion. We have added a\ncorresponding description in the Result section of the manuscript.\nIn the Result section, page 12, line 408-414, we now state:\n“To evaluate the impact of SLC6A20 on glycine-conjugated bile acids, which are essential for\ncholesterol excretion, we measured GCA and GDCA levels in SLC6A20 knockdown HepG2 cells.\nLC-MS analysis revealed significantly reduced levels of GCA and GDCA concentrations in\nsiSLC6A20 cells compared to siControl cells. Following BDE-47 exposure, GCA and GDCA levels\nexhibited a further decrease in siSLC6A20 cells, while siControl cells showed a more moderate\nreduction (Fig. 6hi).”\nFig.6, Boxplots showing GCA (h) and GDCA (i) levels in SLC6A20 knockdown and control HepG2\ncells under non-exposure or BDE-47 exposure.\n6. Integrate transcriptomics to identify downstream effects of PBDE-exposed genotypes.\nResponse: Thank you for your insightful suggestion. In response, we performed RNA-seq in HepG2\ncells exposed to BDE-47 (12.5 μM, 24 h). The analysis identified 1221 significantly differentially\nexpressed genes, including SLC6A20, CETP and CAMK1D. KEGG pathway analysis revealed\nenrichment in cholesterol metabolism and lipoprotein biosynthesis pathways, indicating PBDE\nexposure disrupts lipid homeostasis through these metabolic pathways. Specifically, SLC6A20\nshowed significant downregulation (log₂FC = -2.21, FDR = 0.0032), supporting its involvement in\nresponse to PBDE exposure.\nIn the Result section, page 11, line 364-372, we now state:\n“To investigate the downstream responsive genes of BDE-47 exposure, we performed RNA-seq in\nHepG2 cells following BDE-47 exposure (12.5 μM, 24 h). As a result (Fig S4), we identified 1221\nsignificantly differentially expressed genes (DEG), including SLC6A20, CETP and CAMK1D\n(|log2(FC)| > 1, FDR < 0.05). KEGG analysis revealed that the DEG were primarily enriched in\ncholesterol metabolism and lipoprotein biosynthesis pathways, indicating that PBDE exposure\ndisrupts lipid homeostasis through these key metabolic pathways. Notably, SLC6A20 was among\nthe downregulated genes (log2(FC) = -2.21, FDR = 0.0032), pointing to its role in response to PBDE\nexposure.”\nFig. S4. Transcriptome Analysis following BDE-47 exposure in HepG2 cells.\na, Volcano plot of log2(FC) for Differentially expressed genes (DEG) versus false discovery rate\n(−log10(FDR)). b, KEGG pathway enrichment analysis for DEG; Dot color and size represent P-\nvalues and gene counts, respectively.\nMinor points\n1. Figure 2 Color Consistency Issue:\nIn panel (a), colors distinguish chromosomes, but in panels (b) (stratified risk) and (c) (PRS trends),\nthe same/similar colors represent PBDE exposure groups (low/high). This is confusing. Use distinct\ncolor schemes (e.g., grayscale for chromosomes, blue/orange for exposure groups) or add explicit\nlabels to legends.\nResponse: Thanks for your suggestion. We have revised Fig. 2a by using grayscale to represent\nchromosomes.\n2. Figure 4 (Metabolome Profiling)\nThere is a legend mismatch in Panel (b3) (HDL-C):\nThe legend describes two circle sizes (likely representing effect size or significance), but three sizes\nappear in the plot (small/medium/large). Please align the legend with the visual output.\nOn panel c line widths effectively represent mediation strength, but key metabolites (e.g., glycine,\ncitrate) could be highlighted with bold labels or arrows to emphasize central pathways.\nResponse: Thanks for your reminding. We have standardized the legends by categorizing circle\nsizes into three distinct levels, corresponding to pathway enrichment ratios (low: 5; medium: 10;\nhigh: 15). We have added bold labels to key metabolites, including glycine, citrate, and glycerol 3-\nphosphate, enhancing their visibility. These revisions have been incorporated into the updated Fig.\n4 in the revised manuscript.\n3. Typos and suggestions\nLine 164: unit instead of units\nLine 168: Suggest rather than demonstrate\nLine 1026 Characterized\nLine 1085: H3K4me1 repeated twice\nResponse: Thanks for your reminding. We have corrected the inappropriate description in the\nrevised manuscript.\nReferences\nR1. Jin, Y., et al. High throughput data-based, toxicity pathway-oriented development of a\nquantitative adverse outcome pathway network linking AHR activation to lung damages.\nJournal of hazardous materials 425, 128041 (2022).\nR2. Liu, X., et al. The role of lysosomes in BDE 47-mediated activation of mitochondrial\napoptotic pathway in HepG2 cells. Chemosphere 124, 10-21 (2015).\nR3. Zhu, L., Ma, B. & Hites, R.A. Brominated flame retardants in serum from the general\npopulation in northern China. Environmental science & technology 43, 6963-6968 (2009).\nR4. Vivot, K., et al. CaMK1D signalling in AgRP neurons promotes ghrelin-mediated food\nintake. Nat Metab 5, 1045-1058 (2023).\nR5. Zou, G.Y. Exact confidence interval for Cohen's effect size is readily available. Statistics\nin medicine 26, 3054-3056 (2007).\nReviewer #1 (Comments to the Author):\nThe authors have adequately addressed my concerns, and I have no further comments.\nResponse to Reviewer #1:\nWe thank the reviewer for the positive feedback.\nReviewer #2 (Remarks to the Author):\nThe autjors did an exceptionally work in revising the mansucript. They have not only addressed all\nthe points that I (and R1) raised, but have done it in a manner that has significantly strengthened the\nmanuscript that is now more robust and impactful. I have no mre reserves and recommend the\nmanuscript for acceptance in its present form.\nThe new data and analyses the authors incorporated are impressive. Regarding the prioritization of\nSLC6A20, they provided a detailed clarification of the multi-tiered prioritization strategy. The\nexplanation, based on the strongest G×E interaction signal (rs1392296 locus), followed by fine-\nmapping to identify rs9869609 as the variant with the highest integrated functional potential, is now\ncompletely clear and convincing. The inclusion of functional annotations and biological discussions\nfor other top candidate genes (CETP, CPNE4, CAMK1D) in the results and discussion sections\nprovides a more balanced view. This effectively addresses my initial concern.\nI appreciated the careful revision of the language throughout the manuscript to refer to rs9869609\nas a \"candidate causal variant\".\nRegarding the relevance of in vitro doses. the addition of the benchmark dose (BMD) analysis for\nboth the cellular model and the human cohort is a major strength of the revised manuscript. The\nauthors now show that experimental concentration is multiple of the cellular BMD, and that a\ncomparable multiple of the human BMD falls within a relevant exposure range (the third quartile)\nof the population. This bridges the gap between the experimental system and human environmental\nexposure. This greatly bolsters the translational significance of findings made by the authors.\nThe new LC-MS data quantifying bile acid levels in the knockdown model is an insightful piece of\nevidence. The results elegantly close the mechanistic loop, linking SLC6A20 knockdown to\nimpaired glycine transport to reduced glycine-conjugated bile acids to compromised cholesterol\nexcretion. This is a major improvement to the study's mechanistic narrative.\nThe new RNA-seq data from BDE-47 exposed HepG2 cells is highly valuable. It provides a\nsystems-level view, identifying a network of dysregulated genes (CETP, CAMK1D, etc.) and\nenriched pathways. This successfully addresses my comment and highlights the polygenic nature of\nthe response to PBDE exposure.\nAuthors completely meet the expectations for a revision and I congratulate them for this excellent\nwork and have no further reservation.\nResponse to Reviewer #2:\nanalyses (BMD, LC-MS, RNA-seq) and clarifications have strengthened the manuscript’s\nmechanistic narrative and translational relevance to your satisfaction.\n2"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have adequately addressed my concerns, and I have no further comments.", "reviewer_2": "The autjors did an exceptionally work in revising the mansucript. They have not only addressed all the points that I (and R1 ) raised, but have done it in a manner that has significantly strengthened the manuscript that is now more robust and impactful. I have no mre reserves and recommend the manuscript for acceptance in its present form.\nThe new data and analyses the authors incorporated are impressive. Regarding the prioritization of SLC6A20, they provided a detailed clarification of the multi-tiered prioritization strategy. The explanation, based on the strongest G×E interaction signal (rs1392296 locus), followed by fine-mapping to identify rs9869609 as the variant with the highest integrated functional potential, is now completely clear and convincing. The inclusion of functional annotations and biological discussions for other top candidate genes (CETP, CPNE4, CAMK1D) in the results and discussion sections provides a more balanced view. This effectively addresses my initial concern.\nI appreciated the careful revision of the language throughout the manuscript to refer to rs9869609 as a \"candidate causal variant\".\nRegarding the relevance of in vitro doses. the addition of the benchmark dose (BMD) analysis for both the cellular model and the human cohort is a major strength of the revised manuscript. The authors now show that experimental concentration is multiple of the cellular BMD, and that a comparable multiple of the human BMD falls within a relevant exposure range (the third quartile) of the population. This bridges the gap between the experimental system and human environmental exposure.\nThis greatly bolsters the translational significance of findings made by the authors.\nThe new LC-MS data quantifying bile acid levels in the knockdown model is an insightful piece of evidence. The results elegantly close the mechanistic loop, linking SLC6A20 knockdown to impaired glycine transport to reduced glycine-conjugated bile acids to compromised cholesterol excretion. This is a major improvement to the study's mechanistic narrative.\nThe new RNA-seq data from BDE-47 exposed HepG2 cells is highly valuable. It provides a systems-level view, identifying a network of dysregulated genes (CETP, CAMK1D, etc.) and enriched pathways. This successfully addresses my comment and highlights the polygenic nature of the response to PBDE exposure.\nAuthors completely meet the expectations for a revision and I congratulate them for this excellent work and have no further reservation."}, "author_response": "Response to Reviewer #1:\nWe thank the reviewer for the positive feedback.\n\nResponse to Reviewer #2:\nWe sincerely appreciate your detailed and constructive assessment. We are glad that the additional analyses (BMD, LC-MS, RNA-seq) and clarifications have strengthened the manuscript’s mechanistic narrative and translational relevance to your satisfaction."}]

Decision letter