Life-cycle mediated effects of urbanization on parasite communities in the estuarine fish, Fundulus heteroclitus

This study examined the relationship between urbanization and parasite community structure in the estuarine fish, Fundulus heteroclitus. We measured landscape and physicochemical factors associated with urbanization at 6 sites from 4 collection periods. Concurrently, we quantified the metazoan parasite community in F. heteroclitus collected at those sites, with 105 fish studied per site during the 4 collection periods. Parasite community composition differed among sites. Host size was the most important variable for direct life-cycle parasite assemblages and indirect life-cycle parasites at the individual fish level, while landscape and physicochemical factors determined the structure of indirect life-cycle parasite assemblages at the population scale. Variation in the prevalence and intensity of infection of two indirect life-cycle parasites, Lasiocotus minutus and Glossocercus caribaensis, were the primary parasites that drove differences across sites. Variation in the presence/absence of these indirect life-cycle parasite species was associated with sediment Ni concentrations, patch density, and marsh size. Our data support the hypothesis that urbanization, acting at both landscape and physicochemical scales, can have a significant impact on parasite community structure. This, however, varied by parasite life history: there was little effect of urbanization on the prevalence and intensity of direct life-cycle parasites, but significant variation was detected for indirect life-cycle parasites. This study demonstrates how anthropogenically driven landscape change influences fine-scale population dynamics of parasites.


Introduction
Urbanization significantly alters estuaries and the free-living animal populations within these environments [1,2]. Animal parasites, which are intrinsically connected to their hosts, are also impacted by the consequences of urbanization. Urbanization can directly affect parasites through exposure to an altered physicochemical environment [3,4] or indirectly through changes in host population dynamics [5]. Parasites co-occurring in an individual-the parasite infracommunity [6]-simultaneously use the same host resource, but each parasite species may PLOS  have a different life history strategy which may be differently affected by urbanization. Notably, indirect life-cycle (more than one host species) parasites have life stages that are exposed to the external environment and are dependent upon trophic interactions; whereas direct life-cycle (one host species) parasite transmission largely does not involve trophic interactions. Consequently, urbanization may differentially affect parasite species [7], but data that associates urbanization with parasite communities and life history strategy are sparse (but see [8,9]). Parasite life-cycles can be broadly characterized as either indirect or direct. Indirect lifecycle parasites require different species of hosts to complete their life-cycle and are largely dependent on trophic interactions which occur more frequently in locally stable communities [10]. This essential interaction has been used to demonstrate how parasites can be indicators of free-living diversity [11], and as such, it has been argued that a healthy ecosystem is one that is rich in parasites [12]. In urbanized estuaries, host population dynamics may be less stable, resulting in reduced diversity or less predictable host dynamics, causing a concomitant decrease in the species richness of indirect life-cycle parasites [5]. On the other hand, direct life-cycle parasites require only one species of host to complete their life-cycle and are largely dependent on host population densities, as abundant hosts increase the probability of a parasite encountering a suitable environment to colonize [13,14]. A feature of urbanized estuaries is habitat fragmentation, which can result in increased host densities and a concomitant increase in the species richness of direct life-cycle parasites [15]. Generally, parasite communities are described using terms such as prevalence (the number of hosts infected with a parasite species divided by the number of hosts examined), and intensity of infection (number of individuals of a parasite species in a single host species) [6]. These measures are calculated by quantifying the parasite species in a host individual-the parasite infracommunity [6]-and then aggregating these data to consider all the parasite species-the parasite "component" community-found in a subset of a host species collected from a subset of a sampled abiotic environment [6].
Many of the studies that examine the interaction between host and parasite in urbanized ecosystems have been conducted with human pathogens [16], consider a single host and a single parasite [17], or a single contaminant and landscape factor. Further, linking the regional process of urbanization to local transmission and diversity of parasite communities has yet to reveal consistent patterns [18][19][20]. One potential explanation for the absence of consistent results is that urbanization differentially affects parasite species: the size and direction of effect is contingent upon parasite life history. In this study, we assessed the effect of multiple landscape composition and physicochemical factors on 24 parasite communities of the salt marsh fish, Fundulus heteroclitus. The parasite community in F. heteroclitus includes indirect and direct life-cycle species, and the impact of urbanization on the structure parasite communities should be mediated through different life-cycle strategies. We predicted that urbanized sites would have low diversity and prevalence of indirect life-cycle parasites; but direct life-cycle parasites would increase in diversity and abundance.

