A food web bioaccumulation model for the accumulation of per- and polyfluoroalkyl substances (PFAS) in fish: how important is renal elimination?

Per- and polyfluoroalkyl substances (PFAS) are a large class of highly fluorinated anthropogenic chemicals. Some PFAS bioaccumulate in aquatic food webs, thereby posing risks for seafood consumers. Existing models for persistent organic pollutants (POPs) perform poorly for ionizable PFAS. Here we adapt a well-established food web bioaccumulation model for neutral POPs to predict the bioaccumulation behavior of six perfluoroalkyl acids (PFAAs) and two perfluoroalkyl ether acids (HFPO-DA, 9-Cl-PF3ONS) produced as PFAA replacements. The new model includes sorption to blood plasma proteins and phospholipids, empirically parameterized membrane transport, and renal elimination for PFAAs. Improved performance relative to prior models without these updates is shown by comparing simulations to field and lab measurements. PFAS with eight or more perfluorinated carbons (ηpfc ≥ 8, i.e., C8 perfluorosulfonic acid, C10–C11 perfluorocarboxylic acid, 9-Cl-PF3ONS) are often the most abundant in aquatic food webs. The new model reproduces their observed bioaccumulation potential within a factor of two for >80% of fish species, indicating its readiness to support development of fish consumption advisories for these compounds. Results suggest bioaccumulation of ηpfc ≥ 8 PFAS is primarily driven by phospholipid partitioning, and that renal elimination is negligible for these compounds. However, specific protein binding mechanisms are important for reproducing the observed tissue concentrations of many shorter-chain PFAAs, including protein transporter-mediated renal elimination. Additional data on protein-binding and membrane transport mechanisms for PFAS are needed to better understand the biological behavior of shorter-chain PFAAs and their alternatives.

unitless -∅ Chemical fraction in the dissolved phase Varies by ecosystem, see Tables S4b and S5b unitless varies Organic carbon content in sediment Varies by ecosystem, see Table S5b kg/L varies Organic carbon -water partitioning coefficient Varies by ecosystem, see Table S2 L/kg varies Contaminant concentration in the diet (weighted average between prey items and sediment)    Table S3a. 19 c -Partitioning coefficient assumes a protein density of 1.36 kg/L to convert from units of kg/L to L/L. 20 d -Empirical values from the same or similar study site should be used whenever possible. 21 e -Values were calculated based on calculation from reported uptake rate, assuming a fish weight of 7g and oxygen concentration 22 of 7.23 mg/L. The italicized value was estimated from an uptake rate that was estimated from a log-linear regression between chain 23 length and uptake rate.  Phospholipids --Partitioning to phospholipids is described by the membrane-water partitioning coefficient . In this 39 model, is parameterized using empirical values from a laboratory-based partitioning study using solid-supported 40 phosphatidylcholine lipid bilayer membranes designed to mimic intestinal epithelium 7 . Similar values were measured in a 41 second laboratory study for all PFAAs except the C11 PFCA, for which the measured was lower than the C10 PFCA 4 . However, 42 the second study was inconclusive as to whether this was due to the lower diffusibility of larger compounds or experimental error, 43 and therefore this experimental value is not used in the current study. Empirically measured Dmw values are also similar to those 44 calculated by the mechanism proposed in the ionogenic model, when based on Kow estimates from COSMOtherm (2011), 45 suggesting that these modeled values may be reasonable approximations for compounds without empirical Dmw 46 measurements (see Table S3a).

