Species composition of introduced and natural minnow populations of the Phoxinus cryptic complex in the westernmost part of the Po River Basin (north Italy)

Invasive alien species are a major driver of biodiversity loss, with their impacts potentially more intense when complexes of cryptic species are involved. In freshwaters, the anthropogenic manipulation of fish communities has resulted in altered fish communities, and in Europe has increased the complexity of Phoxinus species assemblages. Here, we investigated the Phoxinus communities of the westernmost part of the Po river basin, where adjacent freshwater ecosystems (Alpine high-altitude lakes and lowland streams) are representative of different management strategies (i.e. manipulated fish communities via stocking in Alpine lakes vs. natural populations in streams). We tested the genetic composition of the cryptic Phoxinus populations inhabiting these waters, as the species are morphologically indistinct. Sequences of the mitochondrial cytochrome oxidase I (COI) were obtained from 239 specimens, with the results indicating that 17 Alpine high-altitude lakes are now populated by a complex of Phoxinus species, comprising P. septimaniae (native to the Mediterranean area of France), P. csikii (native to the Central Balkans) and P. lumaireul (native to the North Adriatic Sea basins). Their introduction resulted from their use as angling live baits. Minnow populations in lowland streams were primarily comprised of native P. lumaireul, with only a single P. csikii specimen detected. While nuclear sequences of the recombination activating gene 1 (RAG1) marker were not useful for tracking the presence of alien alleles in these stream populations, the COI data emphasised the importance of using molecular tools to investigate cryptic species complexes that have been modified by anthropogenic activities.


Introduction
Invasive alien species are an important driver of biodiversity loss in aquatic ecosystems ) and can result in serious economic impacts (Gallardo et al. 2016). Freshwater fish communities can be highly impacted by alien species, with invasions facilitated by the high connectivity of river habitats that enables rapid dispersal (Hermoso et al. 2011). In the Mediterranean region, alien fish account for more than a quarter of the total number of fish species present (Leprieur et al. 2008). For example, in Italy there are 41 established alien fish species present, compared with 48 native species (Nocita et al. 2017); in some rivers, such as the downstream areas of the River Po, the native fish community has largely been substituted by alien fish species .
When a new alien fish introduction has occurred, its effective management requires it to be detected early, ideally before the species establishes, followed by a commensurate intervention, such as eradication for a high risk species (Vander Zanden et al. 2010;Britton et al. 2011). Early detection usually results in fewer resources being required for eradication as the alien species is still of low abundance and spatially constrained (Britton et al. 2011). However, this early detection can be hindered by the unresolved taxonomic status of species, or by the presence of cryptic species (Bickford et al. 2007). This increases the likelihood of a 'hidden' invasion, i.e. one that goes unnoticed until the invader has established and dispersed, and so has achieved a relatively wide distribution before being detected (Morais and Reichard 2018). Indeed, while cryptic alien species are widespread across different taxa (Pfenninger and Schwenk2007;Adams et al. 2014), many are assumed to remain undetected within communities of native species (Morais and Reichard 2018). Consequently, their detection is increasingly reliant on the application of molecular tools (e.g. Smith et al. 2012;Uchii et al. 2016).
An example of a cryptic species complex whose species composition and distributions remain subject to considerable knowledge gaps is that of the Eurasian minnows of the genus Phoxinus Rafinesque 1820. A group of cyprinoid fishes (Mayden and Chen 2010), the Phoxinus genus is characterized by a wide distribution range (basins of Atlantic, North and Baltic Seas, the Arctic and the northern Pacific Ocean), with its species able to populate a wide range of habitats, from large lowland rivers and lakes to the cold, welloxygenated waters of fast-flowing mountain streams (Kottelat and Freyhof 2007). As minnow populations only show subtle morphological differences and phenotypic plastic responses to different habitat typologies (Ramler et al. 2016), only one European species was initially described, P. phoxinus (Kottelat and Freyhof 2007). In the last decade, however, morphological and molecular studies have suggested that the Phoxinus genus comprises of at least 19 valid species (Eschmeyer et al. 2020). Within these species, the Italian minnow P. lumaireul Schinz 1840 has been re-validated as a species (Kottelat 2007) and it is endemic to basins of the north Adriatic Sea (P. lumaireul sensu stricto i.e. clade 1a; Palandačić et al. 2017).
The distribution of European minnows has been altered by introduction practices that have occurred via their use as angling live baits or through being contaminants of stocks of salmonids being used to enhance angling, especially in high-altitude or highlatitude lakes (e.g. Museth et al. 2007;Miró and Ventura 2015). Natural hybridization zones can further increase the complexity in defining the actual distribution ranges of these species (Palandačić et al. 2017;Corral-Lou et al. 2019), despite the relevance of this information for planning suitable strategies for conserving the native lineages and managing the distributions of alien ones. As a result of these processes, there is now minimal information available on the composition of the Phoxinus communities in many areas in southern Europe, including in the River Po basin. Correspondingly, in this study, we investigated how anthropogenic activities have impacted the cryptic species complex of the Phoxinus genus in the adjacent freshwater ecosystems of Alpine high-altitude lakes and lowland streams, focusing on the westernmost area of the River Po basin. In this area, the high-altitude lakes are regarded as being naturally fishless (Tiberti et al. 2014) and the lowland streams are recognised as the natural distribution range of P. lumaireul sensu stricto (i.e. locus typicus; Schinz 1840), with no known minnow stockings. To determine the species present within this cryptic species complex, we applied molecular (mtDNA and nDNA) tools, an approach that also enabled assessment of the genetic integrity of the native populations in the lowland streams, an aspect important for their conservation.

