In-Vitro and In-Silico Assessment of Per- and Polyfluoroalkyl Substances (PFAS) in Aqueous Film-Forming Foam (AFFF) Binding to Human Serum Albumin

Drinking water contaminated by fluorosurfactant-based aqueous film-forming foams (AFFF) is a source of human exposure to poly- and perfluoroalkyl substances (PFAS). However, assessment of bioaccumulation potentials of diverse PFAS in commercial products such as AFFF have been insufficient and challenging, especially due to a lack of analytical standards. Here we explore the value of suspect screening, equilibrium dialysis, and molecular-docking simulations to identify potentially bioaccumulative PFAS. We exposed human serum albumin (HSA) protein to dilutions of a legacy AFFF produced by 3M in 1999 using equilibrium dialysis and screened in-vitro protein-binding affinities using high-resolution mass spectrometry (HRMS). Through suspect screening, we identified 32 PFAS and 18 hydrocarbon surfactants in the AFFF that bound to HSA. Quantification of noncovalent association constants for 26 PFAS standards confirmed that many PFAS, including the short-chain perfluoropropane sulfonic acid (log Ka= 4.1 ± 0.2 M−1), exhibit strong binding affinities with HSA. At least five PFAS in AFFF (including three PFAS with less than five perfluorocarbons) remained bound to the precipitated HSA pellet after extensive solvent washing—an indication of high PFAS binding potential. Three PFAS (PFBS, PFOS, and PFOA) were confirmed in the protein pellet with analytical standards and quantified after acid digestion—this sample fraction accounted for 5 to 20% of each compound mass in the sample. We calculated pseudo-bioconcentration factors (BCFpseudo) for PFAS that suspect screening flagged as noncovalently bound or potentially covalently bound. Most PFAS exhibiting high BCFpseudo, especially those with seven perfluorocarbons, contained a carboxylic acid or a sulfonic acid. Finally, we used molecular docking to simulate HSA binding affinities for 62 ligands (26 PFAS targets, 18 PFAS qualified in AFFF, and 18 hydrocarbon surfactants qualified in AFFF). We found that molecular docking can effectively separate HSA-binding and -nonbinding compounds in AFFF. In-vitro and in-silico approaches described in this study provide replicable, high-throughput workflows for assessing bioaccumulation potentials of diverse PFAS in commercial products.


Introduction
Application of aqueous film-forming foams (AFFF) for fire-suppression at military bases and airports is a cause of drinking water contamination with poly-and perfluoroalkyl substances (PFAS) [1]. Fluorosurfactant-based legacy AFFF formulations include complex mixtures of perfluoroalkyl carboxylic acids (PFCA), perfluoroalkyl sulfonic acids (PFSA), and highly diverse PFAS, including polyfluorinated precursors [2]. Consumption of AFFFcontaminated drinking water can lead to elevated PFAS levels in human blood [3,4]. Occupational exposure to PFAS in AFFF may also present health risks to firefighters [5][6][7]. Human exposure to PFAS has been linked to cancer, cardiovascular disease, kidney disease, liver disease, immune suppression, neurological disease, type II diabetes, osteoarthritis, respiratory disease, among other impacts [8,9]. Given these problems, researchers have gained interest in studying the health impacts of novel PFAS in AFFF [10,11], including compounds with one perfluorinated carbon that sometimes are not classified as PFAS (e.g., fluorinated aromatics) [12]. In particular, constructing physiologically-based pharmacokinetic (PBPK) models for PFAS exposure requires researchers to determine affinities (quantified as partition coefficients or association constants) among different PFAS structures in different biological tissues. The high numbers and structural diversity of existing and emerging PFAS renders this task experimentally infeasible. An alternative approach is to evaluate associations of PFAS mixtures in an AFFF with abundant model proteins (commonly serum albumin and liver fatty acid binding proteins) to identify potentially bioaccumulative PFAS and to yield quantitative relationships between PFAS exposure, bioaccumulation, and tissue distribution [13]. In this study, we assess binding affinities of diverse PFAS in AFFF to human serum albumin (HSA). HSA is the most abundant protein in human blood plasma, presenting in tissues throughout the body, and serves important biological functions (e.g., transportation of fatty acids, drugs, and thyroid hormones) [14]. Data on the binding affinities of AFFF-derived PFAS with HSA will support development of PBPK models for such PFAS.
Accurately incorporating protein binding affinities into PBPK models requires accurate understanding and quantification of the molecular mechanisms at play. Several studies reported that PFCAs and PFSAs bind with serum albumin proteins noncovalently through specific site binding or non-specific surface adsorption [15][16][17][18]. Two studies identified a potential for covalent binding between PFAS and albumin proteins [19,20]. While no studies have investigated ultrastrong noncovalent bindings (K A > 10 9 M −1 ) between PFAS and proteins, ultrastrong binding was observed for PFAS in aqueous supramolecular polymerization [21]. Common serum-extraction protocols are likely to overlook or discard strongly bound ligands, including covalently bound or ultrastrong noncovalently bound PFAS. In organic solvent extraction, for instance, the precipitated protein pellet is disposed after extraction-along with any strongly bound ligands or residual targets. In online or offline solid phase extraction (SPE), covalently bound ligands and denatured proteins are lost on SPE cartridges [22]. Additionally, many studies accounted for matrix effects by spiking calibration standards into blank serum before analysis [23]. Extraction efficiencies are typically reported as satisfactory (e.g., 70-130% spike recoveries) by using such matrix (serum) matched calibration curves for quantification, thereby inadvertently masking strong protein-ligand interactions [24,25].
In this study, we combined experimental and modeling techniques to identify potentially bioaccumulative PFAS present in an AFFF and to investigate multiple binding pathways for diverse PFAS structures. We first utilized HSA as a model protein system to quantify both noncovalently bound and potentially covalently bound PFAS using targeted analysis. This analytical approach facilitates evaluation of the degree to which strongly bound or residual PFAS may be discarded in precipitated protein pellets. We then estimated bioconcentration factors for legacy and novel PFAS using linear (L)-PFOS as a bioaccumulative benchmarking compound. Finally, we predicted protein binding affinities for novel PFAS using molecular docking, a traditional drug design tool that simulates interactions between small molecules and large proteins [26,27]. The results comprised the most comprehensive quantification of relative PFAS-HSA binding affinities to date, providing valuable inputs for bioaccumulation and PBPK models.

