Genetic and morphological analyses reveal a complex biogeographic pattern in the endemic barbel populations of the southern Italian peninsula

Abstract The Italian peninsula is a biodiversity hotspot, with its freshwater fish fauna characterized by high levels of local endemism. Two endemic fluvio‐lacustrine fishes of the genus Barbus (barbel, family Cyprinidae) have allopatric distributions in the Tyrrhenian and Adriatic basins of Italy. Barbus plebejus inhabits the mid‐ to northern Adriatic basins, while B. tyberinus is widespread in all central‐northern basins draining into the Tyrrhenian Sea. For basins in Southern Italy draining into the southern parts of these seas, there remains a knowledge gap on their barbel populations due to no previous genetic and morphological studies, despite their apparent biogeographic isolation. Correspondingly, this study quantified the presence and distribution of barbels in the Adriatic and Tyrrhenian basins of Southern Italy through genetic and morphological analyses of 197 fish sampled across eight populations. Testing of how local isolation has influenced the evolution and persistence of these populations was completed by examining sequence variation at two mitochondrial loci (cytochrome b and D‐loop) and performing geometric morphometric analyses of body shape, plus measuring 11 morphometric and meristic characters. Phylogenetic and morphological analyses revealed the presence of two genetically distinct lineages that differed significantly from adjacent B. tyberinus and B. plebejus populations. These two new taxa, here described as SI1 and SI2 Barbus lineages, are highly structured and reflect a complex mosaic biogeographic pattern that is strongly associated with the underlying hydrographical scenarios of the basins. The geographic isolation of these basins thus has high evolutionary importance that has to be considered for maintaining endemism.


