Interactions between Environmental Contaminants and Gastrointestinal Parasites: Novel Insights from an Integrative Approach in a Marine Predator

Environmental contaminants and parasites are ubiquitous stressors that can affect animal physiology and derive from similar dietary sources (co-exposure). To unravel their interactions in wildlife, it is thus essential to quantify their concurring drivers. Here, the relationship between blood contaminant residues (11 trace elements and 17 perfluoroalkyl substances) and nonlethally quantified gastrointestinal parasite loads was tested while accounting for intrinsic (sex, age, and mass) and extrinsic factors (trophic ecology inferred from stable isotope analyses and biologging) in European shags Phalacrocorax aristotelis. Shags had high mercury (range 0.65–3.21 μg g–1 wet weight, ww) and extremely high perfluorooctanoic acid (PFOA) and perfluorononanoic acid (PFNA) residues (3.46–53 and 4.48–44 ng g–1 ww, respectively). Males had higher concentrations of arsenic, mercury, PFOA, and PFNA than females, while the opposite was true for selenium, perfluorododecanoic acid (PFDoA), and perfluooctane sulfonic acid (PFOS). Individual parasite loads (Contracaecum rudolphii) were higher in males than in females. Females targeted pelagic-feeding prey, while males relied on both pelagic- and benthic-feeding organisms. Parasite loads were not related to trophic ecology in either sex, suggesting no substantial dietary co-exposure with contaminants. In females, parasite loads increased strongly with decreasing selenium:mercury molar ratios. Females may be more susceptible to the interactive effects of contaminants and parasites on physiology, with potential fitness consequences.


INTRODUCTION
Wild animals are concurrently exposed to a suite of natural and anthropogenic stressors that can have additive, synergistic, or antagonistic effects on individuals, with potential consequences on populations. Parasites, environmental contaminants, and their interactions are increasingly a focus of attention because of their ubiquity and their capacity to exert physiological and fitness effects. 1,2 Parasites impose direct costs to their hosts, by competing for nutrients and triggering expensive immune, stress, and behavioral responses, with lethal or sublethal consequences. 3,4 Environmental contaminants such as some essential and nonessential trace elements (TEs) and perfluoroalkyl substances (PFAS) have immunosuppressive and endocrine-disrupting properties and impose energetic and oxidative costs. 5,6 There is some evidence of enhanced susceptibility to infection and parasite pathogenicity linked to TE and PFAS exposure in multiple taxa. 1,7 Yet, our understanding of these interactive effects is still extremely limited, 3 making it difficult to predict the long-term consequences of exposure to multiple stressors for both individuals and populations of wild animals.
Contaminant and parasite interactions and their consequences could be particularly severe for apex predators, especially those with diverse and/or flexible foraging strategies, which can cause them to accumulate a wide range of contaminants and trophically derived parasites. 8,9 Contaminant and parasite burdens in these animals are governed by intrinsic and extrinsic factors that drive variation within and between populations and species. This is partly due to the way individuals accumulate and respond to these stressors, which can be linked to intrinsic factors, such as sex, age, reproductive status, and body mass. Specifically, females have additional elimination routes for contaminants compared to males, through egg, placental, and milk transfer. 10 Males can have a higher susceptibility to parasites than females, e.g., through testosterone-induced immunosuppression. 11 Age can drive differences in tissue concentrations of TEs such as mercury (Hg) and cadmium (Cd) in apex predators, 10 but little is known about age-related trends in tissue PFAS residues. In addition, food ingestion is a significant exposure route to contaminants and parasites in these animals. Diet, feeding habitat, and trophic level are thus key extrinsic drivers of between-individual variation in both contaminant burdens 12,13 and loads of trophically derived parasites. 14 If prey items have high contaminant burdens and high numbers of infective stages of parasites, co-exposure to the two stressors can happen, 15,16 which complicates the detection and understanding of their potential interactive effects on host physiology and fitness.
The marine environment is the final, long-term sink for several TEs and PFAS. 8,10 Anthropogenically modified coastal environments receive a plethora of contaminants from agricultural, industrial, and urban sources that can accumulate in sediments. 10,17 Hg and some PFAS (e.g., perfluooctane sulfonic acid, PFOS) biomagnify strongly up marine food webs, with cumulative quantities transferred at each trophic step. 8,18 This implies a risk of high, chronic dietary exposure in coastal marine predators that remains poorly studied. However, advances in biologging 19 and analytical techniques such as compound-specific stable nitrogen isotope analyses on amino acids (δ 15 N AA ) 20 can now provide detailed, quantitative data on predator feeding strategies. Combining these methods, while taking into account sex and other intrinsic traits, is an emerging and powerful approach to quantify co-exposure to contaminants and parasites and better understand their combined effects on individuals and subsequently populations.
Here, we focus on a population of European shags Phalacrocorax aristotelis (hereafter shags) breeding in the Firth of Forth, a Scottish estuary with over 200 years of contaminant inputs that originate from agrochemical and petrochemical activities and domestic discharges from large urban agglomerations. 21 Studies on marine predators in the region have shown high and slowly declining concentrations of Hg and organochlorine pollutants, 22−24 while exposure to other TEs and PFAS is unknown. The shag is a long-lived, predominantly piscivorous, and pursuit-diving seabird that feeds primarily on the seabed in inshore waters 25−27 and is thus exposed to sediment-derived contamination. 24 At the focal study site, shags are of known sex, age, and breeding history and can be monitored nondestructively for the presence of gastrointestinal parasites through endoscopy. 28 Shags have variable burdens of the nematode Contracaecum rudolphii that negatively impacts adults and nestlings. 4,29 Shags become infected via their fish diet, but drivers of betweenindividual differences in parasite loads are poorly understood. 30 Sex differences in parasite loads and their fitness effects have been reported previously, 29 and sexual size dimorphism is likely to influence shag foraging strategies. 25,31 Therefore, this is a unique system to address critical gaps in the understanding of how intrinsic and extrinsic factors determine individual contaminant and parasite loads and their interactions.
Our aim was to quantify the interactions between contaminant concentrations (11 essential and nonessential TEs and 17 PFAS) and C. rudolphii loads in breeding shags, while accounting for intrinsic (sex, age, and body mass) and extrinsic factors (feeding habitat and trophic level). We characterized sex-and individual-specific feeding ecology using bulk stable carbon (δ 13 C Bulk ) and nitrogen (δ 15 N Bulk ) isotopes and δ 15 N AA analyses supported by biologging. Individuals with a higher trophic level were expected to have higher residues of the biomagnifying Hg and PFOS. 8,18 We expected C. rudolphii loads to increase with increasing Hg and PFAS residues in shags, as a result of immunotoxicity and/or co-exposure. We also considered the alternative hypothesis of lower contaminant residues in the most parasitized individuals, since parasites can act as sinks of contaminants and modify their uptake and metabolism in the hosts, especially for nonessential TEs. 2

