Contaminants of emerging concern in tributaries to the Laurentian Great Lakes: II. Biological consequences of exposure

The Laurentian Great Lakes contain one fifth of the world’s surface freshwater and have been impacted by human activity since the Industrial Revolution. In addition to legacy contaminants, nitrification and invasive species, this aquatic ecosystem is also the recipient of Contaminants of Emerging Concern (CECs) with poorly understood biological consequences. In the current study, we documented the presence, concentrations, and biological effects of CECs across 27 field sites in six Great Lakes tributaries by examining over 2250 resident and caged sunfish (Lepomis ssp.) for a variety of morphological and physiological endpoints and related these results to CEC occurrence. CEC were ubiquitous across studies sites and their presence and concentrations in water and sediment were highest in effluent dominated rivers and downstream of municipal wastewater treatment plant discharges. However, even putative upstream reference sites were not free of CEC presence and fish at these sites exhibited biological effects consistent with CEC exposure. Only the Fox River exhibited consistent adverse biological effects, including increased relative liver size, greater prominence of hepatocyte vacuoles and increased plasma glucose concentrations. Canonical Redundancy Analysis revealed consistent patterns of biological consequences of CEC exposure across all six tributaries. Increasing plasma glucose concentrations, likely as a result of pollutant-induced metabolic stress, were associated with increased relative liver size and greater prominence of hepatocyte vacuoles. These indicators of pollutant exposure were inversely correlated with indicators of reproductive potential including smaller gonad size and less mature gametes. The current study highlights the need for greater integration of chemical and biological studies and suggests that CECs in the Laurentian Great Lakes Basin may adversely affect the reproductive potential of exposed fish populations.


Introduction
The Laurentian Great Lakes contain one fifth of the world's surface freshwater and are home to diverse aquatic and terrestrial habitats that support a large number of fish and wildlife resources. Throughout the Great Lakes, multiple areas have been designated important for resource management, including National Wildlife Refuges (https://www.fws.gov/refuges) and National Estuarine Research Reserves (https://coast.noaa.gov/nerrs). The Great Lakes are also an important economical resource due to recreational and commercial fishing as well as multiple commercial ports. These ports were instrumental in the expansion of manufacturing during the Industrial Revolution in North America. The detrimental environmental legacy of the Industrial Revolution in urban centers around the Great Lakes has been the focus of toxicological studies for many years [1], including recently, studies of Contaminants of Emerging Concern (CECs) [2] which may interact with molecular pathways in exposed organisms [3][4][5][6]. Numerous laboratory and a limited number of field studies have identified consequences of CEC exposures ranging from molecular to organismal and population-level effects [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. These compounds, including pharmaceuticals, personal care products, pesticides and industrial chemicals enter aquatic environments through agricultural and urban storm water runoff, on-site septic system discharge, and outflows from municipal wastewater treatment plants and industrial sources [7,[22][23][24][25][26][27][28][29][30][31][32][33][34][35][36].
The multitude of CEC sources result in complex mixtures [37] that challenge our ability to predict effects to biota based on chemistry alone. Consequently, the evaluation of biological effects of CEC exposure in fish populations requires an integrated approach to acknowledge the complexity of the associated aqueous chemistry and fish biology. The objectives of this study were two-fold. First, we examined the presence and severity of biological effects commonly associated with exposure to CECs. Second, we assessed whether biological effects in resident and caged sunfish correlated with the presence of CECs in the adjacent aquatic environment.
We selected six Great Lakes tributaries for the current study to include the array of land use patterns common in the Great Lakes Basin (Fig 1). These land uses ranged from predominately rural, including forest interspersed with agriculture and limited municipalities, to intensely agricultural with limited urban influence, to exclusively urban. Most rivers contained a range of influences across their watersheds and, whenever land use allowed, a putative reference site (usually in the headwaters) was included. Concurrent resident and caged sunfish (Lepomis ssp.) assessments were used to more fully assess potential effects of CEC exposures to fish [38] (Fig 1).
The resultant multi-dimensional matrix of chemical, biological, and geographic data required an integrative multivariate model to identify relationships between two sets of variables, chemical occurrence and biological consequences. Canonical analyses, a permutation of principal component analyses, provide frequently used tools in community ecology to reduce the multidimensionality of data matrices and identify differences and relationships between groups of samples [39]. For this study, Canonical Redundancy Analysis (RDA) was used with CEC concentration as explanatory variables and biological effects as response variables [40]. Using RDA, the relationship between both sets of variables can be evaluated by determining how much variation in the response variables can be explained by the set of explanatory variables. This is accomplished by regressing multiple response variables onto multiple explanatory variables [41]. The resultant matrix is then analyzed using a principal component analysis. However, while relationships may be observed between the groups, this analysis cannot attribute causality.
We expected that CEC concentrations would be higher downstream of wastewater treatment plant outfalls and lower at upstream sites, with resident and caged fish showing increasing evidence of CEC exposure along a downstream gradient. Concurrently, we expected that an increase in the hepatosomatic index, liver hepatocyte vacuolization, glucose concentrations, and pathologies would be associated with an increase in CEC concentrations. Condition factor, gonadosomatic index and gonadal maturity would decline as CEC concentrations increase. Plasma vtg concentrations are likely to increase with increasing CEC concentrations, however, since many CECs are not estrogenic, this prediction may not accurately reflect the complexity of CEC mixture effects. When combined, these two objectives allow for the largest study to-date examining the biological consequences of CEC exposure on resident and caged fish in the Laurentian Great Lakes tributaries.

Ethics statement
Permission for access to field sites and scientific collectors permits were secured by the U.S. Fish and Wildlife Service. All resident fish collections occurred on public land and did not  Table 1. River basins are outlined on overview map (top). Maps of the six river watersheds are highlighted below. White circles indicate approximate resident fish collection and caging sites in each river. Stars indicate approximate resident fish collection and caging sites in each river that are immediately below wastewater discharge outfalls. Land use is indicated with light-and medium-grey patches of forest and agricultural and dark gray urban areas with linear road features. Map was prepared using ArcGIS 10.3.1 software by Esri.
involve any protected species. This study was carried out in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. The protocol was approved by the Institutional Animal Care and Use Committees of St. Cloud State University (permit number: 0213). Fish were maintained in aerated coolers filled with river water until they could be processed (usually within eight hours). Fish were euthanized by a sharp blow to the head or by immersion in the fish anesthetic MS-222. All efforts were made to minimize distress to the animals.