PFAS Terminology
In this paper, we adopted acronyms and rules established in previous studies for PFAS terminology. We mainly followed acronym naming rules established by Barzen-Hanson et al. (2017) [2]. We followed the rules published by Buck et al. (2011) for naming polyfluoroalkyl phosphate esters (PAPs) and fluorotelomer acrylates [28]. Table S6.2 details our abbreviation rules for compounds that were not previously characterized with abbreviated names. The number of perfluorocarbons (n) in a compound is indicated as C n , which should not be confused with the number of carbons in PFAS. For instance, PFOA is a C 7 PFCA, while PFOS is a C 8 PFSA. Based on categorizations from the Organisation for Economic Co-operation and Development (OECD) [12], we used the term "short-chain" to refer to PFCAs containing fewer than six perfluorocarbons and PFSAs with fewer than five perfluorocarbons. We used the term "long-chain" to refer to PFCAs containing six or more perfluorocarbons and to PFSAs with five or more perfluorocarbons. We applied the short-chain and long-chain convention for PFCAs to all other PFAS evaluated in this study. For analogs of well-known PFSAs identified, we added common contractions to enhance recognizability. For instance, PFOS analogs include chloro (Cl), ketone (K), unsaturated (U), hydrogen (H), linear (L), and branched (br)-substituted PFOS [29].

Study Design and Workflow
The overall workflow consists of the following components. First, we exposed HSA to either (a) dilutions of an AFFF produced by 3M in 1999 or (b) in-house mixtures of 26 PFAS (listed in Table S6.2) through equilibrium dialysis. We used data-dependent (tMS/MS, for target PFAS) and data-independent (all-ion fragmentation, for suspect-screening) mass spectrometry against an in-house library of PFAS and hydrocarbon surfactants to identify compounds that were bound noncovalently to HSA. Second, we acid-digested residual protein pellets that were free of noncovalently bound PFAS. We applied suspect screening to identify residual PFAS in the precipitated protein pellet. Residual PFAS were considered candidates for forming ultrastrong or covalent associations with HSA. We used targeted MS for three PFAS to quantify residual levels in the protein pellet. Third, we quantitatively evaluated protein association constants predicted by molecular docking between 26 target PFAS structures and two HSA crystal structures. Fourth and finally, we used the molecular docking workflow to classify PFAS and non-PFAS surfactants in the AFFF as either HSAbinding or non-binding.

Equilibrium Dialysis
Equilibrium dialysis was performed in a 96-well system (Harvard Apparatus, Holliston, MA, USA) in which each polypropylene cell was separated into two chambers by a 10-kDa regenerated cellulose membrane. One side of each dialysis cell was dosed with 7.97 mg HSA (≤0.02% fatty acids, Sigma-Aldrich, Munich, Germany) in phosphatebuffered saline (PBS, pH 7.4 prepared in HPLC grade water) to a final concentration of 600 µM, which mimics physiological conditions [30]. The other side of each cell was then dosed with 200 µL of an AFFF dilution (4000 to 16,000-fold in PBS) or an in-house mixture of 26 PFAS (prepared with PFAC-24PAR, PFPrS, and br-PFOS from Wellington Laboratories Inc., Guelph, ON, Canada). Seventeen out of the 26 target PFAS in this study are commonly measured in drinking water using EPA method 533 or EPA method 537.1.
An aliquot of legacy AFFF (3M, 1999) was provided by Professor Christopher Higgins at the Colorado School of Mines. Experimental batches consisted of six concentration levels of AFFF or PFAS standard dilutions and were replicated four times. For each batch of experiments, a method blank was prepared with HSA free of AFFF or PFAS. A negative control cell containing AFFF or the PFAS standard mix was prepared without HSA to assess free movement of PFAS through the membrane ( Figure S1.1). The system was incubated at 37 • C while rotating at 30 RPM for 108 h ( Figure S1.2).

