The possible interaction of the PFAS series with the listed protein targets has been estimated by means of molecular docking simulations performed with Glide (Schrodinger 2022–4). The procedure returns, for each possible interaction, a score measured in kcal/mol, which estimates the binding free energy as a sum of contributes, according to Friesner et al. (2006).

Then for specific cases, MD simulations were run, to verify the binding stability, in particular in the case of highly flexible proteins, whose binding site might exhibit different residue orientation than those obtained through docking simulation, and where ligands might assume unexplored or additional poses. MD simulations were also run to compare the stability of legacy PFAS in the considered targets, with respect to novel PFAS, expected to have lower affinity and lower residence time (Bruno et al. 2022).

HSA is known to be the PFAS primary vector in the human body (Forsthuber et al. 2020; Crisalli et al. 2023), thanks to PFAS structure being like that of fatty acids (Gao et al. 2019; Alesio et al. 2022).

Multiple possible pockets, named fatty acid (FA) pockets, have been identified by means of X-ray diffraction and docking simulations, but only few of them are known to bind PFAS with varying affinity (Maso et al. 2021; Perera et al. 2022; Moro et al. 2022). In particular, four pockets have been found occupied by PFAS in PDB ID 7AAI (Maso et al. 2021; Alesio et al. 2022), but further structural analysis and literature search revealed that up to five additional pockets could accommodate molecules as big as myristic acid.

We, thus, decided to expand our docking calculations to all nine pockets (Fig. 1).

a HSA pockets shown as red contours (PDB ID 7AAI). The protein is shown as cartoon. The red labelled pockets are occupied by PFOA. b PFOA docked in FA4 pocket. c cC6O4 docked in FA4 pocket. Ligands and lining residues are shown as capped sticks, H-bonds as yellow dashed lines

As can be observed by the docking scores reported in Fig. S2, narrower pockets FA1, FA2, FA4 and FA5 (Zunszain et al. 2003; Moro et al. 2022) easily accommodate linear PFAS, whilst larger ones also allow bulkier PFAS to bind, even if with overall low affinities. PFOA, known as one of the best binders in FA4, resulted in a minimum of − 13.4 kcal/mol, whilst slightly smaller or larger molecules, like PFHxA and PFNA, have higher scores and therefore lower affinity. A similar trend was observed for sulfonic acids, with a minimum of − 10.8 kcal/mol for PFDS. PFBA constitutes a peculiar case in HSA, and in other targets. Indeed, with respect to longer PFAS likely involved in a higher number of hydrophobic contacts, PFBA showed low scores, that is high affinity, in different pockets, despite the short hydrophobic tail. In general, the highest affinity site for PFOA, i.e. FA4, shows a trend similar to the results reported in literature (Moro et al. 2022; Crisalli et al. 2023): PFOA and PFHxA show the highest affinity, whilst short-chained cC6O4 and GenX present lower affinity than linear, longer PFAS. PFOO, which lacks a functional group, shows affinities of about -7 kcal/mol in different pockets, but very few poses. This is reasonably given to the absence of polar interactions and is indicative of an unfavourable binding.

Overall, the results obtained for the other pockets are in line with those described for FA4, with 8–10 carbon long PFAS showing the most promising interaction, few to no poses for longer molecules, i.e. from 12 to 18 C atoms, and PFBA being present and with good, predicted affinity, despite the absence of experimental indication about a possible interaction with HSA. Interestingly, less hydrophobic, and more branched novel PFAS show, in general, less stable interactions, likely suggesting a lower probability of being bound and transported by HSA. No appreciable differences were observed amongst the analysed PFAS when docked in low affinity sites, as FA6 or cleft.