Study area and design
A total of six rivers, containing 27 study sites, were sampled over two years (Fig 1; Table 1) and represented an array of land uses common in the Great Lakes watershed. For the purposes of the current study, three stream reaches (upstream, middle, downstream) were identified for the collection of resident sunfish (six rivers x three reaches = 18 stream reaches). Furthermore, within each stream reach, we identified at least one and as many as three study sites (for a total of 27 study sites) for the caging of hatchery-reared sunfish and the collection of water and sediment samples (Table 1). Land uses were evaluated using the National Land Cover Database [42]. Land use categories were combined into the more general groups of Urban, Forest, Grass, Crop, Water, and Other [37]. For every river basin, the total number of pixels for each land use category were determined to calculate the area of each land use, the areas of each land use category were summed, and the percentage of each land use category was determined within each basin.
Limitations to migration, including dams or sufficient distance between reaches, were used to identify reaches of likely distinct Lepomis populations. For all rivers, the middle reach was a recipient of treated wastewater from a wastewater treatment plant (however, other river reaches could also be recipients of treated wastewater, see Table 1). Fish were caged at three to six sites in each river across the three reaches to represent each reach as well as specific potential point sources. Where logistically feasible we placed a cage in the discharge zone of major wastewater treatment plants, as well as upstream and downstream of these plants (Table 1, Fig 1). Whenever land use permitted (i.e., upstream reaches with little shoreline development), an upstream putative reference site was included.
The Cuyahoga River (Ohio) and Raquette River (New York) were sampled in 2013 and 2014 to examine annual changes in CEC presence and biological effects. The Clinton River (Michigan) was sampled in 2013, while the Lower Wolf and Upper Fox Rivers (Wisconsinhereafter "Fox River"), the North Shore Channel of the Chicago River (Illinois-hereafter "Chicago River"), and the Little Calumet River (Illinois) were sampled in 2014 (S1 Table).
The Fox River is dominated by cropland and dairy farming with small towns located along the river. The intensity of agricultural land use in this watershed precluded the identification of a suitable reference site. The north branch of the Chicago River is densely urban with study sites upstream, at, and immediately downstream of a major municipal wastewater treatment plant. The Little Calumet River is a commercial shipping waterway located just south of Chicago, IL. Land use in the basin is densely urban and industrial with minor cropland, forest or grassland. The middle reach is located immediately downstream of a major municipal wastewater treatment plant, while the other two sites are upstream and downstream of the plant. The Chicago and Little Calumet Rivers both consist primarily of wastewater effluent from the multiple discharge stations located along the rivers. The Clinton River basin, located northeast of Detroit, contains a mix of forest, grass, and cropland in the northern section of the drainage basin and primarily urban land use farther south. The upper reach is above identified point sources but is surrounded by urban land uses. The middle and downstream are influenced by Sample year, basin drainage area (km 2 ), estimated land use in the drainage basin (as percent of total land use), river reach section, sample site identifier, latitude and longitude, distance (in km) from each site to the next downstream site (including total distance to river mouth at a Great Lake), number of dams between that site and the next site upstream, and the number of identified point sources (e.g., wastewater treatment plants or combined sewer overflows) between that site and the next site upstream. Due to the complex and managed hydrology of the Chicago Area waterways system, distance from Lake Michigan for Little Calumet River sites do not imply an intuitive upstream-downstream flow direction. https://doi.org/10.1371/journal.pone.0184725.t001 multiple municipal treatment plants. No suitable reference sites could be located in the intensively used Chicago, Little Calumet, and Clinton Rivers. The lack of suitable reference sites had been noted previously [4] for the former two rivers. The Cuyahoga River represents a mix of urban, agricultural, and forested land uses. We identified a primarily forested upstream reach as a putative reference site. Its middle reach is influenced by combined sewer overflows located just downstream of the city of Akron and a major municipal treatment plant. The Raquette River is forested with occasional dairy farming and some municipalities along the river. The upstream reach is expected to be minimally impacted and was identified as a putative reference site, while the middle reach is located just downstream of the Tupper Lake municipal wastewater treatment facility. The downstream reach contained sites that bracketed the Potsdam, NY municipal wastewater treatment plant.

Water sample collection
The U.S. Geological Survey (USGS) collected surface water samples at each site twice during each sampling period, near the beginning and end of the caged fish exposures. Sediment samples were collected once from each site. Detailed sampling and analytical methods, along with a summary of quality-control data, are published elsewhere [43,44]. Briefly, water samples were collected in 1-L, baked amber glass bottles using a depth-integrated sampling method with a weighted-bottle sampler. The most recently-deposited bottom sediments (top 10 cm) were sampled using a stainless-steel Ekman dredge. All equipment was stainless steel or other inert material. Water and sediment samples were shipped on ice to the USGS National Water Quality Laboratory.

Chemical assessment
At the USGS National Water Quality Laboratory, water samples were analyzed for over 200 chemicals. For wastewater indicators (e.g., fragrances, flame retardants, plasticizers, and domestic pesticides), unfiltered water samples were analyzed by continuous liquid-liquid extraction and capillary-column gas chromatography/mass spectrometry (GC/MS) following published methods [45]. Steroid hormones, sterols, and bisphenol A were analyzed in unfiltered water samples by solid-phase extraction, derivatization, and gas chromatography with tandem mass spectrometry (GC/TMS) following the methods in Foreman et al. [46]. Filtered water samples were analyzed for pharmaceuticals by direct aqueous injection high-performance liquid chromatography with tandem mass spectrometry following published methods [47]. Sediment samples were analyzed for over 100 chemicals. For wastewater indicators, sediment samples were dried and solvent-extracted. Extracts were cleaned with a solid-phase extraction, and then chemical concentrations sediment samples were dried and solventextracted. Extracts were cleaned with a solid-phase extraction, and then chemicals were determined by capillary-column GC/MS following methods in Burkhardt et al. [48]. For steroid hormones, sterols, and bisphenol A, sediment samples were analyzed by adapting the methods in Burkhardt et al. [48], which are described in Fischer et al. [49]. Accelerated solvent extraction techniques were used to extract pharmaceuticals and antidepressants from sediment samples. High performance liquid chromatography, as described in Lee et al. [50], was then used to determine concentrations of pharmaceuticals and antidepressants. All chemical data, including quality assurance and quality control information are previously published [43,44]. and long term exposed resident fish. The second was differences in continuous exposure (caged fish confined to the exposure area) to periodic exposure (in-stream movement of resident fish). Resident sunfish (genus Lepomis) were collected during the spawning season in spring 2013 and/or spring 2014 using standard boat electro-shocking equipment and transferred to processing areas using aerated holding tanks. An attempt was made to collect 40 sunfish in the genus Lepomis at each river reach [38] as they are native and locally abundant throughout the Great Lakes watershed and maintain site fidelity during the spring when males guard nest sites [51,52]. While sunfish are common, the abundance of each species often differs widely between reaches and collection events. As a consequence, we combined the most commonly found sunfish species, bluegill (Lepomis macrochirus), pumpkinseed sunfish (Lepomis gibbosus), green sunfish (Lepomis cyanellus) and redear sunfish (Lepomis microlophus) for the analysis of the effects of CECs on resident sunfish [4]. This approach is justified by the shared habitat of these species, their phylogenetic closeness [53,54] their frequent hybridization [53,55] and the fertility of hybrid offspring [56].
Sunfish were sacrificed on site or at nearby laboratory facilities within eight hours of their collection. Whole body wet weights (0.01g precision) and length were measured for each fish. Blood was drawn from the caudal vasculature and centrifuged at 5,000 g for 8 min to separate plasma. Gonads and livers from each fish were excised and immediately weighed as wet weight (0.001 g precision, Mettler Toledo AG245, Columbus, OH). Care was taken to weigh organs quickly to avoid desiccation while keeping the organs away from any fluid that would further wet the tissue. Three to five tissue samples, each approximately 3 mm 3 in volume, were collected from each organ and placed into histology specimen cassettes, immersed in 10% neutral buffered formalin, and shipped to St. Cloud State University for histological processing. Plasma samples were shipped on dry ice to St. Cloud State University and stored at -80˚C until analysis.
Similar to other studies using caged fish as sentinels for environmental exposures in aquatic environments [26,57,58,59], at least 40 hatchery-reared sunfish were caged and deployed at multiple sites per river ( Table 1, Fig 1) for 14 days concurrent with resident fish collections. Prior to deployment, a subsample of sunfish was processed to determine baseline conditions and to assure that fish were mature at the time of deployment. Fish were transported in aerated coolers from nearby hatcheries to the field sites, placed in cylindrical cages (1 m length x 0.4 m diameter) made of plastic mesh (1cm mesh size) reinforced with PVC pipe and anchored to the bottom of the waterway, ensuring that the fish cage was in direct contact with the bottom sediment [4]. Two replicate cages were deployed at each site, each containing 25 mixed-sex bluegill sunfish (hybrid bluegill in the Raquette River in 2013 due to a lack of bluegill sunfish). After retrieving the cages, fish were placed in aerated coolers and processed as described for resident sunfish (see above).