MATERIAL AND METHODS
2.1. Study Site and Sampling Procedure. Fieldwork was carried out on the population of shags breeding at the Isle of May National Natural Reserve (56°11′N, 2°33′W), Scotland. This population has been closely monitored over three decades: most birds are either individually marked as chicks and therefore of known age, or assumed to be three years old at first capture as adults, and sex is established from vocalizations. 32 Shags are sexually dimorphic, with females being 15% smaller than males. Forty-four breeding adult shags (26 males, 18 females) were captured during early chick rearing in 2018, with one or both partners sampled per nest. Adults were captured with a crook on the end of a pole or by hand in the early morning before the onset of foraging trips (between 04:00 and 07:00). A blood sample was taken from the brachial vein with heparinized syringes within 3 min of capture. Blood was kept on ice during field procedures before being centrifuged for 10 min at 8000 rpm within 4 h of sampling. Plasma and red blood cells were stored separately in microtubes and frozen at −80°C. After blood sampling, we used endoscopy to quantify the load of gastrointestinal worms in a nondestructive manner in all birds. An endoscope (103 cm long, 9 mm diameter) was gently inserted into the esophagus and stomach, and worms were counted visually using the endoscope's camera images (details in Burthe et al. 28 and Hicks et al. 29 ). All videos were recorded for further count verifications. Endoscopy lasted less than 5 min, and then birds were weighed to the nearest 10 g with a Pesola balance. A subsample of endoscoped birds was fitted with data loggers (see below) before release at the nest. Previous work based on dissections of seven adults and two chicks identified C. rudolphii as being the only nematode species present, 28 and this is in agreement with other studies of shags. 30,33 However, we cannot exclude the possibility that other parasite species are less commonly present in the birds.
2.2. TE and PFAS Determination. Eleven TEs (arsenic, As; Cd; cobalt, Co; chromium, Cr; Cu; Hg; molybdenum, Mo; nickel, Ni; Pb; Se; and Zn) and 17 PFAS (13 perfluoroalkyl carboxylic acids, PFCAs, with carbon chain lengths C 4 −C 18 : perfluorobutanoic acid (PFBA) (C 4 ), perfluoropentanoic acid (PFPeA) (C 5 ), perfluorohexanoic acid (PFHxA) (C 6 ), perfluoroheptanoic acid (PFHpA) (C 7 ), perfluorooctanoic acid (PFOA) (C 8 ), perfluorononanoic acid (PFNA) (C 9 ), perfluorodecanoic acid (PFDA) (C 10 ), perfluoro-n-undecanoic acid (PFUdA) (C 11 ), perfluorododecanoic acid (PFDoA) (C 12 ), perfluorotridecanoic acid (PFTrDA) (C 13 ), perfluorotetradecanoic acid (PFTeDA) (C 14 ), perfluorohexadecanoic acid (PFHxDA) (C 16 ), and perfluoro-n-octadecanoic acid (PFODA) (C 18 ); four perfluoroalkane sulfonic acids, PFSAs, with carbon chain length C 4 −C 10 : perfluorobutane sulfonate Environmental Science & Technology pubs.acs.org/est Article (PFBS) (C 4 ), perfluorohexane-1-sulfonic acid (PFHxS) (C 6 ), PFOS (C 8 ), and perfluorodecane sulfonate (PFDS) (C 10 )) were measured at the UK Centre for Ecology and Hydrology laboratory in red blood cells and plasma, respectively. Analytical methods, quality assurance (QA) and quality control (QC), limits of detection and quantification (LOQ), quantification frequency (QF), and concentrations of all individual PFAS are presented in the Supporting Information (SI) (Section S1 and Table S1). The results are given as mean ± standard deviation (SD) in μg g −1 and in ng g −1 for TEs and PFAS, respectively, both on a wet weight (ww) basis. "Blood" within the text refers either to red blood cells for TE and isotopic values or to plasma for PFAS residues. Total Hg concentrations approximate methyl-Hg concentrations in seabird blood. 34 Hg and other TEs partition preferentially in red blood cells, where the concentrations typically represent recent exposure, but time integration can differ between TEs and physiological conditions. 10 Se has a well-known protective role against Hg toxicity. 35,36 Se:Hg molar ratios of at least one in apex predators' tissues are thought to provide protection against Hg toxicity, but the protective ratio may depend on species and other factors that are currently poorly studied. 35 Nevertheless, Se:Hg molar ratios can provide a more accurate indicator of Hg toxicity than Hg concentrations alone. 37 In vertebrates' blood, PFAS are mainly found in association with plasma proteins. 38 Plasma PFAS residues are thought to represent a large portion of the total body burden in birds, with the liver being another important storage tissue, 39,40 and weak excretion in feathers. 41 Blood PFOS half-life was estimated to be ∼200 days in an experimental study in chickens, but can vary greatly depending on species, and the kinetics of other PFAS are unknown. 42 Thus, plasma PFAS concentrations may reflect a long, unknown temporal period of exposure. Here, birds were sampled before the onset of daily feeding trips. Therefore, potentially confounding effects of very recent dietary intake on blood TE and PFAS concentrations were negligible.
2.3. Characterizing Feeding Ecology. 2.3.1. Biologging and Sediment Data. Twenty-four of the 44 endoscoped birds were equipped with Axytrek data loggers (TechnoSmArt, Italy) and 4 birds simultaneously with IgotU GT-120 GPS (Mobile Action, Taiwan) and D3GT data loggers (Little Leonardo Co., Tokyo, Japan). Loggers were attached on the midline of the mid back using Tesa tape. All 28 equipped birds were recaptured and loggers retrieved after at least 72 h of deployment. Details on data logger specifications, determination of foraging dive locations, and kernel density estimations are given in the SI (Section S2). Dive locations were coupled with seabed sediment data collated from the British Geological Survey through the Marine Digimap download service 43 on QGIS 44 to quantify the percentage of dives carried out in specific habitats (Section S2).
2.3.2. Isotopic Analyses. All samples (N = 44) were analyzed for δ 13 C Bulk and δ 15 N Bulk values. A sex-balanced subset of samples (N = 20) was selected for δ 15 N AA analyses. Isotopic values were determined in freeze-dried red blood cells, which represent dietary information over approximately 2−3 weeks preceding sampling 45 and typically do not require lipid extraction before analyses. Stable isotope analyses were carried out at the Liverpool Isotopes for Environmental Research Laboratory, University of Liverpool, and are reported in δ notation in per mil (‰) relative to Vienna Pee Dee Belemnite for δ 13 C Bulk , and atmospheric N 2 for δ 15 N Bulk and δ 15 N AA . Details of sample preparation, method precision and accuracy, and δ 15 N values of all identified amino acids (alanine, aspartic acid, glutamic acid, glycine, leucine, phenylalanine, and valine) are provided in the SI (Section S3 and Table S2).