Sampling
The sampling of the Phoxinus communities covered the westernmost part of the Po River basin in the Piedmont region of Italy. Among 30 Alpine lakes visited in the sampling period, introduced minnows were detected in 17 lakes that were in the Western Alps at 1500-2500 m.a.s.l. (Figure 1; Online resources Table ESM 1). Of these 17 lakes, four were ponds with depths below 3 m (CFO, ISC, RES, ORG), with the others all being typical Alpine lakes with maximum depths between 3 and 5 m, and of sizes of approximately 0.5 (ISC, RES and ORG) to 158 ha (CER). Most were of natural, glacial and karstic origin, but three (CER, TER and CHI) were manmade. Most had a perennial outflow, although the outflow was intermittent for five lakes (VIS, ROB, TER, ANN, RES). The common feature of these lakes was that they are naturally fishless (Giussani 1997) and their fish community is artificial, reflecting the local practice of species introductions for angling. Indeed, fish releases for angling are a widespread practice in the Alps (Jersabek et al. 2001), as well as in other high-altitude lakes, such as in the Pyrenees in Europe (Miró and Ventura 2015). In the lakes, we also detected the presence of other fish species, mainly salmonids (Online resources Table ESM 2).
In the westernmost lowland streams which were all tributaries of the Po River, the fish assemblages included Phoxinus spp., with other species present being typical of those inhabiting the hyporhitral 'Grayling zone' and the epipotamal 'Barbel zone'. Alien fishes such as Gobio gobio, Salmo trutta and its hybrids with the native S. marmoratus, Oncorhynchus mykiss and Salvelinus fontinalis were present across the sampled sites (i.e. PEL, VAR and PES) ( Table ESM 2).
Between June 2016 and 2019, a total of 157 and 82 Phoxinus specimens were collected from the 17 Alpine high-altitude lakes and the four lowland streams respectively (Table ESM 1). This sampling was conducted using electric fishing in lowland streams, and by line fishing and hand nets in the lakes. Following their capture, all fish were anesthetized (MS-222) and a biopsy of the pelvic fin taken and preserved in ethanol (98%). The fish were then released back to their site of capture.  Table ESM 1 for further details on sampling sites where minnows were collected (i.e. altitude, geographic coordinates, number of fish analysed at both the mtDNA and nDNA)