PFAS Extractions
Aliquots from each dialysis cell were processed to generate three extracts, shown in Figure 1. These were (1) free PFAS from the aqueous fraction, (2) noncovalently bound PFAS associated with the dissolved protein, and (3) residual PFAS in the precipitated protein pellet (candidates for ultrastrong noncovalently bound or potentially covalently bound PFAS). Briefly, post-dialysis aqueous samples (100 µL) from the chemical chambers were equilibrated with 50% methanol. Protein aliquots (100 µL) of the post-dialysis protein cells were extracted with formic acid (FA) acidified acetonitrile (ACN) for protein denaturation and precipitation. Noncovalently bound PFAS were determined by taking the concentration Toxics 2021, 9, 63 4 of 18 difference between protein aliquot fraction and chemical fraction. The residual protein pellets were further washed with 1 mL ACN five times, and the last wash was concentrated (to 200 µL) and saved to verify the absence of PFAS. We then applied a standard acid hydrolysis protocol for amino acid analysis to break peptide bonds and to release any PFAS that were possibly covalently bound to HSA [31][32][33][34][35] Extracts of noncovalently bound PFAS and residual PFAS in the protein pellet were solvent exchanged into 50% methanol to match the final solvent composition of the aqueous extracts. Additional details on the sample preparation can be found in Supporting Information (S1.1). An internal standard mix (ISTD , Table S6.1) in 50% methanol was added into each extract prior to analysis. Details on instrumentation and acquisition settings are provided in detected consistently in the protein pellets (at least 16 out of 18 pellets among 3 trials) after AFFF exposure, but could not be quantified due to lack of available standards.

Analytical Instrumental Set-Up
Quantification of 26 PFAS targets and qualification of diverse PFAS via suspect screening were performed on each AFFF dilution and each PFAS extract from dialysis. Data were acquired on an Agilent 1260 Infinity HPLC system paired with a 6530 QTOF MS. Sample extracts (10 µL) were injected onto a C18 column (ZORBAX RRHD Eclipse Plus C18 column; 2.1 mm × 150 mm, 1.8 µm, Agilent Technologies, Inc., Santa Clara, CA, USA) at a flow rate of 0.40 mL/min, with a total run time of 31.5 min. The aqueous mobile phase (A) was 20 mM ammonium acetate (Fisher Scientific, Pittsburgh, PA, USA) in Optima TM HPLC grade water (Fisher Scientific, Pittsburgh, PA, USA) and the organic mobile phase (B) was 100% Optima TM HPLC grade acetonitrile (Fisher Scientific, Pittsburgh, PA, USA). The mass spectrometer ionized samples in a negative mode using collision energies (CE) of 0, 10, 20, and 40 eV. A quality-control run of the 26-PFAS standard mix was analyzed after every 8 samples to ensure that concentrations of targeted PFAS remained within 30% of known concentrations.

Suspect Screening
For suspect screening, mass-to-charge-ratios (m/z) of 50-1200 were fragmented in the collision cell with CE of 0, 10, 20, and 40 eV in All-Ions acquisition modes. Data were processed using Agilent MassHunter Qualitative Analysis (B.08.00) by applying the "Find by Formula" search against an in-house AFFF Personal Compound Data Library (PCDL). The PCDL contained 3793 PFAS extracted from the Norman Suspect List Exchange (OECDP-FAS) [36], and 727 hydrocarbon surfactants (monoisotopic mass: 100-1200) extracted from the surfactant suspect list curated by   [37]. The PCDL also included 63 MS/MS spectra, of which 31 spectra were acquired from in-house standards, 24 spectra were extracted from MassBank [38], and six spectra were generated with CFM-ID 3.0 [39]. PFAS were considered qualified with level 2-3 confidence as outlined in   [40]. Suspect-screening search settings are listed in Table S5.2. Qualified PFAS were reported only if the following additional criteria were met: (1) the abundance of the qualified ion was greater than three times the experimental blank ion abundance; (2) the ion was qualified in at least 67% of dialysis cells in each batch of experiments; and (3) the PFAS qualified in dialysis cell extracts were also qualified in neat AFFF dilutions. All qualified PFAS that met these additional criteria were confirmed by re-running extracts with data-dependent targeted analysis (see Section 2.6).

Targeted MS/MS
For data-dependent targeted analysis, a list of targeted mass-to-charge ratios (m/z) and corresponding retention times (RT) was compiled based on pre-runs with All-Ions acquisition described in Section 2.5. To avoid overlapping peaks of targeted compounds, duplicate injections were performed for each sample to ensure at least 0.4 min RT difference between peaks of any two targeted compounds within one injection. This approach can minimize false identifications of PFAS due to the instrument's inherent mass error.

PFAS Quantification
Quantification was performed with 19-ISTD dosed ten-point calibration curve (0.5-250 ng/mL). Whole method limits of quantification (LOQ) range from 0.025 to 1 ng/mL. Details of analytical standards and extraction recoveries are available in Table S6.1.

Experimental Determination of PFAS Noncovalent Binding Affinities
The concentrations of 26 PFAS targets that partitioned into the protein chamber, remained in the chemical chamber, and remained associated with the protein pellet were directly determined by HPLC-QTOF-MS. Noncovalent binding affinities, measured as association constants (K A ), were calculated assuming a one site specific binding as shown in Equations (1) and (2).
To enable the assessment of multiple specific binding sites, we also fitted the data to the Langmuir isotherm model [41] with a limited binding sites assumption, following Equations (3) and (4).
In these equations, i refers to the compound of interest, q m is the concentration of total binding sites, and q 0,i is the concentration of empty binding sites. [PFAS i ] is the concentration of free PFAS i measured in the chemical side of the equilibrium dialysis set-up. [HSA − PFAS i ] is calculated by taking the difference between the concentration of PFAS in the protein side and in the chemical side, as shown in Figure 1. Additional Toxics 2021, 9, 63 6 of 18 isotherm models, including linear adsorption and Freundlich adsorption models, were also evaluated ( Figure S1.3).