Estimation of Trophic Position and Feeding
Habitat. δ 13 C Bulk in tissues of marine predators can trace the origin of their food resources, and δ 15 N Bulk can be used to estimate their trophic level. 46 However, δ 15 N Bulk values of consumers are sensitive to δ 15 N values at the base of their food web or baseline, 20 which can hinder accurate trophic level determination. This can be overcome by using δ 15 N AA analyses. Individual amino acids differ in their δ 15 N values, because they are subject to metabolic pathways with different trophic fractionations (reviewed in McMahon and Newsome 47 ). We used phenylalanine, a source amino acid that undergoes minimal trophic fractionation within organisms and food webs, to track the δ 15 N at the base of the food web (δ 15 N Phe 47 ). We subtracted δ 15 N Phe from δ 15 N Bulk values, providing baselinecorrected δ 15 N Bulk (corr-δ 15 N Bulk ) values using a similar approach as in Dolgova et al. 48 (Section S3 and Table S3, SI). The values of δ 13 C Bulk were used to characterize feeding habitat 45 in combination with biologging and sediment data. δ 13 C Bulk and corr-δ 15 N Bulk values were used in statistical models as a proxy of the feeding habitat and relative trophic position, respectively.
2.4. Data Analyses. Data exploration, plotting, and statistical analyses were carried out in R version 3.5.2. 49 2.4.1. Intercontaminant Relationships. A principal component analysis (PCA) was performed on untransformed, scaled TEs and PFAS with QF above 65% (Section S1, SI). For those, concentrations below the LOQ were substituted with randomly generated values between zero and the LOQ.
2.4.2. Characterizing Feeding Ecology. The effect of sex on bulk δ 13 C and δ 15 N and δ 15 N AA values and the logittransformed percentage of dives carried out in rocky seabeds, as well as correlations between isotopic values, were tested through linear models, with model fit being checked by residual analyses. 50 δ 13 C Bulk and δ 15 N Bulk values were strongly correlated ( Figure S1), while there was no correlation between δ 13 C Bulk and corr-δ 15 N Bulk (see Section 3), allowing us to use them as noncollinear covariates in further multifactorial analyses. For those, separate models were fitted for the two sexes (N = 26 males; N = 18 females), given the nonindependence of nest pairs and the strong sexual differences in feeding ecology, contaminant residues, and parasite loads (tested using linear models, see Section 3).

