Effects of pharmaceutically active compounds (PhACs) on fish body and scale shape in natural waters

Ethics statement

This study followed all relevant national and international guidelines concerning the care and welfare of fish (Algers et al., 2009; Johansen et al., 2006). Fish samplings were authorized by the Minister of Agriculture (Permit no.: HHgF/298-1/2016) and fish collection for laboratory examinations was authorized by the Government Office of Pest county (Permit no.: XIV-I-001/2302-4/2012). During sampling, an effort was made to minimize the suffering of fish and all fish were anaesthetized with a lethal dose of clove oil after collection. No endangered species (according to the IUCN Red List of Threatened Species v. 13 (www.iucnredlist.org) and National Law Protected (http://www.termeszetvedelem.hu/)) were caught during this study.

Study area

The study was performed in the suburban area of Budapest, which is the capital and the biggest city in Hungary and in the Carpathian Basin. Altogether, 22 points were sampled for chemical analysis during 2017–2018, and 420 specimens of three species (140 roach Rutilus rutilus, 180 chub Squalius cephalus, 100 gibel carp Carassius gibelio) were collected in 20 sampling points from 10 streams during 29 sampling occasions (Fig. 1). Body-and scale-shape data of 20 specimens/sites were included in the analyses, the number of sampling sites, where the necessary number of specimens were available has been indicated in Table 1.

Water sampling and chemical analysis

Water samples were taken during low water-level periods. General water chemical analysis was performed in the field (Hanna HI 98194 for dissolved O₂, electric conductivity, pH, total dissolved solids, temperature; Macherey-Nagel VisColor PF12 spectrophotometer for NO₂⁻, NO₃⁻, NH₄⁺, PO₄³⁻). For further laboratory analyses (F⁻, Cl⁻, SO₄²⁻, NO₂⁻, NO₃⁻, PO₄³⁻, NH₄⁺, Ca²⁺, Mg²⁺, Na⁺, K⁺) samples were collected in 500-ml borosilicate glass containers. Samples for total organic carbon (TOC) measurements were taken in white, borosilicate containers (50 ml sample with 500 µl 2M hydrochloric acid (VWR International, Monroeville, PA, USA)). For the elemental analysis, a 10-ml water sample was filtered through a 0.45 µm diameter syringe filter, into polypropylene centrifuge pipes free from metal pollutants, and 100 µl NORMATOM nitric acid (VWR International, Monroeville, PA, USA) was added. TOC and total nitrogen (TN) concentrations were measured by using a Multi N/C 3100 TC-TN analyzer (Analytik Jena, Germany). For the determination of anions (F⁻, Cl⁻, SO₄²⁻, Br⁻, NO₃⁻) and cations (NH₄⁺, Ca²⁺, Mg²⁺, Na⁺, K⁺), a Dionex ICS 5000+ dual channel ion chromatograph (Thermo Fisher Scientific, Waltham, MA, USA) was used. PO₄³⁻, NO₂⁻ concentrations, alkalinity as well as total hardness were measured by standard titrimetric and spectrophotometric methods (Eaton et al., 2005). The concentration of heavy metals was determined by using PlasmaQuant MS Elite inductively coupled plasma mass-spectrometer (Analytik Jena, Germany).

For the PhACs measurements, a brown borosilicate glass container with Teflon faced caps (Thermo Fisher Scientific, Waltham, MA, USA) was filled with 2 l water sample, into which 2 ml of HPLC purity formic acid (VWR International, Monroeville, PA, USA) was added. The samples were immediately stored in 4 °C, and transported to the laboratory in a dark cooler box (Dometic CFX40W) within 4 hours, where they were then extracted.

Details of the sample preparation, extraction and analysis process for PhACs have also been described in our earlier papers (Jakab et al., 2020; Kondor et al., 2020; Maasz et al., 2019). Briefly, for sample quantification, the water samples were acidified with formic acid and spiked with the corresponding mass-labelled internal standard (IS). Because of the relatively low concentrations, analytes were isolated by an AutoTrace 280 automatic solid-phase extraction system (Thermo Fisher Scientific, Waltham, MA, USA) using Strata X-CW cartridges (#8B-S035-FCH, Phenomenex, Torrance, CA, USA). To reach the adequate sensitivity, dansyl-chloride was used in the derivatization of steroid agents. A supercritical fluid chromatography (ACQUITY UPC2 system, Waters) coupled with tandem mass spectrometry (MS/MS) (Xevo TQ-S Triple Quadrupole, Waters) was used to analyze and quantify the selected drug residues. Data were recorded by MassLynx software (V4.1 SCN950) in triplicates using TargetLynx XS software for evaluation. The compound separation was performed on an ACQUITY UPC2 BEH analytical column (#186007607, Waters) with 3.0 mm × 100.0 mm, 1.7 μm particle size.

Fish sampling

Fish were caught by electrofishing, and all sampling was undertaken based on the EU Water Framework Directive (EU WFD) (European Commission, 2009) and Hungarian Biodiversity Monitoring System protocols (www.termeszetvedelem.hu). Sampled watercourse sections belonged to River1 (bed width under 5 m, water depth <1 m) and River2 (bed width over 5 m, water depth <2 m) categories, therefore a battery-powered electrofishing device (HANS-GRASSL IG200/2) was used, with a 150-m section length wading in the water upstream. Two watercourses belonged to the River3 (bed width under 30 m, water depth >2 m) category; therefore an aggregator-powered electrofishing device (HANS-GRASSL EL63II) was used, with a 300-m section length leading from a rubber boat going downstream. At every sampling point, 20 specimens comprised of common fish species (not endangered and not protected) were euthanized by using clove oil and stored at −20 °C.

Environmental characterization of sampling sites

The most important environmental variables were recorded at two levels: local level and landscape level (Table 2). The two levels of environmental variables were analyzed separately.

Morphometric analysis

For body morphometrics, after defrosting, a high resolution digital picture was taken of the left side of all specimens using a NIKON D7200 DSLR camera, with a AF-S NIKKOR 35 mm 1:1.8G objective, to avoid variability of side-effects (Takács et al., 2018). Standard length and wet weight were measured with an accuracy of 1 mm and 0.1 g, respectively. Sex was determined by dissection, after the digital photo was captured (Table S1). Five well-developed scales were removed from every individuals’ left side from the flank. Scales were placed between glass slides and scanned using an upper-light scanner (EPSON Perfection V850 Pro) with high resolution (2,400 dpi). One scale per specimen was used for the analysis. Body and scale shape were analyzed using landmark-based geometric morphometry (Zelditch et al., 2004). Ten landmarks were placed on fish body and seven landmarks on fish scales (Fig. 2). For further multivariate analysis, we used the MorphoJ software package (Klingenberg, 2011). To derive shape variables from the raw landmark coordinates, a generalized least-squares Procrustes superimposition was applied to scale, translate and rotate the coordinates (Rohlf, 1990). To eliminate the variances associated with allometric growth, a regression analysis was performed between the logarithm of centroid sizes and the Procrustes coordinates. The regression residuals were used for further analysis (Zelditch et al., 2004). The Procrustes distance (Pd) was used in Canonical Variates Analysis (CVA) for computing group differences, and permutations tests with 1000 iterations were performed to test for significance.

Ecological risk assessment

Ecological risk characterization for PhACs is usually performed by calculating and categorizing a risk quotient (RQ). RQ is a ratio of MEC/PNEC, in which PNEC (predicted no effect concentration) is the estimated highest concentration of an individual PhAC not affecting the aquatic ecosystem, and MEC is the maximum measured environmental concentration in the studied surface water. In general, RQ < 0.01 refers to a negligible risk, 0.01 < RQ < 0.1 denotes a low risk, 0.1 < RQ < 1 indicates a medium risk, while RQ > 1 represents a high risk to the aquatic ecosystem.

Predicted no effect concentration (PNEC) derives from the ratio of available ecotoxicological data (e.g., NOEC, EC50, LC50, HC5) and an assessment factor (AF). When the PNEC value was not available in the literature, we used a selected ecotoxicological data/AF quotient keeping in mind the priorities between the raw data (e.g., applying experimental results instead of extrapolated modelled data, and chronic outcomes in place of acute test results). The magnitude of the AF varies between 1,000 and 5, and it depends on the available ecotoxicological information. The uncertainty (i.e., AF) of the data decreases with expanding of the relevant data set. If PNEC can be calculated only based on acute test results, then AF = 1,000. If PNEC can be derived from chronic data of a species, then AF = 100. Its value further decreases if ecotoxicological chronic test results are available at multiple different trophic levels: AF = 50 (two levels) or AF = 10 (three levels). If PNEC can be determined knowing of hazardous concentration for 5% of species investigated (HC5 based on ecotoxicological results of at least five species), then AF = 5. When data are available for each trophic level, the lowest concentration was selected to determine PNEC since environmental risk assessment is based on the most sensitive elements of the ecosystem (Molnar, Maasz & Pirger, 2020). PNECs with raw ecotoxicological data and AFs are presented in Table S2.

Statistical analysis

Background variables were categorized into four groups: PhAC data, general water chemistry data, local environmental variables data and landscape-scale environmental variables. All variables were numeric and log10 transformed before further analyses. An unconstrained Principal Component Analysis conducted on the shape datasets (x and y coordinates of the regression residuals) was followed by the passive projection of the explanatory variables. The number of permutations in a Monte-Carlo simulation were set to 1,000. In the first model, body shape data, while in the second model, scale shape data, were used with all the environmental variables listed in the dataset. Where forward selection revealed significant effects, variance partitioning was used to assess the relative contribution of the different variable groups (Borcard, Legendre & Dapeau, 1992). Additional Mantel tests were performed on shape-variables (Mahalanobis and Procrustes distances) and PhACs concentrations, to assess the site-specific component of differences.

PhAC data from sampling points

Altogether 54 different types of PhACs were found in the water samples from the sampling points (Table 3). Three compounds were detected in a µg/l concentration range in examined samples, lamotrigine (maxMEC = 14 338.3 ng/l), caffeine (maxMEC = 13 635 ng/l), and diclofenac (maxMEC = 2 201.7 ng/l). The remaining 51 PhACs were measured in a few hundred, a few tens, or a few ng/l concentration ranges each above the limit of detection. A total of 27 PhACs were used in analysis based on their RQ-values, from which eight showed high, eight showed medium and the remaining eleven PhACs received a low risk classification based on the environmental risk assessment (Table 3). To perform the risk assessment using relevant ecotoxicological data, we used the AF and PNEC values of detected PhACs (see Table S2).

Morphometric analysis

Significant differences were found between the average shape of fish stocks in all three species based on both fish body- and scale shape. In the case of roach body-shape, the differences based on stream, as well as in scale shape (Fig. 3), significant differences and Pd-values are shown in Table 4 for body shape and Table 5 for scale shape.

Sampling points of Tápió Stream were discriminated from the others (Szent László Stream, Gerje Stream) along the first axis of CVA, according roach body shape. Significant differences were observed between GERTOS and every other points, based on Hotelling’s t-test (Fig. 3; Table 4). SZEBIC has been differed significantly only from TAPTAP. Scale shape of roach proved to be different in TAPUJS, than most of other sites.

In the case of chub body- and scale shape, there were no clear connections found with the stream (Fig. 4); significant differences and Pd-values are shown in Table 6 for body shapes and Table 7 for scale shapes. Figure 4 suggests negative correlation between the distance from the estuary and CV2 (HOSTOR < HOSKEL < HOSKAM; BUKTOR < BUKSZE < BUKIZB) in case of Hosszúréti Stream and Bükkös Stream also, however CVA-plot for scale shape not support this finding. In the case of gibel carp body shape, all sampling points differed significantly. In the case of gibel carp scale shape, there was a connection with stream, but there are similarities between the sampling points from different streams as well (Fig. 5); significant differences and Pd-values are shown in Table 8 for body shape and Table 9 for scale shape. Interesting pattern of sites could be observed in case of gibel carp body shape, since within-stream difference (GERTOR–GERCEG) seems to be higher than between-stream (GERTOR–SZEBIC; GERTOR–HOSKEL) difference. Regarding gibel carp scale, GERTOR site have not been differed such harshly from others, like in case of body shape. BENBIA proved to be the most different site along CV1.

Significant background variables

Numerous significant background variables were found, which affect fish body shape and scale shape. Local- and landscape-scale environmental variables, water chemistry data and also PhACs were found to be significant. In case of roach scale shape, the significant variables were As (9%) and SO₄²⁻ (3%), and for body shape, TRIM (6%) and CITA (4%) were found to be significant (1% joint effect). In the case of chub scale shape, water chemistry data (significant variables: Mg, As, Ca) was responsible for 5% of the variance, local environmental variables (significant variables: emergent macrophytes, water depth) were responsible for 2% of the variance, while PhACs (significant variable: CODE) were responsible for 1% of the variance. The local environmental variables and CODE had 1% joint effect. In the case of chub body shape, only two variables were significant, Cd as water chemistry data and detritus as a local environmental variable, for 4% and 3% respectively, with 8% joint effect. In the case of gibel carp scale shape, the water chemistry variable Pb (2%) and the landscape scale environmental variable wetland (6%) were significant, with 1% joint effect. For gibel carp body shape, three different type of variables were significant, the PCB PROP, the water chemistry variable Zn, and the landscape-scale environmental variable catchment size, for 6%, 11% and 2% respectively, with 4% joint effect for Zn and catchment size (Table 10).

Mantel tests did not show significant correlation among the site-specific shape variables and the significant background variables, in most of the cases (Table S3). In case of chub scale, Ca shows significant correlation with Procrustes distances, although in case of Mahalanobis distances the correlation was not significant. In case of roach scale, both As and SO₄²⁻ showed significant correlation with Mahalanobis distances, although in case of Procrustes distances the correlation was not significant.

Potential effects of environmental variables on shape

In the case of chub and gibel carp, significant environmental variables were found. The effects of local (section) level variables on chub scale shape could be explained by the life-history characteristics of the species. Different environmental characteristics of the given habitats may cause changes at the population level (Haas, Blum & Heins, 2010). Coverage of emergent macrophytes, water depth and the quantity of detritus were previously found to be connected to the life history parameters of chub (Bolland, Cowx & Lucas, 2008; Ünver & Erk’akan, 2011), therefore these variables might affect the scale and body shape of the fish. In the case of gibel carp, significant environmental variables included landscape-scale variables, wetland (scale shape) and catchment size (body shape). There are several known examples regarding the shape-modification effects of environmental differences in fish. Species of the genus Carassius are characterized by a high level of phenotypic plasticity. In the case of crucian carp (Carassius carassius), the presence or absence of predators and the feeding behavior (zooplankton versus benthic chironomids) have a complex effect on body shape (Andersson, Johansson & Söderlund, 2006).

Potential effects of general water chemistry on scale shape

Water chemistry had a significant impact on roach and chub scale shape. The effects of arsenic (As) on muscle development in fish have already been reported (D’Amico, 2012), and this compound can accumulated in scales (Allen, Awasthi & Rana, 2004) as well, which might affect scale shape itself. Fliedner et al. (2014) studied the water chemistry, especially the heavy metal concentrations in rivers Rhine, Elbe, Danube, Saar, Mulde, Saale and in Lake Belau in Germany. Throughout the study As, Pb, Cu and Hg concentrations were measured from tissue samples of zebra mussel (Dreissena polymorpha) and bream (Abramis brama). Arsenic found to be the only compound, where increase in concentration was detectable while analyzing in bream muscle tissue samples from 1990s to 2014 (Fliedner et al., 2014). Mg²⁺ and Ca²⁺ significantly impacted the scale shape of chub. Ca²⁺ is an essential building component of fish scales (Sankar et al., 2008) while the Mg²⁺ content of water affects calcium uptake in fish (Dabrowska, Meyer-Burgdorff & Gunther, 1991; Van der Velden et al., 1991). Cadmium is a Ca²⁺ uptake inhibiting agent which was also shown to affect chub body shape. The presence of Cd has a negative effect on Ca²⁺ uptake through the gills (Franklin et al., 2005). Lead concentrations are also connected to gibel carp scale shape formation. This heavy metal cannot be excreted physiologically (via the gills or kidneys), and Pb impairs fish scale development to a greater extent than in other organs (Çoban et al., 2013). Zinc also has a significant impact on gibel carp body shape, and is associated with higher (11%) variance. Zinc uptake is related to Ca²⁺ concentrations where high Ca²⁺ concentrations may decrease Zn uptake; excess Zn then accumulates in fish skin, muscle and bones (Hogstrand & Wood, 1996), and therefore might have an effect on body shape.

Potential effect of PhACs on shape

TRIM is a cytoprotective, anti-ischemic agent with a strong antioxidant effect (Sedky et al., 2017). In zebrafish (Danio rerio) TRIM can decrease the ototoxic effects of neomycin on hair-cell loss in the neuromasts (Chang et al., 2013). Phenotypic alterations have not been discussed previously, however, a significant effect was detected on roach body shape in this study. CITA as a SSRI, have also been shown to significantly affect roach body shape. A strong anxiolytic effect has been reported in fish previously (Olsén et al., 2014; Porseryd et al., 2017), and alterations in behavioral patterns might also affect the phenotype as well, because the use of different habitats might alter the phenotype of different species (Faulks et al., 2015). CODE an opiate derivative, is used to treat rheumatic pain (Ytterberg, Mahowald & Woods, 1998), and significantly modulates chub-scale shape. There is evidence of the presence of codeine in fish tissues (Epple et al., 1993; Valdés et al., 2016), however, phenotypic alterations have not been detected. It might be in relation with the inhibition of the expression of receptors for vascular endothelial growth factor, which can affect the early life-stage development of fish (Karaman et al., 2017). PROP, a non-selective β-blocker, affected gibel carp body shape. It is used to treat heart diseases, and has proved to be the cause of decreased testosterone and estradiol levels in zebrafish, and has showed anxiolytic effects, and decreased growth (Mitchell & Moon, 2016). As we discussed in the case of roach and CITA, the anxiolytic effects of drugs might also alter phenotype. Based on RQ-values, CITA was ranked to be high risk, while CODE and PROP were medium risk and TRIM was low risk. These results also suggest that the widely used “traditional” risk assessment may have weaknesses when compared to a “real-life” measured effects.