Comparative genetic analysis of natural and farmed populations of pike-perch (Sander lucioperca)

Pikeperch (Sander lucioperca) is a native fish species in the European and Asian basin, representing high natural and economical value. However, our information on the genetic background and diversity of European populations is still limited, despite of that the production and the number of bread stocks has increased significantly in the last decade. Our aim was to develop new useable species-specific microsatellites and compare genetic diversity of ten pikeperch population from the Danube drainage basin. Thirty-four novel species-specific DNA markers were isolated and seven polymorphic microsatellite markers were selected for the population genetic analysis. The results indicated strong anthropogenic effects among the populations. The FST estimated by the method of Weir and Cockerham (1984) was 0.214, showing moderate genetic difference among the populations. The Structure analysis and the neighbour-joining dendrogram are showing the same result: the ten examined populations aggregate six genetically distinct units. Significant lack of heterozygosity was detected in most of the populations. Presumably, it is indicating the effects of human activities and such as vigorous stocking of natural waters. In addition, the effect was proved that fish stocks involved in pond culture and RAS rearing do not necessarily originated from the geographically closest natural populations.


Introduction
Pikeperch (Sander lucioperca) is a valuable and traditionally important carnivorous fish species in European and Asian Aquaculture. The species is becoming more and more popular among anglers and consumers. However, environmental pollution, loss of habitat, and overfishing reduced the number of pikeperch in natural waters, thus consumer demands cannot be fulfilled (Hedrick and Miller 1992). To overcome the problem, several European research projects were initiated to developing production technologies and techniques (Abdolmalaki and Psuty 2007;Björklund et al. 2007; Bódis and Bercsényi 2009;Menezes et al. 2013;Zakeś et al. 2013;Zarski et al. 2013;Blecha et al. 2016;Lappalainen et al. 2016). Classical intensive breeding and rearing methods can only be used with strong limitations due to the high stress sensitivity of the species and its need for a relatively large territory (especially in the breeding season) which is probably due to the predatory nature of the species (Molnár et al. 2004). If (when) the breeding space is too small, cannibalism begins to take place (Policar et al. 2013;Szkudlarek and Zakęś 2007). In addition pikeperch is highly sensitive to dissolved O 2 levels and to mechanical injuries. There is significant progress in technology development of artificial propagation, incubation, pre-nursing, nursing, and formulated feed-based rearing of pikeperch (Lappalainen et al. 2003;Molnár et al. 2004;Blecha et al. 2016;Miroslav et al. 2016). Information is still scarce on the genetic structure and diversity of both natural water pikeperch populations and farmed stocks. However, it would be essential for the conservation of natural populations. Some of them are strongly influenced by anthropogenic affects, like uncontrolled stocking (Louati et al. 2016;Poulet et al. 2009;Salminen et al. 2012) and admixture (Eschbach et al. 2014) which are reduce the genetic differentiation among the natural populations.
On the other hand, genetic information and molecular genetic marker-based methods (e.g. marker-assisted selection, parentage analyses) are also important for efficient genetic improvement programmes, inbreeding, and selection programmes or for the maintenance of genetic variance of farmed stocks (Vandeputte et al. 2011;Kaczmarczyk and Fopp-Bayat 2013;Fopp-Bayat 2007). Until now, 27 polymorphic and speciesspecific microsatellite DNA markers were isolated from pikeperch (Kohlmann and Kersten 2008;Han et al. 2016) and some more microsatellites were described (Zardoya et al. 1996;Borer et al. 1999;Wirth et al. 1999;Leclerc et al. 2000;Li et al. 2007;Yang et al. 2009;Zhan et al. 2009;Grzybowski et al. 2010;Coykendall et al. 2014) from closely related species (walleye-Sander vitreus, yellow perch-Perca flavescens, perch-Perca fluviatilis). Although some of them (after adaptation) proved to be suitable for population genetic studies on the pikeperch too (Leclerc et al. 2000;Kusishchin et al. 2018), there is a significant need for new, well characterised, highly polymorphic genetic markers that are applicable for gaining new and more detailed information on the genetic variability and unique genetic characteristics of different populations and stocks of the species. In case of species where the genomic sequence is not available and therefore SNP markers are not applicable, microsatellite DNA markers are still widely used for population genetic studies. These codominant markers are usually located in the non-coding regions of the genome and therefore, the mutation rate is high and selection pressure against mutations is relatively low resulting in high variability (Sunnucks 2000). Due to their small size, easy amplification and fragment separation techniques and easy evaluation, microsatellites are ideal tools for conservation biology studies of less studied (fish) species (Chistiakov et al. 2006). With a sufficient, relatively large number of microsatellites, genetic mapping can be carried out (Botstein et al. 1980;Knapik et al. 1998) or breeding programmes can successfully be supported (Beuzen et al. 2000). Moreover, microsatellites can be applied in aquaculture during the selection of spawners based on genotype profiles (Kaczmarczyk and Fopp-Bayat 2013) or in genotyping of cryopreserved fish spermatozoa (Fopp-Bayat and Ciereszko 2012). In aquaculture genome engineering, microsatellites are generally used for identification of maternal or paternal inheritance (Fopp-Bayat 2007) or during the investigation of ploidy level of manipulated fish (Fopp-Bayat and Woznicki 2006).
The aim of the present study was to isolate and characterise novel microsatellite DNA markers from the pikeperch and to describe the genetic diversity, population genetic status, and the effects of breeding and stocking in a few natural populations and farmed stocks in the River Danube drainage system.