47
Following Armitage et al. 2 , the volume fraction of phospholipids ( ) is estimated at 1% in fish. In this study, is 48 approximated at 1% in all other organisms as well, although it is likely that may be more variable in lower trophic level organisms, 49 particularly phytoplankton. 50 51 52 Albumin proteins -PFAA binding to albumin is widely believed to drive high concentrations observed in blood plasma, and 58 has been measured with human and bovine serum albumin. Albumin binding patterns show significant species-specific variation 8 , 59 but PFAA partitioning coefficients have not been measured for any fish-specific proteins. Available binding and partitioning studies 60 using bovine serum albumin (BSA) and human serum albumin (HSA) have shown widely variable results due to a diversity of 61 measurement methods, including poor standardization of ligand (i.e. chemical) and receptor (i.e. protein) concentrations used 62 during measurement 8 . The three available albumin partitioning studies show contrasting patterns for headgroup and PFCA chain 63 length partitioning to BSA. Bischel et al. 9 reported decreasing partitioning strength with increasing chain length for the C8-C12 64 PFCAs, and nearly identical partitioning coefficients for identical chain length PFCAs and PFSAs, whereas Allendorf et al. 10 and Aleiso 65 et al. 11 found increasing partitioning strength with increasing chain length, and greater chain-length dependence for PFCAs 66 compared to PFSAs (Table S3b). Allendorf et al. 10 attributed discrepancies to differences in study design: all studies used equilibrium 67 dialysis, but Bischel et al. 9 conducted their experiment at a higher PFAA:albumin molar ratios, such that over 98% of PFAAs were 68 bound at equilibrium, which can in turn lead to oversaturation of protein binding sites that leads to an underestimation of true 69 partitioning strength. The pattern observed by Allendorf et al. 10 also align more closely with patterns of fish blood protein binding 70 strength measured by Zhong et al. 12 in carp, as well as overall patterns of bioconcentration and biomagnification observed in 71 laboratory studies of fish 3,5 , in which accumulation or partitioning is stronger for PFSAs compared to PFCAs of the same chain length. 72 We use values from Allendorf et al. 10 in this model.