Computational Simulations of Noncovalent PFAS Protein Binding
We used AutoDock Vina (v 1.1.2) [42] to dock 62 ligands (26 PFAS targets, 18 qualified PFAS, and 18 qualified hydrocarbon surfactants in AFFF) to two HSA crystal structures (Protein Data Bank entries 1E7G and 1AO6). 1E7G was chosen as the native structure of HSA, which complexed with tetradecanoic acid (myristic acid) [43,44]. 1AO6 is an unliganded HSA structure and may be more similar in conformation to the HSA we used experimentally, since the protein standard we purchased contained low levels of fatty acids (<0.02%). We followed the workflow outlined by Ng and Hungerbuehler (2015) with several modifications. Specifically, in the ligand-preparation step, we used the "Generating Conformers" function in DataWarrior V5.2.1 [45] to generate 3D structures for all ligands. Then, we optimized ligand structures using the MMFF94s forcefield in Avogadro V1.90.0 [46]. In addition, we used PyMol (v2.3.3) [47,48] for structure visualization, redocking alignment, and crucial residue identification. Simulations were repeated 100 times for 6 binding pockets, and each simulation generated 9 binding modes, yielding 5400 predictions in total for each PFAS. Further details on the docking method as well as simulation precision and accuracy are available in Supporting Information (S3.1 and S3.2).
The simulation method was evaluated by redocking PFOS on an experimentally determined HSA (Protein Data Bank entry 4E99) structure that was originally complexed with two PFOS in fatty-acid binding site (FA) 3/4 and 5. The atomic root-mean-squaredeviation (RMSD) of redocked PFOS on FA 3/4 and FA 5 was determined to be less than two angstroms, indicating successful redocking. The redocking search information and RMSD statistics are available in Table S2.3.