Intrinsic and Extrinsic Drivers of TE and PFAS Residues in Males and
Females. The effect of the explanatory variables body mass, age, feeding habitat (δ 13 C Bulk ), and relative trophic position (corr-δ 15 N Bulk ) was tested on As, Hg, Se, Se:Hg molar ratio, PFOS, PFHxS, PFDA, PFOA, and PFUdA residues. These response variables were selected because of their toxicological relevance or/and their high variable loadings on principal components of the PCA (Section S1 and Table S4). Their concentrations were not transformed because their distributions were close to normal. Model specification was validated via residual analysis of maximal models. 50 Explanatory variables were not significantly collinear (variance inflation factors <3) and were standardized (mean = 0, SD = 1) to facilitate comparison of effect sizes. We adopted an information-theoretic approach through the use of Akaike's information criterion corrected (AICc) for small sample sizes Environmental Science & Technology pubs.acs.org/est Article (R package MuMIn). 51 We compared a list of biologically meaningful candidate models, with the maximal model being contaminant = body mass + age + δ 13 C Bulk + corr-δ 15 N Bulk .
Interactions among the explanatory variables were not included to avoid overfitting and because they were not considered biologically essential. For each specific model, we calculated the AICc, the difference between AICc of the specific model and the best model (ΔAICc), and the AICc weight (normalized weight of evidence in favor of the specific model, relative to the whole set of candidates). 52 If all models had a higher AICc than the null model, the effects of all explanatory variables were considered nonsignificant. When multiple models performed better than the null model and had ΔAICc < 4, we used model averaging to make an inference. This produced averaged parameter estimates (β) of the predictor variables included in those models, weighted using AICc weights. A weight of evidence for each predictor included in model averaging was also produced. Only predictors with weight of evidence >0.90 were considered to have a clear statistical effect. The averaged β ± standard error (SE), 95% confidence interval, and the weight of evidence are reported for each predictor.

Intrinsic and Extrinsic Drivers of C. rudolphii Loads and Their Correlates with TE and PFAS Residues in Males and
Females. The effect of selected contaminant residues and other explanatory variables (body mass, age, feeding habitat (δ 13 C Bulk ), and relative trophic position (corr-δ 15 N Bulk )) was tested on parasite loads. Only Se:Hg molar ratio, PFOA, and PFOS residues were selected among contaminants to minimize collinearity, as inferred from the PCA results (Section S1 and Figure S2, SI). Continuous predictor variables were scaled. Again, we adopted an information-theoretic approach through the use of AICc and model averaging to establish which predictor(s) best explained the variation in parasite loads. Given the small sample size, only models simultaneously including three predictors or fewer were considered.