Molecular analyses
Total genomic DNA was extracted from all individuals using the salting out method based on a proteinase K digestion, followed by sodium chloride extraction and ethanol precipitation (Aljanabi and Martinez 1997). A fragment (635 bp) of the mitochondrial cytochrome oxidase subunit I (COI) locus was amplified from all the specimens, whilst in the lowland populations a fragment (852 bp) of the nuclear recombination activating gene 1 (RAG1) was also amplified. Primer pairs used for each locus are available in Palandačić et al. (2017).
The PCR assays were performed using a Multiplex PCR kit (Qiagen) in 10 ll reaction volume that contained approximately 10 ng of template DNA and 0.25 lM of each primers pair. For the COI locus, thermal cycling was performed as follows: denaturation of 15 min at 95°C, followed by 35 cycles at 94°C for 30 s, 55°C for 30 s, and the extension step at 72°C for 45 s; the final elongation was at 72°C for 10 min. Thermal cycling conditions for RAG1 gene followed those of Palandačić et al. (2017). PCR products were purified using ExoSAP-IT TM (USB, Cleveland, USA) and directly sequenced by MACRO-GEN Inc (Amsterdam, The Netherlands; http://www. macrogen.org) using a 3730XL DNA Sequencer. All haplotypes generated in this study were deposited in the GenBank database (COI: MK984796-MK984822 and MT385939-MT385940; RAG1: MT410529-MT410545).

Phylogenetic data analyses
All the COI and RAG1 sequences obtained were automatically aligned, and then visually checked and adjusted. Identical COI sequences were collapsed into haplotypes in order to facilitate the computational processes and reduce redundancy. For the nuclear RAG1 dataset, the haplotype phases were solved using the PHASE algorithm available in DnaSP v 6 (Rozas et al. 2017). The phylogenetic analyses were exclusively performed on the COI mtDNA, as the RAG1 nDNA fragment had a low signal that could not cluster supported phylogenetic clades (Palandačić et al. 2017;Corral-Lou et al. 2019). The COI dataset was enlarged with 191 evolutionary related sequences (Online  resources Table ESM 3), with trees rooted on P. phoxinus (MF407769) sensu Palandačić et al. (2017).
Phylogenetic tree reconstructions were performed using neighbour joining (NJ), maximum likelihood (ML) and Bayesian inference (BI), through PAUP v 4.0 b10 (Swofford 2002), GARLI v 2.0 (Bazinet et al. 2014) and MrBayes v 3.1.2 (Ronquist et al. 2012) software respectively. The best evolutionary model, identified under the Akaike's information criterion (AIC), was GTR. Statistical support for the NJ and ML phylogenetic tree nodes was estimated as bootstrap probability (btp) values over 1000 replicates and then mapped onto the respective tree with PAUP software (Swofford 2002). Congruence in phylogenetic structure for the BI was performed using four independent Markov Montecarlo Coupled Chain (MCMC) runs of 10 6 generations each to estimate the posterior probability (pp) distribution. Topologies were sampled every 100 generations, and the majority-rule consensus tree was estimated after discarding the first 25% of generations. COI mtDNA uncorrected p-distances were also estimated.
On both the mitochondrial and nuclear datasets, Minimum Spanning Networks (MSNs) were built on haplotypes through a statistical parsimony criterion as available in PopART v1.7 software (Leigh and Bryant 2015). Mitochondrial MSN was built on the sequences produced in the study, while the nuclear data set has been enriched with 32 unphased sequences retrieved from GenBank that belonged to the minnow species detected at the COI locus (Online resource Table ESM  4). For each sampling site and marker, genetic variability indices were also estimated.

mtDNA phylogeny and MSN
A total of 52 haplotypes were identified in the 635 bp length multiple COI alignment that was obtained from the 239 analysed minnows plus the 191 sequences retrieved from GenBank (Online resources Table ESM  5). There were 70 variable nucleotide positions detected, 16 of which were singletons and 54 were parsimony informative sites. Among the sequences produced in this study, 26 haplotypes were recorded of which 17 were detected for the first time ( Fig. 2 (Fig. 2); 14 haplotypes were identified as the endemic minnow (P. lumaireul sensu stricto; clade 1a in Palandačić et al. 2017), 10 as Phoxinus septimaniae Kottelat 2007(clade 12 in Palandačić et al. 2017) and two as Phoxinus csikii Hankó 1992 (clade 5a in Palandačić et al. 2017). Uncorrected p-distance values between the Phoxinus lineages ranged between 2.5 and 6.1%, with the closest distance between P. lumaireul and P. csikii, and the largest between P. lumaireul and P. septimaniae (Fig. 2).
The P. lumaireul, P. csikii and P. septimaniae lineages revealed contrasting distributions. In the four lowland streams, only endemic P. lumaireul were found, except a single individual recognised as P. csikii in the Varaita stream (VAR; Table 1). In the Alpine lakes, however, different Phoxinus cryptic species were detected (Table 2): all three Phoxinus species were present together in two lakes only (ARP and MAD), endemic P. lumaireul were allopatric in four lakes (VIA, AFF, RES and ORG), and sympatric with the aliens P. csikii or P. septimaniae in five and  resource Table ESM 3). Haplotype labels are reported beside each branch where bold labels indicate haplotypes found in the study and asterisks where haplotypes corresponded to sequences already present in GenBank. Clades attribution follows Palandačić et al. (2017). NJ and ML bootstrap (below branches in cursive and above branches in bold respectively) and BI posterior probability (above branches) values are given next to relevant nodes four lakes respectively (FIO, VIS, ROB, ANN and TER with P. csikii and CES, CHI, CAN and ISC with P. septimaniae; Table 2).
The mtDNA COI network analyses also detected the three Phoxinus species, highlighting three different structures (Fig. 3). The 14 mtDNA P. lumaireul haplotypes revealed a haplotype group with the most frequent haplotype (AP01) being central that was linked to three other frequent haplotypes by up to two mutational steps (AP02, AP05 and AP26) (Fig. 3). The remaining 10 haplotypes were shared by only a small number of individuals (one to three) and were linked to AP01 by one to three mutational steps (Fig. 3). AP01 and AP02 were more frequent in lowland streams, with lower frequencies in the Alpine lakes ( Fig. 3; Table ESM 5). Conversely, AP05 and AP26 were more strongly represented in the Alpine lakes, with AP26 only found in these populations ( Fig. 3; Table ESM 5). Of 10 P. septimaniae haplotypes, two were widespread (AP07 and AP53) and were connected by two mutational events, and six were single haplotypes that were separated from the two more frequent ones by between one and six mutational steps (Fig. 3). There were 16 unobserved haplotypes also detected (Fig. 3). P. csikii group comprised of only two haplotypes, with one widespread (AP03) and the other unique (AP77).
In P. lumaireul, there were similar levels of both haplotype (H) and nucleotide genetic diversity indices (p) between the two ecosystem types (H = 0.65 (stream) and 0.70 (lake); p = 0.18% (stream) and 0.23% (lake); Tables 1 and 2 respectively). For genetic diversity levels across the natural and introduced P. lumaireul populations, the highest levels of genetic diversity (H [ 0.60) were from one lowland stream population (VAR; Table 1) and four Alpine lakes (VIA, CES, CHI and ISC) ( Table 2). There was no genetic variability in seven populations (one lowland stream (Table 1) and six Alpine lakes (Table 2)). In the alien minnows, genetic variability was the highest in P. septimaniae and the lowest in P. csikii (Table 2). In P. septimaniae, diversity indices were high (H [ 0.70; p [ 0.40%) in all those populations with more than one specimen present, except for two populations (ARP and ISC; Table 2). In P. csikii, genetic variability was only evident in two out of nine populations (Table 2).

Nuclear DNA variability and MSN
There were 44 haplotypes identified in the 852 bplength multiple alignment of the 194 RAG1 alleles that were obtained by phasing the 65 nDNA sequences produced in this study and the 32 sequences retrieved from GenBank (Online resources Table ESM 5); 45 polymorphic sites were detected of which 12 were singletons and 33 were parsimony informative sites. Thirty-three single nucleotide polymorphisms were detected, whilst no insertions or deletions were observed. Of the 44 haplotypes, 17 haplotypes were detected between the sequences produced in this study and of these, two (Hap2 and Hap16) were already present in GenBank ( Fig. 4; Table ESM 6). They corresponded to P. lumaireul haplotypes (accession N°MF408060 and MF408066, respectively;Palandačić et al. 2017). Indices of genetic variability ranged   Table ESM 1 between 0.21 and 0.42 (H), and between 0.03% and 0.10% (p). PEL and PES populations were the least variable, whilst VAR was the most variable (Table 1).
The network analysis revealed a complex structure in which haplotypes of different mitochondrial lineages did not form a group but were instead spread with varied interconnections. Many of the sequences obtained in the study grouped within the most frequent haplotype (Hap2) that was shared only by P. lumaireul specimens ( Fig. 4; Table ESM6). Three haplotypes were separated by four (Hap6) and five (Hap7 and Hap9) mutation steps from Hap2, a high number that may be representative of a different species, although the weak correspondence between RAG1 haplotypes and mitochondrial lineages prevented further interpretation of this result.

Discussion
The results revealed, for the first time, the presence and co-existence of three European species of the Phoxinus genus in the westernmost part of the Po River basin of Italy: P. lumaireul, P. csikii and P. septimaniae. Whilst the three species have minimal morphological differences (Kottelat and Freyhof 2007), they formed independent and well-supported groups that were also in agreement with recent phylogenetic studies (Palandačić et al. 2017;Corral-Lou et al. 2019). P. lumaireul Schinz (1840) (subclade 1a sensu Palandačić et al. 2017) has been revalidated as an endemic to the north Adriatic basin (Kottelat and Freyhof 2007;Palandačić et al. 2017Palandačić et al. , 2020. Phoxinus csikii (subclade 5a sensu Palandačić et al. 2017) was first described by Hankó 1922 from a karstic brook located in Montenegro and has been recently revalidated, based on both molecular (Palandačić et al. 2017) and morphological evidence (Ramler et al. 2016). Its natural distribution is the lower and middle part of the Danube River watershed of the Balkan Peninsula (Palandačić et al. 2020). P. septimaniae was originally described by Kottelat (2007) as native to the French Mediterranean basins (south-eastern France); it has been also detected in river basins in north-eastern Spain (Corral-Lou et al. 2019;Garcia-Raventós et al. 2020). Although its natural or introduced origin could not be unambiguously discerned in this area of Spain, its geographic distribution resembles the Catalano-Provenzal endemics Squalius laietanus and Barbus meridionalis (Doadrio 2002) that do not naturally populate rivers of the north Adriatic basin. The endangered white clawed crayfish Austropotamobius pallipes occurs both in the Mediterranean French basins and Northern Italian basins, although its presence in the Italian Alps has been suggested as due to human translocations (Stefani et al. 2011). Here, the molecular results revealed P. septimaniae was not recorded in the lowland streams, corroborating that this species is not native to the western Po River basin.
A nuclear marker (RAG 1) was included in analyses, particularly for the lowland stream populations, in order to strengthen the results obtained from the mitochondrial DNA. The RAG 1 results revealed only a few of the least frequent haplotypes that were separated from the main haplotype by several mutational steps (four to five), could possibly belong to other species than the native P. lumaireul. The utility of RAG 1 marker was limited in depicting the phylogenetic relations between minnows that produced an uninformative network where well-defined mitochondrial lineages were not clear, as was evident in other studies (Corral-Lou et al. 2019;Palandačić et al. 2020). Although this could be due to the tendency of this species to hybridize (Corral-Lou et al. 2019), issues including incomplete lineage sorting and artefacts resulting from the gametic phasing or sequencing processes cannot be excluded from the evaluation of these results (Palandačić et al. 2020).
Despite the results from the nuclear DNA were unclear, the results of the mtDNA analyses revealed two important considerations on the distribution and community composition of minnows. Firstly, P. lumaireul was recorded from all sampled sites in the lowland rivers, where they were considered as natural, endemic populations. However, the presence of the exotic lineage P. csikii was detected in one stream, underlying a risk of a 'hidden invasion', even in these non-managed populations. Secondly, in the Alpine lakes, there is now a complex of cryptic Phoxinus species present that have resulted from unregulated introduction events. Although the introduction  resource  Table ESM 5). Circle size is proportional to the observed haplotype frequencies and black points represent unobserved haplotypes and potential intermediates. Colours indicate mtDNA species attribution pathways of these minnows can only be speculated, it is likely that they occurred either unintentionally within contaminated batches of salmonid fishes being released into the lakes to enhance angling, or from their use as live-baits by anglers; both pathways are known sources of invasive minnows elsewhere in Europe (Museth et al. 2007;Miró and Ventura 2015). Patterns of fish introductions in high mountain lakes are often difficult to reconstruct due to a lack of associated information, although evidence suggests that in the Italian Alps at least, angling activities have been popular since the 1960s (Cantonati et al. 2006). The genetic structure and variability between the three cryptic species potentially provide valuable information about their introduction processes. Despite P. septimanae and P. csikii being recorded in a similar proportion of the lakes (eight and nine respectively), P. septimanae had the highest genetic variability, whereas P. csikii had a relatively simple genetic structure. P. lumaireul also showed a substantial gradient of genetic variability (negligible to high). Several studies have revealed invasive populations having much higher diversity than their native range due to local adaptation and/or admixture between different invasive lineages (Lawson Handley et al. 2011;Bock et al. 2015). Distinct genetic patterns observed for the three cryptic minnows here thus suggest multiple and/or older introductions for P. lumaireul and P. septimaniae, in contrast to a recent and/or single introduction event for P. csikii.
Although currently confined to the Alpine lakes, there is an inherent risk that the alien P. csikii and P. septimaniae could be unintentionally dispersed into streams in the Po basin by anthropogenic translocations (e.g. illegal stocking, Fernández et al. 2019), as has occurred in Germany (Knebelsberger et al. 2015), Austria and Croatia (Palandačić et al. 2020). If this happens, it is then probable that they would colonise streams where only endemic P. lumaireul is currently present. Given the phenotypic plasticity and morphological similarity of these minnows, there is then a threat that researchers and practitioners fail to detect their invasion unless genetic studies are completed, and this might only occur once the minnows have established and dispersed, a point when management options to conserve the native communities become inherently difficult (Britton et al. 2011). Indeed, identifying between native and introduced fishes in the field is a problem common to many invaded communities of cryptic cyprinid fishes, including between endemic and alien Barbus fishes in Italy (Zaccara et al. 2019), and native and alien Carassius in some European countries (Hänfling et al. 2005).
In conclusion, these results provide information on the composition of this cryptic complex of Phoxinus species and highlights a concerning conservation scenario for the lowland stream fish assemblages. The potential spread of these alien and cryptic minnows from highly managed ecosystems can contribute to the loss of integrity of freshwater communities that are already highly impacted by other stressors. The resolution of these issues thus requires on-going cooperation and collaboration between researchers, practitioners and regulators to ensure that introductions are prevented and, should this fail, then ensure the released fishes are rapidly identified, and then managed appropriately and according to their risk in the environment (Rytwinski et al. 2019).