Feminization and masculinization of western mosquitofish (Gambusia affinis) observed in rivers impacted by municipal wastewaters

Municipal wastewaters have been known to contain various estrogens and androgens. Little is known about the joint action of these chemicals from wastewaters on fishes in the aquatic environment. The objectives of this study were to investigate the estrogenic and/or androgenic effects in wild mosquitofish (Gambusia affinis) of two effluent-impacted rivers in South China by determining morphological changes and hepatic mRNA expression levels of relevant genes such as vitellogenin (Vtg), estrogen receptor (ERα) and androgen receptors (ARα and ARβ), and to assess the linkages of those morphological changes and hepatic mRNA expression levels to the chemical concentrations measured by in vitro bioassays and chemical analysis. The results showed a significant induction of Vtg and ERα mRNA in the livers of the males and a gonopodium-like anal fin in the females collected at the majority of sites. Redundancy analysis and Pearson correlation analysis showed that the chemical concentrations obtained by in vitro bioassays and chemical analysis had significant correlations with some of the endpoints for the estrogenic and/or androgenic effects in mosquitofish. The findings from this study indicate that the estrogens and androgens present in the two rivers could cause the observed estrogenic and androgenic effects in mosquitofish.

In addition, both E2 equivalent (EEQ) and DHT equivalent (DEQ) in each water sample were calculated using the relative potencies and environmental concentrations of the chemical concentrations based on the addition model ( Table 1). The maximum calculated EEQ and DEQ values (i.e. CEEQ and CDEQ) were 16.0 ng/L and 2280 ng/L, respectively.
Estrogenic and androgenic activities by in vitro bioassays. The estrogenic activity of water samples measured by YES bioassay (MEEQ) ranged from not detected (3% of samples) to 82.6 ng/L ( Table 1). The average estrogenic activity of all water samples was 12.6 ng EEQ/L. The androgenic activity of water samples measured by YAS bioassay (MDEQ) varied in the range of not detected to 46.0 ng DHT/L, with a mean of 9.29 ng DHT/L. However, only 42% of all collected water samples were above the limit of detection for the YAS bioassay (Table 1).
Morphological characteristics of mosquitofish. In general, the changes of P/D ratios (P: perpendicular distance from the point of attachment of the hemal spine to its distal tip; D: perpendicular distance from the distal tip of the hemal spine to the center of the vertebral column) were more obvious than those of L/D ratios (L: total hemal spine length; D: perpendicular distance from the distal tip of the hemal spine to the center of the vertebral column), 3/4W ratios (width ratios of ray 3 and ray 4) and 4/6L ratios (length ratios of ray 4 and ray 6) in the collected mosquitofish from the study area used for the morphological analysis (Figs 2 and 3). Increased P/D ratios in female mosquitofish collected from the Shima River and Danshui River (except for S1) were observed when compared with those in the reference site (S0), and significant differences were found at almost all sampling sites ( Fig. 2A,B). A majority of female mosquitofish from the sampling sites exhibited a masculinized anal fin and hemal spines (Figs 2A,B,E,F and 3A,B), suggesting masculinization of female mosquitofish in the rivers.
The P/D ratios were significantly decreased in male mosquitofish at the majority of sampling sites compared with the reference site (S0) (Fig. 2C,D). A majority of feminized male mosquitofish from the sampling sites exhibited a feminized gonopodium and hemal spines with female characteristics (Figs 2C,D,G,H and 3C,D), indicating feminization of male mosquitofish in the two rivers. In addition, ovotestes were also found in masculinized adult females from some sampling sites such as S2 and S6 by gonadal histopathology (Supplementary Fig. S2G and H).

