Phylogenetic prospecting for cryptic species of the genus Merluccius (Actinopterygii: Merlucciidae)

Hakes of the genus Merluccius include 11 valid species as well a number of rare morphotypes suspected to be “cryptic species”. Concatenated nucDNA ITS1-rDNA and mtDNA cyt b sequences plus nested ITS1Nes sequences allowed to ascribe 14 specimens of nine rare morphotypes from the South Pacific and the South Atlantic to the phylogenetic backbone of this genus. Bayesian analyses pointed to M. bilinearis and M. albidus as the oldest species of the genus and the New World cluster, respectively. The phylogenetic status of M. angustimanus from the upper Gulf of California suggests its hybrid origin between M. gayi and M. productus from about 0.25 MYA, although an ever since confinement of a subset of those species cannot be ruled out. The molecular phylodiagnostic test suggests a common origin of all rare morphotypes and the absence of cryptic hake species in the Southern Cone. The molecular background of the morphotypes distributed between the Western Pacific South of New Zealand and the western Atlantic South of Argentina is compatible with their hybrid origin between M. gayi and both, M. australis or M. hubbsi, respectively.

The genus Merluccius comprises 11 valid species that occur on most temperate and tropical continental shelves except the Asian shores of the Pacific Ocean 1 . Hakes show an anti-tropical distribution in the Atlantic Ocean and the Eastern Pacific and a latitudinal bathymetric overlap between isotherms 7 °C and 23 °C [2][3][4] . Based on osteological data 5,6 it is believed that genus Palaeogadus as ancestor of genus Merluccius, arose near Greenland in the early Eocene (ca. 50 MYA) 7 , dispersed southwards along the North American and Eurasian shelves and entered the Pacific 8 . The earliest known merluccid fossils date back to the Middle Oligocene (ca. 27-33 MYA) in a large inland sea that covered much of central Europe and connected to a temperate Arctic Ocean 9 . It is believed that either an ancestral species of Merluccius or the extant M. bilinearis experienced an evolutionary radiation in two superclusters, i.e. Old World hakes (the Euro-African supercluster) and New World hakes (the American supercluster) [5][6][7][8] . Also, it is hypothesized that a widening rift between Europe and North America plates prompted vicariant speciation and that recurrent dispersal events and adaptation to temperature regimes also played a role in the speciation of this genus 10 . Subsequent geological events such as the closure of the Panama Seaway over 3.5 MYA acted as a geographical barrier between Atlantic and Pacific lineages 8,11 . Successive population fragmentation and expansion due to climatic oscillations during Pleistocene glaciations allowing founder phenomena cannot be ruled out 12 . Such origin and dispersal hypotheses are congruent with the actual phylogeny of the genus worked out after parasite data 8 morphology 3,5,6,13 and genetic data [14][15][16] .
The Old World supercluster comprises five well-defined species and with occasional hybrids between sympatric species, e.g. M. capensis x M. paradoxus 17 . The New World supercluster comprises three clusters of two species each, an Atlantic north cluster that comprises M. bilinearis and M. albidus, a southern cluster that comprises M. hubbsi and M. australis, and a Pacific cluster that groups M. gayi and M. productus. Molecular systematics of Merluccius generally distinguishes those 11 species 18 , however there are still knowledge gaps in hake taxonomy as some specimens found in regions of species overlap show significant morphological divergence from extant species 13 . For instance, the phylogenetic relationships within the New World supercluster have been repeatedly  24 . Several of those morphotypes have been recently taken as synonymous with extant species. For instance, the three stocks of M. productus believed to exist in the northeast Pacific corridor from Washington State to Costa Rica 25,26 likely belong to a single hake species 27,28 . The high variability and overlap of meristic traits between M. hernandezi and M. angustimanus suggests they are synonymous forms 13 and that M. angustimanus is a subpopulation of M. productus confined to the northern Gulf of California 27 . Also, meristic (see Fig. 1 from 23 ) and molecular analyses 29,30 showed that M. patagonicus could be synonymous with M. hubbsi and that M. tasmanicus could be synonymous with M. australis 31 . Taxonomic uncertainties are expected in closelyrelated species that show extensive overlap in morphological and meristic traits. Resolution of those questions has been complicated by the small number of morphotype samples available, the lack of a systematic sampling plan over a reasonable number of localities and the choice of the most appropriate phylogenetic reconstruction algorithm.
We address those issues with the most comprehensive sample collection yet made in this genus, which comprises all 11 valid species 13 as well as M. angustimanus from the Gulf of California 21 and 14 specimens of nine rare morphotypes classified as M. tasmanicus from New Zealand waters 23 , M. polylepis from the Pacific coast of Chile 32 , M. patagonicus from the Atlantic South 23 , uncatalogued specimens of M. hubbsi from the Beagle Channel and Puerto Madryn in Argentina 23 , and rare specimens of M. australis from Chile 23 and New Zealand 24 . An integrative multimarker phylogenetic reconstruction of genus Merluccius is enforced to determine both, the phylogenetic congruence and synergy between mtDNA (cyt b) and nucDNA (ITS1) sequences for inference of phylogenetic relationships in this genus and the genetic prospecting for cryptic species within the New World supercluster, as suspected on the rare morphotypes described so far.