Biological assessment
Morphometric indices. Body weight and total length were used to calculate the condition factor (CF = body weight/[total length] 3 ), a measure of the overall health of the fish [60]. Gonad and whole-body weights were used to calculate the gonadosomatic index (GSI = gonad weight/whole body weight X 100). Liver and whole-body weights were used to calculate the hepatosomatic index (HSI = liver weight/whole body weight X 100).
Histological analysis. For histological assessments, tissue samples were fixed for at least 1 week, dehydrated through a series of ethanol and xylene baths in a Leica automated tissue processor TP 1050 (Leica, Wetzlar, Germany) and embedded in paraffin using a Thermo Scientific Microm EC 350-1 embedding station (Waltham, MA). A previous study [61] demonstrated that a central location and greater number of tissue sections taken from the testes of smallmouth bass was positively correlated with the probability to detect intersex in this organ. Consequently, embedded tissues were sectioned at approximately 1/3 and 2/3 of the depth of the sample (resulting in tissue slices~100 μm apart) using a Reichert-Jung cassette microtome (Leica, Wetzlar, Germany; 5 μm sections).
At least six sections from each organ (gonad, liver) were stained using standard hematoxylin and eosin techniques [62,63] in a Leica Autostainer XL, similar to methods used in other studies [12,64,65]. Histological sections were assessed by an experienced histologist (HLS) and ranked on a semi-quantitative scale (1-4) for vacuolization of liver hepatocytes (1-vacuoles visible in <5% of total area; 2-vacuoles small but throughout image in <25% of area; 3broad presence of large vacuoles in 25%-50% of area; 4 -vacuoles prominent covering more than 50% of the field of view) (S1 Fig). The developmental stage of the gonad (testis or ovary) was also evaluated based on the proportion of cell types visible in the field of view (female: perinuclear oocyte, cortical alveolar, early vitellogenic, late vitellogenic; male: spermatogonia, spermatocyte, spermatid, and spermatozoa) (S2 Fig). The overall maturity of the sample (on a scale of 1 = immature to 4 = only mature sperm present) was calculated as: Ovary ¼ ðð% primary oocyte Þ þ ð% cortical alveolar x 2Þ þ ð% vitellogenic x 3Þ þ ð% mature x 4ÞÞ=100 Seven macroscopic and microscopic pathological observations (1 = pathology present; 0 = pathology absent for each pathological category) were summed for each fish to generate a pathology score (0 = no observed pathologies; 7 = all pathology categories represented in sampled fish). All macroscopic pathologies (lesions, deformities, missing structures, parasites visible without magnification) in livers and gonads were combined into one value as these observations were uncommon. For liver and gonad separately, we recorded the microscopic presence of eosinophilic fluids, parasitic cysts and/or other pathologies (to include intersex, atretic oocytes, etc.) (S3 Fig). All histological assessments were blinded to eliminate observational bias by the assessor of the tissues. As a quality control measure, we re-assessed a subsample of histological sections a second time and compared the resultant maturity values between the initial and the repeat analysis. The calculated mean maturity values obtained for the two analyses differed by <1%.
Hematological analyses. Plasma samples were used to measure glucose and vitellogenin concentrations. We chose glucose as a biomarker for the overall metabolic physiology of the organism as it is related to short-term stress and protein metabolism and may also be indicative of malnutrition. A previous study [39] measured 16 blood related biomarkers, including six metabolic bioindicators, in wild-caught fish (including bluegill sunfish) and found only glucose to exhibit patterns that were consistent with exposure history (ash spill). Other blood parameters, such as triglycerides have also been used in previous ecotoxicological studies as indicators of bioenergetic status in fish (for example [66]), but these authors also noted the variability of this endpoint in resident fish. Given the number of fish assessed in the current study, as well as the ease and reliability of glucose measurements, we focused on glucose as an energetic indicator in the current study. Using 1 μL of plasma, a TRUEbalance Blood Glucose Monitor (Moore Medical, Farmington, CT) was used (detection range of 20 to 600 mg/dL) to determine plasma glucose concentrations, a proxy measure for whole body energetics.
Plasma vitellogenin was measured by antibody-capture competitive enzyme-linked immunosorbent assay (ELISA) incorporating a sunfish-validated anti-vitellogenin polyclonal antibody and purified sunfish vitellogenin as standard [67,68] and followed previously published protocols [4]. R-squared values associated with the standard curve plots were above 0.97, with the majority at 0.99, using at least seven standard concentrations ranging from 4.8 to 0.075 μg/mL. Each plate contained a set of standards for curve generation, and was read precisely at 20 minutes post-TMB addition. The minimum detection limit utilizing this standard curve was 3.75 μg/mL. Analysis of plasma samples was randomized across ELISA plates and analysis days to minimize assay drift. A high degree of homology within monophyletic taxa has been reported for the vitellogenin gene of at least three monophyletic taxa: the order cypriniformes [69], the genus Mugil [70], and the genus Micropterus [69]. Given these data, it is likely that the antibody-capture competitive ELISA incorporating a sunfish-validated anti-vitellogenin antibody and purified sunfish vitellogenin as standard would perform well for the four Lepomis species in the current study. Four plasma samples from each plate were replicated on a second plate to examine cross-plate variability.