RESULTS AND DISCUSSION
3.1. Blood TE and PFAS in Isle of May Shags. Among essential elements, As, Cu, Se, and Zn were quantified in all individuals and Cr in 80% of them (Table 1). Mo (range The reference is females. Values in bold are significant differences. b Only contaminants with a quantification frequency (QF) over 65% are presented. Abbreviations: perfluoroalkyl carboxylic acids, PFCAs; perfluoroalkyl sulfonic acids, PFSAs; limit of quantification, LOQ.
Environmental Science & Technology pubs.acs.org/est Article <LOQ-0.02 μg g −1 ww) and Ni (<LOQ-0.148 μg g −1 ww) were quantified in >50% of the individuals, while Co was not detected in any individual (Table S1). Among nonessential elements, only Hg had a QF of 100% (Table 1), while Pb (<LOQ-1.210 μg g −1 ww) was quantified in 11 individuals and Cd in one male (0.005 μg g −1 ww, Table S1). Concentrations of As, Cu, Cr, and Zn were in the same ranges observed in seabirds worldwide. 53 Hg concentrations were high and comparable to those of double-crested cormorant P. auritus from the polluted North American Great Lakes. 54 Shag Se concentrations and Se:Hg molar ratios were relatively low for seabirds 55 and similar to those of Sternidae species. 53 Among the 17 targeted PFAS, 10 had quantifiable concentrations in at least 65% of individuals. All individuals had quantifiable concentrations of PFHxS, PFOS, PFDA, and PFUdA (Table  1). PFOS had by far the largest contribution of all PFAS in both males and females (mean contribution of 50 and 58%, respectively, Figure S3 and Table S5). PFOS concentrations were particularly high in females (Table 1). Among PFCAs, the main contributor to the total PFAS burden was PFUdA in males and PFDoA in females (9.7 and 12.8%, respectively, Figure S3 and Table S5). Short-chain PFAS (C 4 −C 7 ) were quantified in a few individuals only, with the exception of PFHxS and PFHpA (100 and 44% of individuals, respectively; Table S1). Longer PFCAs (C 14 −C 18 ) were quantified in <10% individuals (Table S1). This is one of the first reports on PFAS exposure in birds from the United Kingdom (see also Leat et al. 56 ), revealing very high concentrations, in particular, for PFOA and long-chained PFCAs, even when compared to avian species from highly contaminated sites in North America and Europe 13,39,57 (see Section S4 (SI) for a detailed comparison with previous avian studies). Our results suggest that Firth of Forth waters, sediments, and organisms are highly contaminated with PFAS and Hg. This calls for an urgent determination of the potential toxic effects of these contaminants in the shags and in other apex predators relying on the waters of this Scottish estuary.

Characterizing Individual Feeding Ecology.
Shags from the Isle of May foraged on a variety of rocky and sandy areas, usually within 12 km from the colony (Figure 1a−c), in agreement with previous findings. 26 However, the biologging and isotopic results show for the first time clear sexual segregation in feeding habitat. Males targeted sites north of the Isle of May, while females traveled further and westward, diving closer to the mainland (Figure 1a−c). Males dived mostly in rocky habitats (mean ± SD, 69 ± 29% of all dives, Table S6). In contrast, females exploited mainly two types of sandy seabeds, the "gravelly sand" and "sand" types (41 ± 36   Table S6). Sexual segregation was clear also from the isotopic results. δ 13 C Bulk values were different between sexes, with the blood of females being depleted in 13 C compared to males ( Figure 1d and Table 1). δ 13 C Bulk values suggest that females fed on pelagic-influenced prey, whereas males fed on a mixture of pelagic-and benthic-influenced prey. At the population level, shags feed mainly on the lesser sandeel Ammodytes marinus (hereafter sandeel), 26,27 with smaller contributions of benthivorous fish of the families Pholidae, Callionymidae, and Gadidae, found in rocky habitats. 27,58 Recently, there has been a marked diversification in shag diet at the population level, with sandeels representing smaller proportions of total prey. 58 Sandeels live part of their life buried in the sediment, but are midtrophic planktivorous schooling fish feeding on pelagic copepods. 59 Sandeels are thus likely depleted in 13 C when compared to benthic-feeding fish found in and around rocky substrates. 46 The 13 C-depleted and small-ranging δ 13 C Bulk values in female shags ( Figure 1d and Table 1), together with a high proportion of dives in sandy areas, strongly suggest that females mainly fed on sandeels. This is further supported by the similarity of δ 13 C Bulk values of female shags and those of sandeels from the North Sea (range −20 to −18‰). 60 In contrast, the large range of δ 13 C Bulk values in male shags (Figure 1d and Table 1) indicates that individual birds fed on various proportions of pelagic-influenced ( 13 Cdepleted) and benthic-influenced organisms ( 13 C-enriched). 46 In addition, there were sex-related differences in δ 15 N Phe values ( Figure 1e and Table 1), indicating that the prey of males and females belonged to different food webs. Specifically, the δ 15 N Phe values of males were similar to the δ 15 N Bulk values of microphytobenthos, a typical primary producer at the base of benthic food webs in coastal environments. 61 Females had lower δ 15 N Bulk values than males (Figure 1e and Table 1). However, there was a clear effect of δ 15 N Phe on δ 15 N Bulk values (Figure 1e), indicating that sexual differences in δ 15 N Bulk values were mainly driven by variation in the δ 15 N baseline and not in the trophic position. Indeed, corr-δ 15 N Bulk values were similar between the sexes (Figure 1d and Table 1), demonstrating that males and females generally fed at the same trophic level despite differences in dietary sources.