The X-ray structure of HSA complexed with four PFOA molecules (PDB ID 7AAI; sites FA4, FA6, FA7, cleft) was then submitted to a 500 ns-long plain MD simulation. To compare PFOA binding and stability in HSA to that of the new generation PFAS cC6O4 (Bruno et al. 2022), we built a construct of HSA complexed with four cC6O4 molecules, obtained by IFD, and submitted it to an equal long MD simulation. The RMSD calculation showed the protein remain stable and compact in both simulations (Figure S3). The stability of the different ligand copies in the four pockets is shown in Fig. S4, where the ligand RMSD, calculated by centring the least square fitting on the pocket residues, and the corresponding moving mean are reported, to better show the ligand movement in and out the pockets.

Both PFOA and cC6O4 appear to be quite stable in FA4, as the RMSD profiles do not present significant variations (Figure S4a). This is also confirmed by the H-bond persistency, showing the maintenance of the interaction with positively charged residues for more than 90% of the simulation time, as shown in detail in Fig. S5. A similar trend can be observed for FA7 pocket, even if, on average, cC6O4 presents slightly larger fluctuations, indicative of a less stable interaction (Figure S4c). On the contrary, the ligands show a completely different behaviour in the solvent-exposed and low-affinity pocket FA6, where PFOA does not undergo significant variations, whilst cC6O4 at ~ 200 ns reports a huge RMSD increment (Figure S4b). This change is given to an abrupt loose of contact with Arg207, and to the subsequent movement of the molecule along the protein surface towards Lys349 (Figure S5b), with which it forms a stable H-bond. The case of the cleft pocket is quite peculiar: it is a largely solvent-exposed pocket and there is only one paper reporting about it being occupied by a ligand (Maso et al. 2021). Both PFOA and cC6O4 are quite unstable in the pocket, suggesting this site can be likely occupied only at very high ligand concentration. The number of contacts between the ligands and any atom within 0.6 nm are reported in Fig. S6. Overall, MD simulations seem to suggest PFOA can be transported by HSA with higher probability than cC6O4. This is in particular highlighted by all performed analyses (Figs. S4–S6), showing that, in all pockets except the cleft, PFOA is more stable than cC6O4. As described in literature (Crisalli et al. 2023), this difference can be mainly attributed to entropic effects arising from different tail length, whilst enthalpic ones, like H-bonds, remain reasonably constant.

These observations have been experimentally confirmed by measuring the percentage of PFOA and of cC6O4 binding to bovine serum albumin (BSA) through dialysis equilibrium experiments. The percentage of PFOA and cC6O4 bound to BSA corresponds to 82% and 62%, respectively (see Table S2). These experiments also highlighted a degree of absorption of compound PFOA (recovery 85%) on the plastic material of dialysis devices, phenomenon not observed in the case of cC6O4 (recovery 99%). These results suggest that novel PFAS should have a lower probability of binding albumin with respect to legacy PFAS.