Data analysis
The richness of the chemical data matrix and its inherent complexity required a priori reductive processes to allow for its integration with the equally rich and complex biological matrix (a more complete treatment of the chemical data is reported in the companion manuscript, Elliott et al., in press). Chemical results were reduced to only one concentration for each chemical per site in each analyzed matrix (sediment and water). This was achieved by identifying the maximum concentration that was detected for each chemical at each site, for both water and sediment. The maximum concentration for each site was used for this screening method because the dataset had an abundance of left-censored data with few samples per site, resulting in insufficient data to appropriately use other estimation techniques, such as maximum likelihood. In addition, the maximum concentration for a given chemical within a given matrix also represents the concentration of greatest toxicological concern, which is relevant when evaluating with biological samples.
Finally, the chemical dataset was further reduced by combining chemicals into 15 classes (S2 Table). The classes used for this analysis includes: Alkylphenols, Fecal Indicators, Flame Retardants, Fragrances, Hormones, Industrial, Insect Repellants, Organohalides, Other, Polycyclic Aromatic Hydrocarbons (PAHs), Pesticides, Pharmaceuticals, Phenolics, Plasticizers, and Sterols. To achieve this, chemical concentrations were summed for each class. Total chemical class concentrations were used to evaluate differences between rivers. The chemical analysis provided by the USGS National Water Quality Laboratory by default included several classes of chemicals not usually considered CECs. These included PAHs, phenolics, sterols, and fecal indicators. Despite the focus of this study on CECs, we elected to include these data in the analyses to provide the most comprehensive assessment of explanatory variables possible.
To more broadly assess the estrogenic potential of the water samples, an estradiol equivalency value (EEQ) was calculated for each study site. EEQs were derived based on those reported in literature [71][72][73][74][75]. The sum EEQ was attained by multiplying individual chemical concentrations by the EEQ multiplier (S2 Table), and then summing the results.
Analysis of biological endpoints was conducted in Graphpad Prism V 6.0a (Oxnard, CA). Gonad maturity values were arcsin-transformed prior to analysis [76]. Plasma vitellogenin and plasma glucose concentrations were log 10 -transformed prior to analysis [4]. To allow for statistical analysis of samples below or above the lower and upper detection limits, respectively, the following convention was applied: vitellogenin or glucose concentrations below the lower assay detection limit were recorded as ½ the calculated lower detection limit of the assay (1.88 μg/mL for plasma vitellogenin, 5mg/dL for glucose). Values above the upper detection limit of each assay were recorded as the nominal upper detection limit (8,000 ug/mL and 600 mg/dL, respectively). Normality of data was tested using a D'Agostine and Pearson omnibus normality test. As many data sets were found to lack normality or were data of ordinal nature (hepatocyte vacuolization, sum of pathologies), all data were analyzed using Kruskal-Wallis with Dunn's post-test. Presence/ absence data for biological endpoints were assessed using a Fisher's Exact Test.
To examine differences between in-river locations and differences between study years (for Cuyahoga and Raquette rivers only), data from male and female fish were analyzed separately using generalized linear mixed models and the GLIMMIX procedure in SAS V9.3 (SAS Institute, Inc., Cary, NC). In the primary analysis, differences between in-river locations (downstream, middle or upstream) were tested using a mixed model in which the collection/caging site was the fixed effect and river and sampling locations within years were random effects. The data from the Raquette and Cuyahoga rivers were used in a secondary analysis to test whether there were differences between fish captured in 2013 and 2014. In these models, year, collection/caging site, and their interaction were included as fixed effects while river and sampling location were considered random effects.
Next, we conducted the canonical redundancy analysis (RDA) on the data matrix using the vegan package in R [77] and obtained a new set of multivariate traits (see S1 File for R code developed for the current study). These multivariate axes remove the effects of correlations between closely aligned variables. The RDA technique was used to relate the matrix of biological response variables (Y) to the matrix of chemical explanatory variables (X; sums of concentrations by class). Explanatory variables were log-transformed to remove considerable skewness in the distributions. At some sites, certain chemical classes contained only censored (non-detect) data. In these cases, we substituted one-half the lowest sum of concentrations for the chemicals that were detected among the remaining sites prior to log-transformation. Although simple substitution techniques may impart bias, the rank order of sums by chemical class was preserved, while allowing the ability to explore general patterns in the data as revealed by the RDA. More roust techniques to estimate censored data, such as maximum likelihood estimation, require more samples than were collected for this study.
A secondary RDA was performed using only pharmaceutical results in water samples as explanatory variables (X, sum of concentrations by subclass). This secondary RDA was performed because the class Pharmaceuticals includes 119 chemicals (S2 Table), and the classification is mostly arbitrary with no common mode-of-action. Their only common feature is their pharmacological origin. Subclasses of pharmaceuticals were classified based on common mode-of-action. Only subclasses with a detection in at least one sample and those that have the potential to impact the biological endpoints measured were used for analyses. Subclasses included: anesthetics, anticonvulsants, antidepressants, antidiabetics, antihistamines, antimicrobials, cardiovascular drugs, muscle relaxants, nonsteroidal anti-inflammatories (NSAID), and opioids (S2 Table).

Water and sediment chemistry and estradiol equivalency values (EEQs)
Contaminants of emerging concern were ubiquitous in water and sediment samples collected from all sites [43,44]. A detailed analysis of water and sediment chemistry is provided in the companion manuscript (Elliott et al., in review). All 15 classes of CECs (S2 Table) were detected in at least one site for both water (Table 2) and sediment (Table 3) samples, although the number of CEC classes detected per site in water samples was always larger than in  sediment (mean 12.9 classes vs. 8.9 classes). Insect repellents, sterols and several other chemicals were documented in water samples from all 27 field sites (Table 2), while pharmaceuticals and pesticides were detected in all sites but one. Fecal indicators and PAHs were the most frequently detected classes in sediment samples ( Table 3). The highest total chemical concentrations in water samples were detected in the Chicago River, all but the most upstream site of the Little Calumet River, and site CUY-3, which is located immediately downstream of a municipal wastewater treatment plant in the Cuyahoga River In water samples, the sum estradiol equivalency (EEQ, Table 2) was the highest at one site in the Raquette River (RAQ-2: 0.028 μg/L EEQ), followed closely by a site in the Chicago River  Table 3). The Little Calumet River did show an increase in sediment EEQ values in a downstream direction, which was congruent with EEQ values in water samples. The highest EEQ value in sediment was detected in the Little Calumet River (CHI-76: 25.86 μg/kg EEQ), followed by the Fox River (LLB-5: 18.59 μg/kg EEQ). The lowest EEQ values (all 0.0000 μg/kg EEQ) were detected in upstream sites in the Clinton River (CLI-5), Cuyahoga River (CUY-5) and Raquette River (RAQ-5). The Little Calumet River had the highest average EEQ value (13.68 μg/kg) as compared to the other rivers. It is interesting to note that while the RAQ-2 site is not the highest value as compared to other sites (11.76 μg/kg), it is much higher than any of the other sites in the Raquette River (next highest was at RAQ-3: 0.29 μg/kg)

Biological endpoints
Biological endpoints (see S2 File for all raw biological data) associated with exposure to CECs exhibited a complex pattern across sites and whether resident or caged fish were assessed (  Table). Our biological sampling design allowed for a nominal total of 224 statistical analyses (eight river sampling events (six rivers plus two years of sampling in two rivers) x two sexes x two fish sources (resident/caged) x seven endpoints (excluding pathologies)). Some data gaps reduced this number to an actual 217 statistical comparisons of biological effects by river, of which 82 yielded statistically significant differences among sites in a river. Overall, 56 significant differences in endpoint expression were identified in resident fish, while only 26 significant differences were present in caged fish (Fisher's exact test, p<0.0001). The frequency of observed effects did not vary significantly between the sexes (Fisher's Exact Test, p>0.05; 30 effects in resident females, 26 in resident male fish; 9 and 17 for caged female and male fish, respectively). Plasma glucose concentrations differed significantly most frequently among fish collected from different sites or reaches within a river (21 instances by sex and resident vs. caged fish), followed by HSI (16 differences), plasma vitellogenin concentrations (13 differences) and prominence of hepatocyte vacuoles (11 differences). It is noteworthy that CF, GSI and gonad maturity varied less frequently among sites (8, 6, 7 differences, respectively) and that most of the variability was observed among resident fish. Caged fish exhibited little variability in these endpoints, likely as a reflection of the homogeneity of hatchery-reared caged fish that were acquired at a similar age and reproductive status.
When statistical analyses for all data sets were compared (both ANOVA and generalized linear mixed model), few patterns were found to be consistent across fish sex and/or resident and caged fish. Only in the Fox River were several biological endpoints consistently altered following a consistent upstream to downstream pattern. Here, CF decreased significantly from upstream to downstream reaches for resident male and female fish ( Figure A in S6 Fig,  Figure A in S7 Fig). This decrease was counter-directional to a consistent increase in HSI and liver hepatocyte vacuolization for both resident and caged fish (except vacuole prominence in resident male fish, Figs 2-5, S4 Table). Interestingly, plasma glucose concentrations showed significant but opposing changes in the Fox River resident female ( Figure D in Fig 2, S4 Table) and male fish ( Figure D in Fig 3, S4 Table), with females having decreasing plasma glucose concentrations and males increasing concentrations with distance downstream. These resident fish results were not paralleled by findings in caged fish from the Fox River.
Fish sampled in other rivers exhibited fewer significant differences in endpoint expression between sites or reaches (Chicago, Little Calumet, Clinton, Cuyahoga 2014) or differences were inconsistent between sites, sexes, and resident or caged fish (Cuyahoga 2013, Raquette 2013). It is noteworthy that in rivers that were sampled twice (Cuyahoga, Raquette), patterns were not consistent between years. Using a generalized linear model, significant differences for CF (p = 0.016) and maturity (p = 0.0016) were found between years for resident male fish (independent of river). Significant differences between years (independent of river) were also found for GSI and maturity in caged female and male fish (p<0.001 for all). In contrast, river location (upstream, middle, downstream reach) was a poor indicator to predict the expression pattern of any endpoint in female or male, resident or caged fish.

