Integration of Genotoxic Biomarkers in Environmental Biomonitoring Analysis Using a Multi-Biomarker Approach in Three-Spined Stickleback (Gasterosteus aculeatus Linnaeus, 1758)

Water is impacted by a variety of increasing pressures, such as contaminants, including genotoxic pollutants. The proposed multi-biomarker approach at a sub-individual level gives a complementary indicator to the chemical and ecological parameters of the Water Framework Directive (WFD, 2000/60/EC). By integrating biomarkers of genotoxicity and erythrocyte necrosis in the sentinel fish species the three-spined stickleback (Gasterosteus aculeatus) through active biomonitoring of six stations of the Artois-Picardie watershed, north France, our work aimed to improve the already existing biomarker approach. Even if fish in all stations had high levels of DNA strand breaks, the multivariate analysis (PCA), followed by hierarchical agglomerative clustering (HAC), improved discrimination among stations by detecting an increase of nuclear DNA content variation (Etaing, St Rémy du Nord, Artres and Biache-St-Vaast) and erythrocyte necrosis (Etaing, St Rémy du Nord). The present work highlighted that the integration of these biomarkers of genotoxicity in a multi-biomarker approach is appropriate to expand physiological parameters which allow the targeting of new potential effects of contaminants.


Introduction
In order to protect, restore and use sustainable water resources of continental hydrosystems, the European Commission set up the Water Framework Directive (WFD, 2000/60/EC) which aims at achieving good water quality in 2027 [1] by promoting the chemical and ecological status of freshwaters. Nevertheless, the chemical approach ignores toxicity, as well as mixture effects of contaminants [2], and ecological assessment is often considered as an a posteriori approach lacking preliminary information on toxicity [3], despite the inclusion of biomarkers in the WFD revision being promoted by scientists [2,[4][5][6]. What is more, some biomarkers are already used in the descriptor 8 (D8C2) of the Marine Strategy Framework Directive (MSFD, 2008/56/EC) to assess the biological effects of contaminants [7,8].
Monitoring biomarkers at the sub-individual level is a way of gaining more knowledge about early-warning signals of potential long-term alterations to population, community and ecosystem health. Therefore, multi-biomarker approaches are increasingly recommended to reflect the main physiological states of an organism [9][10][11]. Notably, this ( Figure 1). Sampled stations were recommended by the Water Agency for the inconsistency between their chemical and ecological status based on the WFD parameters.

Materials and Methods
Each experiment was conducted in accordance with the European directive 2010/63/ UE on the protection of animals used for scientific purposes at INERIS facilities (registration number E60-769-02). For the present work, all fish came from a well-characterized population of the French National Institute for Industrial Environment and Risks (INERIS) (Verneuil-en-Halatte, Oise, France) which were maintained in an outdoor pool with natural vegetation and macro-invertebrate communities.

Experimental Design and Sampling Sites
In Autumn 2019, 30 sticklebacks, with a balanced sex ratio and homogeneous length (47.04 ± 3.76 mm) and weight (1.47 ± 0.4 g) were caged in a 36 L cylindrical structure in six distinct stations in the north of France ( Figure 1, Table 1). The Cligneux river at the Saint-Rémy du Nord station (Nord, France) and the Rhonelle river at the Artres site (Nord, France) were considered as polluted small and natural streams [35,43], poorly urbanised and characterised by a strong pressure of agricultural practices. The Scarpe at the Biache-Saint-Vaast station (Pas-de-Calais, France) was considered highly polluted [35,43] and was located on an artificial canal. The Deule at Courrières (Pas-de-Calais, France) was a highly modified river used as a waterway. The Sensée river was studied on two points near Bouchain, (Sensée Bouchain, Nord, France), where the river was natural, and at the Etaing station (Sensée Etaing, Pas-de-Calais, France), a small stream surrounded by agricultural practices. All these sites were contaminated by polycyclic aromatic hydrocarbons (PAHs) [44]. Temperature, dissolved oxygen, pH and conductivity were monitored at the first day and after 21 days of caging for each station (Table 1). After 21 days of caging, as recommended [21] and already tested by Catteau et al. [11,15,45], sticklebacks were directly anesthetised in situ by balneation with MS222 (tricaine methanesulfonate, 100 mg/L, Tricaine Pharmaq, Overhalla, Norway) to avoid handling stress and then killed by cervical dislocation. Each fish was weighed, measured and sexed and each organ was sampled for further analysis. Briefly, 2 µL of blood was collected and directly diluted in 120 µL of citrate buffer [39]. Then, 50 µL of blood was maintained at 4 • C ± 1 • C until erythrocyte mortality analysis; the rest was stored at -80 • C for genotoxicity analysis. Liver and muscle tissues were collected, weighed and stored at −80 • C for biochemical analysis (approximately 10 mg for muscle). Spleens were removed, pressed through 40 µm sterilized nylon mesh with 1 mL Leibovitz 15 medium modified with penicillin and streptomycin (both at 500 mg/L) and maintained at 4 • C for 12 h before immune analyses to eliminate any bias due to stressful conditions of fish sacrifice [46].

Peripherical Erythrocyte Counts, Mortality and Genotoxic Biomarkers
The 50 µL of blood maintained at 4 • C was diluted with 150 µL of citrate buffer before performance of the flow cytometry analysis (MacsquantX, Miltenyi Biotec, Bergisch Gladbach, Germany). The peripherical erythrocyte counts were detected using forward scatter (FSC-size of cells) and side scatter (SSC-complexity of cells) parameters. Erythrocyte mortality was detected using a double labelling method with YO-PRO ® -1 (1 mM in DMSO, Invitrogen TM Thermo Fisher Scientific, Bothell, WA, USA) for apoptotic cells (B1, green fluorescence) and propidium iodide (1 mg/mL in water, Invitrogen TM Thermo Fisher Scientific, Bothell WA, USA) for necrotic cells (B3, red fluorescence) as preconized for leucocytes by Bado Nilles et al. [46]. After 10 min of incubation on ice and in darkness, samples were analysed after cell excitation by a 488 nm argon laser.
Chromosomal damage was determined for erythrocytes according to the method developed by Vindeløv & Christensen [39] and adapted for turbot by Goanvec et al. [50] and Marchand et al for stickleback [41]. Briefly, erythrocytes were adjusted, under 4 × 10 7 cells/mL, before successive treatments to disrupt plasma membranes to allow access to DNA, disrupt RNA and stabilise amino acids and mark DNA. A sample of stabilized chicken red blood cells (Fitzgerald) was used as a standard and analysed at the same time as the fish blood samples. Each FL3 coefficient of variation (CV) corresponds to a measure of nucleotide size. DNA damage corresponds to the CV of fish erythrocyte samples minus the CV of chicken red blood cell samples. The primary DNA strand break level was assessed by the comet assay on cryopreserved erythrocytes according to the protocol defined by Singh et al. [30] and previously adapted for three-spined sticklebacks by Santos et al. [35,51] with some modifications. All analyses were realized at 4 • C with inactinic light (dim red light) to avoid DNA damage. The blood samples were rapidly defrosted (18 s) in a 37 • C-bath, one-fiftieth (1/50) diluted with PBS buffer to adjust cell density and finally half-diluted with agarose LMP (low melting point, type VII, 1% w/v in PBS, 0.01 M, pH 7.4, without Ca or Mg) and maintained in liquid state at 37 • C. After chilling (15 min on ice), the slides were treated in several baths to disrupt plasma membranes to access nucleoids (lysing solution, Ph 10, 1 h) to obtain single-stranded DNA (denaturing buffer, Ph > 13, 40 min) and DNA strand breaks (SSB, DSB and ALS) (electrophoresis, 24 min, 20 V, 500 mA, 6 W). Then, the slides were washed (neutralization buffer, twice 20 min), dried (absolute ethanol, 10 min) and kept at room temperature and in the dark until their staining with Sybr Green 1 X, according to the supplier's recommendations. Thereafter, slides were scanned with a fluorescence microscope (EVOS TM FL Auto 2, Invitrogen TM Thermo Fisher Scientific, Bothell WA, USA, optical filter for Sybr Green (497 nm/520 nm)). For each slide, 200 comets by fish (corresponding to 100 nucleoids/gel with two gels for one slide for each fish) were scored by an image analysis system (Comet Assay IV version 4.3.1, Instem, Stafforshire, UK). To express DNA strand breaks, the percentage of tail intensity (TI%), corresponding to the percentage of damaged DNA, was chosen, according to Langie et al. [52]. For each gel (100 comets) the median of tail intensity was calculated. Comet assay results for the fish were expressed as the means of the two gel medians.

Statistical Analysis
All statistical analyses were performed with R software version 3.6.1. Before any other statistical tests, normality was checked using Shapiro's tests (p ≤ 0.05) and homogeneity of variance was tested using Levene's test (p ≤ 0.05). When the normality assumption was not fulfilled, data were log-transformed if this improved normality. Then, to assess the difference between stations in function of sex, each biomarker was analysed separately via a two-way ANOVA (package car) followed by a Tukey test (package lsmeans) for parametric data or two Kruskal-Wallis tests (package stats) on site and sex factors separately, followed by a Nemenyi test (package agricolae) for non-parametric data. Finally, hierarchical agglomerative clustering (HAC), including a standardised principal component analysis (PCA) (package FactoMineR), was performed. The purpose of these tests was to explore the relationship between biomarkers and the stations' characteristics and to illustrate the discrimination of stations (package FactoMineR). Information about the biomarkers contributing to the build of each cluster and the percentage of fish from each site gathered in the same cluster were collected. All hypotheses were tested at the level of p ≤ 0.05. The statistical analyses are detailed in the Tables S2-S4.

Results
The number of fish collected after 21 days of exposure and the responses of all biomarkers for each group of sticklebacks in each station are shown in Table 2 (Mean ± SD). These raw values were used for all the analyses presented below.
Twenty-nine fish were recovered in Bouchain and Etaing, 28 in Rémy du Nord, Artres, Biache-Saint-Vaast and only 11 in Courrières. The sex ratio was near 50/50 for each site. For all the fish studied, out of 153 sticklebacks collected, 51.6% were males and 48.4% were females. A slight mortality was observed, with one or two fish per cage for all sites, except for fish caged in Courrières, where a huge individual mortality (63%) was recorded.
Globally, based on measured biomarkers, two types of stations would appear. Courrières and Bouchain seemed to be characterized by low responses of innate immune capacities but with a high respiratory burst index and high cholinesterase activity. The four other stations (St Rémy du Nord, Artres, Biache-Saint-Vaast and Etaing) would be affected by high immune capacities, high antioxidant system responses and genotoxicity. With the aim of improving understanding of the results, a data analysis through a principal component analysis (PCA) followed by hierarchical agglomerative clustering (HAC) was proposed. Table 2. Biomarker responses assessed in sticklebacks after 21 day of caging in the six sites studied. Data are expressed as means ± standard deviation and letters represent differences between sites (p ≤ 0.05) according to a Tukey test or Nemenyi test. n: number of fish dissected for each site; Che: Cholinesterase; EROD: 7-ethoxyresorufin-Odeethylase; GST: glutathione-S-transferase; GSH: total glutathione; GPx: glutathione peroxidase; SOD: superoxide dismutase; CAT: catalase; TBARS: lipid peroxidation; DNA: desoxyribonucleic acid.

Effect of Current Biomarkers on Site Discrimination
The repartition of stations as functions of biometric, biochemical and immune responses was explained by 33.2% of the overall variance in the first two principal components of the PCA (Figure 2, Table S1). The first axis represented 18.2% of variance because of the respiratory burst index (14.80% of contribution), HSI (10.54% of contribution) and leucocyte mortality (10.29% of contribution). The second component explained 15.0% of variance, mainly GPx (20.51% of contribution), GSI (18.54% of contribution), CAT (12.72% of contribution) and granulocyte-macrophage percentages (11.28% of contribution) (Figure 2, Table S1). The fish caging in Courrières and Bouchain formed globally one group in the PCA, which was mostly divided by the first dissimilarity node separating Cluster 1 from the others around 1.6 on the HAC, constituted of 100% of fish from Courrières and 58.6% from Bouchain (Table 3, Figure 2). As shown with the PCA and the HAC (Figure 2, Table 4), the specificity of the fish caged in these two stations concerned their high respiratory burst index (6.95 ± 2.38 U for Courrières, 3.87 ± 1.94 U for Bouchain, Table 2; v-test of 9.9 for Cluster 1, Table 4) and ChE activity (65.34 ± 24.12 U/g protein for Courrières, 43.03 ± 21.66 U/g protein for Bouchain, Table 2; v-test of 5.8 for Cluster 1, Table 4), their weak granulocyte-macrophage percentages (40.36 ± 7.42% for Courrières, 39.52 ± 9.39% for Bouchain, Table 2; v-test of −5.6 for Cluster 1, Table 4) and their lower phagocytosis efficiency (10.95 ± 2.22% for Courrières, for 13.01 ± 4.09% Bouchain, Table 2; v-test of −5.3 for Cluster 1, Table 4). Nevertheless, Bouchain was also under-represented by the other clusters (13.8%, 20.7% and 6.9% of fish from Clusters 2, 3 and 4 respectively, Table 3).

Evolution of Site Discrimination by Adding Genotoxic Biomarkers
The inclusion of genotoxic biomarkers into the analysis improved site discrimination mainly along the first axis of the PCA (Tables 3 and 4, Figure 3) and made a good contribution for chromosomal damage (13.06% of contribution) on the first axis and for erythrocytes necrosis (9.12% of contribution) on the second axis, inducing a little less of the overall variance of the PCA (30.8%). The first axis, which represented 16.8% of variance, was still mainly characterised by respiratory burst index (18.90% of contribution) and by some other immunomarkers that were also prevalent, such as granulocyte-macrophage percentages (13.90% of contribution) and lysosomal presence (11.00% of contribution). The second axis, which explained 14.0% of variance, was still built by GPx (25.67% of contribution), GSI (16.47% of contribution) and CAT (13.91% of contribution) and now by HSI (12.72% of contribution) (Figure 3, Table S1). As shown in the first analysis, caged fish in Courrières and Bouchain formed globally one group better separated in this second PCA by higher correlation with the first axis ( Figure 3, Table S1) and by the first dissimilarity node separating Cluster 1 from the others around 1.3 in the HAC (Figure 3). Always, 100% of the fish from Courrières and now 62.1% from Bouchain (Table 3) fell completely under Cluster 1, as demonstrated by the HAC. This group was still characterised by a higher respiratory burst index (v-test of 9.8, Table 4), higher neurotoxicity (v-test of 5.9, Table 4), weak phagocytosis efficiency (v-test of −5.1, Table 4) and low granulocyte-macrophage percentages (v-test of −5.5, Table 4) and was modified for the most part by the better representation of lysosomal presence (334.09 ± 52.36% for males, 378 ± 66.64% for females for Courrières, and 409.31 ± 90.56% for males, 462.36 ± 78.41% for females for Bouchain Table 2, v-test of 5.4, Table 4) and by the addition of negative CVs of DNA content (-0.60 ± 2.18 for Courrières and −2.34 ± 1.35 for Bouchain, Table 2, with v-test of −5.9, Table 4). Bouchain was still also reflected by three other clusters (24.1%, 3.5% and 10.3% for Clusters 2, 3 and 4 respectively, Table 3).
Another group was still constituted by this second PCA (Figure 3) corresponding to the fish caged in the same sites (St Rémy du Nord, Artres, Biache-St Vaast and Etaing) already highlighted by the first PCA ( Figure 2). A HAC, more detailed due to the formation of a fifth cluster, was obtained with the inclusion of genotoxic biomarkers ( Figure 3, Table 4). Most of the fish caged in Etaing (51.7% of fish), some of those at St Rémy du Nord (25% of fish) and only 3.6% of the fish of Artres and Biache-St-Vaast were gathered in Cluster 5 (Table 3). Cluster 5 was mainly characterised by higher erythrocyte necrosis only at Etaing (11.99 ± 7.09% for males and 6.11 ± 3.58% for females) and St Rémy du Nord (9.21 ± 7.50% for males and 1.56 ± 0.67% for females,  Table 3) and is represented especially by higher leucocyte mortality (10.46 ± 3.77% for leucocyte necrosis and 8.70 ± 3.75% for leucocyte apoptosis for St Rémy du Nord, 12.33 ± 3.54% for leucocyte necrosis and 7.14 ± 3.30% for leucocyte apoptosis for Artres, 10.29 ± 4.24% and 5.27 ± 2.53% for Biache-St-Vaast, Table 2; v-test of 5.7, Table 4). HSI, GSI and SOD activities are well represented in this clustering, but statistical analysis showed non-significant differences between stations for these three parameters (p > 0.05) ( Table 2).
Furthermore, even if the biomarker of DNA strand breaks measured by the comet assay was not well projected in the PCA obtained, this parameter made a lower but significant contribution to Clusters 2 and 3 (v-test of -3.6 and 4.6 respectively,  Table 3) and Bouchain (24.1% of fish, Table 3) were located especially in Cluster 2, characterised by the lowest DNA strand breaks levels (22.36 ± 12.70% TI for Etaing and 18.94 ± 7.85% TI for Bouchain, Table 4).

Discussion
Many xenobiotics can impact DNA integrity, but until today few studies have assessed jointly DNA strand breaks, chromosomal damage and erythrocyte mortality in three-spined sticklebacks. The innovative aspect of this work was to assess genotoxicity at different scales of the genome in peripherical erythrocytes of this sentinel species in the context of active biomonitoring. The first objective was to gain more knowledge of biomarker responses of genotoxicity and erythrocyte mortality. The second was to determine whether these biomarkers could improve the already existing multi-biomarker approach in sticklebacks and their potential contribution to discriminate the studied stations.
As observed by Catteau et al. [11], regardless of the station, a high survival rate of fish was noticed (>93.3%), highlighting a good general health status, except for fish caged in Courrières, which were characterised by high mortality (63%). No significant weight loss was recorded at all stations (Table S4). These observations allow the interpretation of the physiological responses measured by the multi-biomarker approach after 21 days of caging.
For erythrocyte necrosis, fish from Etaing were the most affected. Etaing was also characterised by the highest CV for chromosomal damages (3.87 ± 4.72) and had one of the highest lipid peroxidation levels (TBARS, biomarker of cell integrity), higher than reference values established by Catteau et al. [65]. Fish erythrocyte necrosis was most affected at St Rémy du Nord, Courrières and Etaing, with higher responses in males compared to females. Even if biotic or abiotic factors, including gender, influence the responses of some biomarkers, such as AChE, catalase, EROD, GST and immunomarkers [65][66][67][68][69], no study in the literature seems to have highlighted the relationship between erythrocyte necrosis and gender in fish. This latter relationship was not observed for nuclear DNA content or DNA strand breaks in this work. However, considering that gender is essential to analyse erythrocyte necrosis inside the multi-biomarker approach, future investigations could lead to a better understanding of the links between erythrocyte necrosis and gender. In addition, Courrières and Etaing had the two lowest densities of erythrocytes, which could be correlated with high erythrocyte necrosis. The lifespan of circulating erythrocytes varies in fish between 100 to 500 days [70,71], with a maturation period of 17-23 days [72]. Erythropoiesis is slower than erythrocyte breakdown, which could be impacted by many factors, including environmental toxicants [73]. So, after 21 days of caging, the decline of erythrocytes could be observed before the renewal of blood cells and could be used as an attractive biomarker.
Concerning chromosomal damage, St Rémy du Nord, Artres, Biache-St-Vaast and Etaing were characterized by CVs greater than the CVs obtained in reference conditions (−0.96 ± 0.62) (current work, personal communication). The highest CV obtained by FCM was for Etaing (3.87 ± 4.72), which could be explained by the increase of the variation in DNA content [74], which resulted from chromosomal aberrations (clastogenic effect) [40] or abnormal cell division (aneugenicity) [75,76]. The CVs obtained for Courrières (−0.60 ± 2.18) were the nearest to reference conditions. This station seemed to be less affected by chromosomal damage, but these results came only from the most resistant fish that had survived. One reason for the massive fish mortality recorded in Courrières could be related to heavy genomic abnormality with aneuploidy patterns affecting the cytogenetic quality of cells shown in another species in a different context, the blue mussel (Mytilus spp.), which hardly withstood the environmental conditions [77]. This high fish mortality can be linked to the high erythrocyte mortality and low erythrocyte density (1.56 ± 0.67 cell/mm 3 × 10 5 ) measured in Courrières. Bouchain was characterized by a CV (−2.34 ± 1.35) under the CV obtained in reference conditions. So, for Bouchain, also, the nuclear DNA content variation is high and could be caused by chromosomal loss or nuclear condensation and shrinkage of fish erythrocytes, marking the early process of cellular mortality [78,79].
In this study, all stations had high levels of DNA strand breaks (between 18.94 ± 7.85% TI and 38.78 ± 19.69% TI) which were two to four times above our reference values. The latter is between reference values measured by Santos et al. and Le Guernic et al. using the same cellular model in sticklebacks [18,35]. Fish caged in Artres were the most affected by DNA strand breaks (38.78 ± 19.69% TI) followed by St Rémy du Nord, Biache-St-Vaast, Courrières, Etaing and Bouchain. Chromosomal aberrations are considered irreversible damages, unlike DNA strand breaks [24,25], but erythrocytes' DNA-repair capacities in fish species, and specifically in the three-spined stickleback, are not well known. As no decrease in DNA strand breaks was observed after 20 days of depuration following an exposure of 12 days to a well-known genotoxicant, MMS (methylmethanesulfonate, with 0.5 µM and 5 µM tested), Santos et al. supposed that DNA repair activity in stickleback erythrocytes was low [51]. DNA damage would correspond to the accumulation of primary DNA damage until the renewal of erythrocytes. So, although this biomarker does not differentiate stations in this work, if the erythrocyte lifespan is taken into account, the measurement of primary DNA damage in stickleback erythrocytes could be considered an integrative biomarker of genotoxicity [51]. So, the high DNA strand breaks measured in this work could be translated into an environmental degradation that could be linked to a genotoxic effect.
To investigate in more detail the health of sticklebacks caged in each station at the sub-individual level, biomarker responses were integrated into a multivariate analysis (PCA) followed by hierarchical agglomerative clustering (HAC) with the inclusion or not of erythrocyte necrosis and genotoxicity biomarkers. Gathering biomarker responses into multivariate statistics is advised to take into account all biomarker data and improve their interpretation [80]. The inclusion of genotoxic biomarkers inside the classic multibiomarker approach maintained and affirmed the global separation of the stations. So, both multivariate analyses mainly discriminated two groups of stations, one formed by Courrières and Bouchain and the second by St Rémy du Nord, Artres, Biache-St-Vaast and Etaing. The analysis demonstrated in the first group that fish from Bouchain had a similar pattern of biomarkers compared to Courrières, being represented by 62.1% and 100% of fish classified in Cluster 1, despite the high fish mortality observed in Courrières. Cluster 1, characterised by genotoxicity in terms of negative CVs for DNA content, has allowed us to separate these two stations from the others. In opposition, the new Cluster 5 grouped together half of the fish from Etaing, a quarter from St Rémy du Nord and some from Artres and Biache-St-Vaast, which was characterised by a positive CV of the DNA content and by higher erythrocyte necrosis. Biomarkers of genotoxicity and erythrocyte mortality refined the physiological description of the analysis by the formation of a fifth cluster, only built by nuclear DNA content variation and erythrocyte necrosis, improving discrimination among stations. The high DNA strand break level detected in all stations explained the weak discriminating capacity of this biomarker in this study.
Globally, along the Artois-Picardie watershed, fish caged in Courrières and Bouchain were mainly defined by innate immune responses (high leucocyte respiration, low granulocyte percentages and low phagocytosis efficiency), high cholinesterase activity and negative CVs of nuclear DNA content; St-Rémy du Nord, Artres and Biache-St-Vaast by biomarkers of the antioxidant system (GPx and catalase activities), leucocyte mortality and genotoxicity (increase of CVs, DNA strand breaks) and erythrocyte necrosis; and fish from Etaing were characterized by the same biomarkers, apart from leucocyte mortality, and the lowest cholinesterase activity. These biomarkers measured at sub-individual levels translated into a degradation of the general health of the sticklebacks in each station, which can be linked to different environmental stresses along the Artois-Picardie watershed. Land use and the pollution of rivers could explain this degradation (presented in detail in Table S5). For instance, all stations were chemically downgraded for a contamination of benzo(a)pyrene and fluoranthene belonging to the polycyclic aromatic hydrocarbon (PAH) family. In addition, four phytosanitary substances (aclonifene, cypermethrine, chlorpyrifos, bifenox) were measured, which are classified as priorities by the WFD. Notably, some of these molecules are potentially genotoxic substances, according to the literature [81][82][83][84][85][86][87][88]. Courrières, with high fish mortality, seemed to be the most degraded. This station is subject to all the contamination coming from upstream because Courrières is located at the confluence of two canals (the Deule canal and the Lens canal). High concentrations of organic, trophic and toxic pollutants, which are high and variable in this station, could explain this mortality. In the case of stations showing inconsistency between their chemical status and their ecological parameters under the WFD, as with the six studied stations of the Artois-Picardie watershed, this work has highlighted the relevance of the multi-biomarker approach with the integration of genotoxicity as a complementary indicator to environmental managers in the analysis of the global quality of continental waters.

Conclusions
This work has expanded the existing battery of biomarkers measured in three-spined sticklebacks during an active biomonitoring of genotoxicity, from DNA strand breaks and chromosomal aberration to erythrocyte mortality. The integration of nuclear DNA content variation and erythrocyte necrosis confirmed and improved the discrimination of the different stations, whereas high DNA strand break levels were detected by the comet assay in all stations. These biomarkers measured at sub-individual levels translated into a degradation of the general health of the fish in each station which can be linked to various environmental stresses along the Artois-Picardie watershed. The most degraded station was Courrières, based on the high fish mortality. The multi-biomarker approach with the integration of genotoxicity has allowed us to consider a set of biological effects highlighting the relevance of this approach as a complementary indicator of chemical and ecological parameters of the WFD. These first conclusions were based on studies of six stations along the Artois Picardie watershed. Further research will be carried out on the same stations that will invalidate or confirm these first results. Work is in progress on the three-spined stickleback to define basal values of these biomarkers of genotoxicity and erythrocyte mortality to improve the interpretation of these physiological responses and their reliability in aquatic biomonitoring.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/toxics10030101/s1: Table S1: Contribution (%) for each biomarker to the construction of the two main components of the PCA taking into account genotoxic and/or biometric, biochemical and innate immune biomarkers (Figures 2 and 3). Variables contributing more than 10% to the building of the axis are indicated in bold. HSI: hepatosomatic index; GSI: gonadosomatic index; Che: cholinesterase; EROD: 7-ethoxyresorufin-Odeethylase; GST: glutathione-S-transferase; GSH: total glutathione; GPx: glutathione peroxidase; SOD: superoxide dismutase; CAT: catalase; TBARS: lipid peroxidation; DNA: deoxyribonucleic acid; Table S2: Results of statistical analysis carried out to assess the differences between stations in terms of sex. Each biomarker was analysed separately via a two-way ANOVA followed by a Tukey test for parametric data or two Kruskal-Wallis tests on site and sex factors separately, followed by a Nemenyi test for non-parametric data (*** = p < 0.001, ** = p < 0.01, * = p < 0.05); Table S3: Biometric index (Fulton's condition, weight, standard condition) expressed as means ± standard deviation for each station. Stars represent a significant difference obtained between the biometric index at T0 (the first day of the caging) and at T21 (after 21 days of caging in the field) according to a Student's t-test or a Wilcoxon-Mann-Whitney test (* = p < 0.05); Table S4: Results of a statistical analysis carried out to assess the difference between each biometric index between T0 and T21 according to a Student's t-test or a Wilcoxon-Mann-Whitney test (* = p < 0.05); Table S5  Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found here: [https://www.eau-artois-picardie.fr/donnees-sur-leau/visualiser-et-telecharger-lesdonnees], accessed on 10 January 2022.
Acknowledgments: This work benefitted from the French GDR "Aquatic Ecotoxicology" framework which aims at fostering and stimulating scientific discussions and collaborations for more integrative approaches. The authors are grateful to ©BioRender for the creation of some of the pictures in the graphics. The authors thank Sandrine Joachim for English reviewing of this paper.