Ethics statement
Field collections were conducted under scientific permits issued by the Georgia Department of Natural Resources approved collection under Scientific Collecting Permit number 29-WJH-15-182. Prior to necropsy, fish were maintained briefly in aquaria following animal care protocols approved by The Animal Care and Use Committee at Georgia Southern University (protocol #I14013). Fish euthanasia was conducted in accordance with the 2000 Report of the American Veterinary Medical Association Panel on Euthanasia, approved by The Animal sites representing a range of urbanization and its effects, however, for the landscape composition analyses we used the watershed level as our spatial scale (described below).

Quantifying the landscape and physicochemical environment
We defined the boundaries of the landscape to the lowest level watershed provided by the United States Geological Survey (USGS) Watershed Boundary Dataset (WBD): the hydrological unit code (HUC) 12-digit boundary [24]. We used the 2011 National Land Cover Dataset (NLCD) and National Oceanic and Atmospheric Administration's (NOAA) Coastal Change Analysis Program (C-CAP) for the regional land cover analysis [25,26]. Regional landscape composition and change were measured using ArcGIS (ArcGIS Version 10.3, Environmental Systems Research Institute, Redlands, California). We also calculated a series of regional landscape metrics: patch density, largest patch index, edge density, mean salt marsh size, % salt marsh, % developed land, and total area using FRAGSTATS software [27]. Patch density is an index that expresses the number of continuous habitat types, or patches, per unit area, and is a metric for landscape heterogeneity. Largest patch index is the size of the largest habitat patch relative to the size of the landscape, and is a measure of landscape dominance. Edge density is an index of the amount of habitat edges per unit area, and the contrast between developed land and salt marsh was weighted at 0.5 to evaluate salt marsh habitat fragmentation by urban structures or landscapes. Mean salt marsh size is an index of the size of a salt marsh habitat patch, with larger patches weighted in order to decrease the sensitivity of the mean to small habitat patch sizes. The % salt marsh and developed land is the percentage of the landscape comprised of those categories, and is a metric for habitat size and dominance. Total area refers to the area of the watershed. In addition, we used the National Oceanic and Atmospheric Administration's (NOAA) Coastal Change Analysis Program (C-CAP) to evaluate % landscape change within the HUC-12 digit watersheds for the period of time from 1996-2010.
We collected information on the physicochemical environment from each site at flow and ebb tide in March, June, September 2015 and February 2016. Flow and ebb tide measurements were averaged for each site and time point for analysis. We measured pH, temperature, salinity, and conductivity using a handheld multiparamenter instrument (Yellow Springs Instruments Professional Plus; Yellow Springs Instruments, Yellow Springs, Ohio). In addition, we quantified the trace metal composition, nutrient concentrations, and chlorophyll a concentration during these dates. To extract chlorophyll a, we followed a standard acidification method (EPA 445.0), whereby a volume of 100 ml was filtered using Whatman GF/F filters (47 mm; 0.7 μm nominal pore size). Phytoplankton pigments were then extracted from the filters in 8 ml of 90% acetone at 0˚C for 24 hours. Fluorescence was measured using a Trilogy fluorometer (Turner Designs, San Jose, California). For the nutrient and trace metal quantification, 2 samples from each collection time and site were sent to the University of Georgia Laboratory for Environmental Analysis (Athens, Georgia). These samples were filtered using a microwave-assisted nitric acid method (EPA 3015), and total nitrogen, total phosphorus, and 72 metals were analyzed with inductively coupled plasma mass spectrometry (ELAN 9000; PerkinElmer, Waltham, Massachusetts). In addition to water analyses, we quantified the trace metal composition of sediment by collecting 3-8 cm core samples from randomized locations within each site in September 2015. These samples were sent to the University of Georgia Laboratory for Environmental Analysis (Athens, Georgia), and were prepared using a microwave-assisted nitric acid method (EPA 3052).

Quantifying the parasite community in F. heteroclitus
We used the parasite community of a focal species, the common marsh killifish F. heteroclitus, to assess the impact of human modification of the environment on parasite establishment and richness. The common killifish inhabits coastal waters along the Atlantic coast of North America, typically living in dense shoals that may include thousands of individuals: they are oviparous and spawn from late spring to early fall, depending on latitude [28]. Notably, F. heteroclitus survive large temperature and salinity ranges [29] and appear to tolerate a wide range of contaminants (Burnett, 2007). Additionally, F. heteroclitus are omnivorous and act as a major prey source for fish-eating birds and other predatory fishes [30]. This trophic centrality has been suggested to be a driver of indirect life-cycle parasite persistence [31], and in other parasite community studies they have been found to host at least 22 species of helminth parasites [32]. In Georgia, USA, at least three species of metazoan parasite with indirect life-cycles have been described infecting F. heteroclitus. These include a trematode, Lasiotocus minutus, that has F. heteroclitus as the definitive host, molluscs as the first intermediate host, and fish or crustaceans as the second intermediate host [33]. A cestode, Glossocercus caribaensis, that has fish-eating birds as the definitive host, crustaceans as the first intermediate host, and F. heteroclitus as the second intermediate host [34]. The last indirect life-cycle parasite was a nematode, Contraceacum sp., that can parasitize fish-eating birds and mammals as definitive hosts, and F. heteroclitus as the sole intermediate host.
At each timepoint, 5 baited minnow traps were randomly placed at every site, during the flood tide and collected during the ebb tide. We placed all fish from the 5 traps into a holding tank, and randomly subsampled 30 killifish to obtain a representative population of each site and date combination. This sampling did not take into consideration the size or age of the fish, however prior research has suggested that a sample size of 30 from a host population is appropriate to capture parasite community diversity within a host population [35]. The fish from each collected site were separately housed in the animal facility at Georgia Southern University and necropsied within 1 week (IACUC #I14013).
Each fish was euthanized in a buffered 300 mg L -1 solution of tricaine methanesulfonate (MS-222) until cessation of opercula movement followed by sex determination and measurements of total length and weight. The MS-222 solution, all external surfaces, gills, and opercula were examined for ectoparasites and helminths. The viscera, gonads, spleen, liver, heart, intestine, and bladder were removed and examined for helminths. All helminths from the first collection were heat fixed, stored in 70% ethanol, stained in acetocarmine, and mounted in Permount. Parasites from all collections were identified with keys and primary literature [32,36].

Statistical analyses
To visualize the variation in parasite component community and infracommunity structure, we performed non-metric multi-dimensional scaling (NMDS). In the NMDS for component communities, average abundances were calculated and square root transformed to control for the effect of high abundance parasite taxa. In the NMDS for infracommunities, raw abundances were ln(x + 1) transformed. Resemblance matrices were constructed using Bray-Curtis dissimilarity distance, and stress was calculated using Kruskal's stress formula 1. To determine whether there were significant differences in parasite component community structure, we performed a permutational multivariate ANOVA (PerMANOVA). For this analysis, we used site as a factor and performed 999 permutations. Additionally, we performed a 2-way crossed similarity percentage (SIMPER) analysis to identify the parasite species that contributed to 5% or greater of the variation in parasite community structure among sites. These analyses were conducted in R statistical environment using package "vegan" [37,38].
To determine the roles of landscape composition and physicochemical factors in structuring the parasite infracommunity and component community, we built multivariate random forest (MRF) models using R package randomForestSRC [39]. MRF is an ensemble machine learning method, that builds a suite of bootstrapped multivariate regression trees [40]. MRFs are increasingly being used in ecological studies because they allow multivariate species composition response data, do not assume independence among samples, make no a priori assumptions on the relationships between the response and predictor variables, and can account for imbalanced data [41,42]. Our response variable was the multivariate parasite species data, and our predictor variables included 8 landscape factors (Table 1), 23 physicochemical concentrations (Tables 2 and 3), and host length, weight, and sex. We built 4 MRFs, examining landscape and physicochemical factors effects on direct and indirect parasite communities at the infracommunity and component community scale. MRFs were built using 1000 trees and "Breiman-Cutler" variable importance was calculated through permutation. Variable importance was standardized to average variable effect across species. In these variable importance permutations, each predictor variable is placed in the out-of-sample (out-of-bag for univariate analyses) data for the tree model. The out-of-sample prediction error is then calculated with and without the permutation and averaged over the MRF trees. Greater differences between the permuted and non-permuted out-of-sample prediction error equate to greater variable importance.

Parasite community in F. heteroclitus
Eight taxa of metazoan parasites were identified from 630 fish, representing 24 different parasite communities. The direct life-cycle parasites included the branchiuran Argulus funduli, the copepod Ergasilus funduli, the leech Myzobdella lugubris, and the monogeneans Swingleus ancistrus and Fundulotrema prolongis. The indirect life-cycle parasites included the digenean Lasiotocus minutus, the cestode Glossocercus caribaensis, and the nematode genera Contracaecum sp. These 8 taxa of parasites infected more than 70% of the killifish examined (Tables 4 and 5).

Landscape and physicochemical predictors of parasite community in F. heteroclitus
We identified 12 metals from water analyses and 5 metals from sediment analyses that surpassed concentration thresholds (Tables 2 and 3), and our spatial analysis included 9 landscape metrics that varied across our collection sites (Table 1). We built MRF models to determine important variables associated with direct and indirect life-cycle parasites at the infracommunity (within an individual host) and component community (within a host population at a collection site) scales. The model predicting direct life-cycle parasite community structure at the infracommunity scale explained 40.18% variance, with fish weight and length being the most important variables (Fig 3A). The model predicting indirect life-cycle parasite community structure at the infracommunity scale explained 21.4% variance, with fish weight, fish length, and percent change in developed land being the most important variables (Fig 3B). The model predicting direct life-cycle parasite community structure at the component community scale explained 70.67% variance, with fish weight and length being the most important predictors (Fig 3C). The model predicting indirect life-cycle parasite community structure at the component community scale explained 73.43% variance, with sediment Ni, patch density, and marsh size being the most important predictors (Fig 3D).

Discussion
The goal of this study was to determine the effects of urbanization, mediated through alterations in landscape composition and physicochemical factors, in structuring parasite communities. Our sampling in coastal Georgia, USA provided evidence supporting the role of urbanization in impacting salt marshes via increased metal contamination, altering landscape composition at a fine spatial scale (250 m), and at the watershed level. We determined that parasite community composition differed among sites with different levels of urbanization. Variation in the prevalence and intensity of infection by indirect life-cycle parasites, the digenean L. minutus and the cestode G. caribaensis, were the primary drivers for the differences in parasite communities. Although the prevalence of direct life-cycle parasites differed among sites, this difference did not contribute to overall variation in component community composition. We found that factors driving community structure differed between the infracommunity and component community for indirect life-cycle parasites. Host size was the most important variable for direct life-cycle parasite assemblages and indirect life-cycle parasites at the infracommunity level, while landscape and physicochemical factors determined the structure of indirect life-cycle parasite assemblages at the component community scale. These data support the prediction that in the common marsh killifish, parasite community composition is explained by metrics associated with human modification of the environment, and parasite responses to urbanization are dependent upon life history.
Our results indicate that anthropogenic physicochemical and landscape factors structure indirect life-cycle parasites at the component community scale. The most important physicochemical factor was sediment nickel concentration, and although nickel occurs naturally and  is essential to most organisms as a trace mineral, increased concentrations due to human activities (i.e., burning fuel [44]) can be toxic to salt marsh benthic organisms. For example, Field et al. (2002) built a suite of logistic regression models based on toxicity assays of marine amphipods, and predicted that 50% of the population would show toxicity effects at concentrations of 47 ppm [45]. The concentrations we observed exceed this level, and may affect intermediate host populations that are an integral part of the transmission of the indirect life-cycle parasites we described in our study, the digenean Lasiotocus minutus and the cestode Glossocercus caribaensis. Additionally, these indirect life-cycle parasites have life stages that are external to the host, and the nickel concentrations we measure may be toxic to these free-living parasite life stages further impacting parasite component community diversity. The most important landscape factors for indirect parasites at the component community level were marsh size and patch density. Marsh size has been shown to affect carbon pathways to resident invertebrates, suggesting potential impacts on animal community and resultant trophic structure and interactions that are exploited by indirect life-cycle parasites for transmission [46]. For example, Lowe and Peterson (2014) found marsh size to be positively correlated with abundances of Gulf Menhaden and brown shrimp, though they did not find an association with blue crab abundance in the same sites [1]. Patch density is a measure for landscape heterogeneity. In urbanized estuaries, salt marshes become increasingly fragmented and the landscape increases in heterogeneity [47]. This often has the result of isolating habitat patches, leading to alterations in animal community structure and trophic interactions [48,49]. For example, Rudershausen (2018) found that salt marsh connectivity influenced fish community structure and the probability of F. heteroclitus movement [50]. These studies illustrate that landscape alterations have the potential to alter free-living host community structure; given that indirect life-cycle parasites require stable trophic interactions [5] between animals across trophic levels, altering the structure of the free-living community will necessarily be reflected in the parasite community [51].
Although the infracommunity is a sample within the component community, we found the factors that structure indirect life-cycle parasites differ between these two levels. Parasite transmission is inherently a local phenomenon, and individual risk factors do not necessarily equate to population risk factors. For example, the encounter-dilution effect hypothesizes that individuals can decrease their individual risk of parasitic infection through the formation of groups [52] even as the population level of infection increases. We found that the diversity of parasites in an individual fish was driven by host size (length and weight), but the diversity of parasites in the population was mostly affected by abiotic measures (sediment Ni, patch density, and marsh size). This subtle difference is because parasite community assembly is influenced by processes operating at a range of spatial and temporal scales. Parasite species are found within a regional species pool that is constrained by evolutionary processes, and only a subset of the species from the regional pool will colonize a site depending on ecological factors such as dispersal and the probability of encountering a suitable host. Subsequently, a subset of the species from the regional pool will colonize a particular site depending on dispersal and exposure probability. In essence, our data suggests that the observed parasite community within a host and the factors that determine what parasites you find in that host are dependent on the scale at which the study is conducted because of different abiotic and biotic filters [20,51,53,54].
Finally, our results indicate that host body size was the strongest determinant of direct lifecycle parasites at the infracommunity and component community scales. Individual hosts can be thought as habitat "islands" for parasites [55]. For direct life-cycle parasites, island biogeography theory can be applied [56], wherein larger hosts are larger islands, and consequently provide a larger habitat for a parasite to colonize, allowing for a greater number of parasite individuals and concurrent coexisting species [57,58]. In fish, body growth is largely indeterminate, and is associated with age, which may result in older fish harboring more parasite species and individuals simply because of an increase in the opportunity for acquiring an infection.
Previous studies have shown a relationship between urbanization and parasite community variation. Blanar et al. (2016) demonstrated a correlation between parasite community structure and local contaminants and land use within 5 km [9]. In this study, the relative abundance of the indirect life-cycle digenean was positively associated with an increase in crude oil contaminants, whereas the opposite pattern was observed for the directly transmitted monogenean parasites. This observation follows results of an earlier meta-analysis, where parasite life history (>1 obligatory host in life-cycle vs. 1 obligatory host in life-cycle) and habitat (external vs. internal) informed whether aquatic pollution had a significant effect on parasite population dynamics [59]. Additionally, Calegaro-Marques and Amato (2014) revealed an association between an urban-rural gradient and life-cycle variation in parasite communities of Rufousbellied thrushes [8]. Our results agree with these studies, and support the proposition that urbanization disproportionately affects indirect life-cycle parasites, which is then reflected in the parasite component community structure.
Coastal urbanization and anthropogenic disturbance are rapidly altering free-living and parasite community dynamics. Our results suggest that urbanization and its effects are important in structuring parasite communities of the common marsh killifish, indirect life-cycle parasites are more affected than direct life-cycle parasites, and the degree to which individual parasite species respond to urbanization is determined by their specific life history strategy. Further, our data suggest that urbanization likely affects interactions of free-living hosts, and their environment which is then detected as variation in the parasite community.