Intrinsic and Extrinsic
Drivers of TE and PFAS Residues. 3.3.1. Sex and Feeding Ecology. Sexual differences in contaminant concentrations are not consistent among seabird species, with feeding ecology usually being a key explanatory factor. 12,40 Here, essential TEs such as Cr, Cu, and Zn were similar between sexes and had low between-individual variation ( Table 1), in accordance with homeostatic regulation. 10 In contrast, males had higher As and Hg concentrations than females, while the opposite was true for Se and for the Se:Hg molar ratio (Table 1). In contaminated coastal environments, sediments are long-term sinks for several TEs. 10 Benthic species are known to accumulate TEs more than pelagic species, especially with respect to Hg, because of strong methyl-Hg production in sediments. 62 This is consistent with the higher As and Hg residues in benthic-influenced males than in pelagic-influenced females and with the positive correlation between male δ 13 C Bulk values and Hg residues ( Figure 2). As predicted, Hg concentrations increased with relative trophic position (corr-δ 15 N Bulk ) in males (linear model, estimate ± SE, 0.26 ± 0.08, Figure S4a and Table S7). This highlights Hg biomagnification in the benthic food webs of the Firth of Forth. In contrast, the Se pattern depicted here was unexpected. Contrary to findings in several marine vertebrates, 35,63 blood Hg and Se concentrations were negatively related, irrespective of sex, as shown by the PCA results ( Figure S2). In addition, and contrary to Hg (linear model, 0.21 ± 0.11), δ 13 C Bulk values were negatively related to Se concentrations in males (γ inverse GLM, 0.06 ± 0.03, Figure  2). High Se concentrations were thus associated with pelagicinfluenced diet, as also confirmed by high Se concentrations in females. Se content in shag prey should be investigated to confirm this habitat pattern in Se exposure. Hg and Se residues were not influenced by feeding habitat and trophic level among females (Table S7), likely because of their low betweenindividual variation in feeding ecology or because of other unknown factors. 64 In the closely related double-crested cormorant P. auritus, females have higher Hg demethylation capacities than males, a potential adaptation to reduce Hg transfer to eggs, 65 which could also imply that females have a stronger Se assimilation capacity from their diet.
Despite strong between-individual variation in PFAS concentrations, sex-related differences were evident for several compounds. Females had overall higher total PFAS concentrations than males, mainly because of their higher PFOS residues ( Figure S3 and Table 1). Among PFCAs, males had overall higher concentrations of C 8 −C 10 compounds than females, while the opposite was true for C 11 −C 14 PFCAs, in Only the clear statistical correlation between δ 13 C Bulk values and Se residues in male shags is represented (Table S7).
Environmental Science & Technology pubs.acs.org/est Article particular PFDoA ( Figure S3 and Table 1). PFAS, as amphiphilic molecules, have high water solubility, but can also bind to sediments and dissolved organic matter, depending on their carbon chain length. 66 PFAS patterns in both male and female shags are consistent with feeding on a mixture of pelagic-and benthic-feeding species. 67 Benthic invertebrates, such as amphipods and annelids, were previously found to have exceptionally high concentrations of PFAS, in particular PFOA. 7,17,67 Annelids are a minor food source in shags at the Isle of May, 58 but could be critical in their exposure to PFAS. The high PFOS concentrations in females are suggestive of high PFOS content in their prey, i.e., sandeels, which might assimilate PFOS while buried in sandy sediments. Further investigations on PFAS concentrations in fish species from the Firth of Forth should confirm this. Sex differences in blood chemistry, such as a higher level of total plasma proteins (and thus PFAS binding potential) in females, 68 may also be an explanatory factor of their higher PFAS residues. Contrary to predictions, PFOS and long-chain PFAS did not increase with the trophic level in shags (Table S7). In some fish species, direct uptake of PFAS from the water is significant and can even exceed dietary exposure. 69,70 This complicates our understanding of PFAS biomagnification dynamics, especially at the interface between benthic and pelagic environments. 67 In addition, the potential temporal uncoupling between PFAS and isotopic values (due to PFAS long half-lives in blood, see Section 2) could imply that diet over a long time period should be considered. Specifically, between-individual differences in diet associated with migratory destinations during the nonbreeding period 71 could play a critical role in PFAS exposure 13 in this partially migratory population.