Characterization of PFAS in AFFF
Targeted analysis using 26 PFAS standards was insufficient for characterizing the AFFF sample: less than 9% of the total organic fluorine was quantified as compared to quantitative 19 F NMR (Table S7.1). In addition to the target PFAS, we identified 18 other PFAS and 18 hydrocarbon surfactant structures using suspect screening analysis for initial qualification. The suspected PFAS and hydrocarbon surfactant structures were further confirmed via data-dependent acquisition or library spectrum match. Manual annotation of the MS/MS spectra supported identification of these compounds (S1.2, Figures S4.1-S4.19). Based on structural categorizations conducted by the OECD [12], PFAS qualified in the AFFF sample included: 19 perfluoroalkane sulfonyl compounds, seven perfluoroalkyl carbonyl compounds, four fluorotelomer-related compounds, and two side-chain fluorinated aromatic compounds. Eight PFAS suspects were qualified in the initial screening but eliminated via manual confirmation. A full list of the hydrocarbon surfactants identified in AFFF is provided in Table S6.2. These hydrocarbon surfactants homologous series were detected using EnviHomolog (http://www.envihomolog.eawag.ch accessed on 5 November 2020). A repeating mass-increment of 14.0156 (-CH 2 -) was observed for four Linear Alkylbenzyl Sulfonates (LAS). Repeating mass-increments of 28.0313 (-C 2 H 4 -) and 44.0262 (-C 2 H 4 O-) were observed for 31 Alkyl Ethoxy Sulfates (AES).

Noncovalent and Potentially Covalent Binding of PFAS in AFFF to Human Serum Albumin
Of 32 PFAS identified in the AFFF, 28 PFAS bound noncovalently to HSA in equilibrium dialysis experiments. We confirmed 14 of these PFAS with analytical standards, and the remaining via data-dependent acquisition. Five PFAS were qualified and confirmed in the precipitated and washed protein pellets. Since PFAS released from hydrolyzed HSA pellets could not be extracted with the organic solvent, natural dissociation of this fraction of PFAS was not expected in a reasonable timeframe. Hence, these PFAS were considered as candidates for ultrastrong noncovalent or potentially covalent binding to HSA. Fourteen additional PFAS were qualified (confidence level 4) in the protein pellet, but were not qualified in the AFFF dilutions. We excluded these compounds from further analysis.

Figure 1.
Equilibrium dialysis set-up and observed mass balance for three PFAS with initial dosages of 40-80 ng. PFAS that were free in aqueous solution PFAS, noncovalently bound PFAS, and residual PFAS in the precipitated protein pellet were measured independently. The time required to reach equilibrium (Teq) was previously determined using 26 PFAS standards in Figure S1.2. Detailed mass balance data are available in Table S1.1.

Evaluation of Molecular Docking to Predict the PFAS-HSA Binding Affinities
Molecular docking of PFAS with two HSA crystal structures (1E7G and 1AO6) was used to simulate K A for the 26 PFAS tested experimentally (Table S7.2). Accurate K A predictions using 1E7G were limited to short-chain PFAS. In Figure 3, for the nine short-chain PFAS (PFCA with less than six perfluorinated carbons and PFSA with less than five perfluorinated carbons), a significant positive correlation between the docking-predicted K A and the experimentally determined values was observed (95% CI: slope = 1.02 ± 0.19, r = 0.900). For the 17 long-chain PFAS, significant negative correlation between the docking-predicted Molecular docking of PFAS with two HSA crystal structures (1E7G and 1AO6) was used to simulate KA for the 26 PFAS tested experimentally (Table S7.2). Accurate KA predictions using 1E7G were limited to short-chain PFAS. In Figure 3, for the nine short-chain PFAS (PFCA with less than six perfluorinated carbons and PFSA with less than five perfluorinated carbons), a significant positive correlation between the docking-predicted KA and the experimentally determined values was observed (95% CI: slope = 1.02 ± 0.19, r = 0.900). For the 17 long-chain PFAS, significant negative correlation between the dockingpredicted KA and the experimentally determined values was observed (95% CI: slope = −1.05 ± 0.30, r = 0.7680).

AFFF Formulation
Most PFAS identified in this study were also reported by Houtz

HSA Noncovalent Binding Affinity Relative to Perfluorocarbon Chain Length
Consistent with previous studies of PFAS-protein associations, PFAS were highly bound to HSA [24,25,41]. HSA contains multiple PFAS binding sites with potentially different binding affinities [18,52] such that measured K A values represent a mixture of affinities for different binding sites. A majority of the PFAS exhibited linear binding isotherms, indicating nonspecific noncovalent associations with HSA ( Figure S1.4 and Table S1.1). K A followed an inverted-V trend by which K A increased with perfluorocarbon chain elongation up to C 6 through C 9 and subsequently decreased ( Figure 2). The trend for C 4 through C 6 and C 8 through C 11 PFCAs is consistent with the pattern for bovine serum albumin (BSA)-water distribution coefficients (K PW ) determined by Bischel et al. (2011) [25]. The PFSA trend for C 4 through C 8 is consistent with the BSA-association constants determined by Allendorf et al. (2019) [24]. Our measurements of K A were generally an order of magnitude lower than K A measured by Allendorf et al. (2019), with the exception of PFBA, PFHpA, and C 6 PFAS. We have greater confidence in the physiological relevance of our experimental results as our results were obtained at physiologically relevant molar ratios of PFAS and HSA, and our K A values were determined from isotherm data rather than single-point experimentation. As low levels of AFFF exposures to humans are most common [53], we tested PFAS-HSA association constants from 0.001 to 0.1 PFAS:HSA.
Overall, a trend of increasing K A with perfluorocarbon-chain length was observed for PFCAs and PFSAs up to C 6 . The HSA binding affinities of perfluorohexanesulfonic acid (PFHxS) and perfluoroheptanoic acid (PFHpA) were exceptionally high (Log K A : 4.99 ± 0.44 and 5.53 ± 0.39, respectively). These observations are consistent with long blood plasma elimination half-lives reported for PFHxS in humans [54,55] and for PFHpA and PFHxS in pigs [56].

Residual PFAS in Precipitated Protein Pellets
Consistent detection of three PFAS in the digested protein pellets precipated from AFFF and PFAS standard exposure experiments is presented in Figure 4. Residual protein pellets from HSA exposed to PFAS standards (which did not include C 4 precursors) contained similar amount of PFBS, PFOA, and PFOS as protein pellets from HSA exposed to AFFF dilutions ( Figure 4 and Table S1.1). PFAS release from the protein pellet could be a result of several factors. First, residual PFAS in the pellet could be present as an analytical artifact resulting from high-concentrations of PFAS in AFFF. However, no quantifiable level of PFAS was observed in the protein pellet washes, blank cells, or control cells used in the equilibrium dialysis experiment (See Figure S1.1). Additionally, AFFF was diluted from 2000-fold to 80,000-fold prior to HSA exposure in the most dilute case, and the PFAS targets were still present in the precipitated protein from these tests. Second, PFAS could be retained in the pellet via ultrastrong noncovalently interations, which was observed in an amphiphilic polymerization system [21]. In conjunction with an exterior aqueous environment, large proteins like HSA that have multiple hydrophobic binding sites can provide an amphiphilic environment in which different protein residues interact with the polar headgroups and hydrophobic perfluorotails of PFAS simultaneously [57]. Third, PFAS could be retained in the pellet via covalent interations, the potential for which we evaluate in further detail below. carboxylic and sulfonic acids in the protein pellet. Formation of covalent bonds between carboxylic acids or sulfonic acids containing ligands and protein residues has not been observed under physiological conditions [59,60]. The low mole ratio of residual PFAS detected in the protein pellet to initial HSA levels indicates that not all PFAS-HSA associations resulted in PFAS retention in the protein pellets. While we are unable to disentangle the mechanisms explaining PFAS in the protein pellet, we consider residual PFAS in protein pellets as candidates for ultrastrong noncovalent or potentially covalent binding to HSA. In addition to the perfluoro alkyl acids (PFAAs) described above, we noted that 14 additional PFAS qualified in the dialysis extracts were not further evaluated in this study (i.e., by acquiring targeted MS/MS). It is possible that these PFAS were generated from reactions with HSA or through transformations during the acid-hydrolysis processing step. We would expect the formation of PFAS-HSA covalent bonds and subsequent digestion of PFAS-containing HSA to yield PFAS-peptide complexes. Future analysis should evaluate whether perfluorocarbon moieties are associated with amino acids or peptides following digestion.

Semi-quantification of Bioconcentration Factors of Qualified PFAS
We calculated pseudo-bioconcentration factors (BCFpseudos) for noncovalently bound and potentially covalently bound fractions separately to evaluate patterns related to perfluorocarbon chain length and functional groups ( Figure 5). BCFpseudo serves as a quantitative benchmarking technique to cross-compare bioaccumulation potentials of novel PFAS using qualitative screening data [29]. The BCFpseudo was calculated for each PFAS as follows:  To our knowledge, only two types of PFAS have been reported to bind covalently to proteins: PFCAs [19] and fluorotelomer unsaturated aldehydes (FTUALs) [20,58]. The mechanism of covalent binding between PFAS and proteins remains unclear, but thiol-and nitrogen-containing nucleophilic amino acids in serum albumin proteins are suspected to play a role [19]. In the case of FTUALs, covalent bond formation occurs via Michael addition [58]. This mechanism cannot explain our observation of perfluoro alkyl carboxylic and sulfonic acids in the protein pellet. Formation of covalent bonds between carboxylic acids or sulfonic acids containing ligands and protein residues has not been observed under physiological conditions [59,60]. The low mole ratio of residual PFAS detected in the protein pellet to initial HSA levels indicates that not all PFAS-HSA associations resulted in PFAS retention in the protein pellets.
While we are unable to disentangle the mechanisms explaining PFAS in the protein pellet, we consider residual PFAS in protein pellets as candidates for ultrastrong noncovalent or potentially covalent binding to HSA. In addition to the perfluoro alkyl acids (PFAAs) described above, we noted that 14 additional PFAS qualified in the dialysis extracts were not further evaluated in this study (i.e., by acquiring targeted MS/MS). It is possible that these PFAS were generated from reactions with HSA or through transformations during the acid-hydrolysis processing step. We would expect the formation of PFAS-HSA covalent bonds and subsequent digestion of PFAS-containing HSA to yield PFAS-peptide complexes. Future analysis should evaluate whether perfluorocarbon moieties are associated with amino acids or peptides following digestion.

Semi-quantification of Bioconcentration Factors of Qualified PFAS
We calculated pseudo-bioconcentration factors (BCF pseudos ) for noncovalently bound and potentially covalently bound fractions separately to evaluate patterns related to perfluorocarbon chain length and functional groups ( Figure 5). BCF pseudo serves as a quantitative benchmarking technique to cross-compare bioaccumulation potentials of novel PFAS using qualitative screening data [29]. The BCF pseudo was calculated for each PFAS as follows: where A i,sample is the peak area of compound i detected in the protein aliquot or pellet; A L-PFOS,sample is the peak area of L-PFOS detected in the protein aliquot or pellet; A L-PFOS,AFFF is the peak area of L-PFOS detected in neat AFFF dilutions; and A i,AFFF is the peak area of compound i detected in neat AFFF dilutions. All peak areas were normalized with their corresponding internal standard peak area prior to the calculations. Analysis of BCF pseudo revealed several key findings in Figure 4. First, sulfonic acids (C 4 through C 9 ) and carboxylic acids (C 2 through C 7 ) consistently exhibited high BCF pseudos in the noncovalently bound fraction. Second, C 7 PFAS across different functional groups consistently exhibited high BCF pseudos in the noncovalently bound fraction. We observed a noncovalent binding trend with a turning point at C 7 for all qualified PFAS, which is consistent with observations for targeted PFSAs (Figure 2). Third, PFOA exhibited the highest BCF pseudo for potentially covalently bound fractions. This is despite low levels in the AFFF; PFOA represented less than two percent of the total organic fluorine mass in the AFFF among all PFAS compounds quantified through targeted analysis. Finally, three PFAS with four or fewer perfluorocarbons exhibited higher BCF pseudos in the potentially covalently bound fraction than L-PFOS. The BCF pseudo for the C 2 carboxylic acid in the potentially covalently bound fraction (BCF diTF-IsoBA, pseudo = 13.8 ± 3.14) was an order of magnitude greater than the BCF pseudo for L-PFOS. This finding is concerning, as shortchain PFAS are typically considered less bioaccumulative than long-chain PFAS and exhibit weaker noncovalent interactions with proteins [61]. Our results indicate that short-chain PFAS may in fact be strongly retained in proteins (and in precipitated protein pellets) even when present at low concentrations in serum. Further studies should consider the impacts of these observations on analytical conclusions as well as potential toxicological risks. Toxics 2021, 9, x FOR PEER REVIEW 12 of 19 Figure 5. Pseudo-bioconcentration factors (BCFi,pseudo ) of PFAS in aqueous film-forming foams (AFFF) that were bound noncovalently (yellow) or potentially covalently (red) to HSA. Bubble size represents the natural logarithm of the BCFpseudo (legend BCFpseudo = 1). For 12 qualified PFAS with multiple functional groups, a separate bubble of the same size is displayed for each functional group (e.g., noncovalently bound Cl-PFOS contained chloride and sulfonic acid groups and is represented as two bubbles with seven perfluorocarbons).

Molecular Docking Predictions
Molecular docking simulations provide a rapid and high-throughput strategy to assess protein-ligand interactions. However, the accuracy of AutoDock Vina for the best protein-ligand conformations is only 60 to 80% [62]. Modeling assumptions integrated into AutoDock Vina include, for example, unrealistic rigidity of the protein and the removal of water [63]. The docking scores generated and shown in Figures 3 and S2.3 (for 1AO6) therefore cannot be used directly as protein-ligand binding energies, especially for large compounds. Nevertheless, we expected that docking predictions for short-chain PFAS would be more accurate than for long-chain PFAS, as docking accuracy declines with increasing numbers of rotatable bonds [42]. Indeed, we observed a positive correlation for short-chain PFAS between docking scores generated and our experimental results (Figure 3), while docking scores for long-chain PFAS were negatively correlated with experimental results.
To the best of our knowledge, Ng and Hungerbuehler (2015) [27] conducted the only other PFAS-HSA molecular docking study, predicting binding affinities of 25 PFAS with HSA. As shown in Figure S2.4, the two studies were in good agreement (least-squarelinear correlation predicted KA for 1E7G slope = 0.993, r = 0.997), even though our mean docking scores were calculated based on 5400 conformations (available in Table S7.2) instead of the 54 conformations simulated by Ng and Hungerbuehler. Unlike the previous study, we did not apply a 20-fold correction factor to decrease the predicted KA values, as this did not improve correlations with our experimentally determined KA values. Ng and Hungerbuehler may have required a correction factor due to the diversity of experimental results compiled in their study and unrealistic conditions in available studies. Many studies evaluated PFAS-HSA binding affinities using oversaturated protein conditions. For Figure 5. Pseudo-bioconcentration factors (BCF i,pseudo ) of PFAS in aqueous film-forming foams (AFFF) that were bound noncovalently (yellow) or potentially covalently (red) to HSA. Bubble size represents the natural logarithm of the BCF pseudo (legend BCF pseudo = 1). For 12 qualified PFAS with multiple functional groups, a separate bubble of the same size is displayed for each functional group (e.g., noncovalently bound Cl-PFOS contained chloride and sulfonic acid groups and is represented as two bubbles with seven perfluorocarbons).
To assess the ability of in-vitro binding studies with HSA to represent bioaccumulation potentials of PFAS in animals, we compared our BCF pseudos from the noncovalently bound fraction to pseudo-bioaccumulation factors (BAF pseudos ) calculated in an in-vivo mouse-dosing study that used the same AFFF commercial product [29]. Our calculation of BCF pseudo .

Molecular Docking Predictions
Molecular docking simulations provide a rapid and high-throughput strategy to assess protein-ligand interactions. However, the accuracy of AutoDock Vina for the best proteinligand conformations is only 60 to 80% [62]. Modeling assumptions integrated into AutoDock Vina include, for example, unrealistic rigidity of the protein and the removal of water [63]. The docking scores generated and shown in Figures 3 and S2.3 (for 1AO6) therefore cannot be used directly as protein-ligand binding energies, especially for large compounds. Nevertheless, we expected that docking predictions for short-chain PFAS would be more accurate than for long-chain PFAS, as docking accuracy declines with increasing numbers of rotatable bonds [42]. Indeed, we observed a positive correlation for short-chain PFAS between docking scores generated and our experimental results (Figure 3), while docking scores for long-chain PFAS were negatively correlated with experimental results.
To the best of our knowledge, Ng and Hungerbuehler (2015) [27] conducted the only other PFAS-HSA molecular docking study, predicting binding affinities of 25 PFAS with HSA. As shown in Figure S2.4, the two studies were in good agreement (least-square-linear correlation predicted K A for 1E7G slope = 0.993, r = 0.997), even though our mean docking scores were calculated based on 5400 conformations (available in Table S7.2) instead of the 54 conformations simulated by Ng and Hungerbuehler. Unlike the previous study, we did not apply a 20-fold correction factor to decrease the predicted K A values, as this did not improve correlations with our experimentally determined K A values. Ng and Hungerbuehler may have required a correction factor due to the diversity of experimental results compiled in their study and unrealistic conditions in available studies. Many studies evaluated PFAS-HSA binding affinities using oversaturated protein conditions. For instance, Chen and Guo (2009) measured the HSA association constants of PFAS at ligand:protein ratio ranging from 4 to 40 in their fluorescence quenching study [64]. Han et al. (2003) used microdesalting columns to measure the PFAS-HSA association constants with a ligand:protein ratio ranging from 4 to 36 [65]. When the ligand-to-protein ratio is greater than one, PFAS can shift from primary binding site(s) to secondary binding sites [41]. PFAS associations with secondary binding sites on HSA are not physiologically relevant, as human exposure to PFAS generally occurs at low levels [66]. Our experimental data represent binding affinities at primary binding site(s) because we made direct measurements with high-resolution-mass-spectrometry at low ligand:protein mole ratios.

Qualitative Prediction of HSA-bound vs. Nonbound Compounds
While docking scores are unreliable for quantitative predictions of binding affinities, docking scores can be used for qualitative comparisons between complexes and to identify candidate ligands for the protein of interest. We used AutoDock Vina to identify HSAbinding compounds in AFFF by comparing PFAS and hydrocarbon surfactants in AFFF. Results of docking are presented as violin plots, which display the distribution of simulated docking scores (Figures S3.1-S3.6). For each PFAS, a kernel density plot of docking scores was derived using the 5400 conformations generated from simulations on six HSA binding pockets. All hydrocarbon surfactants identified in AFFF had more than 10 carbons in the backbone, so we selected docking results for PFAS with 10 or more perfluorocarbons in the backbone for comparison ( Figure 6). For both sets of long-chain PFAS and hydrocarbon surfactants, we observed a separation of distributions that appeared to distinguish binding (low docking scores) or non-binding (high docking scores) compounds. The low and high docking scores corresponded, respectively, to PFAS that were observed as bound or unbound to HSA in in-vitro experiments. We performed a Kruskal-Wallis test of docking predicted binding scores for the compounds in Figure 6. Most compounds in AFFF that we observed to bind to HSA in the in-vitro experiments also had predicted docking scores that were significantly different from the unbound compounds (p < 0.05 in Table S7.3). However, the binding score of one experimentally determined unbound PFAS (N-diMAmP-PBSAP) and three experimentally determined bound hydrocarbon surfactants (C10-LAS, C11-LAS, and C12-LAS) were not significantly different from each other (p = 1, in Table S7.3, labeled with black X in Figure 6). Experimental results for suspect screening were thus largely consistent with simulated docking scores when comparing HSA binding affinities within the same class of chemicals. (a) (b) Figure 6. Violin plots of molecular docking simulated 1E7G docking scores for (a) PFAS with greater than 10 perfluorocarbons and (b) hydrocarbon surfactants identified in AFFF that have greater than 10 carbons. The shape of each violin represents a rotated kernel density plot of 5400 HSA-PFAS binding conformations generated from simulations for six binding pockets. Blue and red colors are used to distinguish experimental results. The compounds with significantly greater peak area (after correction with ISTD peak area) in the noncovalently bound fraction of the protein chamber relative to the chemical chamber are shown in blue. The compounds identified experimentally in the protein chamber that were not significantly greater in peak area relative to the chemical chamber are shown were red. Distributions in red were significantly different than distributions in blue Kruskal-Wallis (p < 0.05) except for those distributions marked with a black X.

Conclusions
Our study explored the value of suspect screening and computational simulations to identify potentially bioaccumulative PFAS from a PFAS-containing commercial product, Figure 6. Violin plots of molecular docking simulated 1E7G docking scores for (a) PFAS with greater than 10 perfluorocarbons and (b) hydrocarbon surfactants identified in AFFF that have greater than 10 carbons. The shape of each violin represents a rotated kernel density plot of 5400 HSA-PFAS binding conformations generated from simulations for six binding pockets. Blue and red colors are used to distinguish experimental results. The compounds with significantly greater peak area (after correction with ISTD peak area) in the noncovalently bound fraction of the protein chamber relative to the chemical chamber are shown in blue. The compounds identified experimentally in the protein chamber that were not significantly greater in peak area relative to the chemical chamber are shown were red. Distributions in red were significantly different than distributions in blue Kruskal-Wallis (p < 0.05) except for those distributions marked with a black X.
The shape of the kernel density plots may also provide insights into different binding processes. The kernel density plots for both bound and unbound PFAS (Figures 6a and S2.5-S2.7) are more varied in shape than those for hydrocarbon surfactants (Figure 5b and Figure S3.0). Similar to kernel plots, cluster analyses were commonly used to identify preferential ligand binding sites [67], suggesting site-specific binding between PFAS and HSA [68]. Predicted HSA binding scores for bound hydrocarbons exhibited bimodal or even trimodal (C10-LAS and C11-LAS) distributions. Predicted HSA binding scores for nonbound hydrocarbons converged and centralized for all HSA binding pockets ( Figures S2.8 and S2.9). To further validate and accurately predict PFAS-HSA binding energies, mechanistic studies using molecular dynamics coupled with molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) methods for these PFAS are ongoing in our research group.