73
The presence and abundance of albumin proteins in fish can be highly variable, with albumin making up anywhere from 0 to 74 60% of total blood proteins in different species 13,14 . The volume fraction of albumin binding proteins is estimated to range between 75 0 to 0.4% in fish (SI Section 7) based on literature values for albumin or albumin-like protein concentrations and tissue volumes 76 estimated for an average fish. In this model, the volume fraction of albumin or albumin-like binding proteins is set at 0.3% for fish, 77 0.1% for invertebrates, and 0% in phytoplankton. 78 79 Liver fatty acid binding proteins -PFAAs have also been shown to bind to liver-and other fatty-acid binding proteins, which 83 are considered to be drivers of elevated PFAA concentrations in the liver and kidney. The pattern of human L-FABP binding strength 84 across PFAA structures is similar to the pattern observed for fish blood proteins, with higher binding affinities for PFSAs compared to 85 PFCAs of the same chain length, and increasing binding affinities with chain length for the C4 to C8 PFSAs and C7 to C11 PFCAs 12,15 . 86 No partitioning coefficients have been measured for any L-FABP, and no published studies have measured PFAA binding to L-FABP in 87 fish species, but lower binding strengths have been reported for human L-FABP compared to albumin 8,11,15,16 . The relative 88 abundance of L-FABP is also substantially lower than that of albumin 17 . We estimate that the contribution of L-FABP to total PFAA 89 binding is not larger than uncertainties in albumin binding strength and abundance (Section 7). Therefore, L-FABP and other 90 potential binding proteins are not explicitly included in this model. Branchial uptake and elimination -Membrane transport of ionized compounds is generally significantly reduced compared 97 to that of neutral POPs, with transport dominated by the neutral fraction of the compound. For weakly ionized compounds, diffusive 98 transport by the neutral chemical fraction results in pH-dependent transport rates. However, for highly ionized compounds, the 99 relative importance of other transport mechanisms, such as paracellular transport, protein-mediated or ion channel transport, or 100 mass transfer of the ionic fraction, are hypothesized to increase as the size of the neutral chemical fraction grows negligibly small. 101 Saarikoski et al. 18 found that for highly ionized compounds, as pH increases chemical uptake initially decreases but then plateaus. 102 Armitage et al. 2 updated the Arnot & Gobas 1 submodel for aqueous chemical transfer efficiency ( ) for ionic chemicals to reflect 103 these behaviors by (1) modeling the pH-dependence of diffusive transport and (2) adding a constant to crudely account for non-104 diffusive transport pathways, which increase in importance as diffusive transport rates decrease for highly ionized compounds. 105 In a series of membrane permeability experiments using phosphatidylcholine liposomes, Ebert et al. 4 found in vitro lipid 106 bilayer membrane permeabilities to be consistent with pH-independent passive transport for PFSAs, and only weakly pH-dependent 107 passive transport for PFCAs. For PFSAs, the effective permeability of the ionic fraction is about eight orders of magnitude greater 108 than that of the neutral fraction. The empirically measured permeabilities matched permeabilities calculated from cellular uptake 109 studies quite well, suggesting they reasonably represented real behaviors in cells 4,17 . Model-predicted permeabilities using either 110 correlation with the hexadecane-water partitioning coefficient or COSMOtherm similarly suggest that permeability of the ionic 111 fraction of PFSA is orders of magnitude greater than the effective neutral permeability at biological pH. The authors attributed 112 greater permeability of the ionic fraction of PFSAs compared to PFCAs (SI Tables 3.2 For PFCAs, both the neutral and ionic chemical fractions play a role in transport, with up to an order of magnitude difference 117 between the effective membrane permeability for the neutral and ionic fractions. This suggests that pH-dependent transport of the 118 neutral chemical fraction may play a larger role for PFCAs, but at much lower rates than expected based on models developed for 119 more weakly ionized compounds 2 . While the Armitage et al. 2 model predicts a two order of magnitude difference in permeability for 120 PFCAs between pH 6 and 8, Ebert et al. 4 estimated a variability of only about 30%, based on passive uptake rates measured in vitro 121 in human HEK293 cells 19 . Here we approximate membrane transport of all PFAAs with a pH-independent parameterization of 122 membrane transport at environmentally relevant pH. 123 The parameter is estimated from branchial uptake rates from lab bioconcentration studies based on the modeled 124 relationship between these values and gill ventilation rate (see Table S1). Ew values calculated from the Martin et al. 3 125 bioconcentration study were based on reported uptake rates and gill ventilation rates calculated based on a reported average fish 126 weight of 7.3g and an estimated oxygen concentration of 10 mg/L (92% dissolved oxygen saturation at 12 °C). The uptake rate for 127 the C9 PFCA was not directly measured, and was therefore instead interpolated from a linear regression between uptake rate and 128 PFCA chain length. For the current model, calculated from the laboratory bioconcentration study is directly used in all model 129 applications based on an assumption that aqueous chemical transfer efficiency does not vary substantially with pH and therefore is 130 comparable between all study conditions. In general, the empirically estimated values follow membrane permeability patterns 131 observed in other literature, showing that total permeability is greater for PFSAs than PFCAs, and that permeability increases with 132 chain length (studies cited in Ng & Hungerbuhler 17 , Ebert et al. 4 ).

134
Gut uptake and elimination -Gut membrane chemical transfer efficiencies ( , also referred to as absorption or chemical 135 assimilation efficiencies) are generally not observed to be significantly different for ionogenic chemicals compared to neutral 136 chemicals, due to the long residence time of chemicals in the gut that enables a longer period for transport to occur compared to gill 137 membrane transport 20 . However, the submodel for dietary absorption efficiencies in the bioaccumulation model for neutral organic 138 chemicals is based on an inverse relationship between absorption efficiency and hydrophobicity 1 , whereas absorption efficiency 139 values modeled from uptake rates measured in laboratory biomagnification studies on both adult-and juvenile rainbow trout show 140 dietary absorption efficiency increases with hydrophobicity 5,6 . Thus, empirical values from these laboratory studies are used in this 141 current study. 142 for juvenile rainbow trout was estimated by the experimental study authors by fitting growth-corrected fish 143 concentration data to a kinetic rate equation for constant dietary exposure. Reported values therefore may reflect potential 144 measurement errors for fish and food concentrations, experimental and model assumption errors, and statistical artefacts of 145 averaging across natural variability in experimental values. This may include factors like the assumption of a steady feeding rate, the 146 assumption that 90% of steady state was reached during the experiment, and growth rate-correction for fish concentrations. 147 Reported assimilation efficiencies were greater than 100% (i.e. > 1) for the C8 PFSA and C10-13 PFCAs, which is not physically 148 possible. In this study, values were set at 1 when assimilation efficiencies were greater than 100%. Although these values 149 likely contain some error, we use them as reported in our biomagnification model, because our model relies on many of the same 150 inputs and assumptions as the model used to estimate (i.e. fish concentrations, food concentrations, feeding rate). By using these 151 parameters together, we may correct for co-occuring errors, while still accurately reflecting measured biomagnification and 152 depuration rates that we then attribute to specific tissue partitioning and elimination mechanisms, respectively. 153 Other studies corroborate the overall finding that dietary assimilation efficiencies for PFAAs in fish are high. Assimilation 154 efficiencies for adult trout more physically plausible, but still high (> 0.1) compared to branchial assimilation efficiencies. These 155 values are used in the current study for the model application to field food web data. High dietary assimilation efficiencies compared 156 even to neutral POPs (i.e. > 0.5) have also been assumed in other toxicokinetic models of PFAS uptake, such as in humans 21 3.5x10 -4 3.5x10 -4 3.5x10 -4 3.5x10 -4 3.5x10 -4 3.5x10 -4 1.4x10 -3 1.4x10 -3 1.4x10 -3 1.4x10 -3 1.4x10 -3 1.4x10 -3 182 1 -This parameter represents organic carbon for phytoplankton and NLOM for all other species. Partitioning to these different 183 fractions utilizes different proportionality constants (see Table S1). The NLOM proportionality constant is estimated to be between 184 0.025 and 0.05 2 , while the OC proportionality constant is estimated to be 0.35 1 . 185 186 187 188 189    Table S1) k1 43.6 L/kg/d Calculated (see Table S1) k2 2.42 x 10 -2 d -1 Calculated (see Table S1) kg 2.5 x 10 -2 d -1 Estimated from Figure 4 Table S1) k1 5.93 x 10 -1 L/kg/d Calculated (see Table S1) k2 8.11 x 10 -3 d -1 Calculated (see Table S1) kg 7.24 x 10 -3 d -1 Calculated 1 Bioconcentration factor Log BCF (Model) 0.587 L/kg k1 / (k2 + kg) Log BCF (Obs) 0.591 L/kg 29 (see Table S6c) Model Bias 1.01 207 a -This is the same method that was used for PFAAs. Uptake rates were reported separately for different organs, so a whole-body uptake rate 208 was estimated using a calculation analogous to the one used to calculate a whole-body BCF (Table S6c), using the same organ volume 209 distribution shown in Table S6c.  210  211  212 213 Bioconcentration factors for HFPO-DA were reported for the carcass, fillet, plasma, and liver. A whole-body concentration factor was calculated 214 based on estimated organ sizes, as estimated for the fish modeled by Ng Table S1 for an 8 g fish 219 b -Assumes blood plasma is approximately 55% of whole blood volume 220 c -Calculated as the remaining fraction 221 d -Calculated as a volume-weighted average 222 223 Section 7. Estimation of fish protein content 224 225 Albumin volume % calculation 226 1. Calculate whole blood volume percent in fishgiven the volumes and body weight percent for each organ tissue modeled in Ng 227 & Hungerbuhler 17 (Table S1, assuming a similar density between organs), the volume of a whole fish is approximately 7.97 mL. 228 Given blood, but varies from 0-60% across the several fish species that have been measured 13,14 . In rainbow trout, the species studied 235 in lab-based evaluation dataset in the current study, the albumin concentration has been measured at about 13.8 +/-0.5 g/L, or 236 about 38% of total protein (35.9 +/-1.3 g/L) 31 ; similar albumin:total protein ratios have been observed in a range of species, 237 including rainbow trout, channel catfish, tilapia, striped bass, salmon, three species of carp 14,32,33 . Assuming a similar ratio 238 between total protein content and blood density in humans and fish, we would expect the albumin contribution to whole blood 239 for these common fish species to be about two-thirds of the value estimated for human blood.