Reproduction-related gene expression in mosquitofish.
Hepatic Vtg mRNA expression in female and male mosquitofish was significantly elevated with a coincident significant increase in ERα mRNA expression at the majority of sampling sites compared to the reference site S0 (Fig. 4A-D). In general, induction of Vtg and ERα mRNA expression was more evident in the males than females at most sampling sites and periods. Significant induction in the relative mRNA expression of ARα and ARβ were detected at some sampling sites ( Fig. 4E-H). Moreover, significant decreases in hepatic ARα and ARβ mRNA expression in male mosquitofish were also observed at sites S2 or S3 from the Shima River. Correlations between morphological and genetic responses in mosquitofish and hormonal activities in water samples. Pearson Table S3 and Supplementary Table S4). The results of redundancy analysis (RDA) for various measured data are displayed in Fig. 6. Some estrogenic and androgenic compounds were not included in the RDA analysis because they had high linear correlation with other compounds and showed high variance inflation factors (VIF > 10). The VIFs of the compounds chosen for RDA were reasonably low (from 1.4 to 7.7). The first ordination RDA axis (horizontal) was mainly correlated to MEEQ, CEEQ, 4-NP and 4-t-OP and explained 58.3% of the variation in the estrogenic and androgenic activities (82.8% of their relation to compounds). Vtg and ERα mRNA expression in the females and males as well as the morphological characteristics in the males were found to be strongly influenced by 4-NP and 4-t-OP. The second ordination axis (vertical), which was strongly associated with CDEQ, ADD, 17α -BOL and T, accounted for 6.9% of the variation (9.7% of their relation to compounds). ADD, 17α -BOL and T were found to have a major influence in the anal fin and hemal spines in the female mosquitofish. It should also be noted that both Pearson correlation analysis and RDA analysis showed many of the morphological traits were correlated with each other (Figs 5 and 6).