Sampling
Fin samples were collected from 376 adult pikeperch individuals and stored in 96% ethanol. To reduce animal suffering cloves oil (Syzygium aromaticum) was used as an anaesthetic bath before sampling the fish (15 oil drops/10 l water). After the fish have become immobile, they were sampled (the outermost 0.5-1 cm 2 of caudal fin were clipped) and transferred to clean water until they wake up. Samples were stored in ethanol on − 20°C until utilisation. Samples (Table 1) (Fig. 1). However, the number of the available samples in case of some population (e.g. De or Ti) are too low to represent the full populations, but they are providing important data.

DNA isolation
Standard phenol/chloroform (Rothi-Phenol, Carl Roth; Chloroform, Reanal) extraction method was used to isolate DNA (Blin and Stafford 1976). The quality of extracted DNA was tested by agarose gel electrophoresis and DNA concentration was measured by a spectrophotometer (IMPLEN, NanoPhotometer). Twenty-microgram DNA was separated for genomic library construction while the rest of the samples were diluted to 50 ng/μl concentration for microsatellite analysis.

Genomic library construction
The modified method of Glenn and Schable (2005) was used to construct a repeat-enriched genomic library. The DNA used for library construction was isolated from male individuals since in the examined species, the males are heterogametic and therefore, they carry both sex chromosomes (Rougeot et al. 2002). Twenty-microgram genomic DNA was digested by restriction endonucleases (Rsa I, Fermentas/HpyCH4 V, NewEngland BioLab) creating blunt ends. After separating the products on agarose gel, all fragments with the length of 300-1000 base pairs were re-isolated from the gel using a NucleoSpin Extract II kit (Macherey-Nagel). The DNA concentration was then measured by spectrophotometer (IMPLEN, NanoPhotometer), and following phosphatase treatment (Shrimp Alkaline Phosphatase, Fermentas) of 10-μg fragments, Box I linker (F: 5′-Phos-ATGTCTGAAGGTACCACTGC TGTCCGAAA-3′;R: 5′-CGGACAGCAGTGGTACCTTCAGACAT-3′) was ligated to them. Ligation was checked by a linker-specific PCR in a final volume of 25 μl: 1× Taq DNA polymerase buffer (Fermentas), 0.4 μM Box I reverse primer, 2 mM MgCl 2 (Fermentas), 0.2 mM dNTP (Fermentas), 1 U Taq DNA polymerase (Fermentas), and 4 μl template (adapter connected DNA-fragment). The PCR profile of the reaction was the following: 2 min 94°C, 35× 94°C 20 s, 60°C 30 s, 72°C 3 min, and finally, 72°C 5 min (Mastercycler 5341, Eppendorf). The PCR was followed by the collection of tandem repeats containing DNA fragments. The fragments were hybridised with a (CA) 10 carrying oligonucleotide (3′ biotinylated) in a thermal cycler device using the following protocol: the initial denaturation (92°C 5 min) was followed by the cooling of the reaction mixture from 70 to 50°C with the decreasing of the temperature by 0.2°C/s. The mixture was kept on 50°C for 10 min and then it was cooled to 15°C with the decreasing of the temperature by 0.1°C/s. Hybridization complexes were bound to the surface of streptavidin covered magnetic beads (Dynabeads M-270 Streptavidin, Invitrogen) according to the method Glenn and Schable (2005). The resulting complexes were removed from the solution by a magnet. The repeat-containing DNA fragments were finally eluted by TLE solution (pH = 8.00) on 95°C. The eluted singlestranded DNA was transformed to double-stranded DNA by the previously described PCR (Mastercycler 5341, Eppendorf) using adapter specific primers. Resulting double-stranded fragments were ligated to T-vector (pGEM-T Easy Vector System I, Promega) and transformed into competent Escherichia coli cells (XL10 GOLD, Stratagene) following the protocol described by the manufacturer. Colonies were screened by blue-white screening (Ullmann et al. 1967). Inserts were checked by agarose gel electrophoresis and colony PCR initiated at the T-vector coded M13 primer binding sites. The PCR was carried out in a final volume of 25 μl with the following ingredients: 1× Taq DNA Polymerase Buffer (Fermentas), 0.26 μM-0.26 μM M13 forward and reverse primers (F: 5′ TGTAAAACGACGGCCAGT 3′; R: 5′ CAGGAAACAGCTATGACC 3′), 2 mM MgCl 2 (Fermentas), 0.2 mM dNTP (Fermentas), 1 U Taq DNA Polymerase (Fermentas), and a few cells of bacteria colonies as template. The PCR profile of the reaction was the following: 3× 95°C 2 min, 55°C 1 min, 72°C 2 min, then 41× 95°C 30 s, 55°C 30 s, 72°C 45 s, and finally, 72°C 5 min. Fragments longer than 300 bp were cleaned up (PCR Advanced Clean Up System, Viogene) and the sequence of the inserts were determined (3130 Genetic Analyser, Applied Biosystems) using the SP6-(5′ CATACGATTTAGGTGACACTATAG 3′) and T7-(5′ TAATACGACTCACTATAGGG 3′) primers and 3.1 BigDye kit (Applied Biosystems). Sequences were aligned with MEGA5 software (Tamura et al. 2011) and the ones containing at least 5 dinucleotide repeats were selected and primers were designed to their flanking regions Primer3Plus software (Untergasser et al. 2007). Their optimal reaction parameters were determined.

