Biomarker selection depends on gene function and organ: the case of the cytochrome P450 family genes in freshwater fish exposed to chronic pollution

Sampling sites and sample collection

This research was conducted in the spring of 2016 and sampled four sites in the Maipo River Basin, Chile. Two sites were characterized in previous studies (Vega-Retter et al., 2014; Vega-Retter et al., 2018) as historically non-polluted (San Francisco de Mostazal (SFM) (33°58′19.97″S, 70°42′56.49″W) and Isla de Maipo (IM) (33°44′58″S, 70°53′26″W)) and two as polluted (Melipilla (MEL) (33°42′49,988″S, 71°12′39,13″W) and Pelvin (PEL) (33°36′21″S, 70°54′33″W); Fig. 1). In the previous studies, the non-polluted sites have shown good water quality related to a low-density population and low industrial development nearby. In recent years, increasing urbanization and industrial development close to the IM site has decreased the water quality. In contrast, the polluted sites are downstream of wastewater plants and historical industrial water discharges (Gomez, De La Maza & Melo, 2014).

Twenty-four B. microlepidotus individuals were collected by electrofishing, six per sampling site. Sampling was performed post-reproductive period; thus, individual sex cannot be determined using gonads. The mean fish weight was 5.04 ± 0.90 g for the IM site, 1.90 ± 0.61 g for the MEL site, 5.41 ± 2.52 g for the PEL site, and 16.03 ± 7.97 g for the SFM site. The mean total length was 10.03 ± 0.50 cm for the IM site, 7.16 ± 0.69 cm for the MEL site, 10.27 ± 1.57 cm for the PEL site, and 14.27 ± 2.16 cm for the SFM site. The fish were sacrificed immediately by neck-breaking, and the liver and gills were removed in the field and stored in liquid nitrogen for subsequent RNA extraction.

All the protocols used in this study were approved by the Ethics Committee of the Universidad de Chile and complied with existing laws in Chile (Resolución Exenta No. 3078 Subsecretaria de Pesca).

Physical and chemical characterization of the sampling sites

The surface water and sediment were physically and chemically characterized to determine the current pollution level at each sampling site, focusing on micronutrients, macronutrients, and metals as an approximation to the complex mixture of contaminants in the river. Three samples of surface water (1 L) and sediment (1 kg) were taken at each site. The water samples were collected in the water column using Nalgene vials and stored at 4 °C for at most 48 h before analysis. In addition, physical and chemical parameters such as electrical conductivity (EC), temperature, dissolved solids (DS), and pH were measured three times in situ using a multiparameter device (Hanna Instruments, Woonsocket, RI, USA).

Superficial water characterization

For water samples, nitrite (NO₂⁻) is determined by the formation of a reddish-purple azo dye by diazotization-coupling reaction of sulfanilamide with N-(1-naphthyl)-ethylenediamine dihydrochloride (NED dihydrochloride) at pH between 2.0 and 2.5 and was quantified by spectrophotometry at 543 nm (UV-1700; Shimadzu; Kyoto, Japan) (Sadzawka et al., 2006).

In the case of ammonium (NH₄⁺), it was quantified by Indophenol blue method, which is a colorimetric method (Solórzano, 1969).

Boron (B) was quantified by the technique proposed by Berger and Truog (Sadzawka et al., 2006), measuring the absorbance at a length of 420 nm (UV-1700; Shimadzu; Kyoto, Japan).

Phosphate (PO₄³⁻) was quantified by a colorimetric method that is based on the formation of a hetero-polyacid with the vanado-molybdic reagent, and a subsequent measure of absorbance at a length of 880 nm (UV-1700; Shimadzu; Kyoto, Japan).