Discussion
The present study made good use of both chemical analysis and bioassays to identify toxicants suspected of causing endocrine disrupting effects in fish. When the results of the morphological indices and reproduction-related gene mRNA expression levels in mosquitofish are taken together, there is sufficient evidence to suggest that mosquitofish in the Shima River and Danshui River experienced strong estrogenic and androgenic effects. The results from the present study also demonstrated that estrogens and androgens in the two rivers may be associated with the observed estrogenic and androgenic effects in mosquitofish.
Previous studies have shown shorter gonopodia (feminization) in male mosquitofish in response to exogenous estrogens 20 , and this is accompanied by a general delay in the development of hemal spines 14,19 . In the present study, male mosquitofish from the Shima River and Danshui River had a reduction in the development of gonopodia and hemal spines, as demonstrated by the decreased ratio of 4/6L, 3/4L 14P/D, 15P/D, or 16P/D. It should be noted that many of the morphological traits are correlated with each other. A similar result was also reported in our previous study, which showed estrogenic effects of male mosquitofish in the sewage-contaminated Hanxi River of South China based on the morphology of hemal spines 14 . It has been demonstrated that gonadal secretions are necessary for gonopodium development in male mosquitofish 21 . The observed shorter gonopodia in The length ratio of ray 4 and ray 6 (4/6L) and width ratio of ray 3 and ray 4 (3/4W) were used as the degree of the change of anal fins. Significant differences between sampling sites and reference site (S0) are indicated by an asterisk (p < 0.05).
mosquitofish from two rivers in the present study indicates the alteration of endocrine function. In fact, delayed development of gonopodia and hemal spines in the males was also in accordance with the induction of Vtg mRNA expression in females and males at the study sites.
Previous studies about estrogenic effects have already demonstrated that the mRNA expression of hepatic Vtg can be induced in male fish exposed to estrogens, which is normally present at undetectable or low levels in males 15,22 . The results of the present study demonstrate that mosquitofish at the majority of sampling sites experienced transcriptional feminization as evidenced by the significant induction in Vtg and ERα mRNA expression. More importantly, Vtg mRNA expression levels in the males from some study sites reached or even exceeded those in females. High Vtg mRNA expression levels in mosquitofish can be associated with higher incidences of testicular malformation, higher amounts of oocyte malformation, lower sex steroid concentrations, and kidney  First axis is horizontal, second axis is vertical. F and M in front of the endpoints of mosquitofish represent female and male, respectively; MEEQ and CEEQ represent measured EEQ by YES and calculated EEQ by chemical analysis, respectively; MDEQ and CDEQ represent measured EEQ by YAS and calculated EEQ by chemical analysis, respectively; BPA: bisphenol-A; 4-NP: 4-nonylphenols; ADD: androsta-1,4-diene-3,17-dione; AED: androsterone (ADR), 4-androstene-3,17-dione, EADR: epi-androsterone, T: testosterone. The angles among arrows denote the degree of correlation between the individual variables, and the smaller the angle, the larger the correlation. In addition, positively correlated variables are shown as arrows pointing in the same direction, negatively correlated variables pointing in opposite directions. dysfunction 23 . Given the relationship between Vtg and ERα , a consistent effect on both genes at most study sites of the present study is not surprising.
Concentrations of some estrogens and estrogenic activities (EEQ) measured in the present study showed strong correlations with those morphological and molecular parameters related to estrogenic effects in mosquitofish. Both Pearson correlation and RDA showed that estrogenic compounds such as NP and 4-t-OP in the rivers could have contributed to the induction of Vtg and ERα mRNA expression in mosquitofish. Previous laboratory studies have shown that NP and 4-t-OP have the potential to induce Vtg in males, form intersex gonads in males, and cause other signs of reproductive impairment 24,25 . In addition, NP and 4-t-OP have been reported to be the primary cause of feminized effects downstream of industrial wastewater discharges 26 . However, the contribution of natural and synthetic estrogens cannot be discarded since some of them were present in the sampling sites and others could be under our detection limit but they may still be high enough to contribute to estrogenic effects in fish.
Likewise, the measured in vitro estrogenic activities in the rivers of the present study support the linkage of estrogenic compounds to the observed morphological and genetic responses related to estrogenic effects in male mosquitofish. Therefore, those estrogenic compounds could be contributing to the observed estrogenic effects in mosquitofish in the rivers impacted by municipal wastewater.
The anal fins and hemal spines of female mosquitofish can serve as biomarker for androgen exposure, either in the laboratory or in the environment 13,14 . Female mosquitofish from the Shima River and Danshui River were found to have masculinized anal fin and hemal spines at almost all sampling sites, despite no significant morphological changes in a few sites when compared to the reference site. It seems that the hemal spines in the females were more sensitive as morphological biomarkers than the anal fin rays. In addition, ovotestes were also found in masculinized adult females from some sampling sites. More interestingly, masculinized female mosquitofish showed significant increase in Vtg and ERα mRNA expression at some sampling sites such as S2 andS3, suggesting the simultaneous presence of super-feminization at the transcriptional level in masculinized female mosquitofish. To date, however, the formation mechanism of this phenomenon is still not known. According to the RDA and Pearson correlation analysis, androgens such as ADD, T, and EADR could be responsible for morphological masculinization in the females. The morphological biomarkers were also found to be better correlated to the CDEQ than MDEQ values. This could be attributed to matrix interferences often experienced during in vitro YAS bioassay and other factors as reported by others 17,27 . Thus we should be cautious to use in vitro YAS bioassay data to assess androgenic effects in fish 28 .
The action of androgens is classically mediated by androgen receptors (ARs), which act as ligand-dependent transcription factors, controlling the transcription of target genes. Therefore, any disruption in the signaling of ARs may lead to impairment of genomic pathways and downstream processes 15 . Changes in ARα and ARβ mRNA expression levels in mosquitofish from the study sites did not have a consistent pattern with an increase, decrease, or no change, which may be attributed to a complex feedback of the AR-mediated gene expression in fish 1,29 . The differential expression currently observed in ARα and ARβ mRNA expression in the female mosquitofish within different sampling sites was further reflected by the weak correlations to the calculated and measured DEQ values in the rivers. This may suggest the presence of other androgenic compounds in the rivers; but more likely this may suggest that ARs are not good biomarkers for mosquitofish exposed to androgens as found in previous studies 15,29 . Although it is still unclear whether these transcription-level effects of ARs in mosquitofish as a consequence of exposure to municipal wastewater are translated into morphology-level abnormalities, abnormal AR mRNA expression levels in wild mosquitofish are undeniably a dysfunctional AR signaling.
In summary, the findings in this study showed strong estrogenic and androgenic effects co-occurring in mosquitofish from the Shima River and Danshui River, and good correlations between the analyzed EDCs levels (and hormonal activities) in the rivers and the morphological effects observed in mosquitofish. Since these chemicals are mainly derived from municipal wastewaters, the discharge of municipal wastewaters into the rivers is sufficient to cause estrogenic and androgenic effects in mosquitofish, resulting in feminization of male mosquitofish and masculinization of female mosquitofish. Toxicity identification evaluation (TIE) and effects directed assay (EDA) should be performed in future studies to identify specific toxicants causing endocrine disrupting effects. Further research is also needed to explore the population level effects in mosquitofish in rivers impacted by municipal wastewaters.