Canonical redundancy analysis
The first two axes of the RDA express a mean 43.5% of the variance (34.5% and 9.0% for the first (X) and second (Y) axis, respectively) across both sample matrices (water and sediment)    Table 1. Column graphs indicate mean + standard deviation in panels (A) and and all fish groupings (resident/caged, female/male; Fig 6; S10 and S11 Figs). The explanatory power of the analyses differed little if biological endpoints were compared only against water, sediment or both matrices combined (43.1%; 46.4%; 41.0%, respectively). Among the fish groupings, only resident male fish produced substantially lower first and second axis values (sum mean: 29.2%) when compared to resident female or caged female or male fish (48.6%; 49.2%; 47.2%, respectively). This is consistent with the partitioning of constrained and unconstrained variance of the data set (S5 Table). When constrained variance is greater than unconstrained variance in an RDA analysis, then the results indicate that much of the variation in the response variables (biological data set) is related to the explanatory variables (chemical data set). As a result, more of the variability in biological data in all female and caged male fish is due to the explanatory variables while the origin for the biological variability of resident male sunfish is less certain (Fig 6; S10 and S11 Figs).
In caged female fish, ordination identified pharmaceuticals (p<0.001) and PAHs (p = 0.008) as the significant explanatory variables in water samples, while fragrances were the significant explanatory variable for caged males (p<0.001). In sediment samples, RDA identified sum EEQ (p = 0.035) for resident females, and fragrances for caged females and males (p<0.001 for both) as significant explanatory variables. When sediment and water were combined, ordination identified pharmaceuticals and PAHs for both caged females (p<0.001 and p = 0.006, respectively) and males (p<0.001 and p = 0.008, respectively).
The RDA identified consistent patterns of endpoint expression (response variable) and chemical exposure (explanatory variable) for resident and caged fish and for both sexes (Fig 6;  S10 and S11 Figs). The vector length for each biological response variable indicates the strength of the relationship between the response variable the explanatory variables. Vectors whose angle approaches 90˚have little correlation with each other; vectors approaching 0˚are strongly correlated and vectors approaching 180˚are negatively correlated. Plasma glucose concentrations are consistently driving the differences along the first axis of the RDA plots and are consistently negatively correlated by GSI and gonad maturity (Fig 6; S10 and S11 Figs) suggesting a strong negative correlation between the former and two latter biomarkers. In contrast, CF aligns to the second axis of most RDA plots while HSI and the prominence of hepatocyte vacuoles are commonly aligned with each other. In caged fish, the association between HSI and prominence of hepatocyte vacuoles is weaker than in resident fish. Plasma vitellogenin concentrations often correlated with GSI and maturity in resident sunfish along the first RDA axis but are more commonly associated with the second RDA axis for caged sunfish. The alignment of pathologies with other biological endpoints is variable and explains little of the differences in the data matrices except in the analysis of sediment chemistry (S10 Fig). When examining endpoint expression in the context of CEC classes, plasma glucose appears to align more closely with urban-associated CECs (e.g., flame retardants and plasticizers), while HSI and prominence of hepatocyte vacuoles align more closely with endocrine active CECs (EEQ, steroid hormones) and are opposed by CF and often maturity.
The results of the pharmaceutical subclass RDA paralleled those observed in the other RDAs (S12 Fig). Resident male fish produced substantially lower first and second axis values (22.3%) when compared to resident female or caged female or male fish (52.1%; 58.3%; 52.6%, respectively). In addition, this observation is consistent with the partitioning of the constrained and unconstrained variance of the dataset (S5 Table). For biological responses, (B). Box-and-whisker plots indicate range, 25 th and 75 th percentiles, and mean values in panels (C) and (D). Statistical significance (Kruskal-Wallis with Dunn's post-test; p<0.05) within panels are identified by letters. Pvalues are summarized in S4 Table. Sample size provided in S3 Table. https://doi.org/10.1371/journal.pone.0184725.g002  Table 1. Column graphs indicate mean + standard deviation in panels (A) and (B). Box-and-glucose continued to be the most important driver of the results and was negatively correlated with CF, GSI, and maturity. Hepatocyte vacuolization and HSI often correlated along the second axis. In the redundancy analysis, no pharmaceutical subclasses were identified as significant explanatory variables in resident females (all p>0.05), while only anesthetics were identified for caged females (p<0.001) and opioids for caged males (p<0.001).

Discussion
We applied two lines of evidence gathered from resident and caged sunfish to assess the biological effects of CEC occurrence on fish in Laurentian Great Lakes tributaries. First, we assessed resident and caged fish for biological effects commonly associated with CEC exposure. Second, we evaluated whether the measured biological effects in fish were associated with the concentration of CECs in the surrounding aquatic environment.

CEC presence in Great Lakes tributaries
Contaminants of Emerging Concern were detected in all water and sediment samples collected across the Great Lakes tributaries. CEC presence extended to two putative reference sites we identified a priori based on their limited shoreline development and upstream locations. This omnipresence was expected as the targeted rivers are known to be subject to anthropogenic impacts. Consequently, the results of the current study are congruent with published reports of ubiquitous CEC presence in water samples from aquatic ecosystems in North America [78,79], South America [80], Europe [80] (reviewed in [81,82]), Asia (reviewed in [83]) and Australia [84] leaving little doubt that CECs are a global pollution concern with uncertain consequences for exposed aquatic life.
The results of our water analysis corroborated our predictions that CECs would be most common and at the highest concentrations downstream of municipal wastewater treatment plant outfalls and in effluent dominated rivers. CUY-3, just downstream of an effluent outfall registered some of the highest total contaminant concentration and was matched by most sites in the effluent-dominated Chicago and Little Calumet Rivers. Increasing total CEC concentrations were paralleled with increasing EEQ values, suggesting that estrogenic CECs contribute substantially to the overall CEC concentrations measured. Conversely, upstream sites usually contained fewer CECs and at lower concentrations with less EEQ activity. However, the ubiquitous presence of CECs across all sampling sites may affect study designs as it calls into question the very concept of reference sites. We had identified two putative reference sites (a mostly forested site in the Cuyahoga River and an upstream segment of the Raquette River in the Adirondack Mountains) and already concluded that no reference sites of suitable characteristics existed in the other four rivers. However, even our reference sites were found to contain CECs. Recently published data assessing CEC presence in 50 Minnesota lakes [85] found only three of 50 sampled lakes to be free of CECs. Several lakes without shoreline development, and in one instance within the Minnesota Boundary Waters Canoe Area (an area with strictly limited human access and prohibition of any use of combustion engines), contained detectable concentrations of multiple pharmaceuticals and even illicit drugs. A companion study of 50 Minnesota rivers found only one site free of CECs [86]. Given the population density and agricultural and urban land use intensity around the US portion of the Great Lakes watershed, whisker plots indicate range, 25 th and 75 th percentiles, and mean values in panels (C) and (D). Statistical significance (Kruskal-Wallis with Dunn's post-test; p<0.05) within panels are identified by letters. P-values are summarized in S4 Table. Sample size provided in S3 Table. https://doi.org/10.1371/journal.pone.0184725.g003 Biological effects of Great Lakes contaminants of emerging concern  Table 1. Column graphs indicate mean + standard deviation in panels (A) and (B). Box-andwhisker plots indicate range, 25 th and 75 th percentiles, and mean values in panels (C) and (D). Statistical significance (Kruskal-Wallis with Dunn's post-test; p<0.05) within panels are identified by letters. P-values are summarized in S4 Table. Sample size provided in S3 Table. https://doi.org/10.1371/journal.pone.0184725.g004 Biological effects of Great Lakes contaminants of emerging concern  Table 1. Column graphs indicate mean + standard deviation in panels (A) and (B). Box-andwhisker plots indicate range, 25 th and 75 th percentiles, and mean values in panels (C) and (D). Statistical reference sites practically do not exist. The use of a RDA analysis rectifies the lack of reference sites partially as it determines the degree to which explanatory variables in the data matrix may account for variability in the response variables (biomarkers). The results of a RDA, therefore, represent a "standard curve" against which other sites (for example hypothetical reference sites significance (Kruskal-Wallis with Dunn's post-test; p<0.05) within panels are identified by letters. P-values are summarized in S4 Table. Sample size provided in S3 Table. https://doi.org/10.1371/journal.pone.0184725.g005  Table 1. Sample class information can be found in Table 2. Biological response information can be found in Table 4, Figs 2-5, S6-S9 Figs, and S3 Table. https://doi.org/10.1371/journal.pone.0184725.g006 Biological effects of Great Lakes contaminants of emerging concern with little to no CEC contamination) could be extrapolated. Future studies should take the paucity of reference into consideration when identifying study locations and analysis tools.
Similar to other studies [87][88][89][90], the site-by-site composition of CECs varied, highlighting the complexity of pollutant source, transport and fate across the Great Lakes, as discussed in greater detail in the companion manuscript (Elliott et al., in review). Despite the large chemical data matrix generated in this study, assessing the full complexity of CEC occurrence and concentration variability remains elusive as these compounds vary dramatically over time. For example, Martinovic et al. [89,90] demonstrated in two studies order-of-magnitude changes in estrogenic activity in effluent discharged from several treatment facilities in the Upper Midwest. In one of these studies [90] the authors also identified highly variable biological activity associated with CEC presence in storm water runoff in an urban environment. The size of the data matrix in the current study to some extent compensates for the lack of temporal sampling (as it likely captured both high and low CEC occurrence across sites), but temporal CEC variability remains a concern for the interpretation of most environmentally derived biological data.
The three most urbanized rivers in our data set, the Chicago, Little Calumet, and Cuyahoga Rivers predictably contained the highest chemical concentrations, largest number of chemicals detected, and largest diversity of chemical classes in water. Previous studies have also noted the presence of a large number of CECs in highly urbanized aquatic ecosystems [4,64,90,91,92], particularly those with wastewater treatment plant effluents present. Wastewater treatment plant effluents are a well-studied and consistent source of CECs (reviewed in [93]) and were a prominent source in the current study. For all three rivers, pharmaceuticals factored prominently into the total CEC concentration (S4 Fig), following observations from other studies that detected higher concentrations of pharmaceuticals just downstream of major wastewater outfalls [79,94]. One caveat to this finding is that the pharmaceutical class encompasses a large number of diverse compounds which may contribute to its common presence.
In contrast to the highly urbanized rivers, the Raquette River had lower chemical diversity, but recorded the highest total EEQ value. While the Raquette River has a larger percent forested land than the other sites, it nevertheless contains multiple municipalities along the river that input CECs via their wastewater treatment facilities and also encompasses dairy farms in its watershed, which may have contributed to the high EEQ values. The differences in chemical diversity between rivers receiving wastewater may be due to differences in the composition of the raw influent [79] and its treatment [95] (reviewed in [96]).

Biological observations
The variability in biological effects paralleled and in some instances exceeded the complexity observed in the CEC occurrence matrix. While biological endpoints frequently differed between sites and reaches in the same river, the pattern of effect occurrence and the congruence between biological effects was not always intuitive and did not follow the a-priori prediction of greater biological effects consistent with increasing CEC exposure at downstream sites. However, some broad patterns were noted between sampling sites, fish collection methods, and sex. Effects were more common in resident fish than in caged fish. Direct comparisons of biological responses between resident and caged fish are infrequent in the published literature (however, see [58,97] for examples) but can be informative. Several factors may have contributed to the less pronounced consequences of CEC exposure in caged fish. First, caged fish are only exposed to the site-specific water for two weeks vs. potentially a lifetime for resident fish. In a previous study, Burki and colleagues [97] reported elevated vitellogenin mRNA levels in both resident and caged brown trout downstream, but not upstream, of a wastewater treatment plant in a small Swiss stream. In contrast to the matching vitellogenin mRNA levels in these fish, vitellogenin protein expression was only upregulated in resident fish, but not caged males. The authors hypothesized that this may be related to the differing exposure history of the two study populations [97]. Second, caged fish were in excellent physiological condition at the onset of the study. They were all hatchery reared with ad libitum feed until caging, and had near optimum growing conditions at the hatchery. The better nutritional status of caged sunfish may contribute to the greater plasma vtg concentrations observed in caged fish in comparison to male resident sunfish from the same sites (Table 4). A better nutritional status may allow for greater resource allocation to protein biosynthesis, especially in caged fish that were caged in the effluent plumes of wastewater treatment plant discharges.
When significant biological effects were observed between reaches or sample sites, their frequency was similar for female and male fish, whether in resident or caged fish. Among the seven endpoints that were statistically analyzed (pathologies were sporadically found in all treatments with no strong patterns), two of the three morphologic indices and staples of fisheries biology (CF, GSI) were the least frequently affected. In caged fish, in particular, these two indices provided little evidence of biological effects even when fish livers documented clear evidence of a chemical insult (i.e., changes in relative liver size and prominence of hepatocyte vacuoles), as was evident in the RDA. The muted response of CF and GSI may either reflect adaptive changes in resident fish [98] or the short exposure duration of the caged fish that had been well fed until deployment and may have had sufficient reserves to buffer the exposure stress. Livers, in contrast, serve as the principal detoxifying organ and respond rapidly to a chemical insult, often resulting in more prominent vacuoles in liver hepatocytes [99,100]. Because of the sensitivity and rapidness of change, liver vacuolization has been recommended as a measure of contaminant exposure [101].

Biological consequences of CEC exposure
RDA identified the response variable of plasma glucose concentrations, an indicator of metabolic stress [39], as a sensitive biomarker of CEC exposure in both resident and caged sunfish. Changes to the bioenergetics status as a result of contaminant exposure have been noted in previous studies including bluegill sunfish [39,66]. While Bevelhimer and colleagues [39] noted differences in glucose concentrations between wild-caught fish captured in areas affected by a coal ash spill, Adams and colleagues [66] noted changes in triglyceride concentrations in sunfish across a contaminant gradient. Both measures, glucose and triglyceride, are considered bioindicators of energetic stress and changes to metabolic processing. In the current study, the use of glucose as indicator was chosen due to its previous documented utility [39] and due to the ease of use and reliability of glucose monitors given the large number of processed fish samples. Given the prominence of alterations in blood glucose concentrations in the current study, it is important to note that a previous study documented changes in glucose concentrations for up to two hours after fish were exposed to electroshocking stimuli but not regular handling [102]. In the current study, resident and caged sunfish were processed within eight hours of either electroshocking and handling (resident fish) or handling alone (caged fish). Since positive correlations between glucose concentrations, HSI and hepatocyte vacuoles, and negative correlations with CF and GSI were consistent across resident and caged sunfish, we may assume that the effect of electroshocking resident fish did not bias our results.
Despite the complexity of observed biological effects, the Fox River provided evidence of multi-endpoint effects that are biologically consistent with CEC exposure. In the Fox River, which had a large number of chemical classes detected in water and had the second highest EEQ value in sediment, increased relative liver size (HSI) was consistently observed in resident and caged fish independent of sex. This was paralleled by an increase in the prominence of hepatocyte vacuoles (and a decrease in CF for resident fish) from upstream to downstream. In the 2014 fish samples collected from the Raquette River, prominence of hepatocyte vacuoles decreased in the middle reach of the river in resident and caged fish, independent of sex, and was matched by a similar pattern for plasma vitellogenin and glucose concentrations. The concentrations of both of these molecules in blood are intrinsically connected to the liver as vitellogenin is synthesized by the liver in response to increased estrogen concentrations, while glucose is stored (as glycogen) in this organ and released as a result of metabolic stress. Even with the above described examples of biologically consistent effects across endpoints, the overall association of CEC presence and biological effects was too multi-faceted to be fully comprehended by an analysis of variance of individual effects between sites within a river. Consequently, we used the RDA to identify patterns and associations among our data matrices and to reduce the overall complexity of the data set by eliminating co-variance. The value of this statistical tool becomes apparent when comparing the results between subsets of the overall data set (female vs. male; resident vs. caged fish). Plasma glucose concentrations consistently dominate the first axis (i.e., the axis with the greatest explanatory power) of the RDA, independent of fish origin (resident, caged) or sex, and is mostly consistent between matrices (water only; sediment only; combined water and sediment). Increased plasma glucose concentrations as a result of pollutant stress have been noted in previous studies of pesticide exposure [103,104], estrogenic compound exposure [105] and whole effluent exposure [106]. The consistent opposition of plasma glucose concentration and relative gonad size (GSI) and gonad maturity along the first axis of the RDA suggests that increasing plasma glucose concentrations are inversely correlated with the reproductive potential of fish (Fig 6; S10 and S11 Figs). Reduction in the relative size of gonads is commonly found in fish from polluted riverine sites [107][108][109].
The linkage between energy balance (represented by plasma glucose concentrations in the current study) and reproduction is well established [110] (reviewed in [111]), as the finite amount of energy available to an organism has to be balanced between basic physiological needs and reproductive output. Heightened metabolic stress as a result of pollutant exposure increases core energetic needs at the expense of reproductive output [110,112]. The energetic costs of CEC exposure are further demonstrated by the opposition of CF and HSI/hepatocyte vacuolization on the second axis of the redundancy analysis. Here again, increased energetic expenses associated with stress response are inversely associated with a reduction of CF and the overall nutritional health of the organism.
In contrast to the prominence of the above mentioned biological endpoints, plasma vitellogenin concentrations and tissue pathologies played only a minor and more variable role in the variability of the data matrix. Patterns for these endpoints varied between resident and caged fish and between water and sediment matrices. Given the focus on plasma vitellogenin concentrations in many CEC studies, especially estrogenic CECs, we expected a greater prominence for this endpoint in our analysis, particularly in male fish. However, total estrogenicity (EEQ) seldom reached concentrations beyond 20 ng/L. Only three of the 27 sites had measured EEQs greater than 20 ng/L, and only seven sites were above 10 ng/L (Table 2). Therefore, plasma vitellogenin synthesis would unlikely be induced in male fish based on previous laboratory studies [89,113].
Pathologies were too evenly spread across fish and sites to provide much insight. In the RDA, pathologies only took on a more prominent role when associated with sediment chemistry, perhaps reflecting the effects of hydrophobic compounds bound to sediment. It is, however, noteworthy that in 72% of resident fish at least one pathology (most frequently parasites in the liver) was observed. The effect of existing interactions between CEC exposures and pathologies, particularly parasite loads, represents a critical knowledge gap-especially considering that the liver is the main detoxifying organ of the organism and parasitic infestations are common in fish exposed to estrogenic CECs [114].
RDA identified several significant explanatory variables including pharmaceuticals, fragrances, and PAHs. Chemicals in all three classes have been implicated in previous studies as causing endocrine disruption in exposed fishes [58,115,116,117]. While pharmaceuticals were identified as a significant explanatory variable when all chemicals were evaluated, a secondary RDA evaluating only the contribution of pharmaceutical (grouped into subclasses with similar mode-of-action) did not identify any of the subclasses as having strong explanatory power. This suggests that the totality of pharmaceutical exposures adversely affected fish. Lastly, it is noteworthy that other abiotic and biotic factors that were not enumerated in the current study may have contributed to the observed effects. Higher pollutant loads are frequently associated with other physico-chemical and physical alterations in the environment (i.e., greater nitrification of waters, reduced dissolved oxygen, greater turbidity, declining habitat quality and availability). Beyond the scope of the current study, the interactions of multiple stressors, including CECs, is urgently needed to comprehensively assess the impact of CECs on the sustainability of fish populations [2,90]. Particularly in the Great Lakes Basin, aquatic organisms are exposed to dynamic variations in abiotic factors throughout the year, such as temperature and salinity, in combination with the variability in pollutant loads, which may be related to observed effects.

Experimental design considerations
The current study examined over 2250 sunfish across 27 sites. Despite this effort, it proved difficult to consistently obtain the desired double-digit sample sizes for robust statistical analysis of female and male resident and caged fish [38]. The variability and complexity of aquatic ecosystems remains a challenge for any field study. Although an attempt was made to conduct all field studies in a short time frame to match environmental conditions (for example water temperature) that influence the expression of biological effects in exposed organisms, physicochemical differences between rivers and stream reaches are unavoidable and will add complexity to the data matrix. For example, urban watersheds usually have higher average temperatures than rural watersheds and stream reaches below major wastewater treatment plant outfalls contain greater nutrient loads than putative reference sites in forested upstream habitats. Consequently, the RDA analyses, while consistent in identifying prominent biological response variables, only explained approximately half of the variation in biological responses. Adding to this natural complexity, the individual variability [118] found in every population further complicates meaningful analysis of the resultant data matrices. This variability requires the use of statistical techniques capable of integrating these data sets. It also requires the investigator to use best professional judgment to reduce the complexity of the matrices whenever possible (for example, by combining chemicals into "CEC classes" and by collapsing the genus Lepomis into a single entity) while recognizing that, especially for CECs, these judgement calls remain tenuous. In the current study, RDA provided an integration of two complex matrices that resulted in the identification of patterns that would have easily been lost using less integrative methods (for example general linear models). However, multivariate techniques are limited by their need for large data matrices to draw from. These limitations highlight the need to develop comprehensive experimental designs to improve our understanding of the effects of CECs in aquatic ecosystems. Mesocosm studies, on-site exposures using temporary ad-hoc laboratory spaces, or non-lethal sampling techniques may be required to close some of the existing knowledge gaps.
The current study illustrates a common pattern of CEC presence in tributaries to the Laurentian Great Lakes that is most pronounced in urban influenced aquatic ecosystems but that is also sufficiently widespread to achieve near omnipresence in the Great Lakes watershed. Biological effects were highly variable and frequently failed to follow an upstream-downstream pattern, but carried a consistent energetic cost as suggested by the strong plasma glucose response that may impact the reproductive potential of exposed fish [110] (reviewed in [111]). Although a field-based study of chemical presence and biological effects cannot attribute causality, the omnipresence of CECs and the documented biological effects support further study to predict causality with confidence.
The Laurentian Great Lakes have been a focal point for understanding legacy contaminants and subsequent remediation for decades [1]. More recently, there has been an increased focus on lakes due, in part, to algal blooms from increased nutrient runoff, increases in invasive species, and increased warming of water temperatures. The realization that the same ecosystem is also impacted by the presence of CECs is disturbing, but not unexpected. The Great Lakes Basin is home to diverse fish and wildlife habitats, and supports an economically valuable recreational and commercial fishery. Resource managers need to be aware of the multiple stressors that may be impacting their resources. Our study raises awareness of the risks that CEC may be posing to these resources, especially in areas with greater human influence and urbanization. With the consistently changing environment, and the abundance of stressors that may be present, CECs represent another threat that needs to be managed to protect and conserve the fish and wildlife resources in the Laurentian Great Lakes.

Conclusions
This study correlated an extensive matrix of CECs in six Great Lakes tributaries with an equally extensive matrix of biological responses in resident and caged sunfish. The study confirmed an almost ubiquitous presence of CECs that calls into question the presence of relevant reference sites in the Great Lakes watershed. The complexity of the chemical occurrence data was matched by the complexity of the biological responses, which exhibited few consistent patterns of change when assessed solely between reaches or sites within a river or between rivers. However, full integration of both explanatory chemical variables and biological response variables revealed consistent biological effects of CEC exposure. Chief among them was a change in blood glucose concentration, an indicator of bioenergetics and stress, which was associated with increased relative liver size and greater prominence of hepatocyte vacuoles. These indicators of pollutant exposure were inversely correlated with indicators of reproductive potential, including smaller gonad size and less mature gametes. The current study highlights the need for greater integration of chemical and biological studies and suggests that CECs in the Laurentian Great Lakes Basin may adversely affect the reproductive potential and sustainability of exposed fish populations.  Table 1. Column graphs indicate mean + standard deviation in panels (A), (B), and (D). Box-and-whisker plots indicate range, 25 th and 75 th percentiles, and mean values in panel (C). Statistical significance (Kruskal-Wallis with Dunn's post-test; p<0.05) within panels are identified by letters, with the p-value identified below graphs. Sample size provided in S3 Table. (TIF)  Table 1. Column graphs indicate mean + standard deviation in panels (A), (B), and (D). Box-and-whisker plots indicate range, 25 th and 75 th percentiles, and mean values in panel (C). Statistical significance (Kruskal-Wallis with Dunn's post-test; p<0.05) within panels are identified by letters, with the p-value identified below graphs. Sample size provided in S3 Table. (TIF)  Table 1. Column graphs indicate mean + standard deviation in panels (A), (B), and (D). Box-and-whisker plots indicate range, 25 th and 75 th percentiles, and mean values in panel (C). Statistical significance (Kruskal-Wallis with Dunn's post-test; p<0.05) within panels are identified by letters, with the p-value identified below graphs. Sample size provided in S3 Table. (TIF) S9 Fig. Biological endpoints in caged male sunfish. (A) condition factor; (B) gonadosomatic index; (C) gonad maturity ranked on a scale of 0 to 4; (D) sum of pathologic observations on a scale of 0 to 7. Sample river location located above panels (A) and (B), with columns representing upstream to downstream within each river from left to right. Specific sample site identification can be found in Table 1. Column graphs indicate mean + standard deviation in panels (A), (B), and (D). Box-and-whisker plots indicate range, 25 th and 75 th percentiles, and mean values in panel (C). Statistical significance (Kruskal-Wallis with Dunn's post-test; p<0.05) within panels are identified by letters, with the p-value identified below graphs. Sample size provided in S3 Table. (TIF)  Table 1. Sample class information can be found in Table 3. Biological response information can be found in Table 4.  Table 1. Sample class information can be found in Table 2a and 2b. Biological response information can be found in Table 4 Table 1. Sample subclass information can be found in S2 Table. Biological response information can be found in Table 4. (TIF) S1 Table. Sampling site information. Additional sample site information, including sample site ID, river reach section, latitude, longitude, sample site station name, identification of nearest city or municipality to the sample site, 2014 population and population density (person/km 2 ) in the drainage basin, and dates of field sample events for each station.  Table. Sample size for fish analyzed in the current study. Sample size may vary slightly between biological endpoints due to factors such as insufficient tissues samples or other postdissection effects. (XLSX) S4 Table. P-values for biological data. P-values for biological endpoints with statistical significance between sites. Biological response information can be found in Figs 2-5 and Table 4. All biological data used in this study including: fish number, sample year, sample river, the sample river reach section, identification of caged or resident fish, genus, common species name, total length (mm), total weight (g), gonadal weight (g), liver weight (g), condition factor (CF), gonadal somatic index (GSI), hepatosomatic index (HSI), plasma vitellogenin concentration (μg/mL), log 10

Acknowledgments
For support during fieldwork we thank US Fish and Wildlife Service (USFWS) personnel at field offices across the Laurentian Great Lakes states. We appreciate the help from students in the Aquatic Toxicology Laboratory at St. Cloud State University. Dr. Glenn Johnson at the State University of New York at Potsdam graciously provided laboratory space for fish processing in 2013 and again in 2014. We thank Dr. Michael Penn and others at the USFWS Northeast Fish Health Lab (Lamar, PA) for their assistance in processing fish from the Raquette River in both years. The Cuyahoga National Park Service provided laboratory space at their headquarters in both years of the study. The Ohio EPA and Akron Water Reclamation Facility assisted with fish collections in the Cuyahoga River. Thomas and Helen Pahler, Norwood, NY, are thanked for permission to access the Raquette River through their property. The U.S. Fish and Wildlife Service Alpena Fish and Wildlife Conservation Office Waterford Substation assisted with wild fish collection in the Clinton River. In-kind support from the Metropolitan Water Reclamation District of Greater Chicago is acknowledged. We would also like to thank the Editor and Reviewers for their thoughtful, constructive, and in-depth comments on all