Carbonate (CO₃²⁻), bicarbonate (HCO₃⁻), chloride (Cl⁻), nitrate (NO₃⁻) and sulfate (SO₄²⁻) were determined by Ion Exchange High Performance Liquid Chromatography (IE-HPLC) Isocratic HPLC Pump with a conductivity detector, using anionic (IC-Pak HJC) columns. Isocratic conditions were used with injection volume 50 μL and mobile phase flow 1.2 mL/min. The mobile phase was composed of concentrated borate/gluconate prepared with 34 g boric acid (Merck, Lebanon, NJ, USA), 23.5 mL gluconic acid (Merck, Lebanon, NJ, USA), 8.6 g lithium hydroxide (Merck p. a., Lebanon, NJ, USA) and 250 mL glycerin (Merck, Lebanon, NJ, USA), diluted to 500 mL with MilliQ deionized water. To 20 mL of this solution, we added 120 mL acetonitrile (Merck HPLC grade, Lebanon, NJ, USA) and diluted to 1 L with MilliQ deionized water. Analytes identification was determined by retention time compared with standards (Merck titrisol, Lebanon, NJ, USA). Anion concentrations were estimated using calibration curves generated with standards. Data quality was monitored by measuring element concentration in procedural blanks and synthetic preparations of deionized water (MilliQ) with analytes (Copaja, Núñez & Veliz, 2014).

For, sodium (Na⁺), potassium (K⁺), calcium (Ca²⁺), magnesium (Mg²⁺) were quantified by atomic absorption spectroscopy (AA-6880; Shimadzu, Kyoto, Japan) after filtering the samples using a cellulose nitrate filter with a 0.45 µm pore diameter (Sartorius, Göttingen, Germany) according to Clesceri, Greenberg & Eaton, 2005. Total solids (TSs) were quantified by taking a 100 mL aliquot of water, evaporating the water until dry, and weighing the resulting solid. DSs were quantified by taking a 100 mL aliquot, filtering it through a membrane with a pore diameter of 0.45 µm, evaporating the water until dry, and weighing the resulting solid. Finally, suspended solids (SSs) were calculated as TS minus DS.

Water samples were fixed with 2% HNO₃ (Suprapur Merck, Darmstadt Germany) to quantify iron (Fe), copper (Cu), zinc (Zn), cadmium (Cd), manganese (Mn), nickel (Ni), aluminum (Al), molybdenum (Mo), lead (Pb), mercury (Hg), arsenic (As), and chromium (Cr) concentrations using an atomic absorption spectrophotometer (AA-6880; Shimadzu, Kyoto, Japan), given the reported effects of these metals on biota.

The dissolved oxygen (DO) concentration was estimated by taking three water samples from each sampling site in 200 mL polycarbonate bottles (Nalgene, Rochester, NY, USA), fixing them with MnSO₄ and alkaline iodide, and analyzing them using the Winkler method according to Strickland (1968).

The oxygen used by microorganisms to decompose organic waste was measured as the biochemical oxygen demand (BOD₅). Briefly, three 300 mL flasks of water from each site were incubated at 20 °C for 5 days before oxygen quantification, as described by Clesceri, Greenberg & Eaton, 2005.

Sediment characterization

At each sampling site, three 1 kg sediment samples were collected from the top 10 cm of the surface using a plastic scoop and polyethylene bag (Copaja & Muñoz, 2018). Samples were stored at 4 °C, subsequently dried at room temperature in polyethylene trays, and then sieved into two fractions: a coarse fraction with a particle size of <2 mm for physical and chemical characterization and a fine fraction with a particle size of <0.063 mm for metal quantification.

The physicochemical variables measured in sediment samples were EC, pH, B, PO₄³⁻, N, Ca²⁺, Mg²⁺, and Na⁺. EC and pH were determined by potentiometric methods (sediment:water = 1:2.5). The Berger & Truog (1939) method was used to quantify B with a spectrophotometer (Pharmaspec 1700; Shimadzu, Kyoto, Japan) at 420 nm. The Kjedhal digestion method (Bremner, 1960) was used to quantify total N. Ca²⁺, Mg²⁺, and Na⁺ were quantified by preparing a saturation extract using 50 g of each sample, adding water to saturate the sample in a 1:1 ratio. The extract was left overnight and then centrifuged to collect the supernatant. Quantification used 10 mL of the supernatant using the same procedure as for superficial water. The same metals quantified in the water were analyzed in the sediment: Cu, Fe, Zn, Cd, Mn, Ni, Al, Mo, Pb, Hg, As, and Cr. The total fraction of each metal was obtained by digesting 0.25 g of sediment with 10 mL of nitric acid (Suprapur; Merck, Darmstadt, Germany) in a high-resolution microwave oven (MarsXpress) using the following conditions: power, 800 W; tower, 100%; time, 11 min; temperature, 175 °C; hold, 15 min; cooling, 15 min. This protocol was based on the US Environmental Protection Agency’s method 3051,Washington, D.C., US (Blakemore, Searle & Daly, 1987). Metal concentrations were determined by flame atomic absorption spectrometry using an atomic absorption spectrometer (AA-6800; Shimadzu, Kyoto, Japan) with an air-acetylene flame.