Microsatellite analysis
Microsatellites were tested on a test panel of 8 individuals. After that, microsatellite analysis was carried out on the above mentioned populations by using PCR. At first, directly labelled (5′FAM) microsatellite specific forward primers were used (in case of MS 84 Sl, MS 146 Sl,MS 150 Sl,MS 192 Sl,MS 195 Sl,MS 198 Sl,MS 203 Sl,MS 260 Sl,and MS 268 Sl markers). Later, we chose a more simple method according to Shimizu et al. (2002) to reduce the costs of analysis and labour. The forward primers were elongated with a, 5′, 17 bp long tail sequence (tail; 5′ATTACCGCGGCTGCTGG-microsatellite specific sequence-3′) to provide attachment site for fluorescently labelled (dye: FAM or VIC or NED or PET) Buniversal^oligos (tail-specific oligo; 5′dye-ATTACCGCGGCTGCTGG-3′). One tail-specific oligo was added to the reaction mixtures for fluorescent labelling of amplicons. Table 2 summarises the ingredients of PCRs. The cycling conditions for the  Table 3) 1 min, 72°C 2 min, 35× or 45× (see Table 3) 95°C 15 s, annealing temperature (see Table 3) 20 s, 72°C 40 s, and finally, 72°C 5 min by using Mastercycler 5341 (Eppendorf) thermal cycler.