Age and Body
Mass. In seabirds, blood Hg concentrations are usually higher in adults than in chicks and reach a steady state early in adult life, with no further increases. 12 Age-related trends in other TEs are less clear, but overall there is no support for strong changes of TE blood concentrations in aging adults. 10 The lack of a clear effect of age on TE concentrations in the shags (Table S7), irrespective of sex and despite the large age range (2−20 years), is thus consistent with previous findings and can be explained by the homeostatic balance of essential TEs and efficient detoxification and via deposition mechanisms of nonessential TEs, notably Hg excretion into feathers. 10 On the other hand, PFAS kinetics in wild animals are poorly known, and this is the first study to investigate the effect of adult age on blood PFAS concentrations in an avian species. 42 Previous investigations have shown higher PFAS residues in the tissues of chicks or juveniles with respect to adults, 39,72 but also the opposite 73 or nonsignificant trends. 74 In marine mammals, PFAS residues in internal tissues decrease with adult age, likely as a result of clearance and growth dilution mechanisms. 75 PFAS concentrations were similar among shags of different ages. This suggests that (i) a steady state in PFAS circulating concentrations is reached before breeding maturity in seabirds, similarly to organochlorine pollutants, 76 and that (ii) shags of different ages are continuously exposed to PFAS in this coastal environment. In addition to longitudinal investigations, further studies are needed on the PFAS distribution in internal organs to understand whether PFAS body burdens change with age in birds.
Body mass had no effect on contaminant concentrations in females (Table S7). In contrast, Hg and PFOS residues had, respectively, a negative (linear model, −0.23 ± 0.08, Figure   S4b and Table S7) and positive (31.8 ± 11.7, Figure S5 and Table S7) relationship with body mass. Mass can influence the diving capacities and prey selection of diving birds. 31 Prey species with low Hg and high PFOS burdens targeted specifically by large male shags could thus explain these results.
3.4. Parasite Loads. 3.4.1. Intrinsic and Extrinsic Drivers. Little is known on the drivers of parasite prevalence and loads in free-ranging birds. 3,14,77 As predicted by testosteroneinduced immunosuppression 11 and in agreement with previous results, 29 males had higher C. rudolphii loads than females (Table 1). Feeding habitat (δ 13 C Bulk ) and trophic position (corr-δ 15 N Bulk ) were not related to C. rudolphii loads in either sex ( Figure S6 and Table S8). Therefore, physiological factors related to immunocompetence may outweigh trophic drivers in explaining the parasitic nematode loads in this population. The life cycle of C. rudolphii is complex and poorly known, 78 but adult stages of the parasite are thought to stay attached to the bird host for several weeks to months. 28,30 The carry-over effects of diet on a longer term, especially from migrating grounds, are therefore possible. Biologging and isotopic data do not provide definitive information on the exact species of prey ingested, which could bear different taxa and numbers of infective stages of parasites. To further investigate the trophic drivers of C. rudolphii exposure, the intermediate host species of the parasite and diet composition of male and female shags over both breeding and migrating grounds should be identified. Importantly though, our results indicate that while at their breeding grounds, there is no substantial dietary co-exposure between contaminants and gastrointestinal parasites in shags, as previously found for Hg in other seabirds. 9, 64 3.4.2. Correlates with TE and PFAS Residues. In female shags, there was a strong negative relationship between C. rudolphii loads and Se:Hg molar ratios (γ inverse GLM, 0.03 ± 0.01, Figure 3 and Table S8): females concurrently having high Hg and low Se residues had higher loads of C. rudolphii. Given the absence of substantial co-ingestion between Hg, Se, and C. rudolphii, this association is likely the result of physiological mechanisms and is consistent with the hypothesis of increased parasite infection susceptibility linked to Hg immunotoxicity. 6 In addition, we show here that this relationship could be mitigated by Se. Gastrointestinal parasite loads increased with Hg concentrations in some seabird species, 15,16 while other studies showed no association. 9,79 Contrary to our results, gastrointestinal parasite loads increased with increasing Se concentrations in the glaucous gull Larus hyperboreus, 15 and Se  (Table S8).
Environmental Science & Technology pubs.acs.org/est Article residues were higher in parasitized than in nonparasitized arctic terns Sterna paradisaea. 80 However, those studies did not test the relationship between parasite loads and the Se:Hg molar ratio, which could explain the contrasting results. Se is a wellknown antioxidant and enhances immunocompetence in birds at normal dietary levels, while becoming immunotoxic at higher concentrations. 63 The dose-related relationship of Hg and Se on health parameters can be mutually modified by the two TEs. 35 In agreement with laboratory-based studies in mallards Anas platyrynchos 81 and mice, 37,82 our results suggest that the Se:Hg molar ratio, rather than single elemental effects, should be considered to rule out Hg and Se immunotoxic potential. In addition, our findings highlight that Se might be more important than previously thought in driving the response of wild organisms to parasitic nematode infections. There was no evidence of an association between Se:Hg molar ratios and C. rudolphii loads in males ( Figure 3 and Table S8), despite their higher Hg and lower Se concentrations and higher parasite loads than females. This is suggestive of a different susceptibility to Hg and/or parasites between males and females and of the underestimated complexity of Hg and Se interactions. 81 Correlative and experimental work shows that female shags have (i) increased flight cost with higher parasite loads 83 and (ii) stronger foraging benefits of parasite removal 4 compared to males. This discrepancy in parasite effects between sexes is thought to result from the tighter energetic constraints in females that result from their higher investment in reproduction. 4,83 The higher energetic demands in breeding females could also make them more susceptible to Hg immunotoxicity, especially at younger ages (γ inverse GLM, 0.02 ± 0.01, Table S8). Energetic constraints could be very important in shaping the response of individuals to contaminants and parasites, which opens up a promising avenue for further research. This is the first study to investigate the relationship between PFAS and gastrointestinal parasite loads in a wild animal population. Contrary to our primary hypothesis, C. rudolphii loads did not increase with PFOS residues in either male or female shags ( Figure S6 and Table S8). Instead, the most parasitized females had lower blood PFOS concentrations (γ inverse GLM, 0.02 ± 0.01, Table S8). This supports the alternative hypothesis that C. rudolphii may act as sinks of PFOS, similarly to some TEs. 1,2 This is also consistent with the more parasitized males having lower PFOS residues than females.
In addition to the potential role of C. rudolphii as sinks of PFOS, gastrointestinal parasites could have influenced the blood residues of other TEs and PFAS, by modifying their uptake, distribution, or/and metabolism in the shags. Interestingly, rainbow trouts Oncorhynchus mykiss infested by the parasitic nematode Raphidascaris acus showed a slowed Se accumulation compared to uninfected individuals. 84 Parasitedriven changes in Se accumulation could be an alternative explanation of the relationship between Se:Hg molar ratios and C. rudolphii loads in females and of the low Se concentrations in males. Studies quantifying the effect of gastrointestinal parasites on contaminant kinetics and dynamics in their bird hosts are drastically needed.
By using a multidisciplinary approach, this study makes a major contribution to our understanding of the drivers of contaminant and parasite exposure in wild animals and of the complexity of their interactions. Organochlorine pollutants should also be the object of future investigations, because of their potential interaction with parasites in eliciting negative effects on fitness 85 and their high concentrations in this shag population. 24 Finally, future work on immunological parameters could provide additional insights into parasite load variation and a deeper mechanistic understanding of the sexspecific relationships between contaminant and parasite burdens depicted here.

■ ACKNOWLEDGMENTS
The present research was funded by the European Union through an H2020 Marie-Sklodovska Curie Fellowship (ECOSEA; project number 752714) awarded to A.C. and an Early Career Researcher and Returner fund granted by the University of Liverpool to A.C. and C.d.l.V. UKCEH was funded by the Natural Environment Research Council award number NE/R016429/1 as part of the UK-SCaPE programme delivering National Capability. We thank Scottish Natural Heritage for allowing us to work on the Isle of May and all of the field workers who contributed to data collection for this study, in particular Ruth Dunn and Olivia Hicks. We are grateful to Katsufumi Sato for lending us the Axytrek loggers and to Carlo Catoni from TechnoSmArt for providing excellent technical support. We thank Paul Wigley and Sian Pottenger for providing chicken blood that was used as an internal standard for contaminant analyses. We thank also Ruth Dunn, Alice Trevail, and James Duckworth for R code assistance for spatial and diving analyses and Olivia Hicks for helpful discussions and comments on an early version of the manuscript. We thank the two anonymous reviewers for their very constructive comments.