The sampled sites were categorized and compared with the Vega-Retter et al. (2014) categorization using a principal component analysis (PCA) performed with the fviz_pca_biplot function of the facto-extra package (Langmead & Salzberg, 2017) of the R statistical software (v.4.1.0; R Core Team, 2023). The PCA was performed with all physical and chemical data from the four sampling sites (36 environmental variables) to detect relationships among sites. The PCA results showed that the water quality of the IM site had considerably deteriorated in recent years. Therefore, this site was considered polluted for this study. A second PCA was performed using ten physicochemical parameters (EC, pH, DS, DO, NO₂^-, NH₄⁺, Na⁺, K⁺, Ca²⁺, and Mg²⁺) also included in historical physicochemical surface water data for 2007, 2011, and 2016 for the sampling sites. The second PCA aimed to evaluate the physicochemical stability of the sampling sites.

RNA isolation and RNA-seq

RNA-seq was performed to quantify transcript numbers and determine differential expression between non-polluted (SFM) and polluted (MEL, PEL, and IM) sites. Total RNA was extracted using the PureLink RNA Mini Kit (Ambion; Life Technologies, CA, USA) according to the manufacturer’s instructions and sent to Genoma Mayor Sequencing Services (Santiago, Chile), where it was purified to retain only mRNA. RNA quality and quantity were determined using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Samples with an RNA integrity number >7 (Schroeder et al., 2006) were subjected to 2 × 100 bp sequencing on an Illumina HiSeq 4000 system (San Diego, CA, USA). Four PEL samples were discarded due to RNA integrity: three liver and one gill sample. Twenty-one liver and twenty-three gill samples were sequenced. The raw sequencing data deposited in SRA database of NCBI with the data accessions BioProject ID: PRJNA1033453 and the BioSample accessions: SAMN38033975–SAMN38034018.