Conclusions
Our study explored the value of suspect screening and computational simulations to identify potentially bioaccumulative PFAS from a PFAS-containing commercial product, AFFF. A majority of the PFAS we identified in a legacy AFFF bind to the most abundant protein in human serum, human serum albumin (HSA). At least five PFAS, including two PFAS with less than five perfluorocarbons, were detected in the precipitated and washed protein pellet. The potential health implications of ultrastrong or covalent binding of PFAS are unclear. Covalent modifications of HSA affect the clearance and metabolic destiny of many drugs, and have been hypothesized as the center of toxicity exhibited by many drugs [69][70][71].
Our observation of binding of short-chain PFAS to albumin is concerning and requires further mechanistic assessment. Short-chain PFAS are largely considered less bioaccumulative, with shorter half-lives in organisms, than long-chain PFAS. Our results indicated that some short-chain PFAS may be retained in the blood for much longer-these PFAS remained associated with HSA after extensive solvent washing. Computational simulations for bioaccumulation potential can provide value by decreasing reliance on time-and labor-intensive laboratory experiments. Though predicted binding scores cannot quantitively describe binding affinities with HSA, the scores can be used to qualitatively identify previously uncharacterized PFAS that are likely to bind to HSA. More broadly, this study offers a framework for evaluating bioaccumulation potentials of thousands of PFAS in comparable biological tissues.