| INTRODUC TI ON
The species richness of southern European freshwaters, including the peri-Mediterranean area, is higher than in central and northern Europe, resulting in these freshwaters having high conservation value (De Figueroa, Fenoglio, & Sanchez-Castillo, 2013).
Biogeographically, the region is highly structured with, for example, the freshwater fish diversity between Southern Europe and Northern Africa comprising 23 different ecoregions (Abell et al., 2008;Geiger et al., 2014). Within this, more than 50 native freshwater fish are currently listed as present in the Italian peninsula (Bianco, 2014). The presence of a large number of rare taxa within this relatively small area was strongly influenced by geological and hydrological events during the glacial cycles of the Pleistocene (Bianco, 1995b(Bianco, , 1998Hrbek & Meyer, 2003). These events resulted in the formation of three distinctive ichthyo-geographic districts that are characterized by distinct evolutionary histories in species of the Cyprinidae family (Bianco, 1990(Bianco, , 1995a. To date, fish biogeographic studies in the Italian peninsula have generally focused on the northern and central regions (e.g., Buonerba et al., 2015;Carosi, Ghetti, Forconi, & Lorenzoni, 2015;Livi et al., 2013;Marchetto, Zaccara, Muenzel, & Salzburger, 2010;Meraner et al., 2013;Stefani, Galli, Zaccara, & Crosa, 2004;Zaccara et al., 2019;Zaccara, Stefani, & Delmastro, 2007). These studies have centered on the Padano-Venetian (PV) district of the Italian northeast region, including basins flowing into the upper and middle Adriatic Sea (north of the Vomano River in Abruzzo Region and the Krka River in Croatia), and on the Tuscano-Latium (TL) district of central western region, including all basins draining into the middle Tyrrhenian Sea (Bianco, 1990(Bianco, , 1995a. Conversely, the Apulo-Campano (AC) district of the southern region of Italy, which includes all basins flowing into southern Adriatic, southern Tyrrhenian, and Ionian seas (Bianco, 1990(Bianco, , 1995a Figure 1), has received little research attention. For studies that have been completed, evidence suggests the AC district has long been isolated, and so might have been less influenced by lowered sea levels that occurred during Pleistocene period than basins further north (e.g., Bianco, 2014;Ketmaier et al., 2004), such as the paleo-Po drainage (Bianco, 2014;Buonerba et al., 2015;Livi et al., 2013;Stefani et al., 2004;Zaccara et al., 2019).
Testing the evolutionary effects of the isolation of the southern Italian hydrographic basins, and the potential patterns and processes relating to vicariance events and local dispersal, can be completed using their cyprinid fish communities, as these generally show strong patterns of local endemism (Avise, 2000;Kottelat & Freyhof, 2007;Reyjol et al., 2007;Zardoya & Doadrio, 1999). Cyprinid fishes are widespread throughout all peri-Mediterranean districts, but have limited capability of moving between hydrographic basins due to impassable watershed boundaries, coupled with low saline tolerances that result in coastal areas being effective barriers to their mixing. Among cyprinid fishes, barbels (species of the genus Barbus) have been used widely to study regional biogeography patterns and dynamic changes in continental and inland waters due to their marked diversity, wide distribution, and varied ecology (Buonerba et al., 2015;Gante, 2011). The genus Barbus includes species adapted to a variety of freshwater habitats, ranging from small mountain streams to large and slow-flowing rivers and lakes (Kottelat & Freyhof, 2007).
To fill this considerable knowledge gap on the endemism of barbels in the AC district, the aim here was to test how local hydrographic history has influenced the evolution and persistence of the fluvio-lacustrine barbels in the southern Italian peninsula.
Mitochondrial sequence data and morphological analyses were applied to examine the extent of diversification of the barbels in the AC district compared with barbel populations in northern and central Italy. The results were then used to construct further hypotheses based on biogeographic scenarios that might have influenced patterns of endemism in the southern Adriatic and Tyrrhenian Sea hydrographical networks.

| Sampling
A total of 197 specimens of Barbus spp. were sampled in AC district between 2017 and 2018 with local authority permission. Fish were sampled from three sites in the Tyrrhenian basins and from five sites in the Adriatic basins. The Tyrrhenian sites were the basins Liri-Garigliano (T1) and Volturno (T2), both close to TL district boundary, and Sele (T3) basin, located in the southern part. The Adriatic sites were in the Aterno-Pescara (A1) basin that represents the first Adriatic drainage in AC district, and the Sangro (A2), Biferno (A3), Fortore (A4) up to Ofanto (A5) basins (see Table 1; Figure 1).
Sampling of the fish was completed using electric fishing.
Captured specimens were removed from the water and then held in aerated tanks of water. Under general anesthesia (MS-222), the fish were attributed to a species according to their phenotypic characteristics (e.g., colouration pattern, spot form and size, fin color; Kottelat & Freyhof, 2007;Lorenzoni et al., 2006), enabling recognition of the B. tyberinus phenotype as per Bianco (1995b). Each fish was then measured (fork length, nearest mm), and a biopsy of the anal fin was TA B L E 1 Sampling site locations (expressed with ID code), watershed, river basin, and the number of individuals of each species sampled by site attributed through D-loop mtDNA phylogenetic tree Note: The sampling site position, geographic coordinates, and barbel composition have been indicated also in Figure 1.
taken, preserved in 90% ethanol, and stored at 4°C for subsequent DNA extraction. For morphological analyses, fish were also photographed (left side) using a Nikon D300 camera (24-85 mm lens) positioned by means of a tripod on a table with a millimetric scale. The fish were then placed into another aerated water tank and, following their recovery to normal behavior, were released back into the river.

| Molecular data
Total genomic DNA was extracted from all individuals using a proteinase K digestion, followed by sodium chloride extraction and ethanol precipitation (Aljanabi & Martinez, 1997). Two sets of primers were used to amplify mitochondrial control region (D-loop) and cy- (Livi et al., 2013). D-loop sequences were ob-

| Phylogenetic analyses
Multiple alignments of all sequences were automatically carried out through ClustalW within Bioedit software (Hall, 1999), with polymorphic sites then checked manually. Identical sequences were collapsed into haplotypes in order to facilitate computational processes, as implemented in DnaSP v 5.0 (Librado & Rozas, 2009) software.
Computation of mitochondrial phylogeny was performed independently for each gene on nonredundant haplotypes and on combined cyt b and D-loop fragments dataset. For all phylogenetic analyses, two different phylogenetic inference methods were used as follows: maximum likelihood and Bayesian analyses. The former was conducted in GARLI v 2.0 (Bazinet, Zwickl, & Cummings, 2014;Zwickl, 2006)  GTR + I (Lanave, Preparata, Sacone, & Serio, 1984;Rodriguez, Oliver, Marin, & Medina, 1990) and HKY85 (Hasegawa, Kishino, & Yano, 1985) for cyt b and D-loop, respectively, and HKY85+I+G (Hasegawa et al., 1985) for the combined dataset. The GARLI tree searches were performed under the default settings. Support was assessed with 1,000 bootstrap replicates in GARLI, under the same settings as the best-tree searches. The resulting bootstrap support values were mapped onto the maximum likelihood phylogeny using PAUP software (Swofford, 2002  using PAUP software (Swofford, 2002) and used as a surrogate for levels of species divergence (Doadrio, Carmona, & Machordom, 2002).

| Minimum spanning network, genetic diversity, and demography
A minimum spanning network was created from the multiple D-loop sequences alignment produced in this study using a statistical parsimony criterion as implemented in PopART v 1.7 software (Leigh & Bryant, 2015). The levels of genetic variation within any new endemic lineages were then calculated by estimating nucleotide differences and haplotype diversity using DnaSP v 5.0 software. To visualize their historical demographic trends, mismatch analyses were performed, as implemented in Arlequin v 3.5 (Excoffier & Lischer, 2010) software, testing the sudden demographic expansion model through sum-of-squared deviation values (SSD) in a coalescent algorithm simulation over 1,000 pseudo-replications with statistical significance (p < .05). To test the isolation between populations (within and between Tyrrhenian and Adriatic basins), population genetic differentiation was calculated using the fixation index Φ ST (Weir & Cockerham, 1984) and its significance assessed (p < .05) by permuting haplotypes between populations 3,024 times, as implemented in Arlequin v 3.5.

| Morphological data
The morphology of the barbel specimens was analyzed by measuring seven morphometric and four meristic traits as per Zaccara et al. . Nonshape variation, introduced through variation in position, orientation, and size, was mathematically removed using generalized procrustes analysis, as implemented in MorphoJ software (Klingenberg, 2011). Shape variations were then analyzed by canonical variate analyses (CVA). Mahalanobis distances were calculated using permutation tests (10,000 replicates). Morphometric traits were standardized to the overall mean standard length (Beacham, 1985) to reduce the effects of size and allometry. Pairwise comparison on morphological traits was then recorded between taxa and between populations by performing the analysis of variance (ANOVA) followed by the Tukey post hoc test. These analyses were carried out using PAST software (Hammer, Harper, & Ryan, 2001).

| Multiple alignments and phylogeny
Across the 197 barbels, 26 haplotypes were identified in the 871 bp length of the multiple D-loop alignment, of which 19 were new and deposited in GenBank (under Accession numbers: MK728797-MK728815) as detailed in  Table 2).

| Networks, genetic diversity, and demography of South Italy lineages
In the network analyses of the complete mitochondrial D-loop dataset, the SI1 and SI2 Barbus lineages (N = 181) were linked by more than 13 mutational steps and revealed some distinct patterns. The

| Haplotype distribution and population structure
In the AC district, the SI1 and SI2 Barbus lineages showed an allopatric distribution. The SI1 Barbus lineage was recorded in middle see Figure 1 and Table 1). Genetic differentiation between the SI1 Barbus lineage of the three middle Adriatic populations revealed high genetic structure, with significant ф ST values over 0.39 (p < .01;  Figure 4) and did not show significant differentiation (p > .05; Table S1).

| Morphological pattern among lineages and among populations
The geometric morphometric analyses of the CVA plot revealed there was partial visual separation in body shape morphology in the two SI Barbus lineages ( Figure 5) Figure 1)  (Table S3), inhabiting two basins located at similar latitude but on the opposite sides of the Italian peninsula ( Figure 1).  Table 2). Circle size is proportional to haplotype frequencies. Colors indicate Adriatic (A1 = Aterno-Pescara; A2 = Sangro; A3 = Biferno; A4 = Fortore; A5 = Ofanto) and Tyrrhenian (T1 = Liri-Garigliano; T2 = Volturno; T3 = Sele) populations may be related to local differences in low sea level drainage patterns, although differences in habitats and in biotic interactions might also have been involved.

| D ISCUSS I ON
The results of the population genetic structure have also demonstrated a nonhomogeneous history in the AC basins, showing the presence of unexpected biogeographic boundary that crossed the Apennine watershed. Across the Italian peninsula, the mosaic biogeographic pattern of the genus Barbus was likely to be associated with the differing hydrographic structure of the basins. For example, the SI1 Barbus lineage appeared to originate and only be maintained in basins A1 to A3 (Pescara River up to Biferno River of the middle Adriatic). These basins were not part of the paleo-Po expansion (Bianco, 1990), and so they remained isolated from the widespread dispersion of B. plebejus that occurred in the upper Adriatic basins (c.f. PV district). Within this restricted area, the SI1 Barbus lineage had high levels of genetic variability and was thus highly structured.
These results suggest that climatic, hydrological, and geological factors probably shaped their local isolation and did not result in dispersion events via temporary connections (Forneris, Merati, Pascale, Perosino, & Tribaudino, 2016). Although the hydrogeographic layout of the AC region is congruent with the current topographic and geological pattern, the main distribution of watercourses has also been influenced by its lithological structure from previous geomorphological stages (Amato, Cinque, & Santangelo, 1995 has shown an asymmetric profile of the watershed line, with a retreat of the Tyrrhenian side and progression of the Adriatic side (Brancaccio & Cinque, 1992;Brancaccio et al., 1991). The temporary change in the draining path occurred between Sele (T3) and Ofanto (A5) basins, promoted by temporary river capture events or transitory mountain lakes, that might help explain the actual distribution of the SI2 Barbus lineage in both the southern Tyrrhenian basins (from T1 to T3; i.e., from Liri-Garigliano to Sele basins) and the southern Adriatic basins (A4 and A5; i.e., Fortore and Ofanto basins; Alvarez, 1999), as also reflected by the absence of genetic structure.
Regarding the congruence of the genetic and morphological data, these Italian fluvio-lacustrine barbels, representing a complex of cryptic species, were only partially identifiable by morphology, with their morphological and molecular divergence not always well correlated across the species (Bianco, 1995b;Kottelat & Freyhof, 2007;Livi et al., 2013;Lorenzoni et al., 2006;Zaccara et al., 2019).
Despite this lack of congruence between the genetic and morphological approaches, there was nevertheless some significant correlation between evolutionary lineages and body shape. The two SI

30.9
Note: Data of morphometric traits were transformed according to Beacham (1985) formula. ANOVA results (F) showing differences among the five Barbus groups are also reported; all p-values were <.001.
remain undetected along this complex of cryptic species and will potentially lead to the loss of local endemism. Consequently, these results highlight the necessity for any fish and fishery management programmes in this region to recognize the inherently high conservation value of these endemic barbels and avoid undesirable mixing with other barbels through, for example, fish stocking exercises.

ACK N OWLED G M ENTS
This study was supported by grants from the University Research Fund (FAR-University of Insubria) to SQ, IV, and VDS. We thank the numerous students who assisted in the morphological analysis and participated in the molecular laboratory activities. We also thank the anonymous reviewers for the suggestions and comments proposed that have improved the quality of the research.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R S ' CO NTR I B UTI O N S
SZ, ML, SQ, and AC conceived the ideas and designed the methodology; ML and AC collected the data; SQ analyzed the photographs and performed morphological analysis; IV and VDS worked in laboratory and analyzed the data; and all authors contributed critically to the drafts and gave final approval for publication. F I G U R E 6 Canonical variate analysis (CVA) output of the body shape comparison between the eight populations of Barbus considered in the present study (see Figure 1). Wireframe graphs indicate the shape changes along each axis (from gray to dashed black)