Adapters were removed, and raw reads were filtered using Trim Galore (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/), prinseq-lite.pl (http://prinseq.sourceforge.net/manual.html), and Cutadapt (Martin, 2011). We removed reads with (i) low base quality (qbase ≤ 5), (ii) ≥10% ambiguous bases, (iii) mean Phred score < Q30, and (iv) length < 50 bp.

The de novo assembly of gill and liver transcriptomes was performed with the Bridger software (Chang et al., 2015) with the following parameters: k-mer length = 25, minimum k-mer seed coverage = 2, and minimum k-mer seed entropy = 1.5. Isoforms were not considered, and contigs with lengths > 200 bp were retained. Clustering was performed with the CD-HIT software (http://weizhongli-lab.org/cd-hit/) to generate a set of non-redundant contigs. An identity cut-off threshold of 80% was used. Non-redundant contigs from both assemblies were annotated using the Blast2GO software (Conesa et al., 2005). First, a BLAST search was performed using Blast2GO’s blastx function and a subset of vertebrate sequences in the US National Center for Biotechnology Information’s non-redundant database with a threshold e-value of 1 × 10⁻⁶. In addition, the InterPro database was used for annotation. Second, contigs were mapped to identify the gene ontology (GO) terms associated with them, including biological processes, cellular components, and molecular function. Finally, GO terms with an e-value threshold of <1 × 10⁻⁶ were retained for annotation.

Differential expression

Differentially expressed CYP genes were detected by mapping transcripts from each individual against the de novo assembly of each tissue using the Bowtie2 software (Langmead & Salzberg, 2017) with the following command line flags: -q –phred33 –sensitive –dpad 0 –gbar 99999999 –mp 1,1 –np 1 –score-min L,0, −0.1 −I 1 −× 1000 –no-mixed –no-discordant -p 6 -k 200. BAM files were used to estimate expression levels using the RSEM software (v.1.3.3; Li & Dewey, 2011) with the following command line flags: -p 6 –paired-end –calc-pme –calc-ci –ci-memory 30000 –sort-bam-by-coordinate. Finally, differential expression between polluted (MEL, PEL, and IM) and non-polluted (SFM) sites was quantified using the DESeq2 software (Love, Huber & Anders, 2014) implemented as a package in the R statistical software (v.4.1.0; R Core Team, 2023). Three independent runs were performed for each tissue, comparing all possible pairs of polluted and non-polluted sites (MEL-SFM, PEL-SFM, and IM-SFM) and retaining contigs with a |log fold change (Log₂FC)| ≥ 1 and false discovery rate < 0.05. The statistical power of this experimental design was estimated for each organ using the RNASeqPower package (Therneau, Hart & Kocher, 2023) in the R statistical software (v.4.1.0; R Core Team, 2023).

CYP gene identification and classification

Differentially expressed CYP genes were identified based on Blast2GO’s BLAST results. First, all differentially expressed CYP genes in at least one polluted site and one tissue were considered, and their expression pattern was plotted as Log₂FC with the SFM site as the reference. Next, only CYP genes present in de novo liver and gill assemblies whose general function was known were retained. The retained CYP genes were classified according to their action on endogenous or exogenous compounds (Uno, Ishizuka & Itakura, 2012). To test for differential expression of CYP genes, we used a generalized linear model (GLM) analysis considering differential expression as the dependent variable, sites as the independent variable, and binomial negative as the data distribution. The analysis was run on normalized gene counts using the MASS (v.7.3-60) package (Venables & Ripley, 2002) of the R statistical software (v.4.1.0; R Core Team, 2023). Additionally, Tukey’s test was used to identify significant differences between pairs of sites.

Physical and chemical characterization of the sampling sites

Four sampling sites were chosen for this study based on previous studies (Vega-Retter & Véliz, 2014; Vega-Retter, Vila & Véliz, 2015). Two sites had historically been considered polluted (MEL and PEL), and two non-polluted (SFM and IM; Fig. 1). These sites are inhabited by genetically independent B. microlepidotus populations. Thirty-six parameters were used to assess pollution levels at the sampling sites since the concentrations of certain metals in the superficial water (Cu, Zn, Cd, Mn, Ni, Al, Mo, Pb, Hg, and As) and sediment (Cd, Al, Mo, Pb, Hg, and As) were below the detection limit. The parameters used were pH, EC, DO, BOD₅, TS, DS, SS, NO²⁻, NO³⁻, B, PO₄³⁻, Na⁺, K⁺, Ca²⁺, Mg²⁺, NH₄⁺, HCO₃⁻, CO₃²⁻, Cl⁻, SO₄²⁻, Fe, and Cr in superficial water, and EC, pH, B, PO₄³⁻, N, Ca²⁺, Mg²⁺, Na⁺, Cu, Fe, Zn, Mn, Ni, and Cr in sediments. These data were analyzed using a PCA to compare sampling sites.

PCAs of surface water physical and chemical characteristics for 2007, 2011, and 2016 and surface water and sediment physical and chemical characteristics for 2016 at the studied sites are shown in Fig. 2. For the historical physical and chemical data (Fig. 2A), the first two principal components (PCs) explained 62.09% of the total variance. PC1 (41.54% of total variance; eigenvalue = 4.15) segregated the SFM site from the other sampling sites. Three variables had their highest loading (L) in PC1: Ca²⁺ (L = 0.426), K⁺ (L = 0.408), and Na⁺ (L = 0.405). Three variables had their highest L in PC2 (20.55% of total variance; eigenvalue = 2.55): DS (L = 0.644), pH (L = 0.448), and DO (L = 0.369).

For 2016 data (Fig.2B), the first two PCs explained 55.96% of the total variance. PC1 (35.53% of total variance; eigenvalue = 12.79) clearly segregated the SFM site from the other sampling sites (IM, MEL, and PEL). Three variables had their highest L in PC1: DS (L = 0.26), Ca²⁺ (L = 0.249), and TS (L = 0.242). Three variables had their highest L in PC2 (20.44% of the total variance; eigenvalue = 7.35): Cl⁻ (L = 0.29), Mg²⁺ (L = 0.277), and sed_Ca²⁺ (L = 0.251). Values of the main physical and chemical parameters for the historical and 2016 data are shown in Table 1.

Five parameters had values higher than the standards of Decree 53 (Decreto 53, 2014; Ministerio del Medio Ambiente (MMA)) at the IM, MEL, and PEL sites in the 2016 data, and two parameters had higher values than the standards at the SFM site. In both the historical and current data, SFM was segregated from the other sampling sites, showing the temporal stability of the physical and chemical patterns found in the 2016 data.

Gene expression estimation in the liver and gills of the silverside B. microlepidotus

The livers of 21 individuals (3–6 per site) and gills of 23 individuals (5–6 per site) were subjected to RNA-seq after an RNA integrity test, resulting in 713,922,789 total raw reads, of which 368,882,159 and 328,822,670 were retained for the gill and liver, respectively, after filtering and trimming. After clustering, 63,424 and 105,131 non-redundant contigs >200 bp were obtained for the liver and gills, respectively. Functional annotation of these non-redundant contigs resulted in 20,338 and 28,499 annotated contigs for the liver and gills, respectively (Table 2).

Differential expression

RNASeqPower was used to assess the statistical power of this experimental design. The power was estimated to be 0.99 for the liver samples, with an average coverage of 59× and a variation coefficient of 0.20, and 0.90 for the gill samples, with an average coverage of 30× and a variation coefficient of 0.32.

Based on the PCA analysis, three comparisons were performed to identify differentially expressed genes (DEGs): MEL vs. SFM, PEL vs. SFM, and IM vs. SFM. The MEL vs. SFM comparison identified 2,245 and 1,166 DEGs for the liver and gill, respectively. The PEL vs. SFM comparison identified 1,315 and 2,089 DEGs for the liver and gill, respectively. Finally, the IM vs. SFM comparison identified 299 and 140 DEGs for the liver and gill, respectively. The comparisons with the highest DEGs numbers were MEL vs. SFM (2,245 DEGs) and PEL vs. SFM (2,089 DEGs; Fig. 3). A total of 48 and 37 contigs were identified as CYP genes in liver and gill, respectively. In the case of IM site, three and one contigs identified as CYP genes were differentially expressed in liver and gill, respectively. In the case of MEL site, ten and five contigs identified as CYP genes were differentially expressed in liver and gill, respectively. In the case of PEL site, four and five contigs identified as CYP genes were differentially expressed in liver and gill, respectively. Most contigs identified as differentially expressed CYP genes showed decreased expression in the three polluted sites compared to the SFM non-polluted site, except in the gills at the PEL polluted site where three of five CYP genes were upregulated (Fig. 4). The IM site had the fewest differentially expressed CYP genes, three in liver and one in gill, which could be related to it being a historical reference site, while the MEL site had the most with ten in liver and five in gill (Table 3).

Among the DEGs, seven CYP genes met the criteria of being present in both assembly sets and having a known function. Two (CYP family 4 subfamily B member 1 (CYP4B1) and CYP1A) were differentially expressed in all three polluted sites. Three were classified as endogenous, and four as exogenous based on their function (Table 4). In this group of seven CYP genes, those classified as endogenous showed differential expression in the liver but not in the gills, except for CYP family 26 subfamily B member 1 (CYP26B1). In contrast, those classified as exogenous CYP genes showed differential expression in the gills but not in the liver, except for CYP family 2 subfamily F member 2 (CYP2F2). The seven differentially expressed CYP genes, endogenous or exogenous, showed downregulation in all sites and both organs, except for CYP family 27 subfamily C member 1 (CYP27C1; Fig. 5). Additionally, the GLM confirmed the general downregulation of CYP4B1 in the liver (Fig. S1) and CYP1A in the gill (Fig. S2) at all polluted sites. The trend was unclear for the other CYP genes, except CYP2F2 in the liver, which showed a similar pattern to CYP4B1 but was not detected as differentially expressed in IM site.