Results
Two CA-dinucleotide enriched pikeperch genomic library were constructed with the use of two restriction enzymes (Rsa I, HpyCH4 V). Colony PCR confirmed insert DNA incorporation in case of 115 clones (out of 208 screened clones). Among these, 109 unique sequences were found, and 101 contained microsatellite repeat regions showing that the enrichment was highly effective (93%). New, unique sequences were deposited in GenBank (Table 3). Thirty-four  Microsatellite DNA markers marked with * were applied for the population genetic study sequences were suitable for the development of fully functional microsatellite DNA markers. These markers were tested on a test panel of eight pikeperch individuals to determine their functionality and characteristics. All the 34 functional markers were proved to be polymorphic and the number of detected alleles ranged between 3 and 20 (Table 3). The highest number (20) of amplified alleles was found in case of the MS 260 Sl marker while MS 84 Sl,MS 192 Sl,MS 412 Sl,and MS 424 Sl were also found to be highly polymorphic. Seven of the newly developed microsatellite markers (MS 192 Sl,MS 195 Sl,MS 198 Sl,MS 203 Sl,MS 260 Sl,MS 268 Sl,MS 397 Sl) were applied to analyse ten populations (Ge, Kb, Gy, Ba, Ak, At, Da, Ny, Ti, and De) originating from the drainage system of the River Danube (Fig. 1). According to the MICRO-CHECKER software, there is no large allele drop outs and stuttering errors. Possibility for null alleles was detected for MS 192 Sl in Ny and Da populations,MS 195 Sl in Ak and Da,MS 198 Sl in De,MS 203 Sl in Da,MS 260 Sl in Ny,and MS 397 Sl in At and Kb populations.

MS 84 Sl
The highest number of alleles (17 and 20), the highest mean number of alleles per loci (7 and 6.7), and the highest allelic richness (10.2 and 8.6) values were detected in case of MS 192 Sl and MS 260 Sl markers respectively, while the total average allele number was also relatively high (9.14). The highest average number of alleles per locus was found surprisingly in the artificially established and bred Kisbajcs (Kb: 5.43) and Dalmand (Da: 4.71) stocks. The reason for this phenomenon can be explained with the establishment storey of these stocks. The farmers purchased broodfish from several different fish farms and natural waters at different times. If population size is also considered, the population from the Danube estuary (De) can be classified as a genetically rich population based on allelic richness values (Table 4).
The expected heterozygosity (H E ) values ranged between 0.4518 and 0.5931 while the observed heterozygosity (H O ) values were found between 0.4150 and 0.5665. The differences between H E and H O values were found to be non-significant in all but two case, meaning that only the Timișoara (Ti) and Győr (Gy) populations are in Hardy-Weinberg equilibrium. The highest average genetic diversity was found in the Danube estuary population (De) and in Attala (At) stock (Table 4). In case of the majority of examined populations and stocks, a low but significant lack of heterozygotes was detected while in case of Kisbajcs stock (Kb) the excess of heterozygotes was present (F IS = − 0.038). MNrA, mean number of alleles per population; Ar, allelic richness; H E , expected heterozygosity; H O , observed heterozygosity; HWE, deviation from Hardy-Weinberg Equilibrium. Level of significance is non-significant (ns) or significant *P < 0.05; **P < 0.01;***P < 0.001. AGD, Average gene diversity; F IS , an index describing the variance within the population (Ak, Akasztó; At, Attala; Gy, Győr; Ti, Timișoara; Kb, Kisbajcs; Ny, Nyíregyháza; Da, Dalmand; Ba, Lake Balaton; De, Danube estuary; Ge, Upper Danube) During the examination of the different populations and stocks, it was recognised that although the genetic variability was not very high between the groups, the different populations and stocks possess individual characteristics such as private alleles. The highest number of individual alleles (4) was found in the Danube estuary (De) population while the Kisbajcs (Kb) and Győr (Gy) stocks had 3 and 3 private alleles, respectively. The MS 260 Sl marker showed nine private alleles while we found three private alleles in case of MS 397 Sl. The frequency of these private alleles is low, except two private alleles (171 and 173 bp) in the German (Ge) population on MS 260 Sl locus with a frequency of 0.179 and 0.143, respectively.
The F ST value calculated for all populations and stocks was 0.214 showing a moderate genetic differentiation between these groups. The pairwise F ST values between pairs of populations were also calculated and in some cases, the genetic differentiation was found high. The highest differentiation was found between two natural water populations. The lowest genetic difference was found between the artificial Dalmand and the natural Balaton stocks, confirming that the former originate from the later one. Furthermore, low genetic differentiation was found between the artificially propagated stocks (Table 5). When summarising the results, one can state that the most of population and stock pairs does not represent high genetic differentiation (F ST < 0.25) which can be explained with the fact that broodstocks of intensive systems and fish farms were established with the use of individuals from different origins and fish were restocked to natural waters sometimes randomly from these stocks.
Nei's genetic distances (Da) (Nei et al. 1983) were specified and based on the results a dendrogram (neighbour Joining) was constructed in order to represent the relationship between populations (Fig. 2). The highest genetic distance (0.807) was detected between the German (Ge) and the Győr (Gy) populations where the highest genetic differentiation was also described. Furthermore, high genetic distances were observed between the Timișoara (Ti) and German (Ge; Da = 0.758), the Kisbajcs (Kb) and German (Ge, Da = 0.778), and between the Danube estuary (De) and German (Ge; Da = 0.713) population pairs, respectively. The smallest genetic distance value was found between the Dalmand (Da) and Lake Balaton (Ba; Da = 0.040) stocks. Similarly, low levels were found between the Akasztó-Attala (Ak-At; Da = 0.073) and Kisbajcs-Nyíregyháza (Kb-Ny; Da = 0.085) stock pairs.
The genetic structure of the examined populations was analysed assuming the presence of a different number of clusters based on the microsatellite results using STRUCTURE. The method  Evanno et al. (2005) implemented in STRUCTURE HARVESTER software indicated that the most probable number of clusters was K = 6 ( Fig. 3) showing that the samples collected from ten different locations can be classified into six genetically different groups.