241
Albumin content in fish can be highly variable among species, and or the same species under varying environmental conditions 242 or life stages. It is also possible that some species that do not express albumin can instead express other proteins serving similar 243 functions, and that can also bind PFAAs. For example, Zhong et al. 12 report PFAA binding affinities for "fish blood proteins" in 244 common carp, which have been reported not to express albumin 13 , but do express lipoproteins that serve similar functions. The 245 lipoprotein:protein ratio measured in one carp study was also in the 40% range, although it was noted that the lipoprotein 246 concentration varied seasonally and based on diet 34 . This variability contributes to making albumin volume percent, or binding 247 protein percent more generally, one of the more uncertain parameters in the current model. 248 249 3. Blood albumin % in whole body -Using an approximate maximum albumin content in whole blood of 3% and a blood volume 250 percent of 3.8%, we calculate that albumin in the blood is about 0.11% of total body weight. bioconcentration model for PFAAs, the estimated fraction of albumin in whole blood is 75.5% (see Table S6a). In contrast, based 272 Liver fatty acid binding protein 273 As discussed in Section 3, L-FABP was not included in this model because its estimated total contribution to protein binding is 274 negligible compared to uncertainty in the albumin volume percent parameter. Protein binding studies using mammalian proteins 275 suggest that PFAAs bind more weakly to L-FABP than albumin, but in this calculation we conservatively assume similar partitioning to 276 both proteins. Based on the protein abundance estimated in Table S7a, the addition of L-FABP increases the total protein binding 277 (albumin + L-FABP) pool by approximately 6% (4.55 nmoL / 84.01 nmoL = 0.057). This in turn increases our estimate of the fraction of 278 total body volume made up of binding proteins by 6%, from 0.15%-0.38% (albumin only) to 0.16%-0.40%. Our current estimate of 279 binding protein content in the model (0.3%) remains well within this revised range of albumin + L-FABP content. Therefore, while L-280 FABP is not explicitly included as a separate compartment in this model, it could be considered accounted for within the 281 compartment currently attributed to albumin proteins alone. 290 b -Calculated using uptake rate from a linear regression between PFCA chain length and uptake rate 291 c -Calculated using measured renal-to-total elimination ratio and modeled branchial elimination rates 292 d -Calculated using linear regression (see Figure S7)' 293 294 295 296 297 Figure S8. Estimation of renal elimination rates from linear regression between renal elimination rate and renal clearance-to-298 reabsorption ratio for the C8 PFSA, C8 PFCA and C12 PFCA. Linear regressions were calculated separately for the bioconcentration 299 (circles, solid line, R 2 = 0.98) and biomagnification (squares, dotted line, R 2 = 0.98) datasets. For the remaining PFAAs, renal 300 elimination rates for fish in the bioconcentration and biomagnification studies were estimated from the linear regression 301 relationships based on renal clearance-to-reabsorption ratios calculated from data reported by Weaver et al. 36 Figure S10. Comparison of model-predicted PFAA concentrations for both benthopelagic and fully benthic fish from the Gironde 321 Estuary, France, simulated using observed prey data or no dietary uptake. The solid black line represents a 1:1 line (perfect 322 agreement), while the dotted lines represent a factor of two and a factor of 10 difference between modeled and observed values. 323 Circles show modeled values based on median water/sediment/prey concentrations and error bars represent minimum and maximum 324 reported exposure concentrations. Observed variability likely includes variability not included in the model, such as variability in 325 organism or environmental parameters. Results show that dietary uptake plays an important role in total exposures for all fish, but has 326 greatest influence in the benthic food web. 327