Results
Description of samples and sequences. We examined 1205 specimens from 11 valid hake species and M. angustimanus (Table 1) as well as nine morphotypes ( Table 2). The aligned ITS1 region had a length of 692 bp and comprised the ITS1-rDNA sequence, 53 bp from the 3′-end of the 18S-rDNA gene and 20 bp from the 5′-end of the 5.8S-rDNA gene 33 . A total of 254 variable sites were found among 85 specimens from 12 hake species (including M. angustimanus as putative species) (Supp. Table S1). The ITS1 sequences showed similar %GC content between species and a general low transversion rate (Supp. Table S2). The full dataset of 85 ITS1 sequences comprised 19 variants (Supp. Table S3). The aligned cyt b region had a length of 465 bp and comprised 428 bp from the 5′-end of the cyt b and 37 bp from the 3′-end of the mitochondrial DNA gene tRNA-Glu. A total of 129 variable sites were found among 66 specimens from 12 hake species as including M. angustimanus (Supp. Table S1). The 66 cyt b sequences comprised 29 haplotypes as one to four per species (Supp.  . Table S1) and identified 25 variants (Supp. Figure S4). The aligned ITS1Nes sequences from 39 specimens of reference and 46 clones from 14 specimens of nine morphotypes had a length of 66 bp and 13 variable sites (Supp . Table S1). Those 85 ITS1Nes sequences comprised 11 variants (HakeITS1Nes.1-11) ( Table 3). HakeITS1Nes.2 was shared among the Pacific species and four morphotypes, HakeITS1Nes.5 was shared between the Atlantic North species, HakeITS1Nes.6 was shared among all morphotypes and HakeITS1Nes.9 was shared between two morphotypes (Table 3). All the morphotypes shared at least one ITS1Nes variant among each other or with known species, but also exhibited specific variants, e.g. HakeITS1Nes.7-8-10-11. The recombination parameter (R) detected three pairs of sites in the ITS1Nes region with at least one recombination event 34 Table 3).
Clustering methods. The major groups of species as inferred from the Euclidean divergence of the PCoA correctly identified a) the two major complexes of hakes, Old World versus New World hakes (Supp. Figure S1a, b) Atlantic versus Pacific New World hakes (Supp. Figure S1b, and c) the closeness of hake morphotypes to both, the Pacific group and the Austral group (Supp. Figure S1c). The AMOVA agreed with the partition among the groups identified by PCoA using ITS1 variants, i.e. Old World/ New World and Atlantic New World /Pacific New World /Austral New World (Table 5). Noteworthy, the joint analysis including the morphotypes in the Pacific cluster produced the largest within cluster variation as compared to any other hierarchical level.  Table 1 and are followed either by the number of ITS1 variants per species ( Fig. 1a and Supp. Table S3) or by the number of cyt b haplotypes per species (Fig. 1b  www.nature.com/scientificreports/ Phylogenetic reconstruction. Evaluation of the relationship between the best phylogenetic signal and the most plausible phylogenetic scenario upon previous knowledge on this genus showed that Bootstrap support from Bayesian inference was maximum on concatenated ITS1-cyt b sequences ( Table 6). The most congruent topology from likelihood-based trees was issued from concatenated ITS1-cyt b sequences ( Fig. 2a; Table 6).
ITS1-based phylogeny. The phylogenetic reconstruction of ML (-lnL = 1285.511) performed with PAUP on ITS1 variants recovered a correct supercluster split, a large polytomy within superclusters and an unclear age status of extant species (Supp. Figure S2a). The Bayesian reconstruction from MRBAYES on ITS1 variants recovered a single polytomy comprising all the well-recognized species clusters (Supp. Figure S2b). The Bayesian reconstruction of BEAST on ITS1 variants supported a correct supercluster split, the absence of polytomies and the odd basal placement of M. hubbsi in the New World supercluster (Fig. 1a).  www.nature.com/scientificreports/ Cyt b-based phylogeny. Likelihood computation of dN and dS using HyPhy package 35 on cyt b sequences showed that six out of 142 codons contained nonsynonymous substitutions. Provided that p-values were not significant, the null hypothesis of neutral evolution was accepted. The best -lnL value and bootstrap support were observed on cyt b haplotypes but its topological features were unsatisfactory due to multiple polytomies (Supp. Fig S3a; Table 6). The ML phylogenetic reconstruction on cyt b haplotypes (-lnL = 1137.282) worked out with PAUP showed a poor definition of the Old World supercluster with a high impact of polytomy (Supp. Figure S3a). The best topology among Bayesian-based supported trees was obtained upon inference on cyt b haplotypes (Fig. 1b). The Bayesian phylogenetic reconstruction from MRBAYES on cyt b haplotypes recovered a strong supercluster split and a polytomy within the New World supercluster to which M. bilinearis was ancestor (Supp. Figure S3b). The Bayesian phylogenetic reconstruction of BEAST on cyt b haplotypes recovered a wellsupported both, supercluster split and cladogenesis within superclusters, where M. bilinearis was basal to New World hakes (Fig. 1b).  www.nature.com/scientificreports/ www.nature.com/scientificreports/  Table 1 and are followed by the alphanumeric entry code for each specimen in the authors laboratory, e.g. aust_ma01_2 is the specimen No. 2 of sample ma01 from M. australis (aust). The maps have been modified after 13 . www.nature.com/scientificreports/ Concatenated ITS1-cyt b phylogeny. The phylogenetic reconstruction using concatenated sequences of both genes showed bootstrap values > 90% for species and clusters. ML (Fig. 2a) and Bayesian (Fig. 2b) Figure S4). The phylogenies inferred from the methods ML and Bayesian were topologically similar to each other (Fig. 3). The Bayesian reconstruction failed to resolve the phylogenetic status of M. bilinearis regarding the two superclusters and exhibited a large polytomy within the New World supercluster.
Coalescence. The relative coalescence times inferred from a typical 2% mutation rate of cyt b in fish was ~ 9 MYA between Merluccius and the codfish genus Gadus, ~ 3.5 MYA between the New World and the Old World superclusters, and ~ 2.3 MYA between M. bilinearis and other New World hakes (Fig. 4). Atlantic and Pacific New World hakes would have diverged some ~ 1.4 MYA.
ITS1Nes-based phylodiagnosis. The Neighbor-Joining reconstruction performed with the ITS1Nes fragment included six valid species and M. angustimanus, and comformed the basal phylogenetic backbone of this genus in the New World against which new phylogenetic hypotheses can be tested. The basal support afforded from the ITS1Nes fragment was quite close to that obtained on the full ITS1 sequence (Fig. 1a) but including a polytomy for M. gayi-M. productus (Fig. 5).  (Fig. 6). Two additional minor clusters of morphotypes closely branched either to the Pacific cluster or to the Austral cluster.

