Quaternary geomorphological and climatic changes associated with the diversification of Iberian freshwater fishes: The case of the genus Cobitis (Cypriniformes, Cobitidae)

Ecology and Evolution. 2022;12:e8635.   | 1 of 22 https://doi.org/10.1002/ece3.8635 www.ecolevol.org Received:22December2021 | Revised:20January2022 | Accepted:25January2022 DOI: 10.1002/ece3.8635 RESEARCH ARTICLE Quaternary geomorphological and climatic changes associated with the diversification of Iberian freshwater fishes: The case of the genus Cobitis (Cypriniformes, Cobitidae) Andrea CorralLou1,2 | Silvia Perea1,3 | Anabel Perdices1 | Ignacio Doadrio1 ThisisanopenaccessarticleunderthetermsoftheCreativeCommonsAttributionLicense,whichpermitsuse,distributionandreproductioninanymedium, providedtheoriginalworkisproperlycited. ©2022TheAuthors.Ecology and EvolutionpublishedbyJohnWiley&SonsLtd. 1BiodiversityandEvolutionaryBiology Department,MuseoNacionaldeCiencias Naturales,CSIC,Madrid,Spain 2ConsultoresenBiologíadela ConservaciónS.L.,Madrid,Spain 3InstitutodeBiología,Departamentode Zoología,UniversidadNacionalAutónoma deMéxico,CiudaddeMéxico,México Correspondence AndreaCorral-Lou,Biodiversityand EvolutionaryBiologyDepartment,Museo NacionaldeCienciasNaturales,CSIC, JoséGutiérrezAbascal,2;28006Madrid, Spain. Email:[email protected] Funding information IndustrialDoctorateofCommunityof Madrid,Grant/AwardNumber:IND2017/ AMB-7699;DueroHydrographic Confederation,Grant/AwardNumber: Life13nat/es/000772 Abstract WestudiedthepopulationgeneticstructureofCobitis vettonica,anendangeredfreshwaterfishspeciesendemictotheIberianPeninsula,inordertoproposeabiogeographicmodeloftheresponsesofspeciestothemultiplechangesthatoccurredin theIberianhydrologicalsystemduringtheQuaternaryperiod.Wealsodeciphered therelationshipofC. vettonica with its sister species C. paludica,particularlyinsympatricareas,andprovidegeneticinformationforconservationpurposes.Toachieve thisgoal,weanalyzedbothmitochondrialandnucleardata(thecytochrome band the nuclear recombination activating 1genes)andabatteryofsingle-nucleotidepolymorphisms(SNPs)of248individualsofC. vettonica or C. paludicafrom38localities, includingsomesympatricones,coveringtheentiredistributionareaofC. vettonica. Wehighlighttheimportantroleplayedbythehydrogeomorphologicalprocessesand climaticchangesthatoccurredintheIberianPeninsuladuringtheQuaternaryonboth thepopulationstructureofC. vettonicaanditsrelationshipwithitssisterspeciesC. paludica.Ourresultssupportthegeneticintrogressionofpopulationsattheeastern limitofthedistributionofC. vettonica.Furthermore,wepostulategeneticintrogressioninsympatricareas.Finally,weproposetheestablishmentorexpansionoffour OperationalConservationUnits(OCUs)forC. vettonica,andhighlightthethreatfaced byitspopulationsduetothelowlevelofgeneticdiversitydetectedforsomeofits populationsandgeneticintrogressionwithC. paludica,whichcouldeventuallydisplaceC. vettonica,resultinginalossofdiversityinthisspecies. KEYWORDS Cobitis,hybridization,IberianPeninsula,populationstructureandconservation,quaternary TAXONOMY CLASSIFICATION Populationgenetics 2 of 22 | CORRALLOU et AL. 1 | INTRODUCTION The Quaternary is a geological period characterized by glacial– interglacial cycles that have dominated the global climate since 2.58Myatothepresentera(Gibbardetal.,2010;Pillans&Naish, 2004),whichhavehaddrasticconsequencesontheevolutionofthe biotaofmanyregions.Duetoadvancesinthefieldofphylogeography,wehaveabetterunderstandingoftheresponsesoforganismstoPleistoceneevents(Hewitt,2004;Weiss&Ferrand,2007). ClassicalstudiessupporttheroleofMediterraneanpeninsulasas refugesforfaunaduringtheQuaternary,whichprovidedthestock forrecolonizationsofnorthernandcentralEurope(Hewitt,1996; Taberlet&Bouvet,1994).TheIberianPeninsulaisconsideredoneof themostimportantrefugesduringthisperiod,asshownbyseveral studies,mainlyofterrestrialfauna(Querejetaetal.,2017;Valdiosera etal.,2008). Primaryfreshwaterfishes(i.e.,fishesstrictlyintolerantofsaltwater; Myers, 1966) have limited dispersal abilities and are often restrictedtospecifichydrographicbasins.Dispersalbetweenbasins during the Quaternary was possible mainly through downstream connectionscausedbythedecreaseinsealevelduringcoldperiods orbyupstreampiracy(Corral-Louetal.,2019;Mesquitaetal.,2005; Perea&Doadrio,2015).TheformationoftheIberianfluvialnetwork alsoculminatedduringthisperiod,whichaffectedtheregion'shydrogeomorphologyandthereforethecurrentevolutionarypatterns foundinprimaryfreshwaterfauna(Alonsoetal.,2007;Corral-Lou etal.,2021;Paisetal.,2012;Pereaetal.,2016).Likewise,thecurrent populationstructureofmanyIberianfreshwaterfisheshasbeenattributedtotheinteractionofvariousnaturalfactorsthatoccurred duringtheQuaternarysuchasthedryingupofbodiesofwater,sea- levelfluctuations,tsunamis,streampiracy,isolationofbasins,hydrogeomorphologicalchanges,andsecondarycontactoftwodifferent basins(Casal-Lópezetal.,2017;Corral-Louetal.,2019;Ganteetal., 2009;Pereaetal.,2016). Despiteeffortsmadeinthelastseveraldecades,thereisstill muchtobelearnedabouttheimpactofQuaternarychangesonthe evolutionary processes and patterns of diversification of Iberian freshwater fish populations. Especially those species with a restricteddistributionrangesincetheywereprobablymoreaffected byclimaticandgeologicalchangesduringtheQuaternarythanspecieswithlargerranges.Inordertoaddresssomeopenevolutionary questions, such as the role of hydrogeomorphologic changes and genetic introgression in species evolution, we analyzed populationsofCobitis vettonicaasacasestudy.Thisspeciesisanendangered Iberian freshwater fish whose distribution is restricted toafewriversintheTagusandDuerobasinsinthemid-western IberianPeninsula(Figure1).ItinhabitsriverswithlowpHandwater hardness levels, andgravel androckybottoms,andismorecommonlyfoundintheheadwatersoftheserivers(Carmonaetal.,1999; Collares-Pereiraetal.,2021;Doadrioetal.,2011;Perdices&Coelho, 2020).Incontrast,thesisterspeciesofC. vettonica,Cobitis paludica (Doadrio&Perdices,2005;Perdices&Doadrio,2001),isageneralistthatinhabitsmostoftheIberiandrainagesincludingTagusand Douro,withapreferenceforstreamsclosetothemainchannelwith highsuspendedsolids,highwaterhardness,lowtransparency,low currentvelocityandmuddybottoms(Carmonaetal.,1999;Doadrio etal.,2011).SympatriczonesofC. vettonicaandC. paludicahave beenreportedinthewesternlimitsofthe distributionareaofC. vettonicabutgeneticintrogressionbetweenthemhasnotbeenreported(Perdices&Coelho,2020).However,insomeeasternpopulationstherehasbeenmentionofgeneticintrogressionbetweenC. vettonicaandC. paludica,butnomoreinformationhasbeendetailed (Doadrio et al., 2011, 2021). Sister species of Iberian freshwater fishesgenerallyhaveanallopatricdistributionthatwasestablished mainlybeforetheQuaternary(Doadrio,1988;Sousa-Santosetal., 2019).Somespeciesnowshowpatternsofsympatryincertainareas as a consequence of secondary contacts during the Quaternary, whichhasledtogeneticintrogression,asisthecaseofsomespecieswithinthegeneraofLuciobarbusandPhoxinus(Corral-Louetal., 2019;Denysetal.,2013;Machordometal.,1990).Inthecaseof C. vettonicaandC. paludica,phylogeneticstudieshaveshownthey diverged during the Plio-Pleistocene (Doadrio & Perdices, 2005; Sousa-Santosetal.,2014).However,boththeoriginofthesympatric zonesbetweenC. vettonicaandC. paludicaandtheplausiblegenetic introgressionbetweentheminthewesterndistributionareaofC. vettonicaareunknown.Allthesepointsmakefisheswithrestricted distributionareasasexcellentmodelstodeciphertheQuaternary effectsonthediversificationoffreshwater-restrictedtaxawithina glacialrefugeastheIberianPeninsula(Ganteetal.,2009). PhylogeographicstudiesofIberianprimaryfreshwaterfishes have mainly used mitochondrial genes combined with nuclear markerssuchasmicrosatellitesorintrons(Casal-López&Doadrio, 2018;Corral-Louetal.,2019;Gonzalezetal.,2014,2018;Perea &Doadrio,2015).Inrecentyears,phylogeographicandpopulation genetics studies have taken advantage of next-generation sequencing(NGS)technologiestobroadlyscreenthegenomeat highresolution,yetsomechallengesremainintheanalysisofNGS FIGURE 1 Samplinglocalitiesincludedinthisstudy.Theorange shadingdelimitstheDueroBasin,andthegreenshadingdelimits theTagusBasin.NumberscorrespondtothoselistedinTable1 | 3 of 22 CORRALLOU et AL. datasuchasthehandlingoflarge-scaleandcomplexdata,theupstreamprocess of pipeline,andthelimitednumber ofavailable reference genomes (McCormack et al., 2013; Tan et al., 2019). Combined analyses of single-nucleotide polymorphisms (SNPs) with mitochondrial and nuclear markers have proven useful to providemorecompleteandmoredetailedphylogeographicand biogeographic models of the past and present relationships of populationsofvariousspecies (Corral-Louetal.,2021;Mendes etal.,2019;Zarraonaindiaetal.,2012). The main aim of this study is to provide a robust biogeographicmodelforthespeciesC. vettonicathroughoutitsentire distributionasawitnessoftheevolutionoftheIberianPeninsula throughouttheQuaternary.Wealsodeciphertherelationshipof C. vettonica with its sister species C. paludicaatthelimitsofits distributionarea.Inaddition,weassessthepopulationstructure andgeneticdiversityofthespeciesacrossitsdistributiontorevise the Operational Conservation Units (OCUs; Doadrio et al., 1996)previouslyestablishedforC. vettonica(Doadrioetal.,2021) ortoestablishnewones.MoredetailedandaccurateOCUsand decipheringtherelationshipwithitssisterspeciesthroughoutits entirerangearekeyfortheeffectivemanagementandconservationoftheendangeredC. vettonica.Toachievethesegoals,we analyzed the mitochondrial marker cytochrome b(MT-CYB), the nuclearmarkerrecombination activating gene 1(RAG1),andaset ofSNPsobtainedthroughnext-generationsequencingofpopulationsthroughouttheentiredistributionrangeofC. vettonicaand someofC. paludicafromadjacentsub-basins. 2 | MATERIALS AND METHODS 2.1 | Sampling, DNA extraction, amplification and sequencing Wesampled 204individualsofC. vettonicaand60ofC. paludica fromatotalof38localities(Figure1;Table1).The27localities(23in theTagusand4intheDuero)inwhichC. vettonicawasfoundcover theentireknowndistributionareaofthespecies,andincludesits typelocality(ÁrragoRiver,TagusBasin).ThesamplesofC. paludica werecollectedfrom13localities,includingitstypelocality(Tiétar River,TagusBasin):ofthese,11wereeitheradjacenttothoseofC. vettonicaorinsympatriclocalitiesintheTagusBasin,andtwowere intheDueroBasin(Figure1;Table1).Tissuesampleswereobtained from the DNA and Tissue Collection at the National Museum of NaturalSciencesofMadrid(MNCN–CSIC;TableS1).Sequencesof MT-CYBofC. vettonicaandC. paludicaavailableinGenBankwere alsoincludedinthestudy(33ofC. vettonicaand7ofC. paludica; TableS1). For each individual, DNA was extracted from ventral fin tissue using the Qiagen DNeasy® Blood and Tissue Kit (Qiagen, Inc., Valencia, CA, USA), following the manufacturer's protocol. Polymerasechainreaction(PCR)wasusedtoamplify1140bpofMT- CYBfrom208individuals(164ofC. vettonicaand44ofC. paludica), and1500bpofRAG1from90individuals(66ofC. vettonicaand24 ofC. paludica)(TableS1).Amplificationswereperformedfollowing theprotocoldescribedbyDoadrioandPerdices(2005)forMT-CYB, andCorral-Louetal.(2021)forRAG1. All sequences (the new ones and the ones downloaded from GenBank;TableS1)forMT-CYB(248sequences,240fromthepresentstudyand40fromGenBank;TableS1)andRAG1(90sequences fromthepresentstudy;TableS1)werealignedusingMAFFT(Katoh & Standley, 2013), as implemented in Geneious 10.1.3 (https:// www.geneious.com;Kearseetal.,2012),andthenmanuallyexamined. Alleles of RAG1were separated using the PHASEalgorithm (Stephens&Donnelly,2003),asimplementedinDnaSPv.6.10.01 (Rozasetal.,2017). 2.2 | Genotyping and SNP filtering FortheSNPstudy,weselectedeightpopulationsofeachC. vettonicaandC. paludica(Table1)basedonthestructureobserved inapreviousstudyusingmitochondrialdataandavailableDNA (Doadrioetal.,2021).Atotalof106individuals(60ofC. vettonica and46ofC. paludica;TableS1)wereusedtopreparelibrariesfor doubledigestrestrictionsite-associatedDNAsequencing(ddRAD- seq)followingtheprotocoldescribedbyKessetal.(2016).AllDNA sampleswereusedasaninputforacustomlibrarypreparation protocol.Librariesweredual-indexedforpostsequencingdemultiplexing.ThesampleswereruninaNovaSeq6000PE150lane. Trimmomatic0.36(Bolgeretal.,2014)wasusedtoremoveadapters (ILLUMINACLIP option). Using the process_radtags program inSTACKS2.4(Catchenetal.,2013),allreadsweretruncatedto thesamelengthof95bp,andlow-qualityreadswereremoved usingthe-qparameteraccordingtophred33system.Adenovo_ mapanalysiswasalsoperformedinSTACKS,inwhichdifferent programswereruntoassemblelociineachindividual(ustacks), buildacatalogue(cstacks),matchallgenericsamplesagainstthe catalogue(sstacks),andreconstructlociusingR2readsandidentifySNPsusingthemetapopulationinformation(gstacks).Priorto runningtheSTACKSmodules,severaltestswereusedtoidentify whichparametersmaximizedthenumberofSNPsobtainedinat least 80% of the individuals (r80 rule; Paris et al., 2017). Since SNPswereobtainedforbothspecies,weexploredparametersin aconservative wayinthede_novomapsmoduleofstacks.The parameter m (i.e., the minimum depth of coverage required to createastack)wassetto5duetothedepthofcoveragevalues obtained.TheMparameter(i.e.,numberofmismatchesallowed betweenstackswithinindividuals)hastobecarefullyset.Ifitis settoohigh,paralogousornonhomologouslocicanbeincorrectly mergedintothesamelocusandifitissettoolow,homologous locicanbelost.Forthisreason,theparameterMwasexplored between3and5.FollowingtheindicationsofParisetal.(2017), thevaluesofn(i.e.,thenumberofmismatchesallowedbetween stacksbetweenindividuals)wereexploredfrom2to4(TableS2). Finally,theselectedparameterswerem =5,M =3,andn = 2. The 4 of 22 | CORRALLOU et AL. TABLE 1 Informationonthesamplinglocalitiesincludedinthepresentstudy No Species River Basin Subbasin Locality Province mtDNA nDNA SNPs 1Sympatricarea Alfrividas Tagus Ponsul Alfrividas Portugal 4/1 – – 2C. vettonica Aravil Tagus Aravil Alcafozes Portugal 1– – 3C. vettonica Aravil Tagus Aravil CegonhasNovas Portugal 4 – – 4C. vettonica Erjas Tagus UpperErjas Valverdedel Fresno Spain 12 9? 5 5C. vettonica SanMartin Tagus UpperErjas SanMartínde Trevejo Spain 6 1 – 6C. vettonica Trevejana Tagus UpperErjas Cilleros Spain 4 4 – 7C. vettonica Erjas Tagus Middle Erjas Termasde Monfortinho Portugal 6– – 8C. vettonica Arades Tagus Middle Erjas Salvaterrado Extremo Portugal 2– – 9C. vettonica Erjas Tagus LowerErjas Salvaterrado Extremo Portugal 2– – 10 Sympatricarea Erjas Tagus LowerErjas CabeçaQueimado Portugal 1/1 – – 11 Sympatricarea Erjas Tagus LowerErjas AzenhadoRoque Portugal 4/1 – – 12 Sympatricarea Erjas Tagus LowerErjas Serrinha Portugal 2/1 – – 13 C. vettonica Árrago Tagus Western Alagón CadalsodeGata Spain 15 3? – 14 C. vettonica Gata Tagus Western Alagón Gata Spain 19 5 5 15 C. vettonica SanBlas Tagus Western Alagón Gata Spain 1– – 16 C. vettonica Acebo Tagus Western Alagón Hoyos Spain 3 1 – 17 C. vettonica Mayas Duero Águeda Descargamaría Spain 2– – 18 C. vettonica Mayas Duero Águeda ElSahugo Spain 20 410 19 C. vettonica Águeda Duero Águeda LaHerguijuela Spain 1– – 20 C. vettonica Turones Duero Águeda LaBouza Spain 1 1 – 21 C. vettonica Alagón Tagus Eastern Alagón SantiIbañezdela Sierra Spain 14 9 4 22 C. vettonica Francia Tagus Eastern Alagón NavadeFrancia Spain 14 5 8 23 C. vettonica Cuerpode Hombre Tagus Eastern Alagón Sotoserrano Spain 19 910 24 C. vettonica Ladrillar Tagus Eastern Alagón LaRebollosa Spain 1– – 25 C. vettonica Hurdano Tagus Eastern Alagón VegasdeCoria Spain 1 1 – 26 C. vettonica Caparro Tagus Eastern Alagón Zarzade Granadilla Spain 23 78 27 C. vettonica Jerte Tagus Eastern Alagón Navaconcejo Spain 15 510 28 C. paludica Alvito Tagus Ocreza MonteGordo Portugal 1– – 29 C. paludica Aurela Tagus Salor Santiagode Alcántara Spain 5 5 9 30 C. paludica Toula Tagus Aravil Alcafozes Portugal 1– – 31 C. paludica Aravil Tagus Aravil CegonhasNovas Portugal 1– – 32 C. paludica Salor Tagus Salor Membrío Spain 2– 2 33 C. paludica Tamuja Tagus Almonte Trujillo Spain 4 – 4 | 5 of 22 CORRALLOU et AL. programpopulations,alsoinSTACKS,wasusedtofiltertheSNPs usingthefollowingparameters:-p16-r0.80,– min_maf =0.05, – max_obs_het =0.75,and– write_random_snp. TwomatricesfortheSNPswerebuilt.Theallowedpercentage ofmissingdataperlocusandperindividualwas30%(TableS3).The firstmatrixincludedthepopulationsofC. vettonicaanddatafrom60 individualsfrom8localitiesand4000polymorphicloci.Thesecond includedpopulationsofbothC. vettonicaandC. paludica,referred toasC. vettonica + C. paludica,anddatafrom106individualsfrom 16localitiesand4538polymorphicloci. 2.3 | Phylogeny Aphylogenetictreewasconstructedtoassessrelationshipsbased onthecollapsedMT-CYBhaplotypesmatrix.TheselectedsubstitutionmodelswereSYM,HKY+I,andGTR+Gforthefirst,second,andthirdposition,respectively,basedontheresultsobtained inPartitionFinder2(Lanfearetal.,2017)usingtheAkaikeinformationcriterion(AIC;Akaike,1974).Theanalysiswasimplementedin MrBayesv3.2(Ronquistetal.,2012),withtwosimultaneousindependentrunseachwithfourMarkovchainMonteCarlo(MCMC), whichwererunfor5× 107generations.Thefirst25%ofgenerationswereremovedasburn-in.Posteriorprobability(pp)values wereusedtoassessthereliabilityofthephylogenetichypothesis. TwosequencesofotherspeciesofCobitiswereusedasoutgroups: C. bilineata(EF605321.1)andC. zanandreai(AF263089.1).Thegenetic divergencebetweenthelineagesobtained intheMT-CYB analysis was evaluated through uncorrected p-distances with 1000permutationsusingMEGAv.6.0(Tamuraetal.,2013).For the SNP data, phylogenetic relationships were evaluated using thecompleteSNPmatrix(C. vettonica + C. paludica)throughthe MaximumLikelihood(ML)methodintheRaxMLprogramimplementedinCIPRESScienceGatewayv3.3.Theevolutionarymodel selectedwasASC_GTRGAMMA,asrecommendedintheprogram manual for data that contain only variable sites, and the Lewis ascertainmentbiascorrectionwasused(Stamatakis,2016).The robustness of the tree was evaluated with 1000 bootstrap (b) replicates. 2.4 | Genetic structure To examine the genetic structure of all studied populations of C. vettonica and C. paludica, haplotype networks for the two genes wereconstructedusingthemedian-joiningalgorithm(Bandeltetal., 1999)implementedinPopART(Leigh&Bryant,2015).Analysisof molecularvariance(AMOVA)wasusedtodeterminethesourceof thegeneticvariationinMT-CYBforthepopulationsidentifiedasC. vettonica(excludingV4subset,seebelow)usingdifferentgroupings basedonbasin,sub-basin,andevolutionarylineage.TheseanalyseswereimplementedinArlequinv.3.11(Excoffier&Lischer,2010) with 10,000 permutations. Population differentiation in terms of Φst(Hudsonetal.,1992)betweenallpopulationpairs(exceptthose representedbyonlyoneindividual)forthecompletedataset(C. vettonicaandC. paludica)wasalsocalculatedinArlequinv.3.11.Forthe localitiesinwhichmitochondrialhaplotypesofbothC. vettonicaand C. paludicawerefound,theΦstwascalculatedbytreatingtheindividualsofeachspeciesasseparatepopulations. FortheSNPsdata,onlyneutrallociweretakenintoaccountto studythepopulationstructureinordertoavoidthebiasthatcould becausedbycandidatelociforselection.Forthisreason,weevaluated the presence of loci under selection with BayeScan v.2.0 (Foll & Gaggiotti, 2008), using the default parameters except for thepriorodds(priorodds=100;Lotterhos&Whitlock,2014).We removedthelociunderselectionandconstructedtwonewmatricesfortheneutralloci,keeping3998lociintheC. vettonicamatrix and4500lociintheC. vettonica + C. paludicamatrix.STRUCTURE (Pritchard et al., 2000) was used to assess the genetic structure basedonthetwoSNPmatrices.Themostprobablenumberofsubpopulations(K)foreachanalysiswasequaltothenumberofpopulationsstudied+1(K =9andK =17,respectively).Weperformed10 independentsimulationsforeachKwithaburn-inlengthof50,000 and50,000MCMCrepetitionsaftertheburn-in.Themostprobable number of subpopulations for each matrix was estimated by takingintoaccounttheresultsofboththeΔK(Evannoetal.,2005) andthePuechmaillemethod(Puechmaille,2016).Thewebserver StructureSelector(Li&Liu,2018)wasusedtomaketheseestimations.WeselectedKbasedonbothmethodsbecauseagreaterprobabilityofK =2existswhenthestructureisanalyzedwiththeΔK No Species River Basin Subbasin Locality Province mtDNA nDNA SNPs 34 C. paludica Almonte Tagus Almonte Jaraicejo Spain 5 5 5 35 C. paludica Ibor Tagus Lower Tagus CastañardeIbor Spain 3– 3 36 C. paludica Tiétar Tagus Tiétar LaIglesueladel Tiétar Spain 11 5 6 37 C. paludica Yeltes Duero Huebra- Yeltes SanMartínde Yeltes Spain 5 5 8 38 C. paludica Huebra Duero Huebra- Yeltes Elcubodedon Sancho Spain 9 4 9 Notes: Alsodetailedarethenumberofindividualsstudiedforeachofthemarkertypes:mtDNA(MT-CYB),nDNA(RAG1),andSNPs.Underthe mtDNAcolumn,foreachsympatricarea,thetwonumbersindicatethenumberofindividualsanalyzedforC. vettonicaandC. paludica,respectively. TABLE 1 (Continued) 6 of 22 | CORRALLOU et AL. method(Janesetal.,2017),andsimilarresultsfordifferentmethods indicateaclearsignal(Puechmaille,2016). 2.5 | Genetic diversity and demography DifferentgeneticdiversityparameterswereestimatedforbothMT- CYBandRAG1inDnaSPforeachriver,sub-basin,species,evolutionarylineage,andOCUs.Thestudiedgeneticdiversityparameters werehaplotypeornuclearallelesnumber(h or aforMT-CYBand RAG1,respectively),haplotypeorallelicdiversity(Hd or adforMT- CYBandRAG1,respectively),nucleotidediversity(π),andnumberof polymorphicsites(S).ForMT-CYB,thepopulationswithsharedhaplotypesbetweenC. vettonicaandC. paludicawereanalyzedconsideringonlytheindividualsidentifiedasC. vettonica.FortheSNPdata, adifferentsetofgeneticdiversityparameterswasestimatedusing the populationsprograminSTACKS2.4(Catchenetal.,2013).These parameterswereevaluatedforthetwomatricesindependently(C. vettonicaandC. vettonica + C. paludica);however,onlytheresultsfor C. vettonicaareshown.Inaddition,athirdmatrixwasstudiedtaking intoaccountthepopulationsofbothspecies,exceptthepopulations oftheAurelaandSalorriversofC. paludica,duetothecloserelationshipthattheselocalitiesshowedwiththelocalitiesoftheErjas sub-basinofC. vettonica. ToevaluatethedemographyofpopulationsofC. vettonica,deviationsfromamodelofmutation-driftequilibriumforMT-CYB were testedusingFu'sFs(Fu,1997)andTajima'sD(Tajima,1989),asimplementedinArlequinv.3.11. 2.6 | Divergence times and niche modeling DivergencetimesamongthepopulationsofC. vettonica were estimatedusingarelaxedlognormalclockandacoalescentmodelonthe collapsedMT-CYBhaplotypesmatrixinBEASTv1.8.4(Drummond etal.,2012).Themolecularclockwascalibratedusing,asanormal prior,anevolutionaryrateof0.34%divergenceperlineagepermillionyears(Doadrio&Perdices,2005),andthesubstitutionmodel usedwasGTR+ I. ToidentifythepotentialnicheofC. vettonicainbothpastand presentscenarios,ecologicalnichemodels,asappliedinthesoftwareMaxEntv.3.4.1(Phillipsetal.,2017),wereusedtopredicta meanhabitatsuitabilityvalueforeachlocality.Toselectthemost appropriateofthe19accessiblebioclimaticvariablesdownloaded from the WorldClim dataset (Fick & Hijmans, 2017) for the final analysis,correlationandcollinearityanalysesandpreliminaryanalyseswerefirstperformed.Threevariablesbasedontheecological characteristics of the species were chosen for the final analysis: Isothermality (BIO3), Precipitation of Driest Month (BIO14), and PrecipitationSeasonality(BIO15).Thedatamatrixonthepresence ofthespeciesincludedonlythestudiedpopulationsofC. vettonica thatdonothaveanyputativegeneticintrogressionwithC. paludica. Thenichemodelingwasmadeforthreeperiodsoftime:Present(the years1970–2000),LastGlacialMaximum(LGM;~22,000yearsbeforepresent,BP),andLastInterglacial(LIG;~120,000–140,000BP). ThequalityofthemodelwasevaluatedbytheAreaUndertheCurve (AUC)derivedfromtheReceiverOperatingCharacteristic(ROC). 2.7 | Ancestral area reconstructions TheancestralareasofC. vettonicawerereconstructedusingboth theS-DIVA(Yuetal.,2010)andS-DEC(Beaulieuetal.,2013;Ree & Smith, 2008) methods implemented in RASP v.4.2 (Yu et al., 2015,2020).TheestablishedareaswerebasedonthecurrenthydrographicsubdivisionoftheIberianPeninsula.Duetothegenetic differentiationdetectedforMT-CYBin the Alagón sub-basin,this areawasdividedintoaneasternandawesternpart.Thecodefor theareaswasasfollows:A:Águeda,B:WesternAlagón,C:Eastern Alagón,D:Aravil,E:Erjas,F:Ponsul,G:Salor,H:Ocreza,I:Aurela,J: HuebraandYeltes,K:TamujaandAlmonte,L:Ibor,andM:Tiétar.We usedtheBEAST-derivedtreesastheinputfilesforRASP.Weeliminated50%ofthetotalinitialtrees,andused100randomtreesfor theanalysis.Ancestralrangeswerelimitedtoincludenomorethan threeadjacentareas.OnlythereconstructedareasforC. vettonica areshownasthislineageisthefocusofthepresentstudy. 3 | RESULTS 3.1 | Phylogeny ThephylogeneticanalysisofMT-CYBrecoveredtwomainlineages correspondingtothespeciesC. vettonica(S1)andC. paludica(S2) withhighsupport(pp=1;Figure2).ThespeciesC. vettonica(S1) was not monophyletic since many individuals identified morphologicallyasC. vettonicafromtheAlagónsub-basinpopulationswere foundinbothlineages(i.e.,Jerte,Francia,Alagón,Caparro,Cuerpo deHombre[Cdh],andGatarivers).ThelineageofC. paludica(S2) wasalsonotamonophyleticgroupsinceoneindividualidentified morphologicallyasC. paludicafromtheSalorsub-basinwasincluded intheC. vettonicalineage(S1).WithinC. vettonica(S1),therewere twowell-differentiatedlineages(Figure2).ThefirstgroupedpopulationsfromtheAlagónsub-basinintheTagusBasin(Árrago,Gata, SanBlas,Acebo,Alagón,FranciaCdh,Ladrillar,Hurdano,Caparro, andJerterivers)andthosefromtheÁguedasub-basinintheDuero Basin(Turones,Mayas,andÁguedarivers).Withinthislineage,these populationsgroupedintotwodistinctlineages,V1andV2.ThepopulationsfromtheeasternAlagónsub-basin(Alagón,Francia,Cdh, Ladrillar,Hurdano,Caparro,andJerterivers)formedtheV2lineage (pp=0.99),whilethosefromthewesternAlagón(Árrago,Gata,San Blas,andAceborivers)andÁguedasub-basins(Turones,Mayas,and Águeda rivers),and one individual fromJerteRiver(fromeastern Alagónsub-basin)thatwasnotincludedinV2,formedtheV1lineage.ThephylogeneticrelationshipsofthesequencesincludedinV1 werenotresolved.However,withinV2,wedetectedanotherhighly | 7 of 22 CORRALLOU et AL. supportedmonophyleticgrouping(pp=1)inwhichonlythesingle individual studied from Ladrillar River was excluded. The second well-differentiatedlineagewithinC. vettonica(V3;pp=1)included populations from the Ponsul, Aravil, and Erjas sub-basins located alongthewesternlimitofthespecies’distributionrange,andone individualfromSalorRiver(Figure2). GeneticdistancesamongthethreelineagesdetectedinC. vettonica(V1,V2,andV3)rangedfrom0.3%to1.2%(Table2). TheunrootedphylogenetictreeofC. vettonicaandC. paludica basedonthe4538SNPswaslargelycongruentwiththeresultsobtainedinthephylogeneticreconstructionforMT-CYB,althoughwith someslightdifferences(Figure3).Theindividualsofthepopulations analyzed for SNPs belonging to the V1 lineage (i.e., Mayas from Águedasub-basin,DueroBasinandGatafromwesternAlagónsub- basin,TagusBasin)werewelldifferentiatednotonlyfromtherest oftheindividualsanalyzed(b>90)butalsofromeachother(b>90 and75< b >90,respectively).ThosebelongingtotheV3lineage (ErjasRiver),whichwereseparatedfromtherestofthepopulations ofC. vettonica,weremorerelatedtoindividualsofC. paludicafrom AurelaRiver,atributaryoftheTagusBasinthatisgeographically closetotheErjassub-basin.Likewise,individualsofC. vettonicafrom theCaparroriverpopulation(belongingtoV2andS2)weremore related to those of C. paludica from the populations in the Salor, Tamuja,Ibor,Huebra,andYeltesrivers.Theremainingpopulations ofC. vettonicaincludedintheSNPanalysisbelongingtotheV2and S2lineages(Alagón,Cdh,Jerte,andFranciarivers)didnotshowa clearphylogeneticstructurewithseveralintermediateindividuals. 3.2 | Genetic structure ThehaplotypenetworkforMT-CYBwasconsistentwiththephylogenetic reconstruction for MT-CYB. The network showed two main groups separated by 20 mutational steps corresponding to the species C. vettonica and C. paludica (S1 and S2, respectively; Figure4).ForC. vettonica(S1),thenetworkwasmoreinformative FIGURE 2 PhylogenetictreebasedonBayesianinferenceforthemitochondrialmarkerMT-CYB.Posteriorprobabilityvaluesare indicatedabovebranches.S1:Cobitis vettonica.Theverticalcoloredbarsrepresentthethreemainlineages(V1–V3)ofthespeciesandtheir localities(seelegend).S2:Cobitis paludica.Also,themapisshownindicatingthesamplinglocationsreferencedinTable1.Thecolorofthe samplingpointindicatesthedetectedlineageineachpopulationbasedontheresultsofthephylogenetictree TABLE 2 Uncorrectedabsolutegeneticdistancesbetweenthe mainlineagesofC. vettonica(V1,V2,andV3)andC. paludica(S1). Thestandarddeviationisshownontheupperdiagonal V1 V2 V3 S2 V1 0.002 0.003 0.005 V2 0.003 0.003 0.004 V3 0.011 0.012 0.005 S2 0.022 0.025 0.026 8 of 22 | CORRALLOU et AL. thanthephylogeneticreconstructionsincetherewerethreewell- differentiatedhaplogroups(V1,V2,andV3).Thefirsthaplogroup (V1)clusteredpopulationsfromtheÁguedasub-basinintheDuero BasinandintheTagusBasin,thosefromthewesternpartofAlagón sub-basin, plus one individual belonging to Jerte River (eastern Alagónsub-basin)(Figure4).Thesepopulationscorrespondedtothe FIGURE 3 UnrootedMLphylogenetic treebasedonthe4538SNPloci. Circlesandsquaresindicateindividuals morphologicallyidentifiedasC. vettonica andC. paludica,respectively.Thenumbers refertothoseassignedtothepopulations inTable1.Thebootstrapvalues(b)are representedwithalineperpendicular tothebranchfor75< b >90,andtwo linesforb>90.Thearrowshighlight thepopulationsofC. vettonicafromthe Caparro(→)andupperErjas(←)rivers. Abbreviation:R.forriver.Also,themapis shownindicatingthesamplinglocations referencedinTable1 FIGURE 4 HaplotypenetworkforthemitochondrialmarkerMT-CYB.Thelocalitiesareindicatedbydifferentcolorsandnumbers(see Table1).Mutationalstepsarerepresentedasfollows:oneshortlinefortwosteps,twoshortlinesforthreesteps,oracirclewiththe numberofstepsindicatedforfourormoresteps.V4*refersonlytothoseindividualsidentifiedasC. vettonicaandfoundintheS2lineage (individualsindicatedbypinkarenotincludedinV4*).Inthelegend,populationsidentifiedasC. vettonicaarerepresentedbycircles,and thoseidentifiedasC. paludica,bysquares.Populationsthatcontainindividualswithcontrastingmorphologicalandgeneticidentifications areunderlined.Abbreviations:R.forriver;S.forsub-basinandB.forbasin | 9 of 22 CORRALLOU et AL. nondifferentiatedlineageobtainedinthephylogeneticreconstruction(i.e.,V1inFigure4).Thesecondhaplogroup(V2)wascomposed ofalltheeasternAlagónsub-basinpopulationsintheTagusBasin, correspondingtotheV2lineageobtainedinthephylogeneticreconstruction.Likewise,thethirdhaplogroup(V3),correspondingto theV3lineageinthephylogeneticreconstruction,includedpopulationsfromthewesternlimitofthedistributionrangeofC. vettonica intheTagusBasin,namelythosefromthePonsul,Aravil,andErjas sub-basins,plusoneindividualfromSalorRiver.Theindividualfrom SalorRiver(identifiedasC. paludica)presentedthemostcommon haplotypeofthehaplogroup.HaplogroupsV1andV2wereseparatedbythree/fourmutationalsteps,andV3wasseparatedfromV1 by10mutationalsteps.Inthegroupcomprisinghaplotypesfound mainly in C. paludica (S2), a subset was composed of individuals identifiedasC. vettonicafromtheAlagónandErjassub-basins(V4; referringtoonlythoseindividualsidentifiedmorphologicallyasC. vettonicabutthataregeneticallyclosertoC. paludicaaccordingto themitochondrialdata)alongwithotherindividualsofC. paludica fromotherTagusandDuerotributaries.However,noevidenceofa geographicstructurewasfoundwithinS2.Thehaplotypenetwork forthenuclear marker(RAG1) showeda singlegroup inwhichall studiedpopulationsofbothspecieswererepresentedandtherefore wasoflowresolution(Figure5). AccordingtotheAMOVAswiththeMT-CYBdata,thehighest percentageofvarianceamonggroupswaswhenthegroupingswere comprisedofthe lineagesfoundinthepreviousanalysis (V1,V2, V3,andV4;Table3).Thepercentageofvarianceexplainedbetween basins(DueroandTagus)wasnotsignificant,withmostofthevariation explained by differences between and within populations. Populationdifferentiationbetweenlocalitieswithineachlineageof C. vettonicaforMT-CYBwaslow(Φst <0.34),andnosignificantdifferenceswerefound,exceptfortheAceboRiver,whichshoweda highlevelofdifferentiationwiththeÁrragoriver(Φst >0.74;Table S4).HighvaluesofΦstweredetectedbetweenpopulationsfromdifferentlineages(Φst >0.75). TheresultsoftheSNPstructureanalysesshowedthatthemost probablenumberofsubpopulations(geneticgroups)forC. vettonicamatrixwasK =4andK =2forPuechmailleandΔKmethods, respectively(TableS5).ForK =4,onlythreepopulationsshowed verylittleadmixture:MayasRiverintheDueroBasin(firstgroup,in yellow)andCaparroandErjasriversintheTagusBasin(secondand thirdgroups,inblueandpurple,respectively)(Figure6).Theother populations(Alagón,Cdh,Jerte,andFrancia)showeddifferentdegreesofadmixtureofthesecond(blue)andfourth(orange)groups. GataRiverhadasimilarlevelofadmixtureofthefirst(yellow)and fourth(orange)groupswithasmallcontributionfromthethird(purple).TheresultsforK =2wereconsistentwiththoseobtainedfor K =4.MayasandGatabelongedtothefirstgroup(orange),Caparro tothesecondgroup(blue),andtherestofthepopulationspresented amixtureofbothgroups.Whenweanalyzedthepopulationsofboth species,themostprobablenumberofsubpopulationswasK =6and K =2forPuechmailleandΔKmethods,respectively(TableS5).For K =6,MayasandGatacomprisedthefirstgroup(orange),withGata havingasmallcontributionfromthethirdgroup(purple)inallindividuals(Figure6).Alagón,Cdh,Jerte,Francia,andCaparroshowed differentdegreesofadmixtureofthefirst(orange)andsecond(blue) groups,withCdhandFranciaalsohavingasmallproportionofthe thirdgroup(purple)andJertehavingasmallproportionofthefourth group(pink).Thesecondgroup(blue)waspresentinallpopulations ofC. paludica.TheupperErjasRiverwasassignedtoasinglegenetic group(Group3;purple),whichalsocomprisedapproximately29% and15%ofthetotalcompositionoftheAurelaandSalorrivers,respectively,withverysmallcontributionsfromthesecond(blue)and fourth(pink)groupsinsomeindividuals.TheresultsforK = 2 were consistentwiththoseobtainedforK =6.MayasandGatabelonged tothefirstgroup(orange),Ibor,Tamuja,Huebra,andYeltestothe secondgroup(blue),andtherestofthepopulationspresentedamixtureofbothgroups. 3.3 | Genetic diversity and demography Overallgeneticdiversityparametersforallindividualsidentifiedas C. vettonicabasedonMT-CYBwerehigh(Table4).ExcludingpopulationsofC. vettonicainwhichonlyoneindividualwasstudied,Hd rangedfrom0to1,andhrangedfrom1to7.Thepopulationswith thelowestvaluesofgeneticdiversitywereMayas,Acebo,Arades, andCaparrowithonlyonehaplotypeeach(Table4).Themostvariablepopulationsforallofthediversityparameterswerethelower andmiddleErjas,Alagón,andGata.Thesub-basinswiththehighest FIGURE 5 HaplotypenetworkforthenuclearmarkerRAG1.Thecolorscorrespondtothedifferentlocalities,exceptpink,whichisused foralllocalitiesofC. paludica.Mutationalstepsarerepresentedasfollows:nolineforonestepandoneshortlinefortwosteps.Inthe legend,thepopulationsidentifiedasC. vettonicaarerepresentedbycirclesandthoseidentifiedasC. paludicaarerepresentedbysquares. NumberscorrespondtothoseusedinTable1 16 of 22 | CORRALLOU et AL. smallnumberofindividualsfromtheriversofthewesternAlagón sub-basin(assuggestedinourancestralareareconstructions). 4.1.2 | (V2)EasternAlagónsub-basin Thisgroup,composedofthepopulationsfromtheeasternAlagón sub-basin,wasnotmonophyleticbasedonMT-CYB,andintheSNP analysis,presentedahighlevelofadmixturewithC. paludica. These resultssuggestextensivegeneticintrogressionofC. vettonica with its sister species C. paludica,exceptinHurdanoandLadrillarrivers, whereonlyoneindividualfromeachwasstudied.Therefore,regardlessofthedifferentlevelsofgeneticintrogression,weproposethat thesepopulationsallbeassignedtotheOCUthatwasrecentlyestablishedbyDoadrioetal.(2021)(OCUIII;Figure10). 4.1.3 | (V3)Ponsul,Aravil,andErjassub-basins According to the MT-CYB analyses, the Ponsul, Aravil, and Erjas sub-basin populations constituted a well-differentiated independentlineageofC. vettonica,supportingevidencefromotherstudies thatsuggesttheuniquenessoftheErjaspopulations(Doadrioetal., 2011,2021).TheupperErjaswasalsoanindependentgroup,accordingtotheSNPanalyses,indicatingitspresent-dayisolationwith respecttotherestofthepopulationsintheAlagónsub-basinand DueroBasin.Forthesereasons,weincludetheupperErjas,Ponsul, andAravilsub-basinsintheOCUpreviouslyestablishedforonlythe Erjassub-basin(Doadrioetal.,2021),therebyexpandingthisOCU’s areaofcoverage(OCUIV;Figure10). 4.2 | Hybridization Hybridizationduetointrogressionasaresultofbothanthropicand naturalfactorshasbeendetectedforotherMediterraneancyprinid and cobitid species (Almodóvar et al., 2012; Bohlen & Ráb, 2001; Choleva et al., 2014; Corral-Lou et al., 2019; Cunha et al., 2004; Denysetal.,2013;Jankoetal.,2005;Pereaetal.,2016;Tancioni etal.,2013).WithintheIberianPeninsula,themajorityofthenaturalhybridpopulationsoriginatedduringtheQuaternaryperiodasa consequenceofhydrographicalrearrangementsand/orPleistocene glaciation–deglaciationcycles(Almodóvaretal.,2012;Gante,2009). Consistentwiththis,wepresentstrongevidenceconfirmingtheintrogressionpreviouslyreportedintheeasternAlagónpopulationsof C. vettonicabyDoadrioetal.(2011),Doadrioetal.(2021),andalso postulatethegeneticintrogressionofpopulationsinsympatricareas (Ponsul,Aravil,andlowErjas;seebelow).Supportingthishypothesis isthelocationofthesepopulationsoutsidethepotentialrangeofthe speciesbasedontheecologicalnichemodels.Theintrogressionof evenafewlociduetohybridizationcanleadtoadaptivedivergenceif thelociarebetteradaptedtoothernichesandpersistinthegenome (Abbottetal.,2013). 4.2.1 | (V1)WesternAlagónsub-basinand DueroBasin TheGatapopulationwithinthewesternAlagónsub-basinwasthe onlyonethatpresentedamitochondrialhaplotypeofC. paludica. However,intheSNPanalysis,noneoftheindividualsshowedadmixturewithC. paludica,andtheyalsoresolvedasaphylogenetically isolatedgroup.Therefore,atpresent,thereisnoevidencetosupportgeneticintrogressionwithC. paludicainthisarea. 4.2.2 | (V2)EasternAlagónsub-basin WefoundstrongevidenceofgeneticintrogressionbetweenC. paludicaandC. vettonicaintheeasternAlagónsub-basin,inlinewith previous studies (Doadrio et al., 2011, 2021). These populations notonlyshowedahighlevelofadmixturewithC. paludicabutalso highvaluesofgeneticdiversityintheSNPanalyses.Indeed,they hadthehighestvaluesofallthepopulations,aswouldbeexpected forhybridspopulations.Thehighdiversityvaluesobtainedforthe SNPdatawerenotcongruentwiththelowonesobtainedforMT- CYB.However,thisobservationisnotsurprisingashybridstypically cannotbedetectedbymitochondrialmarkersduetotheirmaternal inheritance. ThecaseofthepopulationfromCaparroRiverisparticularlyinteresting.ItwastheonlyeasternAlagónsub-basinpopulationthat showed a closer phylogenetic relationship with populations of C. paludicathanwiththoseofC. vettonicaandahigherdegreeofgeneticintrogression.Thisfactcouldbeexplainedbytwohypotheses. ThefirstisthatthesourcepopulationwasC. paludicaand,through afoundereffect,afewindividualsofC. vettonicaestablishedthemselvesinthislocality,givingrisetoapopulationofC. paludicaintrogressedbyC. vettonica.Thesecondhypothesisisthatthesource populationwasC. vettonica,butitwasdisplacedbyC. paludicaafter thisspecies’arrivalduetoitsbetteradaptedgenome(Abbottetal., 2013). 4.2.3 | (V3)Ponsul,Aravil,andErjassub-basins Accordingtothemitochondrialanalyses,thePonsul,Aravil,andErjas sub-basinpopulationsofC. vettonicawererelatedtoapopulation ofC. paludicafromSalorRiverwhosemostcommonhaplotypewas alsosharedbyallofthesepopulationsofC. vettonica.Inaddition,the resultsoftheSNPanalysisrevealedastrongrelationshipbetween theErjaspopulation(upperErjas)withthoseintheSalorandAurela sub-basins,particularlythelatterwithwhichitwasmostcloselyrelatedphylogenetically.Therelationshipbetweenthesepopulations wasmoreevidentwhentheAurelaandSalorpopulationswereremovedfromtheSNPgeneticdiversityanalysis,whichresultedinan increaseinthenumberofprivateallelesfromeightto141.However, intheSNPstructureanalysis,onlyoneindividualbelongingtothe Erjaspopulationhadaverylowcontributionofthegeneticgroup | 17 of 22 CORRALLOU et AL. assignedtoC. paludica.Altogether,theresultsdonotsupportthe presenceofgeneticintrogressionintheupperErjas,wherenomitochondrialhaplotypesofC. paludicawerefound.Moreover,higher valuesofgeneticdiversitywouldbeexpectedfortheSNPdataif therewereintrogressions,asoccursintheeasternAlagónsub-basin populations. However,wepostulategeneticintrogressioninthePonsul,Aravil, andlowerErjassub-basin,wherebothspeciesarepresent(Perdices &Coelho,2020).Thispostulationisbasedontheirlocationoutside theareaindicatedbytheecologicalnichemodelsandthefactthat introgression by hybridization is a common phenomenon that has beendetectedfornotonlyC. vettonicaandC. paludicabutalsoother speciesofEuropeanCobitis(Bohlen&Ráb,2001;Jankoetal.,2007; Ráb&Slavík,1996).Furthermore,thisgeneralpatternhasoftenbeen observedintheIberianPeninsulafortwosympatricspeciesasaconsequenceofsecondarycontactafterdiversification(Corral-Louetal., 2019;Doadrioetal.,2021;Pereaetal.,2016,2021). Withinthisgroup,hybridlocalitiesseemtoberestrictedtothe lowerreachesofrivers,whiletheupperreachesareinhabitedby onlyC. vettonica(Perdices&Coelho,2020).Thehabitatpreferences ofthetwospecies,C. paludicainthemiddleandlowerreaches,and C. vettonicaintheupperreaches(Doadrioetal.,2011),likelyexplains therestricteddistributionofthehybrids. 4.3 | Biogeography Accordingtoourdivergencetimeestimations,thetimingofthese divergencesagreeswithacommonperiodofpopulationdifferentiationformanyIberianfreshwaterfishspecies,whichinturncoincideswiththeculminationoftheformationoftheIberiandrainage networkandthePleistoceneglacialcycles(Casal-Lópezetal.,2017; Corral-Louetal.,2019,2021;Paisetal.,2012;Pereaetal.,2016). Althoughgeologicalinformationforthestudyareaisscarce,several tectonicandclimaticeventsareknowntohaveoccurredduringthe Quaternary near the study area, which influenced the hydrogeomorphologyoftheaquaticnetworksuchasthefluvialcapturephenomena,thepresenceofpaleoglaciersduringthelastglacialcycle, andchangesinthetrajectoryofsomerivers(Benitoetal.,2003; Carrascoetal.,2013,2015;Goyetal.,2020).Therefore,wepropose severalhypothesesbasedonthebiogeographyofthespecies,which mayrevealthehydrogeomorphologicalevolutionofthestudyarea withintheIberianPeninsula. TheearliestdivergenceinC. vettonicawasthePonsul-Aravil-Erjas lineage(V3),whichsplitfromthetwoAlagónsub-basinlineages(V1- V2)duringthePleistocene(Calabrianperiod,~250,000Mya).Based ontheancestralareareconstructions,thediversificationofthese lineagesoccurredintheareacurrentlycomprisingboththeErjasand Alagónsub-basins.ThisimpliesthattheErjasacquireditscurrent configurationduringthisperiodoftime,interruptinganygeneflow thatmayhavebeenestablishedbyconnectionsbetweentributaries ontherightbankofErjasRiverandtheleftbankofÁrragoRiver. Subsequently,followingtheisolationanddiversificationoftheErjas sub-basinpopulationfromtherestofthepopulations,individuals from this sub-basin dispersed to downstream areas of the Tagus Basin(reachingPonsul,AravilandSalorsub-basins).Thesedispersal eventscouldhavetakenplaceasaconsequenceofclimaticchanges associatedwiththeLGM,assupportedbyournichemodelswhich indicated the Ocreza, Nisa, Sever, Ponsul, Aravil, and Salor sub- basinsasmoresuitablehabitatsthantherestoftheriversinthe currentdistributionrangeofthespecies.Theseeventsgaveriseto boththecolonizationofnewareastothewest(PonsulandAravil sub-basins)andthesympatrywithC. paludica,whichledtointrogressionevents. Thedivergenceofthetwootherlineagescomprisingthewestern andtheeasternAlagónsub-basinpopulations(V1andV2)alsooccurredduringthePleistocene(Chibanianperiod,~80,000Mya).The complexgeologyoftheAlagónsub-basinisevidentbytheanomalouscurrenttrajectoryofsomeofitsrivers(Carrasco&Pedraza, 1991;DíezHerrero,2003;Goyetal.,2020;Jiménez,1994;Schnabel &Gutiérrez,2014).Forinstance,thetwobanksoftheAlagónsub- basin,representedbythetwomitochondriallineages,arecurrently joinedbythemouthoftheÁrragoRiverintheAlagónRiver,that is,themainriverofthewesternandtheeasternAlagónsub-basin, respectively.Theunionbetweenthetwobanksoccurredrelatively recentlyasaconsequenceoftheupriveractionoftheAlagónRiver, whichcapturedtheCaparroRiverthateventuallyflowedintothe Tagus through the Fresnedosa riverbank (Schnabel & Gutiérrez, 2014).Giventhiscontext,weproposetwoscenariosforthedivergenceofV1andV2.Inthefirst,thetwolineagesdivergedpriorto theunionofthetwobanksandlatergeneticflowbetweenpopulationswasnotpossibleduetothepoorhabitatconditionsforC. vettonicainthelowerAlagónsub-basin.Inthesecond,theydiverged aftertheunionofthebanksduetoclimaticchangesduringtheLGM. According to our niche models, the optimal ecological conditions forC. vettonicawerelocatedfurtherdownstreamand,aspredicted byphylogeographytheory,thelatitudinaldistributionoforganisms retractedsouthwards(Ehrichetal.,2007;Hewitt,1996;Rodríguez etal.,2011).Onceoptimalconditionsreturnedtotheriverheadwaters,thespecieswentbackupstreamtowardboththeeastandwest marginsoftheAlagónsub-basin,andsubsequentlydiversified. AfterthedivergenceofthetwobanksoftheAlagón,thepopulations in the Duero and western Alagón sub-basin diverged. SeveralstudieshavereportedonthegeneticstructureoffreshwaterfaunaassociatedwithconnectionsbetweentheTagusBasin and adjacent basins (e.g., Alagón, Alberche and Lozoya rivers; Carmonaetal.,2000;Casas-Sainz&DeVicente,2009;Doadrio, 1988;Pérez-González,1980;Sousa-Santosetal.,2007).Although recent connections may be explained as a consequence of the tectonic activity in the region (Goy et al., 2020), there are no geologicalstudiesthatsupportarecentconnectionbetweenthe westernAlagónsub-basin(TagusBasin)andtheÁguedasub-basin (DueroBasin),despitethecloseproximityofsomeoftheirriver headwaters(insomecases,aslittleas150moflineardistance). 18 of 22 | CORRALLOU et AL. A connection, however, has been hypothesized to explain the co-occurrence of individuals with the mitochondrial genome of Squalius carolitertii(distributedintheDueroBasin)andthoseof S. pyrenaicus(distributedintheTagusBasin)intheÁrragoRiver (Pereaetal.,2021).ThecloseproximityofMayas(DueroBasin) andÁrrago(TagusBasin)rivers,whosecoursesareonlyseparated by~150m,isparticularlynotable.PiracyeventsthatoccurredbetweenthesetworiversinthePleistocenemayexplaintheclose mitochondrialrelationshipofthepopulationsinhabitingthem,as suggestedbythedispersaleventestimatedinourancestralarea reconstructionanalysis.LaterisolationoftheDueroandTagushydrologicalbasinsthenledtothedivergenceofthesepopulations, ascorroboratedbytheSNPanalysis. 4.4 | Conservation ThehabitatsofC. vettonicahavebeenandcontinuetobethreatenedbythemaincausesofbiodiversityloss(e.g.,overexploitation, waterpollution,flowmodification,habitatdestructionanddegradation,andtheintroductionofinvasivespecies)(Doadrioetal.,2011; Dudgeon,2019;Sousa-Santosetal.,2014).Giventhesethreats,the populationsinOCUI(Águedasub-basin)areparticularlyvulnerable astheyhadthelowestgeneticdiversityvaluesforbothMT-CYBand theSNPs,andonlyoneMT-CYBhaplotypesharedbyallthepopulations.Thislowlevelofgeneticdiversitymakesthemmoresensitivetoextrinsicchangesandthereforeatgreaterriskofextinction (Frankhametal.,2002).Inaddition,populationsinOCUIareshrinkingasaresultoftheformationofdamsalongtheÁguedasub-basin, andinsomelocalities(i.e.,Turones),theyhaveevendisappeared. Thespecies’distributioninthisareaisnowrestrictedtoafewtributariesofMayasRiver.AlthoughthegeneticdiversityvaluesofpopulationsinOCUIIwerehigherthanthoseofOCUI,thepossibilityof abottleneckfollowedbyapopulationexpansioncannotberejected forthesepopulations,despiteadecreaseinthenumberofindividualsobservedinrecentyears(Doadrioetal.,2011).Exoticspecies, suchasLepomis gibbosusandMicropterus salmoides,havecolonized theupperpartsoftheriverscoveredbyOCUII,whichmaybeone ofthemaincausesforthedeclineinthenumberofC. vettonicain theserivers.Otherpotentialcausesofthisdeclineremainunknown. For OCU IV, there was no genetic distinction between the non-introgressed(upperErjas)andputativelyintrogressed(Ponsul, Aravil,andlowerErjas)populations.Despitethis,anymanagement planmusttakeintoconsiderationthispotentialintrogressionasthe indiscriminate mixing of individuals from these two groups could leadtogeneticintrogressionwithC. paludicathroughouttheentire areaofOCUIV.Cobitis paludica,whichismoreofageneralistthan C. vettonica,iswidelydistributedthroughouttheIberianPeninsula, occupyingagreatvarietyofecologicalniches(Doadrioetal.,2011). Thus, genetic introgression by C. paludica in sympatric localities would likelyprove disadvantageous forC. vettonicaingeneral,as twoofthefourOCUs(IIIandIV)wouldbeaffectedbythisphenomenon.ItcouldleadtotheextinctionofC. vettonicaintheseareasor toadaptivevariations,resultinginamajorlossofthegeneticdiversityofthisspecies. ACKNOWLEDGMENTS WethankP.Garzón,I.DoadrioJr,J.L.González,andG.Gonzálezfor collectingalmosteverysamplefrom2009to2010undertheproject “BasisforthemonitoringofSpanishFreshwaterFishes”,supportedby theMinisteriodeMedioAmbienteyMedioRuralyMarino(MARM), andC.Marcos,J.CVelaco,andG.Gonzálezfortheirhelpintheproject Life13nat/es/000772.WegreatlythankL.Alcarazforlaboratoryassistance.WearegratefulforthevaluablesuggestionsandEnglishediting performedbyM.Modrell.ThisstudywasalsosupportedbyamanagementcommissionfundedbytheDueroHydrographicConfederationin Spain(Análisisyestudiosfilogeográficosydevariabilidadgenéticade endemismosibéricosdeciprínidos:Life13nat/es/000772“Actuaciones paralaprotecciónyconservacióndeciprínidosibéricosdeinteréscomunitario”)andbygrantIND2017/AMB-7699toAC-L,fundedbythe Community of Madrid (Spain) through its Industrial Doctorate program. We acknowledge the Santander Supercomputacion support groupattheUniversityofCantabriaforprovidingaccesstothesupercomputerAltamiraattheInstituteofPhysicsofCantabria(IFCA-CSIC), amemberoftheSpanishSupercomputingNetwork,whichwasusedto performsimulations/analyses. CONFLICT OF INTEREST Theauthorsdeclarenoconflictofinterest. AUTHOR CONTRIBUTIONS Andrea CorralLou: Conceptualization (lead);Datacuration (lead); Formal analysis (lead); Methodology (equal); Resources (lead); Software (lead); Validation (lead); Visualization (lead); Writing – originaldraft(lead);Writing–review&editing(lead).Silvia Perea: Conceptualization (equal); Data curation (equal); Formal analysis (supporting); Methodology (supporting); Supervision (supporting); Writing–review&editing(equal).Anabel Perdices:Resources(supporting);Supervision(supporting);Writing–review&editing(supporting). Ignacio Doadrio: Conceptualization (lead); Data curation (equal);Fundingacquisition(lead);Investigation(lead);Methodology (lead); Resources (equal); Supervision (lead); Validation (lead); Visualization(equal);Writing–review&editing(equal). DATA AVAILABILITY STATEMENT The new sequences of the mitochondrial (MT-CYTB; OM234794- OM235001) and nuclear (RAG1; OM235002-OM235091) markers obtainedfromthisstudyareavailableinGenBank.TheSNPdatahave beendepositedinVCFformatinFigshare(https://doi.org/10.6084/ m9.figshare.18778148.v1). ORCID Andrea CorralLou https://orcid.org/0000-0003-0895-1305 Silvia Perea https://orcid.org/0000-0003-0436-8577 Anabel Perdices https://orcid.org/0000-0003-0980-215X Ignacio Doadrio https://orcid.org/0000-0003-4863-9711 | 19 of 22 CORRALLOU et AL. REFERENCES Abbott,R.,Albach,D.,Ansell,S.,Arntzen,J.W.,Baird,S.J.E.,Bierne,N., Boughman,J.,Brelsford,A.,Buerkle,C.A.,Buggs,R.,Butlin,R.K., Dieckmann,U.,Eroukhmanoff,F.,Grill,A.,Cahan,S.H.,Hermansen, J. S., Hewitt, G., Hudson, A. G., Jiggins, C., … Zinner, D. (2013). Hybridizationandspeciation.Journal of Evolutionary Biology,26(2), 229–246.https://doi.org/10.1111/j.1420-9101.2012.02599.x Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control,19(6),716–723.https://doi. org/10.1109/TAC.1974.1100705 Almodóvar,A.,Nicola,G.G.,Leal,S.,Torralva,M.,&Elvira,B.(2012). NaturalhybridizationwithinvasivebleakAlburnus alburnusthreatens thesurvivalofIberianendemiccalandinoSqualiusalburnoides complex and Southern Iberian chub Squalius pyrenaicus. Biological Invasions, 14(11), 2237–2242. https://doi.org/10.1109/ TAC.1974.1100705 Alonso, J. L., Pulgar, F. J. Á., & Pedreira, D. (2007). The relief of the CantabrianMountains.Teaching Earth Sciences [In Spanish],15(2), 151–163.https://doi.org/10.1109/TAC.1974.1100705 Bandelt,H.J.,Forster,P.,&Röhl,A.(1999).Median-joiningnetworksfor inferringintraspecificphylogenies.Molecular Biology and Evolution, 16(1), 37–48. https://doi.org/10.1093/oxfordjournals.molbev. a026036 Beaulieu, J. M., Tank, D. C., & Donoghue, M. J. (2013). A Southern Hemisphere origin for campanulid angiosperms, with traces of thebreak-upofGondwana.BMC Evolutionary Biology,13(1),1–17. h t t p s : / / d o i . o r g / 1 0 . 1 1 8 6 / 1 4 7 1 - 2 1 4 8 - 1 3 - 8 0 Benito, G., Sopeña, A., Sánchez-Moya, Y., Machado, M. J., & Pérez- González,A.(2003).PalaeofloodrecordoftheTagusRiver(centralSpain)duringtheLatePleistoceneandHolocene.Quaternary Science Reviews, 22(15–17), 1737–1756. https://doi.org/10.1016/ S 0 2 7 7 - 3 7 9 1 ( 0 3 ) 0 0 1 3 3 - 1 Bohlen, J., & Ráb, P. (2001). Species and hybrid richness in spined loachesofthegenusCobitis(Teleostei:Cobitidae),withachecklistofEuropeanformsandsuggestionsforconservation.Journal of Fish Biology,59,75–89.https://doi.org/10.1111/j.1095-8649.2001. t b 0 1 3 8 0 . x Bolger, A. M.,Lohse,M.,& Usadel,B. (2014).Trimmomatic: aflexible trimmerforIlluminasequencedata.Bioinformatics,30(15),2114– 2120.https://doi.org/10.1093/bioinformatics/btu170 Buj, I., Marčić, Z., Ćaleta, M., Šanda, R., Geiger, M. F., Freyhof, J., Machordom, A., & Vukić, J. (2017). Ancient connections among theEuropeanriversandwatershedsrevealedfromtheevolutionaryhistoryofthegenusTelestes(Actinopterygii;Cypriniformes). PLoS One, 12(12), e0187366. https://doi.org/10.1371/journ al.pone.0187366 Carmona, J. A., Doadrio, I., Márquez, A. L., Real, R., Hugueny, B., & Vargas,J.M.(1999).DistributionpatternsofindigenousfreshwaterfishesintheTagusRiverbasin,Spain.Environmental Biology of Fishes,54(4),371–387.https://doi.org/10.1023/A:1007535200837 Carmona, J. A., DomÍnguez, J., & Doadrio, I. (2000). Congruence betweenallozymeandcytochromebgenesequencedatainassessing geneticdifferentiation withintheIberianendemicChondrostoma lemmingii(Pisces:Cyprinidae).Heredity,84(6),721–732.https://doi. org/10.1046/j.1365-2540.2000.00720.x Carrasco, R. M., & Pedraza, J. (1991). In the Jerte valley [In Spanish]. Cuaternario Y Geomorfología, 5, 15–25. https://doi. org/10.1046/j.1365-2540.2000.00720.x Carrasco,R.M.,Pedraza,J.,Domínguez-Villar,D.,Villa,J.,&Willenbring, J. K. (2013). The plateau glacier in the Sierra de Béjar (Iberian Central System) during its maximum extent. Reconstruction and chronology.Geomorphology,196, 83–93.https://doi.org/10.1109/ TAC.1974.1100705 Carrasco, R. M., Pedraza, J., Dominguez-Villar, D., Willenbring, J. K., & Villa, J. (2015). Sequence and chronology of the Cuerpo de Hombrepaleoglacier(IberianCentralSystem)duringthelastglacial cycle. Quaternary Science Reviews, 129, 163–177. https://doi. org/10.1016/j.quascirev.2015.09.021 Casal-López,M.,&Doadrio,I.(2018).TheMessinianimprintonthe evolutionoffreshwaterfishesofthegenusLuciobarbusHeckel, 1843 (Teleostei, Cyprinidae) in the western Mediterranean. Journal of Biogeography, 45(7), 1593–1603. https://doi. org/10.1111/jbi.13244 Casal-López,M.,Perea,S.,Sousa-Santos,C.,Robalo,J.I.,Torralva,M., Oliva-Paterna, F. J. & Doadrio, I. (2017). Paleobiogeography of anIberianendemicspecies,Luciobarbussclateri(Günther,1868) (Actinopterygii,Cyprinidae),inferredfrommitochondrialandnuclearmarkers.JournalofZoologicalSystematicsandEvolutionary Research,56(2),127-147.https://doi.org/10.1111/jzs.12192 Casas-Sainz,A.M.,&DeVicente,G.(2009).Onthetectonicoriginof Iberiantopography.Tectonophysics,474(1–2),214–235.https://doi. org/10.1016/j.tecto.2009.01.030 Catchen,J.M.,Hohenlohe,P.A.,Bernatchez,L.,Funk,W.C.,Andrews,K. R.,&Allendorf,F.W.(2013).Stacks:Ananalysistoolsetforpopulationgenomics.Molecular Ecology,22(11),3124–3140.https://doi. org/10.1111/mec.12354 Choleva,L.,Musilova,Z.,Kohoutova-Sediva,A.,Paces,J.,Rab,P.,&Janko, K.(2014).Distinguishingbetweenincompletelineagesortingand genomicintrogressions:completefixationofallospecificmitochondrialDNAinasexuallyreproducingfish(Cobitis;Teleostei),despite clonalreproductionofhybrids.PLoS One,9(6),e80641.https://doi. org/10.1371/journal.pone.0080641 Collares-Pereira,M.J.,Alves,M.J.,Ribeiro,F.,Domingos,I.,Almeida,P. R.,DaCosta,L.,Gante,H.,Filipe,A.F.,Aboim,A.,Rodrigues,P.M., &Magalhães,M.F.(2021).Guide to freshwater fish and migrants from mainland Portugal [In Portuguese](p.292).EdiçõesAfrontamento. Corral-Lou,A.,Perea,S.,Aparicio,E.,&Doadrio,I.(2019).Phylogeography andspeciesdelineationofthegenusPhoxinusRafinesque,1820 (Actinopterygii: Leuciscidae) in the Iberian Peninsula. Journal of Zoological Systematics and Evolutionary Research, 57(4), 926–941. https://doi.org/10.1111/jzs.12320 Corral-Lou,A.,Perea,S.,&Doadrio,I.(2021).HighgeneticdifferentiationintheendemicandendangeredfreshwaterfishAchondrostoma salmantinumDoadrioandElvira,2007fromSpain,asrevealedby mitochondrialandSNPmarkers.Conservation Genetics,22(4),585– 6 0 0 .  h t t p s : / / d o i . o r g / 1 0 . 1 0 0 7 / s 1 0 5 9 2 - 0 2 1 - 0 1 3 8 1 - y Cunha, C., Coelho, M. M., Carmona, J. A., & Doadrio, I. (2004). PhylogeographicalinsightsintotheoriginsoftheSqualius alburnoides complexviamultiplehybridizationevents.Molecular Ecology,13(9), 2807–2817.https://doi.org/10.1111/j.1365-294X.2004.02283.x Denys,G.P.J.,Dettai,A.,Persat,H.,Doadrio,I.,Cruaud,C.,&Keith,P. (2013).StatusoftheCatalanchubSqualius laietanus(Actinopterygii, Cyprinidae) in France: Input from morphological and molecular data.Knowledge and Management of Aquatic Ecosystems,408,04. https://doi.org/10.1051/kmae/2013039 Díez Herrero, A. (2003). Geomorfología e Hidrología fluvial del río Alberche. Modelos y SIG para la gestión de riberas. Serie Tesis Doctorales2. Publicacionesdel Instituto Geológico yMinerode España,Madrid.587 pp. Doadrio,I.(1988).DelimitationofareasintheIberianPeninsulaonthe basisoffreshwaterfishes.Bonner Zoologische Beiträge,39,  1 1 3 – 1 2 8 . Doadrio,I.,Carmona,J.A.,&Machordom,A.(2002).Haplotypediversity andphylogeneticrelationshipsamongtheIberianbarbels(Barbus, Cyprinidae) revealtwoevolutionarylineages.Journal of Heredity, 93(2),140–147.https://doi.org/10.1093/jhered/93.2.140 Doadrio,I.,Perdices,A.,&Machordom,A.(1996).Allozymicvariationof theendangeredkillifishAphanius iberusanditsapplicationtoconservation.Environmental Biology of Fishes,45(3),259–271. Doadrio,I.,&Perdices,A.(2005).Phylogeneticrelationshipsamongthe Ibero-Africancobitids(Cobitis,Cobitidae)basedoncytochromeb 20 of 22 | CORRALLOU et AL. sequencedata.Molecular Phylogenetics and Evolution,37(2),484– 493.https://doi.org/10.1016/j.ympev.2005.07.009 Doadrio,I.,Perea,S.,&Corral-Lou, A. (2021).Genetic characterization and identification of Conservation Units in endemic fish of community interest[InSpanish].Retrievedfromhttp://hdl.handle.net/10261/ 241564 Doadrio,I.,Perea,S.,Paloma,G.H.,&González,J.L.(2011).Ictiofauna continental española: bases para su seguimiento, Ministerio de Medio Ambiente y Medio Rural y Marino. Centro de Publicaciones.DGMedio NaturalyPolíticaForestal,MARM. Drummond,A.J.,Suchard,M.A.,Xie,D.,&Rambaut,A.(2012).Bayesian phylogeneticswithBEAUtiandtheBEAST1.7.Molecular Biology and Evolution,29(8),1969–1973.https://doi.org/10.1093/molbev/ mss075 Dudgeon,D.(2019).Multiplethreatsimperilfreshwaterbiodiversityin theAnthropocene.Current Biology,29(19),R960–R967.https://doi. org/10.1016/j.cub.2019.08.002 Ehrich, D., Gaudeul, M., Assefa, A., Koch, M. A., Mummenhoff, K., Nemomissa, S., Consortium, I., & Brochmann, C. (2007). Genetic consequences of Pleistocene range shifts: Contrast between the Arctic, the Alps and the East African mountains. Molecular Ecology, 16(12), 2542–2559. https://doi. org/10.1111/j.1365-294X.2007.03299.x Elith,J.,H.Graham,C.,P.Anderson,R.,Dudík,M.,Ferrier,S.,Guisan,A., J.Hijmans,R.,Huettmann,F.,R.Leathwick,J.,Lehmann,A.,Li,J.,G. Lohmann,L.,A.Loiselle,B.,Manion,G.,Moritz,C.,Nakamura,M., Nakazawa,Y.,McC.M.Overton,J.,TownsendPeterson,A.,…E. Zimmermann,N.(2006).Novelmethodsimprovepredictionofspecies’distributionsfromoccurrencedata.Ecography,29(2),129–151. https://doi.org/10.1111/j.2006.0906-7590.04596.x Evanno,G.,Regnaut,S.,&Goudet,J.(2005).Detectingthenumberof clusters of individuals using the software STRUCTURE: a simulation study. Molecular Ecology, 14(8), 2611–2620. https://doi. org/10.1111/j.1365-294X.2005.02553.x Excoffier,L.,&Lischer,H.E.(2010).Arlequinsuitever3.5:Anewseries ofprogramstoperformpopulationgeneticsanalysesunderLinux andWindows.Molecular Ecology Resources,10(3),564–567.https:// doi.org/10.1111/j.1755-0998.2010.02847.x Fick,S.E.,&Hijmans,R.J.(2017).WorldClim2:New1-kmspatialresolutionclimatesurfacesforgloballandareas.International Journal of Climatology,37(12),4302–4315.https://doi.org/10.1002/joc.5086 Foll,M.,&Gaggiotti,O.(2008).Agenome-scanmethodtoidentifyselectedlociappropriateforbothdominantandcodominantmarkers:ABayesianperspective.Genetics,180(2),977–993.https://doi. org/10.1534/genetics.108.092221 Frankham,R.,Ballou,S.E.J.D.,Briscoe,D.A.,&Ballou,J.D.(2002). Introduction to conservation genetics.CambridgeUniversityPress. Fu,Y.X.(1997).Statisticaltestsofneutralityofmutationsagainstpopulation growth, hitchhiking and background selection. Genetics, 147(2),915–925.https://doi.org/10.1093/genetics/147.2.915 Gante,H.M.F.(2009).A role for introgressive hybridization in the evolution of freshwater fishes: The casestudy of Iberian Barbus (Teleostei, Cyprinidae).ArizonaStateUniversity. Gante, H. F., Micael, J., Oliva-Paterna, F. J., Doadrio, I., Dowling, T. E., & Alves, M. J. (2009). Diversification within glacial refugia: tempo and mode of evolution of the polytypic fish Barbus sclateri. Molecular Ecology, 18(15), 3240–3255. https://doi. org/10.1111/j.1365-294X.2009.04264.x Gibbard,P.L.,Head,M.J.,Walker,M.J.,&SubcommissiononQuaternary Stratigraphy(2010).FormalratificationoftheQuaternarySystem/ Period and the Pleistocene Series/Epoch with a base at 2.58 Ma. Journal of Quaternary Science, 25(2), 96–102. https://doi. org/10.1002/jqs.1338 Gonzalez,E.G.,Cunha,C.,Ghanavi,H.R.,Oliva-Paterna,F.J.,Torralva, M.,&Doadrio,I.(2018).Phylogeographyandpopulationgenetic analyses in the Iberian toothcarp (Aphanius iberus Valenciennes, 1846)atdifferenttimescales.Journal of Heredity,109(3),253–263. https://doi.org/10.1093/jhered/esx076 Gonzalez,E.G.,Pedraza-Lara,C.,&Doadrio,I.(2014).Geneticdiversity andpopulationhistoryoftheendangeredkillifishAphanius baeticus. Journal of Heredity,105(5),597–610.https://doi.org/10.1093/ jhered/esu034 Goy,J.L.,Cruz,R.,Martínez-Graña,A.,Valdés,V.,&Yenes,M.(2020). Geomorphological map and quaternary landscape evolution of theMonfragüePark(Cáceres,Spain).Sustainability,12(23),10099. https://doi.org/10.3390/su122310099 Hewitt,G.M.(1996).Somegeneticconsequencesoficeages,andtheir roleindivergenceandspeciation.Biological Journal of the Linnean Society,58(3),247–276.https://doi.org/10.1111/j.1095-8312.1996. t b 0 1 4 3 4 . x Hewitt,G.M.(2004).Geneticconsequencesofclimaticoscillationsin the Quaternary. Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences,359(1442),183–195.https:// doi.org/10.1098/rstb.2003.1388 Hudson,R.R.,Boos,D.D.,&Kaplan,N.L.(1992).Astatisticaltestfordetectinggeographicsubdivision.Molecular Biology and Evolution,9(1), 138–151.https://doi.org/10.1093/oxfordjournals.molbev.a040703 Janes,J.K.,Miller,J.M.,Dupuis,J.R.,Malenfant,R.M.,Gorrell,J.C., Cullingham, C. I., &Andrew,R.L. (2017).TheK=2conundrum. Molecular Ecology, 26(14), 3594–3602. https://doi.org/10.1111/ mec.14187 Janko, K.,Culling, M. A.,Rab, P., &Kotlik,P.(2005). Iceage cloning– comparisonoftheQuaternaryevolutionaryhistoriesofsexualand clonalformsofspinyloaches(Cobitis;Teleostei)usingtheanalysis ofmitochondrialDNAvariation.Molecular Ecology,14(10),2991– 3004.https://doi.org/10.1111/j.1365-294X.2005.02583.x Janko,K.,Flajšhans,M.,Choleva,L.,Bohlen,J.,Šlechtová,V.,Rábová,M., Lajbner,Z.,Šlechta,V.,Ivanova,P.,Dobrovolov,I.,Culling,M.,Persat, H.,Kotusz,J.,&Ráb,P.(2007).DiversityofEuropeanspinedloaches (genusCobitisL.):anupdateofthegeographicdistributionofthe Cobitis taeniahybridcomplexwithadescriptionofnewmolecular toolsforspeciesandhybriddetermination.Journal of Fish Biology, 71,387–408.https://doi.org/10.1111/j.1095-8649.2007.01663.x Jiménez,A.S.D.C.(1994).Themountain-lowersurfacecontactinthe southern slope of the Sierra de Gata (Spanish Central System, Cáceres)[InSpanish].Cuaternario Y Geomorfología,8(1–2),45–55. Katoh, K., & Standley, D. M. (2013). MAFFT multiple sequence alignmentsoftwareversion7:Improvementsinperformanceandusability.Molecular Biology and Evolution,30(4),772–780.https://doi. org/10.1093/molbev/mst010 Kearse,M.,Moir,R.,Wilson,A.,Stones-Havas,S.,Cheung,M.,Sturrock, S., Buxton, S., Cooper, A., Markowitz, S., Duran, C., Thierer, T., Ashton,B.,Meintjes,P.,&Drummond,A.(2012).GeneiousBasic: Anintegratedandextendabledesktopsoftwareplatformforthe organizationandanalysisofsequencedata.Bioinformatics,28(12), 1647–1649.https://doi.org/10.1093/bioinformatics/bts199 Kess,T.,Gross,J.,Harper,F.,&Boulding,E.G.(2016).Low-costddRAD methodofSNPdiscoveryandgenotypingappliedtotheperiwinkleLittorina saxatilis. Journal of Molluscan Studies,82(1),104–109. https://doi.org/10.1093/mollus/eyv042 Lanfear, R., Frandsen, P. B., Wright, A. M., Senfeld, T., & Calcott, B. (2017).PartitionFinder2:Newmethodsforselectingpartitioned modelsofevolutionformolecularandmorphologicalphylogenetic analyses.Molecular Biology and Evolution,34(3),772–773.https:// doi.org/10.1093/molbev/msw260 Leigh,J.W.,&Bryant,D.(2015).popart:full-featuresoftwareforhaplotypenetworkconstruction.Methods in Ecology and Evolution,6(9), 1110–1116.https://doi.org/10.1111/2041-210X.12410 Li,Y.L.,&Liu,J.X.(2018).StructureSelector:Aweb-basedsoftwareto selectandvisualizetheoptimalnumberofclustersusingmultiple methods. Molecular Ecology Resources,18(1),176–177.https://doi. org/10.1111/1755-0998.12719 | 21 of 22 CORRALLOU et AL. Lotterhos,K.E.,&Whitlock,M.C.(2014).Evaluationofdemographic historyandneutralparameterizationontheperformanceofFST outlier tests. Molecular Ecology, 23(9), 2178–2192. https://doi. org/10.1111/mec.12725 Machordom, A., Berrebi, P., & Doadrio, I. (1990). Spanish barbel hybridization detected using enzymatic markers: Barbus meridionalis Risso × Barbus haasi Mertens (Osteichthyes, Cyprinidae). Aquatic Living Resources, 3(4), 295–303. https://doi.org/10.1051/ alr:1990030 Marić, S., Stanković, D., Šanda, R., Ćaleta, M., Čolić, S., Šukalo, G., & Snoj, A. (2019). Genetic characterisation of European mudminnow (Umbra krameri) populations from the Sava River system. Knowledge & Management of Aquatic Ecosystems,420,46.https:// doi.org/10.1051/kmae/2019035 Marić,S.,Stanković,D.,Wanzenböck,J.,Šanda,R.,Erős,T.,Takács, P.,Specziár,A.,Sekulić,N.,Bănăduc,D.,Ćaleta,M.,Trombitsky, I.,Galambos,L.,Sipos,S.,&Snoj,A.(2017).Phylogeographyand populationgeneticsoftheEuropeanmudminnow(Umbra krameri) with a time-calibrated phylogeny for the family Umbridae. Hydrobiologia, 792(1), 151–168. https://doi.org/10.1007/s1075 0 - 0 1 6 - 3 0 5 1 - 9 McCormack,J.E.,Hird,S.M.,Zellmer,A.J.,Carstens,B.C.,&Brumfield, R.T.(2013).Applicationsofnext-generationsequencingtophylogeographyandphylogenetics.Molecular Phylogenetics and Evolution, 66(2),526–538.https://doi.org/10.1016/j.ympev.2011.12.007 Mendes,S.L.,Coelho,M.M.,&Sousa,V.C.(2019).PopulationgenomicsdatasupportsintrogressionbetweenWesternIberianSqualius freshwaterfishspeciesfromdifferentdrainages.bioRxiv,585687. Mesquita, N., Hänfling, B., Carvalho, G. R., & Coelho, M. M. (2005). Phylogeography of the cyprinid Squalius aradensis and implicationsforconservationoftheendemicfreshwaterfaunaofsouthern Portugal. Molecular Ecology, 14(7), 1939–1954. https://doi. org/10.1111/j.1365-294X.2005.02569.x Myers,G.S.(1966).DerivationofthefreshwaterfishfaunaofCentral America.Copeia,766–773.https://doi.org/10.2307/1441405 Pais, J., Cunha, P. P., Pereira, D., Legoinha, P., Dias, R., Moura, D., González-Delgado, J. A. (2012). The Paleogene and Neogene ofwesternIberia(Portugal):aCenozoicrecordintheEuropean Atlantic domain. In The Paleogene and Neogene of Western Iberia (Portugal)(pp.1–138).SpringerBriefsinEarthSciences, Springer,SeriesID:8897,158pp.https://doi.org/10.1007/978-3- 6 4 2 - 2 2 4 0 1 - 0 _ 1 Paris,J.R., Stevens, J.R., &Catchen, J. M.(2017).Lost in parameter space: A road map for stacks. Methods in Ecology and Evolution, 8(10),1360–1373.https://doi.org/10.1111/2041-210X.12775 Perdices, A., & Coelho, M. M. (2020). First records of the Vettonian spinedloachCobitis vettonicainPortugalwithupdateonitsIberian distribution.Fishes in Mediterranean Environments,2020(002),15. https://doi.org/10.29094/FiSHMED.2020.002 Perdices,A.,&Doadrio,I.(2001).ThemolecularsystematicsandbiogeographyoftheEuropeancobitidsbasedonmitochondrialDNA sequences.Molecular Phylogenetics and Evolution,19(3),468–478. https://doi.org/10.1006/mpev.2000.0900 Perea,S.,&Doadrio,I.(2015).Phylogeography,historicaldemography and habitat suitability modelling of freshwater fishes inhabiting seasonallyfluctuatingMediterraneanriversystems:acasestudy using the Iberian cyprinid Squalius valentinus. Molecular Ecology, 24(14),3706–3722. Perea,S.,Cobo-Simon,M.,&Doadrio,I.(2016).Cenozoictectonicand climatic events in southern Iberian Peninsula: Implications for theevolutionaryhistoryoffreshwaterfishofthegenusSqualius (Actinopterygii,Cyprinidae).Molecular Phylogenetics and Evolution, 97,155–169.https://doi.org/10.1016/j.ympev.2016.01.007 Perea,S., Sousa-Santos, C., Robalo, J., &Doadrio, I. (2021).Historical biogeographyoftheIberianPeninsula:Multilocusphylogenyand ancestralareareconstructionforthefreshwaterfishgenusSqualius (Actinopterygii,Leuciscidae).Journal of Zoological Systematics and Evolutionary Research, 59(4), 858–886. https://doi.org/10.1111/ jzs.12464 Pérez-González, A. (1980). Geology and stratigraphy of the Aridos deposits in the alluvial plain of Arganda (Madrid) [In Spanish]. In M. Santonja, N. López Martínez, & A. Pérez-González (Eds.), Ocupaciones achelenses en el valle del Jarama (ArgandaMadrid)(pp. 49–61).PublicacionesdelaExcelentísimaDiputaciónProvincialde Madrid. Phillips,S.J.,Anderson,R.P.,Dudík,M.,Schapire,R.E.,&Blair,M.E. (2017).Openingtheblackbox:Anopen-sourcereleaseofMaxent. Ecography,40(7),887–893.https://doi.org/10.1111/ecog.03049 Pillans, B., & Naish, T. (2004). Defining the quaternary. Quaternary Science Reviews,23(23–24),2271–2282.https://doi.org/10.1016/j. quascirev.2004.07.006 Pritchard,J.K.,Stephens,M.,&Donnelly,P.(2000).Inferenceofpopulationstructureusingmultilocusgenotypedata.Genetics,155(2), 945–959.https://doi.org/10.1093/genetics/155.2.945 Puechmaille, S. J. (2016). The program structure does not reliably recover the correct population structure when sampling is uneven: Subsampling and new estimators alleviate the problem. Molecular Ecology Resources, 16(3), 608–627. https://doi. org/10.1111/1755-0998.12512 Querejeta, M., Fernández-González, A., Romero, R., & Castresana, J. (2017). Postglacial dispersal patterns and mitochondrial genetic structure of the Pyrenean desman (Galemys pyrenaicus) in the northwesternregionoftheIberianPeninsula.Ecology and Evolution, 7(12),4486–4495.https://doi.org/10.1002/ece3.3034 Ráb, P., & Slavík, O. (1996). Diploid-triploid-tetraploid complex of the spinedloach,genusCobitisinPšovkaCreek:Thefirstevidenceof newspeciesofCobitisintheichthyofaunaoftheCzechRepublic. Acta Universitatis Carolinae. Biologica,39,201–214. Ree,R.H.,&Smith,S.A.(2008).Maximumlikelihoodinferenceofgeographicrangeevolutionbydispersal,localextinction,andcladogenesis. Systematic Biology,57(1),4–14.https://doi.org/10.1080/10635 150701883881 Rodríguez,R.,Ramirez,O.,Valdiosera,C.E.,Garcia,N.,Alda,F.,Madurell- Malapeira,J.,&Dalen,L.(2011).50,000yearsofgeneticuniformity inthecriticallyendangeredIberianlynx.Molecular Ecology,20(18), 3785–3795.https://doi.org/10.1111/j.1365-294X.2011.05231.x Ronquist, F., Teslenko, M., van der Mark, P., Ayres, D. L., Darling, A., Höhna,S.,Larget,B.,Liu,L.,Suchard,M.A.,&Huelsenbeck,J.P. (2012).MrBayes3.2:efficientBayesianphylogeneticinferenceand modelchoiceacrossalargemodelspace.Systematic Biology,61(3), 539–542.https://doi.org/10.1093/sysbio/sys029 Rozas, J., Ferrer-Mata, A., Sánchez-DelBarrio, J. C., Guirao-Rico, S., Librado, P., Ramos-Onsins, S. E., & Sánchez-Gracia, A. (2017). DnaSP6:DNAsequencepolymorphismanalysisoflargedatasets. Molecular Biology and Evolution, 34(12), 3299–3302. https://doi. org/10.1093/molbev/msx248 Schnabel,S.,&Gutiérrez,Á.G.(2014).AdvancesinGeomorphologyin Spain2012-2014.InS.Schnabel&Á.G.Gutiérrez(Eds.),Actas de la XIII Reunión Nacional de Geomorfología. Cáceres, Universidad de Extremadura– Sociedad Española de Geomorfología. Seifertová,M.,Bryja,J.,Vyskočilová,M.,Martínková,N.,&Šimková,A. (2012).MultiplePleistocenerefugiaandpost-glacialcolonizationin theEuropeanchub(Squalius cephalus)revealedbycombineduseof nuclearandmitochondrialmarkers.Journal of Biogeography,39(6), 1024–1040.https://doi.org/10.1111/j.1365-2699.2011.02661.x Sousa-Santos,C.,Collares-Pereira,M.J.,&Almada,V.(2007).Reading the history of a hybrid fish complex from its molecular record. Molecular Phylogenetics and Evolution,45(3),981–996.https://doi. org/10.1016/j.ympev.2007.05.011 Sousa-Santos, C., Jesus, T. F., Fernandes, C., Robalo, J. I., & Coelho, M.M.(2019). Fishdiversification atthepaceof geomorphological changes: Evolutionary history of western Iberian Leuciscinae 22 of 22 | CORRALLOU et AL. (Teleostei: Leuciscidae) inferred from multilocus sequence data. Molecular Phylogenetics and Evolution, 133, 263–285. https://doi. org/10.1016/j.ympev.2018.12.020 Sousa-Santos, C., Robalo, J. I., & Pereira, A. (2014). Threatened fishes of the world: Cobitis vettonica Doadrio & Perdices, 1997 (Cobitidae).Croatian Journal of Fisheries,72(4),174–175.https://doi. org/10.1109/TAC.1974.1100705 Stamatakis,A. (2016). The RAxML v8. 2. X Manual. Heidleberg Institute for Theoretical Studies. Retrieved from https://cmeh-itsorg/exeli xis/resource/download/NewManualpdf https://doi.org/10.1109/ TAC.1974.1100705 Stephens,M.,&Donnelly,P.(2003).Acomparisonofbayesianmethods forhaplotypereconstructionfrompopulationgenotypedata.The American Journal of Human Genetics,73(5),1162–1169.https://doi. org/10.1086/379378 Swets, J. A. (1988). Measuring the accuracy of diagnostic systems. Science, 240(4857), 1285–1293. https://doi.org/10.1109/ TAC.1974.1100705 Taberlet,P.,&Bouvet,J.(1994).MitochondrialDNApolymorphism,phylogeography,andconservationgeneticsofthebrownbearUrsus arctosinEurope.Proceedings of the Royal Society of London. Series B: Biological Sciences, 255(1344), 195–200. https://doi.org/10.1098/ rspb.1994.0028 Tajima,F.(1989).StatisticalmethodfortestingtheneutralmutationhypothesisbyDNApolymorphism.Genetics,123(3),585–595.https:// doi.org/10.1093/genetics/123.3.585 Tamura, K., Stecher, G., Peterson, D., Filipski, A., & Kumar, S. (2013). MEGA6: Molecular evolutionary genetics analysis version 6.0. Molecular Biology and Evolution, 30(12), 2725–2729. https://doi. org/10.1093/molbev/mst197 Tan,M.P.,Wong,L.L.,Razali, S.A.,Afiqah-Aleng,N.,MohdNor,S.A., Sung,Y.Y.,VandePeer,Y.,Sorgeloos,P.,&Danish-Daniel,M.(2019). Applicationsofnext-generationsequencingtechnologiesandcomputational tools in molecular evolution and aquatic animals conservation studies: A short review. Evolutionary Bioinformatics, 15, 1176934319892284.https://doi.org/10.1177/1176934319892284 Tancioni,L.,Russo,T.,Cataudella,S.,Milana,V.,Hett,A.K.,Corsi,E., &Rossi,A.R.(2013).TestingspeciesdelimitationsinfourItalian sympatricleuciscinefishesintheTiberRiver:Acombinedmorphologicalandmolecularapproach.PLoS One,8(4),e60392.https://doi. org/10.1371/journal.pone.0060392 Valdiosera,C.E.,García-Garitagoitia,J.L.,Garcia,N.,Doadrio,I.,Thomas, M.G.,Hänni,C.,&Götherström,A.(2008).Surprisingmigration andpopulationsizedynamicsinancientIberianbrownbears(Ursus arctos).Proceedings of the National Academy of Sciences of the United States of America, 105(13), 5123–5128. https://doi.org/10.1073/ pnas.0712223105 Weiss,S.,&Ferrand,N.(2007).Phylogeography of southern European refugia(pp.341–357).Springer. Wetjen,M.,Hübner,D.,Seehausen,O.,&Schulz,R.(2020).GeneticdiversityofendangeredChondrostomanasusintheRiverRhinesystem:Conservationgeneticsconsiderationsonstockingandreintroduction.Knowledge & Management of Aquatic Ecosystems,421,25. https://doi.org/10.1051/kmae/2020016 Yu, Y., Blair, C., & He, X. (2020). RASP 4: ancestral state reconstruction tool for multiple genes and characters. Molecular Biology and Evolution, 37(2), 604–606. https://doi.org/10.1093/molbev/ msz257 Yu,Y.,Harris,A.J.,Blair,C.,&He,X.(2015).RASP(ReconstructAncestral StateinPhylogenies):atoolforhistoricalbiogeography.Molecular Phylogenetics and Evolution, 87, 46–49. https://doi.org/10.1016/j. ympev.2015.03.008 Yu, Y., Harris, A. J., & He, X. (2010). S-DIVA (Statistical Dispersal- Vicariance Analysis):A tool forinferring biogeographic histories. Molecular Phylogenetics and Evolution,56(2),848–850.https://doi. org/10.1016/j.ympev.2010.04.011 Zarraonaindia, I., Iriondo, M., Albaina, A., Pardo, M. A., Manzano, C., Grant,W.S.,Irigoien,X.,&Estonba,A.(2012).MultipleSNPmarkersrevealfine-scalepopulationanddeepphylogeographicstructureinEuropeananchovy(Engraulis encrasicolusL.).PLoS One,7(7), e42201.https://doi.org/10.1371/journal.pone.0042201 SUPPORTING INFORMATION Additional supporting information may be found in the online versionofthearticleatthepublisher’swebsite. How to cite this article:Corral-Lou,A.,Perea,S.,Perdices,A., &Doadrio,I.(2022).Quaternarygeomorphologicaland climaticchangesassociatedwiththediversificationofIberian freshwaterfishes:ThecaseofthegenusCobitis (Cypriniformes,Cobitidae).Ecology and Evolution,12,e8635. https://doi.org/10.1002/ece3.8635