Methods
Study area. Shima River and Danshui River were selected as the study areas since they are located in a rapidly urbanized Pearl River Delta region of South China (Fig. 1). Among the 10 selected sampling sites, 5 sites were on the Shima River (S1-S5), and 5 sites on the Danshui River (S6-S10), with S1 used as the reference site as it is located in the upstream section of the Shima River with little human activity. Another reference site (S0) was also included, and it is located in the upstream of the Liuxi River in the region with little human activity. Both reference sites are located at the outlet of reservoirs.
Fish sampling. Mosquitofish (Gambusia affinis) were captured by electrofishing and/or netting from the selected sites in July 2012 (wet season) and December 2012 (dry season). Surface water samples were also collected simultaneously with the fish. During the sampling, water quality parameters were measured according to standard methods 30 , and these data can be found in Supplementary Table S1. All the collected fish were kept alive in plastic storage boxes filled with river water from their respective collection sites and aerated to ensure sufficient oxygen supply. Once in the laboratory, fish were sorted into sexually mature females, immature females, gravid females, mature males, and immature males. The criteria to determine whether fish are mature were based on previous studies 13,14 . In order to avoid sampling bias, only mature females and males from each site were randomly selected for further measurement. A total of 1157 female mosquitofish and 1143 male mosquitofish were used for standard length and body weight. The numbers, mean standard length and mean body weight of mosquitofish in each site can be found in Supplementary Table S2. The bodies of a total of 501 females and 496 males from all sites (20)(21)(22)(23)(24)(25)(26)(27) of each sex at each site) were preserved in 70% ethanol for determination of morphology. The livers of 15 females and 15 males from each site were preserved in RNAlater (Sigma) for determination of target gene expression. All procedures with fish described in this study were performed in accordance with the 2004 Guidelines (Use of Fishes in Research 2004) 31 . They were also approved by the Animal Care and Use Committee of South China Agricultural University.
Water sample collection and extraction. Surface water samples were collected in 1 L amber glass bottles from each site. The collected water samples were transported in coolers to the laboratory and stored at 4 °C, and then processed within 48 h. Six replicates of the surface water samples were collected from each site. Three replicates were used for YES and YAS bioassays, whereas the other three replicates were spiked with the internal standards for chemical analysis of estrogens and androgens. The water samples for YES and YAS bioassays were extracted by HLB cartridges (Waters Oasis 6 mL, 500 mg) according to the method reported by Zhao et al. 32 . The target estrogenic and androgenic compounds in water samples were extracted using solid phase extraction according to our previous method 5 . A pentafluorobenzoyl chloride derivatization method was applied to the quantification of estrogens in the extracts prior to analysis by gas chromatography-mass spectrometry with negative chemical ionization (Agilent 6890N/5975B, USA) as described by Zhao et al. 32 . The androgens in the extracts of water samples were measured by rapid resolution liquid chromatography-electrospray ionization tandem mass spectrometry (Agilent 1200 LC-Agilent 6460, USA), and its detailed operating method is detailed in the previous study 5 . All data obtained from the analysis were subject to strict quality assurance and quality control (QA/QC) procedures. The limits of detection and quantification are given in Table 1.
In vitro bioassays. The YES and YAS bioassays were based on the methods as described by Zhao et al. 33 .
Methanol was used as the blank control; E2 was used as the positive control of YES with an initial concentration of 0.2 μ M; and DHT was used as the positive control of YAS with an initial concentration of 2 μ M. The levels of estrogenic and androgenic activities were then expressed as E2 equivalent (EEQ), and DHT equivalent (DEQ), respectively.
Characterization of basic fish parameters and morphology. Fish were killed in an ice bath, and their standard length (from the snout to beginning of the caudal fin) and body weight were measured. The method for preparation of the anal fin and skeletons is described by Xie et al. 14 . Briefly, preserved mosquitofish were soaked in purified water to rehydrate and then placed in 1% potassium hydroxide to remove the soft tissue. The anal fin and skeletons of mosquitofish were photographed using a camera mounted on a stereomicroscope.
The length ratio of ray 4 and ray 6 (4/6L) and width ratio of ray 3 and ray 4 (3/4W) were used to assess impacts to anal fin development. Morphological measurements of hemal spines were performed according to the method developed by Rawson et al. 19 . These measurements included total spine length (L), perpendicular distance from the point of attachment of the hemal spine to its distal tip (P, positive where the spines were anteriorly directed and negative where the spines were posteriorly directed), and perpendicular distance from the distal tip of the hemal spine to the center of the vertebral column (D) on the 14th, 15th and 16th hemal spines (labeled as 14L, 15L and 16L; 14P, 15P and 16P; and 14D, 15D and 16D) (Supplementary Fig. S1). In the present study, the ratios of 14P/D, 15P/D, 16P/D, 14L/D, 15L/D and 16L/D were used to assess the elongation and the anterior bending of the hemal spines on the 14th, 15th and 16th vertebrae in mosquitofish.
The mRNA expression of reproduction-related genes in mosquitofish. Total RNA was extracted from the pooled livers of five female or male fish at each sampling site (three replicate pooled liver samples per site) using Trizol reagent (Invitrogen) according to the manufacturer's instructions. First-strand cDNA synthesis was performed with 1 μ g of total RNA using ReverTra Ace ® qPCR RT Master Mix with gDNA Remover (Toyobo, Japan) in a total volume of 20 μ L according to the manufacturer's instructions.