Discussion
Genetic homogeneity is the null hypothesis in widely-distributed marine taxa with large population size and sometimes morphological divergence may be a first clue to hypothesize on the existence of genetically divergent units 13 . Demonstration of intraspecific genetic differentiation requires comprehensive spatio-temporal sample designs, the identification of suitable genetic markers at the resolution level concerned, and the choice of appropriate phylogenetic algorithms 18 . The nuclear ITS1-rDNA region has been successfully applied in the identification of closely related taxa 36 , in fish phylogeography 37 , in phylogenetic inference 38 as well as in forensic authentication of species 33 . Also, the cyt b gene is a well-known mtDNA gene in structure and function 39 and is useful in phylogenetic reconstruction at many taxonomic levels, including congeneric species and confamiliar genera 40 but see exceptions 41 . The synergy afforded from concatenated analyses of the nuclear ITS1-rDNA and the mitochondrial DNA cyt b allows comparison of interspecific levels of divergence and to achieving phylogenetic scenarios unaffordable from single-markers approaches 42 . The higher interspecific variation observed between New World and Old World superclusters (e.g. ITS1, D xy ≈ 0.090) than within superclusters (e.g. ITS1,  www.nature.com/scientificreports/ D xy ≈ 0.040) is expected if previous calibrations were robust 33 . In consequence a higher net evolutionary divergence (d) existed among Old World species (d ≈ 0.080) than among New World species (d ≈ 0.030). In contrast, the minimum interspecific divergences (d and D xy ) observed among Pacific New World hakes (M. gayi-M. productus-M. angustimanus) represented the most recent evolutionary scenario in this genus, as reported from morphometric data 5 , allozyme data 43 and mtDNA plus microsatellite data 27 . The PCoA performed on ITS1 sequences agreed with AMOVA partitions that separated New World hakes into Atlantic, Austral, and Pacific groups. Those exploratory analyses also allocated morphotypes to either the Pacific or the Austral groups. The phylogenetic support for those groups tested between methods (Bayesian inference vs. ML) and markers (ITS1, cyt b) with non-parametric bootstrapping produced the best-supported Bayesian-based reconstruction on cyt b haplotypes. Such reconstruction was characterized by a well-defined supercluster split and the placement of M. bilinearis at the base of the New World supercluster (BEAST, Fig. 1b). However, the best Bayesian tree and the best ML tree were not fully congruent with previous studies, a handicap if congruence among marker reconstructions is a relevant asset, e.g. the better ML-value on cyt b does not grant any better species tree than that from ML-ITS1 because the topology of the latter is more congruent with previous studies. Such scenario highlights the insufficiency of the ML method on cyt b variation as compared to the Bayesian approach (but see 44 ). The Bayesian phylogeny performed on both genes with BEAST was the unique method recovering both, a well-supported supercluster split and a full intracluster resolution, including the full definition of the Austral cluster. Discrepancies between markers using BEAST consisted on both, the branching of M. hubbsi as the oldest taxon among New World hakes with ITS1 and the branching of M. albidus within the Austral cluster with cyt b. Nonetheless, caution is needed on rejecting those odd positions regarding previous studies, since they could result from differential evolutionary rates among the markers applied.
It is expected that deep phylogenetic rooting of ancient hake lineages could be better afforded from conservative mtDNA haplotypes 45 . However, more recent evolutionary processes such as hybridization and drift could be better unveiled by highly recombinant nuclear DNA markers. Therefore, speciation histories based on mtDNA alone can be extensively misleading and large phylogenetic discrepancies have been reported between nuclear DNA and mitochondrial DNA 46 . Successful reconstructions have been achieved on concatenated sequences of the ITS region and COI 47 or on nuclear genes 48 . Also, current data showed that concatenated data from ITS1-rDNA and cyt b at reconstructing hakes phylogeny using Bayesian inference have dramatically improved  Fig. 2a). Those unexpected scenarios in view of previous studies can be explained as a synergical advantage of data concatenation as genes add up to produce a more balanced signal of the interspecific evolutionary divergence than that afforded from single genes 47 . Since each gene responds to distinct evolutionary dynamics it is expected that concatenated gene reconstructions approach the average evolutionary signal of their common history 50 . The major phylogenetic split in Merluccius comprises two monophyletic superclusters, the Old World one comprises five species and the New World one comprises six species, as shown with parasites 5,8,51 , morphology and meristic traits 6 , allozymes 14,52 and nucDNA 33 . This genus is believed to emerge in the Cretaceous after the opening of the southern Atlantic Ocean basin between South America and Africa 7 . However, its evolutionary bifurcation into New World and Old World superclusters is believed to have begun in the Oligocene 53 Table 1 and are followed by the alphanumeric entry code for each specimen in the authors laboratory, e.g. aust_ma01_3 is the specimen No. 3 of sample ma01 from M. australis (aust).   Table 1 and are followed by the alphanumeric entry code for each specimen in the authors laboratory, e.g. aust_ma01_3 is the specimen No. 3 of sample ma01 from M. australis (aust). Codes for morphotypes are given in Table 2 and are followed by an alphanumeric entry from the authors laboratories (Supp. www.nature.com/scientificreports/ benthonic species closer to the New World supercluster, they have been proposed both, either to descend from an early Old World Merluccius cluster 52 or to represent a distinct speciation process in the Eastern Atlantic 59 . The variance partition, the multivariate approach, and the phylogenetic inference, all them support three clusters within the New World supercluster, namely the Atlantic, the Austral and the Pacific. In the Atlantic cluster, the evolutionary status of M. bilinearis is disaggregated in classical hake phylogenies. Generally, this species has been placed as ancient to New World species with various methods and markers, including likelihood 52 and parsimony 14 on allozymes; NJ, ML, UPGMA on the D-loop 15 and current Bayesian analyses on cyt b haplotypes. However, new scenarios appear after concatenated ITS1-cyt b sequences which indefectibly place M. bilinearis as the oldest among extant hake species (see Fig. 2).
Most previous studies agree on that hake originated in the North Atlantic and entered the Pacific, but disagree on the origin of M. hubbsi, i.e. as diverted either from an eastern south Pacific stock 6,60 or from a western north Atlantic stock 5,8 . A novelty in the current concatenated approach is a well-supported Austral cluster that is congruent with parasite studies in which M. hubbsi and M. australis are closely related taxa of Austral origin 6 43 . Recently, it has been suggested that M. gayi 68 could be the single pan-Pacific species 28 although the North-Eastern Pacific hake has to be named Merluccius productus by priority 28 . Despite being considered by FAO as a variant of M. angustimanus 3 and reproductively isolated from M. productus and M. angustimanus 14 , Merluccius hernandezi 22 from Sinaloa state would likely be a dwarf morphotype of M. productus 28 . After its first description as a "tropical deep-water species off the Pacific coast of Central America from Mexico to Panama which is sandwiched between M. productus to the north and M. gayi to the south" 14 , M. angustimanus has been a more recognized taxon than M. hernandezi 69 . In the first genetic description of M. angustimanus 33 it was named M. hernandezi provided its northern California origin. Following recovering of its trawling data from the seventies 70 as well as the examination of its large scales (Mathews, personal communication) such sample was properly renamed as M. angustimanus 18 . That genetic analysis of M. angustimanus showed that the PCR-RFLP pattern on the ITS1 spacer was very similar to those of M. productus and M. gayi 18,33 . Also, current data showed that M. angustimanus shared its two unique ITS1 variants with M. productus. That result is congruent with the within species concerted evolution of the ITS1-rDNA gene family 37 and also with the recently proposed confinement and drift of M. productus in the Gulf of California giving rise to its present divergent population 27,27 . Present data also suggest the putative hybrid origin of M. angustimanus by means of a trans-equatorial incursion of M. gayi 5,6,14,22 in the territories of M. productus in the Gulf of California (see the bootstrap weakness within the Pacific cluster in all concatenated reconstructions, Fig. 2). Moreover, the six specimens of M. angustimanus examined had no evolutionary novelties in their cyt b region but rather were a chimeric genome from extant neighboring species plus an A-139 residue shared with M. albidus (see Supp. Table S4). In summary, M. angustimanus seems to be a hake population trapped in the northern Gulf of California 27 whose origin dates back some 0.25 MYA either from hybridization between M. productus and M. gayi or from a confinement of a subset of those species therein. Such hypotheses would explain the weak morphological and molecular differentiation of M. angustimanus from its neighboring hake species 22 . Whether that confined population can be self-sustainable as a relic of an ancient hybridization or if its viability depends on ongoing genetic contribution from the surrounding M. productus, are unexplored questions.
The morphotypes under test other than M. angustimanus were analyzed with the nested fragment ITS1Nes which offered a more conservative view of the interspecific diversity than the full ITS1 sequence, i.e. 11 HakeITS1Nes variants versus 19 HakeITS1 variants. In addition to morphotype-specific variants, all morphotypes shared the specific HakeITS1Nes.2 sequence from Pacific hakes, what suggests their common origin. Such commonality does not preclude the presence of specific sequences within the ITS1 spacer family that trace back the origin of their carriers. For instance, although some ITS1Nes variants of a given morphotype fully grouped into the M. hubbsi clade (e.g. HUBBM131) other intraindividual variants branched intermediately between the Pacific cluster and the Austral cluster. That is a typical scenario for interspecific hybrids which can be ascertained by inspection of the dramatic decrease of the basal tree bootstrap support of the Austral and the Pacific clusters in the phylogenetic backbone (see Figs. 5 and 6). The significant No. of recombination events detected in those sequences also comes to support the working hypothesis of a hybrid origin for those morphotypes.  Table S9). The Atlantic cod (Gadus morhua) was used as a phylogenetic outgroup. The world-wide sampling was performed on board of factory ships and research vessels between 2000 and 2003 33 when 9-277 specimens per species were preserved in 95% ethanol or frozen upon collection and their GPS recorded on board (Table 1). All specimens were identified using speciesspecific morphological traits such as otoliths, shape of abdominal vertebrae (parapophysis), cranial shape and pectoral fin length, following classification criteria previously established 6 . Whole specimens were boiled to recover those structures and to facilitate bone cleaning. Six specimens of M. angustimanus were sampled as preserved frozen from a bottom trawl campaign (80-523 m in depth) carried out in 1973 by Instituto Nacional de la Pesca, Mexico 70 . Fourteen additional specimens as holotypes, paratypes and rare morphotypes were either sampled in situ by the authors or taken from museum specimens as preserved in formalin. Two specimens of M. tasmanicus (DM5566 and DM3963) and four specimens of M. patagonicus were sampled and described in 23,24 , respectively. The paratype USNM 157765 of M. polylepis from Puerto Montt (Chile) was described by 23 and the M. polylepis holotype off Chiloe (Chilean Pacific) was described by 32 . Three rare morphotypes of M. hubbsi from the Beagle Channel (IIPB92/1987) and from Puerto Madryn were described by 23 . Finally, two rare morphotypes of M. australis from New Zealand (MOVI 27492-27493, formerly NMNZ P.13122) were described by 24 and two additional ones from the Chilean Pacific (off the Aysén Region) were described in 23 .

