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

Background In recent years, there are growing concerns about pharmaceutically active compounds (PhACs) in natural ecosystems. These compounds have been found in natural waters and in fish tissues worldwide. Regarding their growing distribution and abundance, it is becoming clear that traditionally used risk assessment methodologies and ecotoxicological studies have limitations in several respects. In our study a new, combined approach of environmental impact assesment of PhACs has been used. Methods In this study, the constant watercourses of the suburban region of the Hungarian capital (Budapest) were sampled, and the body shape and scale shape of three fish species (roach Rutilus rutilus, chub Squalius cephalus, gibel carp Carassius gibelio) found in these waters were analyzed, based on landmark-based geometric morphometric methods. Possible connections were made between the differences in body shape and scale shape, and abiotic environmental variables (local- and landscape-scale) and measured PhACs. Results Significant connections were found between shape and PhACs concentrations in several cases. Despite the relatively large number of compounds (54) detected, citalopram, propranolol, codeine and trimetazidine significantly affected only fish body and scale shape, based on their concentrations. These four PhACs were shown to be high (citalopram), medium (propranolol and codeine), and low (trimetazidine) risk levels during the environmental risk assessment, which were based on Risk Quotient calculation. Furthermore, seven PhACs (diclofenac, Estrone (E1), tramadol, caffeine 17α-Ethinylestradiol (EE2), 17α-Estradiol (aE2), Estriol (E3)) were also categorized with a high risk level. However, our morphological studies indicated that only citalopram was found to affect fish phenotype amongst the PhACs posing high risk. Therefore, our results revealed that the output of (traditional) environmental/ecological risk assessment based on ecotoxicological data of different aquatic organisms not necessarily show consistency with a “real-life” situation; furthermore, the morphological investigations may also be a good sub-lethal endpoint in ecotoxicological assessments.


INTRODUCTION
The first detection of pharmaceutically active compounds (PhACs) in aquatic ecosystems and drinking water dates back to the 1980s (Richardson & Bowron, 1985;Watts et al., 1983). Since then, an emerging number of studies have reported the distribution and the potential threat posed by these compounds (Boxall et al., 2012;Datel & Hrabankova, 2020;Dietrich, Webb & Petry, 2002). Selective Serotonin Reuptake Inhibitors (SSRIs), β-blockers and anti-inflammatories are considered to be the most abundant drug residuals occurring in surface waters (Boxall et al., 2012). These compounds can be released into natural waters via several ways. The main sources of pollution are Wastewater Treatment Plants (WWTPs) (after the excretion of human waste) (Subedi et al., 2012), the pharmaceutical industries and the excretion of drugs from animals used in agriculture (Boxall et al., 2012). The recent technologies of WWTPs cannot eliminate these compounds fully from wastewater (Golet et al., 2001;Ternes et al., 1998;Tsui et al., 2014;Yang et al., 2020). To minimize the potential environmental risk posed by PhACs, several regulations for ecotoxicological testing have been enacted (EMEA, 2006). In recent years, several weaknesses of these regulations have been reported in scientific articles (Ankley et al., 2007;Boxall et al., 2012) such as: (1) official tests usually use lethal endpoints, (2) little attention is paid to metabolites, (3) different regulations for human and for veterinary drugs, (4) tests for unique agents, (5) calculating the degradation of compounds and (6) overabundant compounds (over 4.000 drug substances) to test all of them. These weaknesses and the resulting shortcomings in risk assessment procedures may cause uncertainties regarding their validity. If these points are not addressed and alternative, more adequate risk assessment techniques would not added to the regulations, then a false illusion of low risk may result in many cases. Therefore, the current shortcomings need to be examined in detail in order to better understand the problem. It is a well-known fact that several biotic and abiotic factors can influence the body shape of fish, such as food availability (Currens et al., 1989;Marcil, Swain & Hutchings, 2006;Park et al., 2001), food type (Day, Pritchard & Schluter, 1994), temperature (Beacham, 1990;Šumer et al., 2005), and the presence or absence of predators (Brönmark & Miner, 1992). In addition, it has also been proven that environmental parameters can affect the shape of fish scales (Ibáñez, 2015;Staszny et al., 2013;Takács et al., 2016). The effect of basic chemical parameters (e.g., ion concentrations) of the water may also affect the phenotype of fish, however their effect on the shape (body or scale) is unclarified (Çoban et al., 2013;Franklin et al., 2005;Schlenk & Benson, 2001). Due to the chronic, multigenerational exposure of fishes to PhACs, phenotypic alterations are possible, and there is evidence that progestogen contaminations can affect somatic indices (Maasz et al., 2017). Therefore, the aim of this study was (1) to find connections between the PhACs measured in small watercourses and the body and scale shape of selected fish species; and (2) to describe which type of PhACs or abiotic environmental factors are responsible for anatomical differences.

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 2 , electric conductivity, pH, total dissolved solids, temperature; Macherey-Nagel VisColor PF12 spectrophotometer for NO − 2 , NO − 3 , NH þ 4 , PO 3− 4 ). For further laboratory analyses (F − , Cl − , SO 2− 4 , NO − 2 , NO − 3 , PO 3− 4 , NH þ 4 , Ca 2+ , Mg 2+ , 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  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 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 mm 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

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 2− 4 (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 2− 4 showed significant correlation with Mahalanobis distances, although in case of Procrustes distances the correlation was not significant.

DISCUSSION
Our results indicated that PhACs can influence fish body shape and scale shape in natural environment and habitats. There are several studies that showed shape differences between fish stocks in natural waters (Ibáñez & Jawad, 2018;Takács et al., 2016). These studies usually explain the variations by different genetic background (Lõhmus et al., 2010;Staszny et al., 2013), phenotypic plasticity (Vasconcellos et al., 2008), or some basic environmental differences, such as food availability (Currens et al., 1989 Hutchings, 2006; Park et al., 2001), temperature (Lõhmus et al., 2010;Šumer et al., 2005), flow-regime (Haas, Blum & Heins, 2010). These effects, and their combination have also affected the phenotype of fish included this study. Moreover, the observed impact of PhACs on shape is considered to relatively small, however it should be taken into consideration during the studies, carried out in natural waters. In addition, the results of this study suggest that the mixtures of PhACs that occur in natural waters have different effects on different species and phenotypes such as body and scale.

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  , 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 2+ and Ca 2+ significantly impacted the scale shape of chub. Ca 2+ is an essential building component of fish scales (Sankar et al., 2008) while the Mg 2+ content of water affects calcium uptake in fish (Dabrowska, Meyer-Burgdorff & Gunther, 1991;Van der Velden et al., 1991). Cadmium is a Ca 2+ uptake inhibiting agent which was also shown to affect chub body shape. The presence of Cd has a negative effect on Ca 2+ 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 2+ concentrations where high Ca 2+ 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.

CONCLUSIONS
In summary, our results suggest that PhACs in natural waters can affect the phenotypic characteristics of fish species. Although a relatively large number of PhACs (54 compounds) were found in the water samples, only 4 compounds were found to have significant effects on phenotype. This study did not aim to find clear cause and effect relationships between the given compounds, or to reveal the mode-of-actions; however, the individual-scale effect of PhACs was identified. The results of this study showed that differences in phenotype can be detected, therefore the morphometric analysis was suitable for an alternative, sub-lethal endpoint of environment-level toxicological investigation. However, in order to get a more accurate picture of the actual phenotypic effect of PhACs in the environment, a more detailed study with a larger sample size is needed. Since the effects of PhACs on scale shape have been observed, scale sampling may be a suitable, effective and ethically acceptable tool to extend studies on different river systems.