Statistical analyses.
Data on morphological and genetic responses in mosquitofish are presented as mean ± standard deviation (SD) in each site. All physiological endpoints and gene expression patterns were assessed for normality and homogeneity of variances using Kolmogorov-Smirnov and Levene's tests, respectively.
Raw data were log transformed to meet the assumptions of one-way analysis of variance (ANOVA). Statistical differences between sampling sites for morphological and genetic responses in mosquitofish were determined by ANOVA, followed by Dunnett (homogeneous variance) or Dunnett's T3 (heterogeneous variance). Differences were considered statistically significant when p < 0.05. Correlations between morphological or genetic biomarkers in mosquitofish and hormonal activities in water were analyzed using Pearson correlation and redundancy analysis (RDA). For Pearson correlation analysis, the data were log transformed to render them more nearly normal. RDA analysis was selected according to an initial detrended correspondence analysis (DCA). DCA revealed that the data of the estrogenic and androgenic responses exhibited a linear, rather than a unimodal, response to estrogens and androgens, thus RDA as a linear constrained ordination method was chosen for data analysis. In RDA, the estrogenic and androgenic responses were used as "response variables", and the ordination axes were constrained to be linear combinations of the estrogenic and androgenic compounds (explanatory variables). All data applied in RDA were log transformed. The Monte Carlo permutation test (499 permutations) was performed to determine the significance of the estrogenic and androgenic compounds in accounting for the variance of the estrogenic and androgenic responses by assessing p-values. All data analysis was performed with software SPSS 13.0, Origin 7.5 and Canoco 4.5 for Windows.