Molecular data.
Genomic DNA was extracted with FENOSALT 71 including a preliminary 24 h hydration step for formalin-fixed samples. Two synergic tools for hakes identification were applied, one based on the ITS1-rDNA spacer 33 and another based on the mtDNA cyt b gene 18 . Those targets are suitable DNA regions for the taxonomic classification of closely related taxa and their combinatory power has proven to be a robust approach to reduce the authentication error from hake-based commercial products 72 . Amplification conditions for the ITS1-rDNA spacer followed previous developments in this genus 33,73 . Electropherograms were revised with Chromas software (Technelysium, Tewantin, Australia). A total of 85 specimens were sequenced as averaging ~ 6 specimens per species (Supp. www.nature.com/scientificreports/ using DnaSP v5.10.1 77 . The average number of nucleotide substitutions per site (D xy ) and estimates of the Net Evolutionary Divergence between species (d) were compared among taxa using the confidence interval calculated as CI = − x±1.96*(S.D./ √ n ). The Principal Coordinates Analysis (PCoA) built with GenAlEx v6.503 78 on the molecular variation of ITS1 variants was used to check for the correct divergence signal between superclusters (New World and Old World) as well as between clusters and between species. Genetic variation within and between clusters were contrasted with an Analysis of Molecular Variance (AMOVA) as implemented in Arlequin v3.5.2.2 79 .
Comparative phylogeny and synthetic phylogeny. The phylogenetic reconstructions using the ITS1 region, the ITS1Nes fragment and the cyt b region included reference samples used at establishing the phylogenetic backbone of the genus 18,33 . The phylogenetic test enforced to allocate morphotypes was based on the phylogenetic backbone of New World hakes built with the ITS1Nes fragment. The most suitable nucleotide substitution models were selected using jModeltest v2.1.10 80 . The HKY85 + I + G model 81 was the best fitted to the ITS1-rDNA variation and the GTR + G model 82 was the best fitted to that of the cyt b gene. Likelihood computations of nonsynonymous (dN) and synonymous (dS) substitution rates were conducted using HyPhy software package 35 . The phylogenetic reconstruction was carried out using full sequences (ITS1, ITS1Nes and cyt b) as well as sequence variants (ITS1 and ITS1Nes) and cyt b haplotypes, as obtained with program DnaSP v5.10.1 77 . That software was also used to estimate the Recombination parameter (R) 34 as (R = 4Nr) where N is the population size and r is the recombination rate per sequence. The rate R allowed to assess whether recombination fingerprints can be detected among ITS1Nes sequences as expected in highly recombinant nuclear DNA markers 37 .
Character-based phylogenies were first constructed using Maximum Likelihood (ML) hypothesis testing after PAUP v4.0 83 , a gamma distribution coefficient (gamma-ITS1 or gamma-ITS1Nes = 0.269; gamma-cyt b = 0.420), a transition/transversion rate (T i /T v -ITS1 and T i /T v -ITS1Nes = 1.761; T i /T v -cyt b = 7.143) and the nucleotide substitution rate matrix observed for each marker. An algorithmic Neighbor-joining tree was built on ITS1Nes using PHYLIP v3.6 upon 5000 bootstrap replicates 84 . Bayesian phylogenetic hypotheses on ITS1 variants as well as on cyt b haplotypes were tested using the programs BEAST v1.8.4 85 and MRBAYES v3.2.6 86 and taking into account the best fitted evolutionary models. The comparative phylogeny consisted on contrasting tree support values and topological consistence per marker type (ITS1 and cyt b) among reconstruction methods (ML and Bayesian). The synthetic phylogeny of the genus consisted on a ML reconstruction made upon concatenated data sets of both genes using IQ-TREE (http://iqtre e.cibiv .univi e.ac.at/ 87 ) which allows partitioning the analysis to simultaneously conjugate two evolutionary models and also enforces the Shimodaira-Hasegawa test (SH) for testing gene trees. The null hypothesis (H 0 ) stated that all the trees tested would be equally performant at explaining the observed data, while the alternate hypothesis (H 1 ) stated that only one among several trees was a better proxy to real data. The Bayesian models implemented in MRBAYES and BEAST were employed in the phylogenetic reconstruction of concatenated sequences. The Net Evolutionary Divergence (d) among sequences (ITS1 variants, ITS1Nes variants, cyt b haplotypes) was calculated upon divergence between OTUs (species, clusters and superclusters) as (P AB ) and corrected for within-OTU variation (P A and P B ) as: www.nature.com/scientificreports/ 87. Trifinopoulos, J., Nguyen, L. T., von Haeseler, A. & Minh, B. Q. W-IQ-TREE: a fast online phylogenetic tool for maximum likelihood analysis. Nucleic Acids Res. 44(W1), W232-W235. https ://doi.org/10.1093/nar/gkw25 6 (2016).