Identifying spatially concordant evolutionary significant units across multiple species through DNA

18 Pages·2017·0.71 MB·English

GlobalEcologyandConservation12(2017)170e187
ContentslistsavailableatScienceDirect
Global Ecology and Conservation
journal homepage: http://www.elsevier.com/locate/gecco
Original Research Article
Identifying spatially concordant evolutionary significant units
across multiple species through DNA barcodes: Application to
the conservation genetics of the freshwater fishes of Java and
Bali
Aditya Hutama a, Hadi Dahruddin b, Fre(cid:1)de(cid:1)ric Busson c, Sopian Sauri b,
Philippe Keith c, Renny Kurnia Hadiaty b, Robert Hannerd,
Bambang Suryobroto a, Nicolas Huberte,*
aBogorAgriculturalUniversity,FacultyofMathematicsandNaturalScience,AnimalBioscience,Jl.RayaDarmaga,Bogor16680,
Indonesia
bLIPI,ResearchCenterforBiology,ZoologyDivision,MZB,GedungWidyasatwaloka,Jl.RayaJakartaBogorKm.46,Cibinong-Bogor
16911,Indonesia
cMus(cid:1)eumnationald’Histoirenaturelle,UMR7208(MNHN-CNRS-UPMC-IRD-UCBN),CP026,43,rueCuvier,F-75231ParisCedex05,
France
dBiodiversityInstituteofOntarioandDepartmentofIntegrativeBiology,UniversityofGuelph,Guelph,ON,Canada
eInstitutdeRecherchepourleD(cid:1)eveloppement,UMR226ISEM(UM2-CNRS-IRD),PlaceEug(cid:3)eneBataillon,CC065,F-34095Montpellier
cedex05,France
a r t i c l e i n f o a b s t r a c t
Articlehistory: Delineating Evolutionary Significant Units for conservationpurposes is a crucial step in
Received3July2017 conservation.Acrossadistributionrange,speciesfrequentlydisplaypopulationstructure
Receivedinrevisedform16November2017 that drives the distribution of geneticdiversity.These patterns of geneticstructureand
Accepted16November2017
diversity result from intricate interactions between biogeographic history and de-
Availableonline1December2017
mographic dynamics. Prior biogeographic knowledge, however, is scarcely available, a
trend particularly pronounced in the tropics where the taxonomic impediment is
Keywords:
hampering biogeographic studies and conservation efforts. DNA barcoding has been
Phylogeography
initiallyproposedtofostertaxonomicstudiesthroughthedevelopmentofanautomated
Biodiversityhotspot
Pleistocene molecular system of species identification. While its utility for species identification is
SoutheastAsia increasinglyacknowledged,itsusefulnessforfast andlarge-scaledelineationofESUre-
Coalescenttheory mainstobeexplored.Ifprovedtobeusefulforthatpurpose,DNAbarcodingmayalsoopen
Crypticdiversity new perspectives in conservation by quickly providing preliminary information about
population conservation status. The present study aims at assessing the utility of DNA
barcodingforthedelineationofESUsamongthemostcommonfreshwaterfishspeciesof
JavaandBalithroughthecomparisonofpopulationgeneticstructuresanddiversification
patterns across multiple species. Substantial levels of cryptic diversity are discovered
among the three widely distributed freshwater fish species analyzed with a total of 21
evolutionaryindependent mitochondriallineages(BINs) observedinBarbodesbinotatus,
Channa gachua and Glyptothorax platypogon. The maximum genetic distance for each
coalescenttreerangesfrom6.78to7.76K2PgeneticdistancesforC.gachuaandG.pla-
typogon,respectively.Diversificationandpopulationgeneticanalysessupportascenarioof
allopatricdifferentiation.TheanalysisoftheBINsspatialdistributionindicatesconcordant
* Correspondingauthor.
E-mailaddress:[email protected](N.Hubert).
https://doi.org/10.1016/j.gecco.2017.11.005
2351-9894/©2017TheAuthors.PublishedbyElsevierB.V.ThisisanopenaccessarticleundertheCCBY-NC-NDlicense(http://creativecommons.org/
licenses/by-nc-nd/4.0/).
A.Hutamaetal./GlobalEcologyandConservation12(2017)170e187 171
distributionpatternsamongthethreespeciesthatallowidentifying18ESUs.Implications
fortheconservationgeneticsofthesespeciesarediscussedatthelightofthehistoryofthe
region.
©2017TheAuthors.PublishedbyElsevierB.V.ThisisanopenaccessarticleundertheCC
BY-NC-NDlicense(http://creativecommons.org/licenses/by-nc-nd/4.0/).
1. Introduction
Conservationaimsatpreservingspeciesevolutionarypotentialtosustaintheiradaptiveabilitiesinfluctuatingenviron-
mentsandfuelevolutiononalong-termperspective.Thenatureofthebiologicalunitstobetargetedforachievingthisgoal,
however,hasbeenacontentiousissuesincethe1990s(Crandalletal.,2000).Withtheobjectivetogetaroundtaxonomic
confusion,Ryder(1986)proposedtheconceptofEvolutionarySignificantUnits(ESU)thathedefinedas‘subsetofthemore
inclusiveentityspecies,whichpossessgeneticattributessignificantforthepresentandfuturegenerationsofthespecies’.He
proposedtodelineateESUsbasedonconcordantevidencesfromecological,physiologicalandgeneticperspectives(Ryder,
1986). Following this initial formulation, several recognition criteriawere proposed with varying emphasis on reproduc-
tive isolation (Waples,1991,1995), historical trends in population structure (Avise,1989; Moritz,1994), shared character
states(VoglerandDeSalle,1994)andgeneticorecologicalexchangeability(Crandalletal.,2000).Fromageneticperspective,
mostofthecriteriafocusondetectingtheimprintofdisruptedgeneflowinmitochondrialandnucleargenomes,eitherfroma
short-termperspectivewhenESUsaredelineatedbasedondifferencesinallelefrequencies(Waples,1991,1995;Dizonetal.,
1992; Waples,1995) or long-term perspective when ESUs are based on reciprocal monophyly of their genes genealogies
(Avise,1989;Moritz,1994).
FraserandBernatchez(2001)highlightedthattheproceduresfordelineatingESUsarebasedoncriteria,notmandatory
properties, which applicability depends on biogeographical and ecological context. Since then, major improvements
happened in terms of genome analysis and statistical tools to link population genetic data with landscape ecology
(HoldereggerandWagner,2006;Andersonetal.,2010;SorkandWaits,2010)orlandscapehistory(Beheregaray,2008).These
methodological developments opened new perspectives in conservation, in particular, landscape genetics and niche
modelingshedanewlightatthewaywelookatdispersalandgeneflow(Endoetal.,2014;Gutie(cid:1)rrez-TapiaandPalma,2016).
Usuallyimplementedatsmallspatialscales,thisprocess-basedapproachrequiresapriorknowledgeofpopulationstructure
attheregionalscaleandidentifyingESUsisstillapreliminarystepinconservation(FraserandBernatchez,2001;Pearseand
Crandall,2004).
DelineatingESUiscontextdependentasEarthbiotasoriginatedfromdiversebiogeographicalhistoriesresultinginspecies
ofvaryingage,distributionrangeandpopulationstructure(WiensandDonoghue,2004;Mittelbachetal.,2007;Weirand
Schluter, 2007;Hua andWiens,2010). Thus, prior biogeographicknowledge should provide a useful frameworktoguide
thedelineationofESUs(Avise,2000;FraserandBernatchez,2001).Thisknowledge,however,isnotsufficientlydetailedfor
manyregionsoftheworld,asituationfrequentlyobservedinthetropicswherehighspeciesrichnesshaslargelyamplifiedthe
problem (Beheregaray, 2008; Hubert and Hanner, 2015). DNA barcoding, the use of the cytochromeoxidase I gene as an
internal species tag for molecular identifications, has opened new perspectives on collecting genomic resources across
multiplespecies(Hebertetal.,2004;HebertandGregory,2005;Janzenetal.,2005;Smithetal.,2005;Steinkeetal.,2009;
Hubertetal.,2012).AquickscanofpopulationstructureacrossmultiplespeciesthroughDNAbarcodemayprovideapre-
liminaryandinsightfulapproachinconservation,priortoamorecomprehensiveassessmentusingbothmitochondrialand
nuclearmarkers,by(i)delineatingmitochondriallineagesthatdepartfrompopulationdynamicsanddisplayindependent
mutation/driftdynamics(Hajibabaeietal.,2007;Vernooyetal.,2010;RatnasinghamandHebert,2013;KekkonenandHebert,
2014),(ii)providinggenomicresourcesforquickspeciesandESUsmolecularidentification(Hajibabaeietal.,2011;Gibson
etal.,2014),(iii)easingthecharacterizationofESUsrangedistributionandtheirspatialmatchacrossmultiplespeciesto
identifyconservationpriorities.
Suchstrategymaybeparticularlyrelevantinareasfacingmassiveanthropogenicthreatsandwhereconservationstra-
tegiesareimpededbythelackofappropriateknowledgeonspeciesevolutionarydynamicsandtaxonomicconfusion,atrend
furtheramplifiedbytherecentrarefactionoftaxonomistsworldwide(i.e.taxonomicimpediment).Thissituationiscurrently
observedinSouth-EastAsiawherethefourbiodiversityhotspotsidentifiedareamongthemostthreatenedtodate(Myers
etal.,2000;Lamoreuxetal.,2006;Hoffmanetal.,2010).ThisisparticularlyevidentfortheinsularhotspotsoftheIndo-
nesianarchipelago(i.e.SundalandandWallacea),wheretheimpactofanthropogenicactivitiesisamplifiedbytheirrefugial
state, particularly so in the Sundaland hotspot (Kottelat, 1989; Woodruff, 2010; Lohman et al., 2011). Considering that
IndonesiaincreaseditsGrossDomesticProduct(GDP)andcarbonemissionsby1000%and400%duringthelasttwodecades
(WorldBank),respectively,withapopulationreaching260Millionspeople,itbecomesevidentthatanthropogenicthreats
haveincreasedseverely.ThepresentstudyfocusesondelineatingESUsandtheirspatialconcordancethroughthedevel-
opmentofaDNAbarcodereferencelibraryforthewidespreadfreshwaterfishesoftheislandsofJavaandBali,twooftheless
exploredislandsoftheSundalandhotspot(Hubertetal.,2015b;Dahruddinetal.,2017).Ourobjectiveistoprovidea3-step
generalframeworkforthedelineationofESUsthroughtheanalysisofthespatialandtemporalpopulationstructuresamong
172 A.Hutamaetal./GlobalEcologyandConservation12(2017)170e187
multiplespeciesbasedonmitochondrialcoalescenttreesattheCOIgene(Fig.1).Finally,thisframeworkisusedtoprovide
recommendationsforevidence-basedconservationstrategiesinthearea.
2. Materialsandmethods
2.1. Samplingstrategyandcollectionmanagement
AlargescaleDNAbarcodingcampaignwasconductedbytheauthorsbetweenNovember2012andMay2015across95
sitesinJavaandBaliislands(Dahruddinetal.,2017).Atotalof3310specimens,including162species,110generaand53
familieswerecollected,providingacomprehensiveassessmentoftheJavaneseandBalineseichthyofauna.Amongthose162
species,threespecies(Barbodesbinotatus,GlyptothoraxplatypogonandChannagachua)weresampledinmorethan50%ofthe
sitesvisited,displayedunusuallyhighmaximumwithin-speciesgeneticdistancesbasedonapreviousassessmentusinga
restrictedsetofindividuals(Dahruddinetal.,2017),andassuchrepresentedsuitablecandidatesconsideringtheobjectiveof
thepresentstudies.Eachofthethreespeciesbelongtodifferentordersdisplayingvariedecologicalpreferences(Froeseand
Pauly,2011),andassuch,theyensurethattheconcordanceoftheirESUsspatialpatternsisunlikelytoresultsfromthehistory
of a restrictedsetof aquatic habitats butoriginated through common evolutionary dynamics of broad impacton aquatic
ecosystems(Avise,2000;Aviseetal.,2016).
Thethreespecieswerecollectedacross51sitesdistributedacrosstheislandsofJavaandBali(Fig.2,SupplementaryTable
S1).Specimenswerecapturedusingvariousgearsincludingelectrofishing,seinenets,castnetsandgillnets.Specimenswere
photographed,individuallylabeledandvoucherspecimenswerepreservedina5%formalinsolution.Afincliporamuscle
biopsywastakenforeachspecimenandfixedina96%ethanolsolutionforfurthergeneticanalyses.Bothtissuesandvoucher
specimensweredepositedatthenationalcollectionsattheMuseumZoologicumBogoriense(MZB)intheResearchCentrefor
Biology(RCB)fromtheIndonesianInstituteofSciences(LIPI).
2.2. Sequencingandinternationalrepositories
GenomicDNAwasextractedusingaQiagenDNeasy96tissueextractionkitfollowingthemanufacturer'sspecifications.A
651-bpsegmentfromthe50regionofthecytochromeoxidaseIgene(COI)wasamplifiedusingprimerscocktailsC_FishF1t1/
C_FishR1t1includingaM13tails(Ivanovaetal.,2007).PCRamplificationsweredoneonaVeriti96-wellFast(ABI-Applied
Biosystems)thermocyclerwithafinalvolumeof10.0mlcontaining5.0mlBuffer2X,3.3mlultrapurewater,1.0mleachprimer
(10mM),0.2mlenzymePhire®HotStartIIDNApolymerase(5U)and0.5mlofDNAtemplate(~50ng).Amplificationswere
conductedasfollow:initialdenaturationat98(cid:2)Cfor5minfollowedby30cyclesdenaturationat98(cid:2)Cfor5s,annealingat
56(cid:2)Cfor20sandextensionat72(cid:2)Cfor30s,followedbyafinalextensionstepat72(cid:2)Cfor5min.ThePCRproductswere
purifiedwithExoSap-IT®(USBCorporation,Cleveland,OH,USA)andsequencedinbothdirections.Sequencingreactionswere
performed using the “BigDye® Terminator v3.1 Cycle Sequencing Ready Reaction” and sequencing was performed on the
automatic sequencer ABI 3130 DNA Analyzer (Applied Biosystems). The sequences and collateral information have been
depositedinBOLD(RatnasinghamandHebert,2007)intheprojectsBIFGAandBIFGinthecontainer‘BarcodingIndonesian
Fishes’ofthe‘BarcodingFish(FishBOL)’campaignandDNAsequencesweresubmittedtoGenBank(accessionnumbersare
accessibledirectlyattheindividualrecordsinBOLD).
2.3. Populationgeneticstructure(Fig.1,step1a)
Weexaminedthedistributionofmolecularvariancethroughanadditivepartitioningacrossincreasingspatialscalesas
implementedinAMOVA(Excoffieretal.,1992;ExcoffierandSmouse,1994).WeoptedfortheSAMOVAversionthatboth
definesgroups ofpopulationswithout apriori partitioningschemeandestimatesmolecularvariancewithinpopulations,
amongpopulationswithingroupsandamonggroups,basedonasimulatedannealingapproach(Dupanloupetal.,2002).
SAMOVAisnotadecision-basedmethodpointingtotheoptimalnumberofgroups,wehenceappliedanempiricalthreshold
andstoppedincreasingthenumberofgroupsoncethepartitioningschemeproducedatleastonegroupconsistingofasingle
sampling site. The SAMOVA were performed for each species using the software SAMOVA 2.0 (Dupanloup et al., 2002).
SAMOVApartitioningschemeswerefurthercomparedtotheresultsofahierarchicalclusteringbasedongeneticdistances.
MeangeneticK2PdistanceswerecomputedamongsamplingsitesusingMEGA6(Tamuraetal.,2013)andusedtoproduce
hierarchicalclustersderivedfromthecompletelinkagealgorithmasimplementedinhclustfunctionoftheRStatsver.3.1.2
package(R_Core_Team,2014).Finally,traditionalparametersofpopulationgeneticdiversityincludingthenumberofhap-
lotypes,nucleotidicdiversityandmeanpairwiseK2PdistanceforeachgroupswerecomputedusingARLEQUIN3.5(Excoffier
etal.,2005).
2.4. Delineatingmitochondriallineagesandinferringtheirdiversification(Fig.1,step1b)
Severalalternativemethodshavebeenproposedfordelineatingmolecularlineages(SchlossandHandelsman,2005;Pons
etal.,2006;Puillandreetal.,2012).Theyhaveallincommontodetecttransitionzonesinbranchingpatternsresultingfrom
differentsegmentsofthegenegenealogiesthatoriginatedfromphylogeneticdiversification(speciationandextinction)or
Fig.1. ConceptualframeworkdevelopedinthepresentstudyforthedelineationofESUs.Step1a,detectionofgroupsofpopulationsdifferentiatedbytheir
haplotypesfrequencies.Step1b,detectionofmolecularlineageswithindependentevolutionarydynamics(e.g.lowerconnectance).Step2,comparingpopulation
groupswithmolecularlineages.Ifpopulationgeneticstructureresultsfromancientfragmentationofthepopulations,acorrelationbetweengeneticgroupsand
BINsisexpected.Conversely,thelackofcorrelationwouldindicatethatpopulationgroupsoriginatedrecentlyandeithershareancientpolymorphismorhave
beenconnectedbygeneflowinarecentpast(Fu,1999;NielsenandWakeley,2001;Wakeley,2001,2003).Alongthesameline,severalBINsmaybedelineated
withinageneticgroupasaconsequenceofthestochasticnatureofthecoalescentbutnotdisruptedgeneflow(Hudson,1982;Kingman,1982;Tajima,1983).Step
3,DelineationofESUsforindividualspecies.ESUsaredefinedbasedoneithervariationofhaplotypefrequenciesorindependentmitochondriallineagesafter
comparisonswiththephylogroupsdefinedasgroupsofpopulationsharingsimilarsetsofmitochondriallineages.
1
7
4
A
.
H
u
ta
m
a
e
t
a
l.
/
G
lo
b
a
l
E
co
lo
g
y
a
n
d
C
o
n
se
rv
a
tio
n
1
2
(2
0
1
7
)
17
0
Fig.2. Locationofthe51collectionsitesforthesamplesanalyzedinthisstudy. 1e
8
7
A.Hutamaetal./GlobalEcologyandConservation12(2017)170e187 175
Table1
Partitioningofthemolecularvarianceatvariousspatialscales,fixationindexesandnumberofpopulationgroupsasinferredfromSAMOVA.Thepercentage
columnindicatestheamountoftotalvarianceexplainedbyeachofthehierarchicallevelsaccordingtothenumberofgroupsofpopulation.F-statistics
estimatethecorrelationamonghaplotypesateachofthehierarchicallevelsexaminedandtheirsignificantdeparturefromarandomdistributionofthe
haplotypewastestedthroughrandomizationacross1000permutations.
NumberofGroup C.gachua G.platypogon B.binotatus
5 7 6
Variancecomponent Variance%of P F-statisticsVariance%of P F-statisticsVariance%of P F-statistics
total total total
Amonggroups 8.788 70.56 <0.001FCT¼0.7068.042 78.75 <0.001FCT¼0.7886.958 69.90 <0.001FCT¼0.702
Amongpopulationswithin 3.071 24.69 <0.001FSC¼0.8391.516 11.34 <0.001FSC¼0.5342.094 10.75 <0.001FSC¼0.361
groups
Withinpopulations 0.591 4.75 <0.001FST¼0.953 1.015 9.91 <0.001FST¼0.9010.903 19.05 <0.001FST¼0.890
coalescent dynamics (mutation and drift). We opted for the Refined Single Linkage (RESL) algorithm that considers the
number of connections of each sequences in a network estimated through the silhouette index (Rousseuw, 1987) as
implemented in BOLD. Sequence connectivity is explored through random walks and optimal partitioning schemes are
identifiedthroughMarkovclustering.Attheend,eachclusterofsequenceisassignedtoaBarcodeIndexNumber(BIN)in
BOLD(RatnasinghamandHebert,2013).
Once BINs were delineated (Table S2), their timing of diversification was explored through the Bayesian approach
implementedinBEAST1.8.1(Drummondetal.,2012).Inordertoestablishrobustprior,thebest-fitsubstitutionmodelwas
selectedthroughtheBayesianInformationCriterion(BIC) asimplemented inJMODELTEST2.1.7 (Darribaetal., 2012)and
furtherusedasapriorforthejointreconstructionoftreetopologyanddivergencetimes.Theinitialtreetopologywasob-
tainedwithanUPGMAstartingtreeandaCoalescentmodelwasusedasatreeprior(Kingman,1982).Weusedthecanonical
fishsubstitutionrateof1.2%ofgeneticdivergenceperMillionsyearsformitochondrialproteincodinggene(Bermingham
etal.,1997)andappliedittothemaximumK2PdistancewithineachspeciestoestimatetheageoftheMostRecentCom-
monAncestor (MRCA) tobe used asa priorfor Bayesiananalyses. The priorupperand lower boundsfor theMRCAs age
intervalswerederivedfromtheknownhighestandsmallestsubstitutionratesforfishmitochondrialgenomes(Hardmanand
Lundberg,2006;Readetal.,2006).WeranoneMCMCof10(cid:3)106steplong,sampledevery1000stateswithaburn-inperiod
of10000.ThemaximumcredibilitytreewasobtainedwithTreeAnnotator1.8.1afteranadditionalburn-inperiodof10000.
Mediannodeagesand95%highestposteriordensity(HPD)intervalswereplottedinthechronogram.Duplicatedsequences
wereremovedfortheseanalyses.WefurtherexploredthepropertiesofBINsdiversificationthroughLineageThroughTime
(LTT)plots(Harveyetal.,1994)andtheGeneralizedSkylinePlots(GSP)(StrimmerandPybus,2001)methods.TheBayesian
LTTwasconductedthroughsimilarMCMCparametersasofthetreeanalysesandduplicatedsequenceswereremoved.The
BayesianGSPanalyseswereconductedontheentiredataset,includingduplicatedsequences,aHKYsubstitutionmodeland
MCMCchainsof50(cid:3)106stepslongsampledasdescribedforLTTandtreereconstructionanalyses.
2.5. DelineatingESU(Fig.1,steps2&3)
WeexaminedtherelationshipbetweenthegeneticgroupsdelineatedbySAMOVAandtheBINs(Fig.1,Step2)bytesting
theindependenceofeachpartitioningschemeusingthePearsonchi-squaretestofindependenceasimplementedintheR
Statsver.3.1.2package(R_Core_Team,2014).TheobjectivewastoexaminetowhatextenttheSAMOVAgroupsofpopulations
weredeterminedbytheBINsdelineatedbytheRESLalgorithm(Fig.1,Step2).Wefurtherexaminedtheconcordanceofthe
BINsspatialdistributionamongspeciesbyperformingageneralhierarchicalclusteranalysisbasedontheaveragephylo-
geneticdistanceoftheBINsamongsitesasimplementedinhclustfunctionoftheRStatsver.3.1.2package(R_Core_Team,
2014). The matrixof phylogenetic distance among sites was computed based on the composite chronogram of the three
graftedmaximumcredibilitytreesofBarbodesbinotatus,ChannagachuaandGlyptothoraxplatypogonandBINsoccurrence
data using the R package PICANTE (Kembel et al., 2010). The composite chronogramwas constructed byadding internal
branchessetto0.001Millionsyearsbetween(i)theMRCAofB.binotatusandtheMRCAofG.platypogon,(ii)theMRCAofboth
B.binotatusandG.platypogonandtheMRCAofthethreespecies.TheinternalbranchbetweentheMRCAofC.gachuaandthe
MRCAofthethreespecieswasfurtheradjustedtoproduceanultrametrictree.BINswereconsideredasabsentfromsitesout
oftheirdistributionrangebutBINsabsencewascodedasmissingdataatsiteswithintheirdistributionrangeinorderto
accountforsamplinguncertaintyi.e.aBINpresentbutnotsampled(TableS3).Wefinallyproduceddistributionmapsofthe
populationgroups(Fig.1,Step1a),theBINs(Fig.1,Step1b)andthegeneralhierarchicalcluster(Fig.1,Step3).Thegeneral
hierarchicalclusterwasfurtheredusedtodefinephylogroupsi.e.groupsofsiteshostingphylogeneticallyrelatedBINs.The
distributionofthephylogroupswasusedasatemplatetospotphylogeographicbreaks,asexemplifiedbythephylogroup
geographicboundaries,andcomparedtothedistributionofSAMOVAgroupsandBINStoproposeESUs(Fig.1,Step3).
1
7
6
A
.
H
u
ta
m
a
e
t
a
l.
/
G
lo
b
a
l
E
co
lo
g
y
a
n
d
C
o
n
se
rv
a
tio
n
1
2
(2
0
1
7
)
17
0
e
1
8
Fig.3. PopulationgeneticstructureofBarbodesbinotatus(a),Channagachua(b)andGlyptothoraxplatypogon(c)asinferredfromthehierarchicalclusteranalysisandSAMOVA. 7
A.Hutamaetal./GlobalEcologyandConservation12(2017)170e187 177
3. Results
A total of 317 new sequences were successfully generated for the three species including 108 sequences of Barbodes
binotatus,99sequencesofChannagachuaand110sequencesofGlyptothoraxplatypogon.Togetherwiththe37sequences
previouslypublishedforthesespecies(Dahruddinetal.,2017),atotalof354sequenceswereanalyzedincluding122se-
quencesofBarbodesbinotatus,109sequencesofChannagachuaand123sequencesofGlyptothoraxplatypogon(TableS1).All
the sequences were above 500 bp of length and no codon stops were detected, suggesting that the sequences collected
representfunctionalcodingregions.The354sequenceswerecollectedfrom39sitesforB.binotatus,28sitesforC.gachuaand
25sitesforG.platypogon.
3.1. Populationgeneticstructureandmolecularvariance
A total of 5, 7 and 6 groups of populations were delineated for Channa gachua, Glyptothorax platypogon and Barbodes
binotatusrespectively(Table1,TableS1).Forthethreespecies,mostofthemolecularvarianceisexplainedbydifferences
amonggroupsofpopulationswithapercentageofthetotalvariancerangingfrom69.9percentinB.binotatusto78.75percent
inG.platypogon.‘Amongpopulationswithingroups’isthesecondlevelintermsofthetotalvarianceexplainedwith24.69
percentforC.gachuaand11.34percentforG.platypogon.B.binotatusdiffersfromthetwootherspeciesinhavingaslightly
higherproportionofthetotalvarianceexplainedbydifferenceswithinpopulations,with19.05percent,thanamongpop-
ulations within groups, with 10.75 percent. The fixation indexes are significant at all spatial scales supporting spatially
structured populations in the three species. This result was further confirmed by the hierarchical cluster analysis of the
populationsthatindicatesasubstantialgeneticdifferentiationofeachgroupwithK2Pgeneticdistancesamonggroupsabove
0.04,0.02and0.01forC.gachua,G.platypogonandB.binotatus,respectively(Fig.3).Withingroupsofpopulations,haplotype
diversityishighonaveragewithhabove0.5exceptingforpopulationgroupsII,VandVIIofG.platypogonwithhbelow0.2
(Table2).Thishighhaplotypediversityisopposedtothelownucleotidediversityandthelowmean-pairwisedifferences
withingroupsforthethreespecies.
3.2. DelineationanddiversificationofBINs
Atotalof9,7and3molecularlineages(i.e.BINs)weredelineatedbytheRESLalgorithmonBOLDforChannagachua,
GlyptothoraxplatypogonandBarbodesbinotatus,respectively(TableS1).BothmaximumandaverageK2Pdistancesarelowon
averagewithinBINsandcontrastwiththehighmaximumK2Pdistancesobservedforeachofthethreecoalescenttrees(Table
3).Amongthe19BINsdelineated,fourwererepresentedbysingletonsinC.gachua(ACQ3951,ACQ6941,ACQ6939,ACQ6940).
TheBayesiananalysesyieldedthreemaximumcredibilitytrees(Fig.4a,b,4C)thatconfirmedthecontrastbetweenthedeep
divergence among BINs, ranging from 0.47 to 2.75 Millionyears ago (Ma), and the shallow coalescent depth of the BINs
rangingfrom0.08to0.89Ma(Table3).TheageoftheMRCAwasverysimilaramongthethreespeciesrangingfrom2.71Ma
forB.binotatusto3.14MaforG.platypogon(Table3).PlottingtheBINscoalescentdepthontheLTTcurvesfurtherconfirmed
thatthemoleculardiversitywithinBINsaccumulatedveryrecentlywithinthethreespecies(Fig.4d,e,4f).TheBayesianGSP
Table2
Summarystatisticsofthegeneticdiversityincludingthesamplingsize(N),haplotypediversity(h),nucleotidediversity(p)andmeannumberofpairwise
differencesamonghaplotypes(mean-pairwisedifferences)foreachofthepopulationgroups.
N h p mean-pairwisedifferences
Barbodesbinotatus
I 4 1 0.019 12.17
II 22 0.87 0.002 1.41
III 33 0.83 0.002 1.37
IV 12 0.55 0.005 2.73
V 36 0.82 0.004 2.35
VI 15 0.89 0.011 6.88
Channagachua
I 14 0.88 0.006 3.91
II 34 0.93 0.010 6.07
III 11 0.51 0.008 0.51
IV 36 0.81 0.009 5.85
V 14 0.58 0.011 6.97
Glyptothoraxplatypogon
I 16 0.58 0.001 0.75
II 34 0.28 0.004 2.29
III 6 0.60 0.001 0.60
IV 39 0.75 0.006 3.62
V 10 0.20 0.002 1.20
VI 7 0.91 0.032 20.80
VII 11 0.18 0.001 0.18
178 A.Hutamaetal./GlobalEcologyandConservation12(2017)170e187
Table3
SummarystatisticsofthegeneticK2Pdistancesandageestimates.ThemaximumandaverageK2PdistancesareprovidedforeachBINsandthemaximum
K2Pdistancesareprovidedfortheentirecoalescenttrees.TheageestimateoftheMRCAisprovidedbasedonthreealternativehypothesesofmolecularclock
including 0.005 genetic divergence per millionyears (Hardman and Lundberg, 2006; Read et al., 2006), 0.012 genetic divergence per millionyears
(Berminghametal.,1997)and0.02geneticdivergencepermillionyears(Readetal.,2006).DivergenceestimatesarederivedfromtheBayesiananalyses
basedontheH2calibrationwithupperandlowerboundsofthepriorageintervalsdefinedbythemolecularclockhypothesesH1andH3.
K2Pdistance(%) MRCAH1a(Myr) MRCAH2b(Myr) MRCAH3c(Myr) Divergenceestimates(Myr) Coalescentdepth(Myr)
Maximum Average
Channagachua
Coalescent 6.78 6.78 2.83 1.70 2.75(6.02-1.70)
ACQ0292 2.1 0.87 2.04(4.56-0.82) 0.66(1.60-0.26)
ACQ0290 1.4 0.61 2.04(4.56-0.82) 0.61(1.49-0.16)
ACQ0291 0.15 0.08 2.75(6.02-1.70) 0.08(0.31-0.01)
ACQ3951 - - 1.23(2.87-0.41) -
ACQ6941 - - 0.73(1.75-0.24) -
ACQ6940 - - 0.47(1.17-0.24) -
ACQ6939 0.00 0.00 0.47(1.17-0.24) -
ACQ3952 1.4 0.31 0.87(2.08-0.31) 0.14(1.40-0.14)
ACQ3950 0.93 0.32 0.87(2.08-0.31) 0.33(0.84-0.08)
Glyptothoraxplatypogon
Coalescent 7.76 7.76 3.24 1.94 3.14(6.74-1.94)
ACQ5850 - - 1.05(2.68-0.21) 0.13(0.45-0.01)
AAY1028 0.00 0.00 1.05(2.68-0.21) 0.08(0.33-0.01)
ACQ5898 0.93 0.23 1.75(4.22-0.51) 0.39(1.13-0.05)
ACP6225 0.31 0.12 1.98(4.71-0.73) 0.23(0.72-0.03)
ACQ6223 1.24 0.55 1.09(2.63-0.32) 0.50(1.35-0.14)
ACP6117 0.15 0.09 0.69(1.75-0.18) 0.08(0.34-0.01)
ACQ6224 0.46 0.05 0.69(1.75-0.18) 0.21(0.56-0.05)
Barbodesbinotatus
Coalescent 6.80 6.8 2.84 1.70 2.71(5.96-1.70)
ACP6290 0.01 0.00 2.71(5.96-1.70) 0.12(0.41-0.01)
ACP6025 0.57 1.09 1.74(4.26-0.48) 0.38(1.12-0.08)
ACP5712 1.94 0.73 1.74(4.26-0.48) 0.89(2.33-0.48)
a 0.005.
b 0.012.
c 0.02.
yieldedverysimilardemographictrajectoriesforthethreespecieswithaglobaltrendofconstantpopulationsizeoverthe
last2millionsyearsandasteepdeclineofpopulationsizeduringthelast100.000years(Fig.4g,h,4i).
3.3. IdentifyingphylogroupsanddelineatingESUs
ThenullhypothesisofindependencebetweenpopulationgroupsdelineatedbySAMOVAandBINswasrejectedbythe
Pearsonchi-squaretest(Х2¼4486.7,p-value<0.001)suggestingthatthedelineationofpopulationgroupsislargelysup-
ported by alternative distribution of BINs among populations. The general hierarchical cluster constructed using the
phylogeneticdistanceamongsitesresultedinthedelineationoffivemajorphylogroups,divergingbymorethan5Million
yearsonaverage,andorderedintotwomainclusters(Fig.5).ThegeographicrangeofphylogroupsIandIIshowednooverlap
andallowedidentifyingthreephylogeographicbreaks(Fig.5b).PhylogroupIisrestrictedtothewesternpartofJavawhile
phylogroupIIshowsanalternativedistributionineasternJavaandBali.Thefirstphylogeographicbreakidentifiedislocated
intheJava-BalistraightbasedonthesegregationofphylogroupsIandIIoneachsideofthestraight(Fig.5,break1).The
secondphylogeographicbreakidentifiedislocatedincentralJava,inbetweenthewesternandcentralvolcanicarches(Fig.5,
break2).Thethirdphylogeographicbreakislocatedonthewesternmostslopeofthewesternvolcanicarchbasedonthe
westernlimitofphylogroupII(Fig.5,break3).ThedistributionrangesofphylogroupsIII,IVandVaremoreintricatedbut
showsomesimilaritieswithrangedistributionofphylogroupsIandII.TheJava-Balistraightisalsoidentifiedasaphylo-
geographicbreaksegregatingphylogroupsIVandV,eachbeingalternativelydistributedinJavaorBali(Fig.5,break1).The
centralJavaphylogeographicbreakisalsoassociatedwiththesegregationofphylogroupsIIIandV(Fig.5,break2)buttwo
additionalphylogeographicbreaksareidentifiedincentraljava,inbetweenthewesternandcentralvolcanicarches,thatare
associatedwiththesegregationofphylogroupIIIfromthephylogroupV(Fig.5,break4)andthesegregationofanisolatedsite
ofphylogroupVfromtheeasternsites(Fig.5,break5).Thewesternmostphylogeographicbreakisalsoassociatedwiththe
segregationofphylogroupsIVandVintheirwesternrangebutwithsomeoverlaps(Fig.5,break3).
TheboundariesofpopulationgroupsandBINsrangedistributionaregenerallyassociatedwiththefivephylogeographic
breaksforthethreespecieswithsomeexceptions(Fig.6).Threecasesoftransdistributionacrossphylogeographicbreaksare
detectedinC.gachua(Fig.6aandb).Thefirstcaseislocatedatthephylogeographicbreak3withthepopulationgroupV
(Fig.6a,red).ThepopulationgroupV,however,isassociatedwithfourindependentBINs.Thesecondcaseisassociatedwith
the phylogeographic break 2 with the population group IV (Fig. 6a, yellow). Finally, the population groups I and II are
A
.
H
u
ta
m
a
e
t
a
l.
/
G
lo
b
a
l
E
co
lo
g
y
a
n
d
C
o
n
se
rv
a
tio
n
1
2
(2
0
1
7
)
1
7
0
e
1
8
7
Fig.4. BINsdiversificationpatternsofChannagachua(a,dandg),Glyptothoraxplatypogon(b,eandh)andBarbodesbinotatus(c,fandi)includingMaximumCredibilityTrees(a,b,c);LineageThroughTimeplot(d,e,f)
andBayesianGeneralizedSkylinePlots(g,h,i).Solidblacklineind,eandfisthemediandiversityaccumulationcurve,blueshadedareaisthe95%highestposteriordensityintervalsanddottedlinesrepresenttheBIN
coalescentdepthinMillionyears.Thesolidblacklineing,handirepresentthemedianeffectivepopulationsize,blueshadedareaisthe95%highestposteriordensityintervalsandPopulationsizescalar(in
millions)¼effectivepopulationsizexgenerationtime.CalibrationsarederivedfrompreviouslypublishedmolecularclockhypothesesappliedtothemaximumK2Pdistanceforeachspecies-i.e.ageoftheMRCA(see
Table3).(Forinterpretationofthereferencestocolourinthisfigurelegend,thereaderisreferredtothewebversionofthisarticle.)
1
7
9