TTR presents two identical and symmetrical binding sites (AC and BD), in which the natural substrate thyroxine (T4) is stabilised by hydrophobic contacts with most of the residues lining the pockets, and by electrostatic interaction with Lys15 and Lys15’ (PDB ID 2ROX, (Wojtczak et al. 1996). Similar contacts are made by PFOA (PDB ID 5JID (Zhang et al. 2016)). According to Purkey et al. 2001, and Tomar et al. 2012, TTR binds thyroxine T4 and other ligands with a stoichiometry lower than 2 in physiological conditions. Considering this indication, the negative cooperativity suggested by Neumann et al. (Neumann et al. 2001), and the fact that TTR is present in the serum at a concentration much higher than PFAS (Ingenbleek and Bernstein 2015; Göckener et al. 2020), we hypothesised that in standard conditions, only one PFAS molecule at a time would be able to bound TTR. We, thus, docked the PFAS series in only one of the two pockets (Figure S1b), assuming the same results could be expected for the second site, considering the two are basically identical (Dharpure et al. 2023).

When re-docked in TTR, PFOA superimposed very well to the co-crystallised ligand (Figure S7), thus further validating the docking procedure, and showed a score about -7 kcal/mol. The latter corresponds to one of the lowest (Figure S8), even if there are not extremely significant differences amongst the considered PFAS. In general, the profound hydrophobic channel of TTR easily accommodates PFAS carbon chain, whilst the two gating lysines, Lys15 and Lys15’, stabilise the negatively charged moiety by means of salt bridges. Specifically, medium-length PFAS, (6–10 C atoms) show a mediocre affinity for TTR, with PFDA resulting in a minimum of -7.3 kcal/mol, whilst PFNA and PFOA differ by less than 0.5 kcal/mol (− 7.1 and − 7.0 kcal/mol, respectively). Smaller or bigger molecules have slightly higher scores (e.g. − 5.7 kcal/mol for PFBA), and therefore slightly minor affinities for this target. cC6O4 returned a score of − 5.6 kcal/mol, more similar to low affinity molecules than medium ones. These results agree with those of (Ren et al. 2016), who reported higher affinity for PFOA (IC50 equal to 378 nM) and lower for PFDA and PFNA (1623 and 1977 nM, respectively).

MD simulations were carried out with the intent of investigating legacy and novel PFAS stability in TTR, again taking as reference PFOA and cC6O4. Specifically, 200 ns-long MD was performed in triplicate, for a total of 600 ns. As shown by RMSD calculation (Figure S9), the protein is stable throughout the whole dynamics, only showing higher flexibility in loops and terminal regions for each monomer. Interestingly, whilst PFOA maintains its pose favourably interacting with the hydrophobic residues lining the pocket and H-bonding the two gating Lys15 and Lys15’, cC6O4 shows higher fluctuations, meaning a reduced stability in the pocket. Moreover, one in three replicas shows displacement of cC6O4 exiting the funnel-shaped binding pocket (Figure S10). The lower number of contacts that cC6O4 forms with apolar residues, compared to PFOA, is mainly due to its shorter and cyclic chain. In addition, the stronger interactions between the ligand carboxylate and the positively charged amino groups of Lys15 and Lys15’ may help to explain the ligand movement towards the exterior of the pocket. T4, which is in pocket A (Figure S1b), shows different stability depending on whether the protein is complexed with cC6O4 or PFOA. In the former case, T4 RMSD values are low (less than 0.1 A) and conserved in the three replicates. In addition, T4 favourably H-bonds to Thr119, a residue located in the innermost portion of the pocket, thus explaining its higher stability in the binding site. On the other hand, T4 in the TTR:PFOA/T4 complex shows higher fluctuations and a reduced H-bond occupancy with Thr119, due to the increased distance between the ligand and the residue (Figure S11), meaning a lower stability in the pocket. The observation of imbalance between the higher affinity of PFOA for site B and the lower stability of T4 in pocket A, as well as the worse affinity of cC6O4 and the better stability of T4, could be explained by negative cooperativity, (Ferguson et al. 1975; Neumann et al. 2001; Tomar et al. 2012; Haupt et al. 2014).

Finally, additional unbinding simulations were carried out on both complexes to rank and classify the ligands according to their residence time in the binding site. Simulations resulted in an unbinding time of 10.8 (± 1.82) ns for PFOA, and 8.5 (± 1.69) ns for cC6O4. Due to the proportional correlation between the residence time and the binding affinity, the lower values for cC6O4 may assess a lower affinity for TTR compared to PFOA.

TBG has an even lower serum concentration than TTR (1.1–2.1 mg/dL), but a much higher affinity for T4 and T3 (Ka equal to 1010 and 2·108 (M−1), respectively (Refetoff 2000)). It presents a single and quite large binding site, characterised by the presence of two arginines. Arg378 forms a salt bridge with the substrate carboxylic moiety, whilst Arg381 forms a π–cation interaction with the substrate aromatic ring (Ren et al. 2016; Figure S12). Even if no crystallographic information is available about PFAS interacting with TBG, binding experiments showed a distinct drop in affinity after performing site mutagenesis on Arg378 or Arg381. We, thus, guided the docking pose selection according to the capability of PFAS carboxylic, sulfonic or hydroxyl group to H-bond to one of the arginine residues.

Molecular docking returned results quite similar to those reported for TTR, with a general poor affinity and no significant variations, apart for the longer molecules, which did not retrieve any pose (Figure S13). In particular, the lowest values were obtained for the R stereoisomer of ADONA and PFUnA (− 6.7 kcal/mol and − 6.6 kcal/mol, respectively). Medium-length molecules gave slightly lower affinities, with PFOS giving the minimum for the sulfonic acid category (− 6.1 kcal/mol). Less favourable values were obtained also for GenX (-5.0 kcal/mol) and for cC6O4, giving a minimum of − 5.7 kcal/mol, but with only a handful of poses. The docking pose of PFOA and cC6O4 in TBG is reported in Fig. S14.

However, the few available experiments showed a different trend, with the longest molecules, i.e. PFODA, telomers, being the only with a significant affinity towards TBG (Ren et al. 2016).

MD simulations of the TBG:PFOA complex showed conserved interactions between the carboxyl group and only one of the two arginines of interest, Arg381, establishing additional H-bonds with Ser23, Ser24 and Lys270, whilst the hydrophobic end of PFOA faced towards Arg378. The docking pose appeared to be conserved over the three dynamics, with main fluctuations related to the fluorinated carbon chain probably due to its reduced size and, thus, inadequate filling of T4 pocket. In the case of the TBG:cC6O4 complex, two out of three replicas reported a similar H-bond pattern to that shown for PFOA. However, the cC6O4 pose appeared to be less conserved since the molecule adopted a different orientation in the third replica, with the carboxylic group exposed towards both Arg378 and Arg381 and the terminal per-fluorinated carbon directed on the other side of the binding site. Analogously to PFOA, cC6O4 showed a discrete degree of movement due to its small and compact size not completely filling the binding pocket (Figures S15, S16).

PFAS are overall good ligands for FABPs (Cheng et al. 2021; Fenton et al. 2021; Crisalli et al. 2023). Given their similarity to fatty acids, they assume a similar binding mode in FABP binding sites, forming hydrophobic interaction with the many lipophilic residues lining the pocket and electrostatic contacts with conserved arginine in all isoforms. However, the different pocket size in FABPs plays a role in docking results. As a reference, we considered and reported in Fig. S17 the crystallographic pose of palmitic acid in L-FABP (PDB ID 3STK). The different FABP pocket size is shown in Fig. S1d-h. In particular, amongst the FABPs included in the present study (L-FABP, I-FABP, H-FABP, A-FABP, PmP2), L-FABP has the biggest pocket (Smathers and Petersen 2011) and is able to accommodate PFAS with up to 12 carbon atoms (Zhang et al. 2013a). This observation has been confirmed by our simulations, showing the lowest value (− 13.0 kcal/mol), and the highest affinity, for PFDoDA in L-FABP (Figure S18a), and higher values for shorter compounds. PFPeA and PFBS showed the highest score values for carboxylic and sulfonic acids, respectively.

A similar trend can be observed for I-FABP, for which PFDoDA gave, again, the lowest score of − 13.6 kcal/mol (Figure S18b). Very good scores, in this case, were obtained also for PFNA and, in general, for shorter molecules as PFOA, in agreement with the reduced dimension of I-FABP binding site, likely able to properly host even smaller compounds. PFBA and PFHpS showed the highest scores for carboxylic and sulfonic acids, respectively. For both L-FABP and I-FABP, ADONA and GenX presented similar score values, with the minimum ranging in the − 12/− 10 kcal/mol range, and numerous poses, indicative of a favourable interaction. cC6O4, on the other hand, reported higher scores: − 9.4 kcal/mol in L-FABP and − 8.1 kcal/mol in I-FABP, with fewer poses when compared to the other two perfluoroethers. Similar trends have also been also reported by (Cheng et al. 2021).

Even if with higher values, also H-FABP seem to properly host linear carboxylic and sulfonic acids (Figure S18c). PFDoDA presented a minimum score of − 11.8 kcal/mol and was the biggest molecule able to dock inside this pocket, Smaller molecules showed increasingly higher scores, that is, lower affinities. Sulfonic acids demonstrated a downward trend with PFOS reaching − 8.0 kcal/mol docking score, whilst fluorotelomers gave comparable results. Amongst novel PFAS, ADONA reached the best affinity of − 10.2 kcal/mol, GenX around − 8.7 kcal/mol and cC6O4 of − 9.2 kcal/mol. It is also important to note that both isoforms SS and RR of cC6O4 gave a single docking pose, whilst other novel PFAS returned multiple poses.

A distinguishing case in our study is represented by A-FABP (Figure S18d). Being the protein pocket smaller with respect to other isoforms, many carboxylic acids were not processed by the docking algorithm, or returned a very few poses, as previously mentioned, indicative of unfavourable binding (Smathers and Petersen 2011). A downward trend is clearly observable for sulfonic acids and hydroxyl compounds, whilst novel PFAS gave good affinities, with cC6O4 showing the lowest minimum (− 9.4 kcal/mol) and numerous docking poses. Indeed, A-FABP is the only target in our study to show a higher affinity for cC6O4 than for other PFAS.

On the other hand, PmP2 is much more in line with L-FABP, I-FABP and H-FABP (Figure S18e). Linear PFAS bind with ease in these pockets and have a docking pose similar to fatty acids, with a minimum energy found for PFNA equal to − 11.4 kcal/mol. For GenX, no docking pose has been identified for this target, whilst ADONA and cC6O4 showed minimum values about − 9 kcal/mol (− 9.4 and − 9.2 kcal/mol, respectively).

No poses were found for PFOO and PFTeCO in any of the analysed FABPs. Again, this can be reasonably attributed to the impossibility of forming effective interactions with the positively charged arginines lining the binding site.

MD simulations were run for PFOA and cC6O4 complexed to L-FABP, being, by far, the most studied FABP, at least concerning PFAS binding (Zhang et al. 2013a; Sheng et al. 2016, 2018; Cheng and Ng 2018; Gao et al. 2019). In particular, similarly to TTR, we run two 200 ns-long replicas for both the L-FABP:PFOA and L-FABP:cC6O4 complexes. The RMSD calculated for the protein shows the conformation is stable in both cases (Figure S19). Differently, whilst the RMSD for PFOA, with fitting on the pocket residues, shows a stable behaviour after a conformational change at 70 ns in both replicas, the RMSD calculated for cC6O4 shows larger variability, suggesting a less stable binding (Figure S20).

As the major excretion pathway of PFAS is through urine (Ng and Hungerbuehler, 2015; Fu et al. 2016; Göckener et al. 2020; Fenton et al. 2021; He et al. 2023), whilst secondary pathways are more gender specific (Zhang et al. 2013b), we focused on the possible interaction of PFAS with kidney transporters, i.e. OATs, having uric acid as main physiological substrate. OATs are expressed on both sides of epithelial kidney cells, with similar targets but opposite function: OAT1 and OAT3 (amongst others) mediate chemical excretion from the blood into the urine, whilst OAT4 and URAT1 reuptake the same ligands, from the urine back into the bloodstream (Fig. 2).

Cellular localisation of uric acid transporters

These antiporters display multiple conformations, each having a specific role in cell ligand homeostasis. In particular, the outward form corresponds to the transporter (on either cell membrane, that is, apical or basolateral) open towards the extracellular side, whilst the inward state presents the transporter open towards the cytosol. According to the specific transporter function and localization, both outward and inward conformations can indicate either PFAS excretion or reabsorption. We, thus, considered both states in our simulations. Whilst an occluded conformation also exists, we did not include it in this paper.

OAT1 show a marked difference between outward and inward states. In the outward one, a minimum of − 13.2 kcal/mol is found for PFDA, and less favourable interactions for compound belonging to other classes than carboxylic acids (Figure S21a1,a2). Novel PFAS show affinities comparable to that of sulfonic acids (minimum of − 10.3 kcal/mol for ADONA). cC6O4 offers a few poses, but still with good affinity. In the inward conformation, instead, the score distribution is much flatter, with a downward trend centred on PFHxA, and novel PFAS showing lower affinities up to − 4.2 kcal/mol for cC6O4 RS.

OAT3 (Truong et al. 2008; Li et al. 2019; Janaszkiewicz et al. 2023) struggles to accept molecules bigger than PFNA in the outward conformation, but still shows good affinities (maximum affinity of − 10.2 kcal/mol for PFPeA; Figure S21b1,b2). Perfluoroethers show a similar trend as for OAT1, but with slightly less affinity for this target, particularly cC6O4, which shows a minimum of − 7.0 kcal/mol for the SS stereoisomer. The inward conformation, on the other hand, shows higher affinities across the entire PFAS series, in particular for carboxylic ethers. Whilst ADONA R shows a minimum for the corresponding group (− 10.7 kcal/mol), GenX R and cC6O4 SS have lower but still comparable affinity (− 9.6 and − 9.4 kcal/mol, respectively). This leads to the hypothesis that whilst both OAT1 and OAT3 excrete PFAS with similar lengths and functional groups, OAT1 releases them in the cytosol with much more facility.

The first reuptake transporter we analysed is OAT4. In the outward conformation, a trend similar to that of OAT1 and OAT3 can be observed for all PFAS classes, even with a downward trend for both carboxylic and sulfonic acids (PFOA shows a minimum of − 11.3 kcal/mol). Amongst novel PFAS, as usual, cC6O4 shows less affinity than the others and than legacy PFAS (RS, − 7.6 kcal/mol). This point is made more remarkable by the overall low scores of the other stereoisomers (Figure S21 c1, c2). In the inward state of OAT4, instead, more poses are found for all classes of PFAS, and the results of PFOA (− 12.0 kcal/mol) are comparable to those of novel PFAS as, specifically, GenX (− 10.2 kcal/mol for S stereoisomer).

In the outward conformation, URAT1 has much more sparse results, with a flatter downward trend for carboxylic acids. Carboxylic molecules having more than six carbons display only a couple of poses (minimum of − 10.9 kcal/mol for PFDA), whilst sulfonic ones of similar length are much more numerous (Figure S21d1,d2). GenX and cC6O4 display significantly lower affinity for this conformation when compared to all other PFAS. Lastly, URAT1 showed a U-shaped trend in the inward facing conformation, centred on PFHpA (− 11.5 kcal/mol; Figure S21d2). Novel PFAS showed again lower affinities, ranging from − 8.8 to − 3.0 kcal/mol (GenX S and cC6O4 SR, respectively).

Our results are in line with those reported by (Bruno et al. 2022), who found a remarkably lack of interaction for cC6O4 with OATs, which translates to easier excretion through urine for this PFAS, compared to known ones like PFOA and PFOS, which also gave some of the highest affinities for these targets. Indeed, PFOA and PFOS, known to last for years in the human body (Olsen et al. 2007; Zhang et al. 2013b; Fu et al. 2016; Xu et al. 2020; Fenton et al. 2021; Fustinoni et al. 2023), give lower docking scores, and therefore better interactions with these targets. On the contrary, smaller, bulkier molecules like cC6O4 are in striking contrast and returned less favourable interactions with OAT1 and OAT3. In Figures S22 and S23, we report some docking poses for PFOA and cC6O4 in OAT1 and URAT1, respectively.