Discussion
In the present study, natural water populations from the Danube catchment area and pikeperch stocks reared in extensive fish farms and in intensive systems in the basin of the river were examined. The newly developed species-specific microsatellite DNA markers could complement previously developed marker sets (some of them adapted from closely related species) and could provide a highly useful tool for gaining more information on the genetic background of the species. Two microsatellite-enriched genomic libraries were constructed using the method (slightly modified) of Glenn and Schable (2005). Ninety-three percent of the specified sequences contained microsatellites and 33.6% of these were successfully developed into genetic markers. Thirty-four microsatellites were tested altogether for functionality and polymorphism. The method of enrichment proved to be very efficient in contrast to traditional isolation methods (Rassmann et al. 1991;Wu et al. 1994;Lench et al. 1996;Hunt et al. 1999) Compared to other enrichment-based microsatellite isolation methods, the efficiency results are still acceptable (Ostrander et al. 1992;Armour et al. 1994;Kandpal et al. 1994).
The newly developed markers were proved to be polymorphic; the use of the ones with higher polymorphism levels and allelic richness is recommended in further studies even in case of closely related species.
Seven out of the new markers (MS 192 Sl,MS 195 Sl,MS 198 Sl,MS 203 Sl,MS 260 Sl,MS 268 Sl,and MS 397 Sl) were used to describe the genetic structure, variability, and relatedness of ten populations and stocks originating from the Danube catchment area. We looked for the effects of increased stocking and intensive pikeperch farming on the natural populations and found that in case of the majority of markers and in all examined populations and stocks, the expected and observed heterozygosity differed significantly from each other. Some analyses reported populations with higher genetic variance. For example, Kusishchin et al. (2018) described that the pikeperch populations in Akhtuba and the Volga rivers have higher heterozygosity (0.792 and 0.750) than in some examined Baltic (0.59-0.71) (Khurshut and Kohlmann 2009;Säisä et al. 2010;) and Western European populations (< 0.67) (Eschbach et al. 2014) populations. However, Bprecise comparison is hardly possible,^because the experiments are used different marker sets. In a similar comparison, the heterozygosity of our examined populations are similar to other European populations (Table 4.), together with a small but significant lack of heterozygotes (F IS = 0.039). This phenomenon is probably due to overfishing and the low number of spawners (broodfish) used for artificial propagations. Heterozygosity could be maintained or increased with the implementation of long term, well-planned breeding programmes (based on parentage analyses, genotypes of the breeders, crossing rare and/or distant genotypes of the stock, or introducing new genotypes).
The analysis of the difference among the population did not show outstandingly high genetic difference (F ST = 0.214); however, pairwise tests showed notable separation between the German (Ge) population and the Győr (Gy), Kisbajcs (Kb), Nyíregyháza (Ny), and Timișoara (Ti) stocks. The highest genetic differences were found between the population originating from the Upper Danube in Germany and all other populations. The highest number of private alleles was also described in the German (Ge) population (instead of the low number of analysed individuals). The segregation of this population is probably the consequence of geographic distance. The highest genetic diversity was found in Kisbajcs (Kb) intensively reared stock, which is most likely the result of recent broodfish import from different external resources or hybridization with pikeperch from unsteadied area. This is also supported by the dendrogram based on the Nei's genetic distances between population pairs, where the wild living populations can be found on the tree in relation to their geographic location (The Lake Balaton (Ba) population can be found between the Danube estuary (De) and Upper Danube (Ge) populations). The stocks reared in extensive fishponds and in intensive RAS cannot be fitted into the above mentioned system. The origin of these stocks can be clearly described and the (sometimes identical) genetic basis can be always identified. For example, the Kisbajcs (Kb) and Nyíregyháza (Ny) stocks (geographic distance3 10 km) are standing genetically closer to each other than the geographically close Kisbajc (Kb) and Győr (Gy) stocks (geographical distance~7.2 km). The latter is an intensively selected and reared, artificial feed consuming stock that shows the highest genetic similarity to the Attala (At) stock (geographical distance~145 km) and Akasztó (Ak) stock (geographical distance~305 km), because Győr (Gy) stock was derived from these two (At, Ak) pond hatcheries (pers. comm.), while the Dalmand (Da) stock is standing genetically close to the Lake Balaton population, since the farm was established with broodfish from Lake Balaton. The study on the genetic structure of the populations showed the same results, the ten stocks can be clustered into six genetically different groups. The Dalmand stock originates from Lake Balaton, while Győr stock was established by brood fish from Akasztó and Attala and the Kisbajcs stock was selected from Nyíregyháza stock. The linkage between the detected alleles and the sex of the examined fish was also tested but no correlation was found. These results proved the presence of anthropogenic effects and draw attention to the increasing requirement of genetic conservation of natural pikeperch populations.
The results suggest that there is a continuous and direct connection between fishpond/ intensive stocks and natural water populations. The broodstocks of intensive systems and fish farms are established by mixing individuals of different origin and fish are restocked to natural waters, sometimes randomly from these stocks.
The breeding programmes (Abdolmalaki and Psuty 2007;Björklund et al. 2007;Menezes et al. 2013;Lappalainen et al. 2016) and technological developments (Zakeś et al. 2013;Zarski et al. 2013;Blecha et al. 2016;Miroslav et al. 2016) have occurred in the last 10 years had fortunately only a minor effect on the natural populations of pikeperch in the Danube basin so far. But since there is a deviation from Hardy-Weinberg equilibrium due to overfishing and poorly planned restocking, the natural populations may face the loss of their original genetic variability. The intensification of breeding and selection programmes together with high levels of stocking can degrade the genetic basis of the species and can lead to the disappearance of varieties and to the uniformization of the different populations as it happened in case of the long time ago domesticated common carp (Hulak et al. 2010) and brown trout (Ward 2006). Recognising this phenomenon suggests that the conservation of pikeperch genetic resources became an actual issue.