Morphometric and molecular differentiation of Pimelodus grosskopfii and Pimelodus yuma (Siluriformes: Pimelodidae)

Abstract Pimelodus grosskopfii and Pimelodus yuma, two species endemic to the Magdalena-Cauca basin in Colombia, overlap in the ranges of some of their diagnostic characters, which hampers their correct morphological identification. Aiming to help discriminate these species, this study conducted an integrative analysis using traditional and geometric morphometrics, phylogenetic analysis based on partial sequences of the mitochondrial cytochrome c oxidase subunit 1 gene (COI, cox1) and the identification of diagnostic Single Nucleotide Polymorphism markers (SNP). The species differ significantly in body geometry, allowing 100% discrimination, which was reinforced by a phylogenetic analysis that recovered well-supported monophyly of each species (posterior probability > 0.95). Additionally, the traditional morphometric results corroborated some previously reported diagnostic traits for both species and let us describe one non-overlapping ratio related to the adipose fin length. Three of five SNP markers had reciprocally exclusive alleles suitable for identifying each species. The morphometric and molecular methods conducted in this study constitute alternative tools for the correct discrimination of P. grosskopfii and P. yuma in the wild and in captive populations used for aquaculture.


INTRODUCTION
Pimelodidae is one of the most diverse Neotropical families within the order Siluriformes, comprising 30 genera and 116 species (Fricke et al., 2022a) distributed from Panama to South America (Nelson et al., 2016).Members of this family are economically important and highly demanded in commercial fisheries and in the ornamental fish trade (Cala-Cala, 2019).Within Pimelodidae, Pimelodus Lacepède, 1803 includes 36 valid species (Fricke et al., 2022b), three of which are endemic to the colombian Magdalena-Cauca basin: 1) Pimelodus crypticus Villa-Navarro & Cala, 2017, restricted to the upper Cauca River drainage; 2) Pimelodus yuma Villa-Navarro & Acero, 2017, distributed in the Magdalena, Cauca and Sinú river drainages, and 3) Pimelodus grosskopfii Steindachner, 1879 which inhabits the Magdalena-Cauca basin.Pimelodus crypticus and P. yuma are recently recognized species that had been subsumed within the concept of Pimelodus blochii Valenciennes, 1840(Villa-Navarro et al., 2017).
Pimelodus grosskopfii and P. yuma occur sympatrically in the lower Cauca and Magdalena River basins, where they migrate up to 500 km (Zapata, Usma, 2013).Both species reproduce externally and are omnivorous with insectivorous preferences (Lasso et al., 2011;Ramírez, Pinilla, 2012), although Jiménez-Segura et al. (2020) considered P. yuma to be carnivorous.Both congeners are among the top five species of fishery interest (Hernández-Barrero et al., 2020) and exhibit two genetic stocks in the lower section of the Cauca River (Joya et al., 2021;Restrepo-Escobar et al., 2021).The International Union for Conservation of Nature -IUCN lists P. grosskopfii as Critically Endangered (CR), whereas P. yuma has not yet been included or evaluated (Villa-Navarro et al., 2016).
Overlapping ranges in some diagnostic characters hamper the correct morphological identification of P. grosskopfii and P. yuma.Additionally, phenotypic plasticity influences some traits that putatively discriminate P. grosskopfii from P. yuma, which could hinder the correct identification of both species.The potentially plastic traits include the dark spots along the entire length of the body, longer adipose fin, and shorter distance between the dorsal and adipose fins in P. grosskopfii (Villa-Navarro et al., 2017).For instance, body coloration in fishes can change in response to light (Ninwichian et al., 2018;Pinto et al., 2020) or other external environmental stimuli (Sugimoto, 2002;Van Der Salm et al., 2004), or to enhance reproduction or survival (Rodgers et al., 2010).Biotic and abiotic conditions and metabolic requirements can affect the size and presence of the adipose fin in some catfish species (Reimchen, Temple, 2004;Temple, Reimchen, 2008).Variation in water velocity is known to induce body shape variation in P. grosskopfii (Hincapié-Cruz, Márquez, 2021).Thus, the ability of morphological characteristics to diagnose these similar species may be questioned.
In this context, the present study assessed the utility of four approaches to discriminate between the catfish species P. grosskopfii and P. yuma: traditional and geometric morphometrics, phylogenetic analysis of a portion of the mitochondrial cox1 gene, and identification of diagnostic single nucleotide polymorphisms (SNPs) in the nuclear genome.Traditional morphometrics quantify a series of morphological characters (length, depth, and width; or their expression as means, ranges, or ratios) and use multivariate methods that allow describing patterns of variation (Marcus, 1990;Rohlf, Marcus, 1993).Together with meristic characters (e.g., Londoño-Burbano et al., 2011;Vanegas-Ríos et al., 2011;Arcila et al., 2013;Ferrer, Malabarba, 2013), and geometric morphometric data (e.g., Maderbacher et al., 2008;Sidlauskas et al., 2011;Gupta et al., 2018;Garavello et al., 2021), the morphological characters are widely used for identifying neotropical fish species.Geometric morphometrics allows the analysis of biological shape variation and its covariance with other variables and thus can find patterns of morphological variation among and within species (Bookstein, 1992).In fishes, this method can discriminate morphologically similar species (García-Alzate et al., 2011), assess sexual dimorphism (Herler et al., 2010), test the covariation of shape with environmental conditions (Restrepo-Escobar et al., 2016a,b;Hincapié-Cruz, Márquez, 2021) and support studies of genetic differentiation (Kocovsky et al., 2013).
Partial sequences of the mitochondrial gene cox1 are often used as a barcode to differentiate animal species (Hebert et al., 2003) due to their wide availability in public databases (Heller et al., 2018) and other advantages of mitochondrial DNA (mtDNA) such as the simple structure, high copy number, uniparental inheritance, and rapid evolutionary rate (Rubinoff, Holland, 2005).This gene has successfully discriminated freshwater fish species from Canada (Hubert et al., 2008), Central America (Valdez-Moreno et al., 2009), andTaiwan (Bingpeng et al., 2008).It has also performed well in several studies of Pimelodus in South America (Rosso et al., 2012;Becker et al., 2015;Frantine-Silva et al., 2015;Díaz et al., 2016;Guimarães-Costa et al., 2019).
However, cox1 sequences do not diagnose all species reliably; given the uniparental inheritance of mtDNA, this gene cannot elucidate hybridization or introgression events (Hubert et al., 2008).Thus, this study also explored the utility of SNPs in nuclear genes Pimelodus grosskopfii and P. yuma discrimination by taking advantage of a RAD-seq library of previously cox1 barcoded individuals of Pimelodus (Martínez et al., 2022).Although closely related species will not have diverged at all genomic regions (Duranton et al., 2018), the occurrence of genomic islands, intraspecific selection (Duranton et al., 2018), and random fixation of alleles by mutation-drift equilibrium (Martin, Jiggins, 2017) usually result in rapid divergence and fixation of diagnostic alleles at some loci.Such exclusive regions are the target for the development of SNPs markers employed to discriminate P. grosskopfii from P. yuma.
The present study explores the utility of morphometric and molecular tools to discriminate these catfish species.Reliable diagnoses would aid in identifying individuals from wild or farmed populations to the correct species and help to scaffold further interspecific genetic analyses.

MATERIAL AND METHODS
Sampling.A total of 159 cox1 barcoded individuals of the genus Pimelodus from the Magdalena and Cauca rivers were analysed (Fig. 1, Tabs.S1 and S2).A subgroup of 125 samples from the middle and lower sector of the Cauca River and lower sector of the Magdalena River (Tab.S1) was evaluated by phylogenetic and morphometric analyses, while the remaining 34 samples from a previously sequenced RAD-seq library (Martínez et al., 2022) were used for genome association studies.This last subgroup also included samples from the upper sectors of Cauca and Magdalena rivers (Tab.S2).
Traditional and geometric morphometrics.A total of 125 individuals (79 females and 46 males) were photographed in left lateral view within a diffuse light box at 0.47 m focal length using a dual pixel 12 MP Samsung digital camera.The sex identification was based on the usual method of applying external abdominal pressure to obtain gametes in ripe adults and visual inspection of the urogenital papillae (Ross, 1984).Digitization of landmarks was conducted using the online application XYOM v. 2.0 (Dujardin, Dujardin, 2019; https://xyom.io/).For each image, nine landmarks were digitalized (Fig. 2): eight type I (juxtaposition of tissues) and one type II (point of maximum curvature in the head; Bookstein, 1992).
To ensure maximum precision in landmark assignment, digitization was repeated twice, and the repeatability index was quantified using XYOM v. 2.0.Centroid size (CS) was calculated based on landmark coordinates (x, y), and then a Generalized Procrustes Analysis was conducted to eliminate shape variation resulting from differences in size, position, and rotation.Finally, the partial (partial warps) and global (uniform components) shape deformations of individuals with respect to the consensus (Goodall, 1991;Warheit et al., 1992) were used as shape variables.
Differences in body shape between species and sexes were explored using a Principal Components Analysis (PCA), and the statistical significance of Euclidean distances among groups was obtained by 1,000 permutations.In addition, the allometric effect was calculated using a multivariate regression between shape variables and centroid size.When significant allometry was detected, a size correction of shape was conducted if allometric slopes fit a common model, which was assessed using a Multivariate Analysis of Covariance (MANCOVA).Additionally, the diagnosability of each group suggested  by PCA was examined via assignment tests using Simple Reclassification and Verified Cross-Classification.All calculations were conducted in XYOM v. 2.0.
To assess differences in centroid size, we conducted a pairwise Tukey's test among compared groups and a two-way ANOVA including groups and sex as main factors.This analysis was performed after verifying that the centroid size met the assumptions of normality (Shapiro-Wilks test) and homoscedasticity (Levene's test with the Brown-Forsythe modification) using, respectively, the packages nortest v. 1.0.4 and stats v. 3.6.0, of the statistical program R (R Development Core Team, 2013).
Finally, we calculated distances between landmarks to obtain traditional morphological data from all 125 photographs and visualized differences among individuals by PCA.The numbers of pixels were converted into centimetres using a reference scale within each image, then transformed using the Neperian logarithm.Because standard length was not directly measured in our landmark set, we approximated that measure as the mean of the distances between the tip of the snout and the dorsal and ventral origins of the caudal fin.All other measurements were then expressed as percentage of the approximated standard length.The measurements that displayed a greater contribution to the total variation were evaluated by parametric (t tests, if assumptions of normality and homoscedasticity were met) or non-parametric tests (Welch t-test, Kruskal-Wallis, and Mann-Whitney test) using the stats package v. 3.6.3implemented in R. The ratio distance/length of the adipose fin was also calculated to explore whether that ratio can discriminate the species.To assess possible allometric variation in that ratio, the ratio of the distance/length of the adipose fin was regressed on standard length using standardized major axis regression in smatr v. 3.4-8 of the statistical program R (R Development Core Team, 2013).

Phylogenetic analyses.
DNA extraction was performed from 125 fin or muscle tissues previously preserved in 96% ethanol using the GeneJET Genomic DNA Purification kit, following the manufacturer's instructions.A partial region of the mitochondrial cox1 gene was amplified using a cocktail of four primers: VF2_t1 5'-TCAACCAACCACAAAGACATTGGCAC-3', FishF2_t1 5'-TCGACTAATCATAAAGATATCGGCAC-3', FishR2_t1  Hall, 1999).The software Geneious v. 11.0.4 was also used to translate consensus sequences into the corresponding amino acids to confirm the absence of stop codons or frame shift errors within the sequence.
The number of haplotypes, nucleotide diversity (π) and haplotype diversity (Hd) were calculated in DNAsp v. 6 (Rozas et al., 2017).The best fit model of molecular evolution was selected using jModelTest v.2.1.10 (Guindon, Gascuel, 2003;Darriba et al., 2012).Haplotypes obtained in this study were compared with known cox1 sequences from Pimelodus crypticus (OK206769; ON596099), P. grosskopfii (MK748281) and P. yuma (NC_061908).In addition, a previously published GenBank sequence of Pseudoplatystoma magdaleniatum Buitrago-Suárez & Burr, 2007 (KP090204) was used as outgroup.Phylogenetic analysis was conducted in BEAST v. 2.6.4 (Bouckaert et al., 2019) implemented on the CIPRES Science Gateway (available at http://www.phylo.org/).To this end, an xml file was created in BEAUti, setting HKY + G as the best-fit model, a strict molecular clock with a clock rate of 1 and a Yule model as a prior on the speciation process; remaining parameters were left at their default settings.Analyses were performed using 900 million Markov Monte Carlo Chains (MCMC) sampled every 90,000 generations, and discarding the initial 10% as burn-in.Convergence (ESS value > 200) was assessed using Tracer v. 1.7 (Rambaut et al., 2018) and a consensus tree was generated using Treeannotator v. 1.10.4.Finally, the obtained phylogenetic tree was visualized using Fig Tree (http://tree.bio.ed.ac.uk/).Following Hubert et al. (2008) recommendations for the DNA barcoding of freshwater fishes, the haplotypic divergences between and within species was estimated by calculating the Kimura 2 parameter model in Mega v. 11.0.13 (Tamura et al., 2021).For species discrimination based on PCR methods, two primer pairs were designed using Primer-Blast (which in turn uses Primer3, Untergasser et al., 2012) to design PCR primers and BLAST and global alignment algorithm to screen non-specific amplifications of primer pairs against a user-selected database (Ye et al., 2012).The same approach was used for designing primers for diagnostic SNP loci identified as described below.v. 3 (Ravindran et al., 2019) to demultiplexing by 16-nt barcodes, filter the samples by quality and construct de novo loci.The script used the following parameters: "-a" psweep mode for loci assembling optimization, Phred score threshold >20, minimum sequencing depth of 3 (-m), maximum allowed distance (in nucleotides) among stacks of 4 (-M), maximum allowed distance among putative loci of 1 (-n), maximum stacks per locus of 2 (-x) and minimum percentage of samples of 0.9 (-S).

Identification of diagnostic
The SNP output files for population analyses (STRUCTURE and VCF) were generated by STACKS v. 2.53 (Catchen et al., 2013) from the RADproc catalogues, running the modules SSTACKS, TSV2BAM, GSTACKS, and POPULATIONS in pipeline mode.In the STACKS map file, each individual was considered as a population; in the POPULATIONS module, the minimum number of populations was set to 100%, the minimum percentage of individuals in a population required to process a locus for that population was set to 90% (-r = 0.9), with a minor allele frequency (MAF) of 0.05, and remaining parameters were left by default.
To corroborate the existence of P. yuma and P. grosskopfii as biological entities for subsequent population analyses, a coancestry analysis was performed in STRUCTURE v. 2.3.4 (Pritchard et al., 2000) using 20 independent runs from K = 1 to K = 4, 1,000,000 Markov Monte Carlo Chains (MCMC) and a burn-in of 100,000 iterations; STRUCTURESELECTOR (Li, Liu, 2018) was used to calculate the best estimate of K.A summary of structure runs was plotted on a coancestry histogram of all individuals using the integrated software CLUMPAK (Kopelman et al., 2015).
Once the biological assignation of the 34 individuals was corroborated, an interspecific analysis was performed in STACKS (POPULATION module) to generate SNPs statistics between the two species.Specifically, we obtained LOD scores (likelihood odds ratio logarithm > 3), which estimates linkage disequilibrium between a SNP and the biological entities, and the Fisher's probability (P) of SNPs fixation index (Fst).Finally, the package "qqman" v. 0.1.4(https://CRAN.Rproject.org/package=qqman) in the R software environment v. 3.4.1 (R Development Core Team, 2014) was used to construct a Manhattan plot using Fisher's P values to identify SNP alleles significantly associated with each species.In addition, SNPs identity was assessed using the BLAST algorithm in Genbank (https://www.ncbi.nlm.nih.gov/genbank/), using the genome of Ictalurus punctatus (Rafinesque, 1818) (Genbank accession number: LBML00000000) as a reference.
The Simple Assignment analysis (Pg: 60/60, 100%; Py: 65/65, 100%) and Verified Cross-Classification (Pg: 60/60, 100%; Py: 65/65, 100%) confirmed that individuals separated by the first principal component in the PCA belong to the differentiated groups.That result demonstrates the existence of two distinct groups corresponding to the two species.That major difference does not result from sexual dimorphism.
Average body shape comparison of each group (Fig. 4) revealed allometric differences showing that the body of P. grosskopfii tends to be more streamlined and hydrodynamic than P. yuma, which seems to have a slightly taller and more robust body shape.In addition, within P. yuma, females (PyF) tend to be slightly taller than males (PyM), while in P. grosskopfii these differences were not significant.Regarding size, statistically significant differences were found between species (F 1, 122 = 163.632,P < 0.000), however, there were no significant effects between sexes (F 1, 122 = 0.361, P = 0.549) nor in the Species x Sex interaction (F 1, 121 = 2.564, P = 0.112).Specifically, P. grosskopfii males and females were larger than those of P. yuma (Fig. 5).Pairwise comparisons between PyF-PgF (P = 0.000), PyM-PgF (P = 0.000), PyF-PgM (P = 0.000), and PyM-PgM (P = 0.000), showed significant differences between the means of each group.
Large individuals of these species were found to differ in body pigmentation.Specifically, larger individuals of P. grosskopfii possess separated and well-defined black spots distributed towards the anterior end or over the whole body.Smaller individuals of P. grosskopfii resemble some individuals of P. yuma by displaying spots with smaller diameters than those present in large individuals of P. grosskopfii group (Fig. 8).

Pimelodus grosskopfii
Pimelodus yuma    Pimelodus grosskopfii and P. yuma discrimination Phylogenetic analyses.After sequence edition, we obtained an alignment 489 bp in length of the cox1 gene (GenBank accession numbers: ON596020 -ON596098; ON596100 -ON596145).Translation into amino acids revealed no indels or stop codons, showing that the alignment corresponds to functional cox1 gene sequences.A total of 17 haplotypes of P. grosskopffi made up two haplogroups that clustered within the reference P. grosskopfii clade (posterior probability = 1), while 11 haplotypes of P. yuma females made up two haplogroups that clustered with reference P. yuma (posterior probability = 1), showing a 100% concordance between with morphometric assignation of individuals to species (Fig. 9).

DISCUSSION
This study conducted an integrative analysis combining traditional and geometric morphometrics, phylogenetic analysis using cox1 gene sequences, and identification of diagnostic SNP markers to improve the discrimination of the catfish species Pimelodus grosskopfii and P. yuma.Our results show that traditional and geometric morphometrics can reliably discriminate these congeners by body shape, as could phylogenetic analysis of information contained in cox1 sequence data.For the first time, genomic SNP markers with exclusive alleles were found to also discriminate these species.
Traditional, geometric morphometrics and phylogenetic analyses.Principal components analyses based on geometric and traditional data revealed the existence of two groups with differences in body shape.Assignment analyses classified all individuals to one of these groups with high confidence.These results were corroborated by cox1 gene sequence analyses, which accurately clustered the haplotypes belonging to each morphological group with those of known sequences of P. grosskopfii and P. yuma.The phylogenetic approach also corroborates that the species described by Villa-Navarro et al., (2017) in the Magdalena-Cauca basin represent monophyletic groups clearly distinguishable by cox1 sequences (Martínez et al., 2022), as evidenced by the high posterior probability values that supports the observed clustering.Based on the divergence of cox1 sequences between species, two pairs of primers were designed to discriminate them, although its experimental evaluation and consequent usefulness remains to be explored.
Although two haplogroups were found in each species, the Kimura 2 parameter genetic distances between those haplogroups are less than 1%, suggesting that they do not represent cryptic species (Ward et al., 2009).However, further integrative studies would be required to fully test a hypothesis of further unrecognized diversity in this clade.The observed differences in average body shape of females and males established that P. grosskopfii has a more streamlined body, and hence may be more hydrodynamic than P. yuma.These results accord with the different environments that these species inhabit.Pimelodus yuma is restricted to the lower sector of the Cauca River (Villa-Navarro et al., 2017) and experiences relatively stable hydrodynamic conditions, while P. grosskopfii migrates to the upper sector (Restrepo-Escobar et al., 2021), and experiences many different flow regimes.It was previously reported that differences between lotic and lentic habitats can induce body shape differences in this species (Hincapié-Cruz, Márquez, 2021).
Exploring intrinsic sources of phenotypic variation, this study detected no sexual dimorphism in body shape or size in P. grosskopfii, which concords with Cala (1997), though Valbuena-Villarreal et al. (2010) suggested dimorphism in size for this species.In contrast, we detected differences in body shape between sexes of P. yuma, especially in body height between the origin of the dorsal fin and the origin of the adipose fin and at the origin of the anal fin.Further studies should explore these results in individuals with similar states of gonadal development (Vazzoler, 1996), since reproductive maturity and gravidity may affect the body shape of the individuals analysed.
Previous traditional morphometric analysis suggested that the size of the adipose fin and the predorsal distance are the measurements that best differentiate P. grosskopfii from P. yuma (Villa-Navarro et al., 2017).However, we found the range of these variables to overlap in the two species, prompting a search for additional traits that allow discrimination.Herein we examined four measurements related to the length of the adipose fin, which differ from those reported by Villa-Navarro et al. (2017).The ranges of three of these measurements also overlap among the species, but the ranges of the ratio L3-L8/L3-L4 do not.This ratio juxtaposes the distance between the adiposefin origin and pelvic-fin origin with the length of the adipose-fin base and provides a measurement that allows easier differentiation of these species.Likewise, the differences in the allometric trajectories of P. grosskopfii and P. yuma represent an additional point of discriminating both species; and their visualization in the same plot (Fig. 7) provides a reference for the diagnosis of those species.
Additionally, this study found variations in body pigmentation patterns in both species.Although all the large-sized individuals of P. grosskopfii showed dark spots on their bodies, smaller individuals showed smaller spots like those seen in some P. yuma specimens, which shows that the skin pigmentation is a variable trait and its validity as a distinctive characteristic is limited.Intraspecific variations in the skin coloration pattern have been described for species of Pseudoplatystoma Bleeker, 1862(Buitrago-Suárez, Burr, 2007;Scarabotti et al., 2020), Pseudopimelodus Bleeker, 1858(Restrepo-Gómez et al., 2020), Pagrus pagrus (Linnaeus, 1758) (Van Der Salm et al., 2004), Trichogaster pectoralis Regan, 1910(Ninwichian et al., 2018), among others (Sugimoto, 2002).
Identification of diagnostic SNPs.This study identified for first time a set of diagnostic SNPs that can be used to discriminate P. grosskopfii from P. yuma, which were previously separated genetically with only mitochondrial cox1 sequences.Five loci exhibited strong differential association with each species, and three of these loci had exclusive alleles in homozygous state in both species.Therefore, 2.5% of the total 120 SNPs perfectly separate these species.A study of four tilapia species found a slightly lower detection efficiency of 1.8% in a panel of 1371 SNPs obtained from partial genome sequencing (ddRADseqs) (Syaifudin et al., 2019).Although SNPs allowed the differentiation by species, we recommended a further validation study using more samples to assess the efficacy of these markers in discriminating both species via PCRbased methods.For this proposal, the primers designed for the loci 5961 and 12198 are highly recommended, whereas the primers for the locus 1997 may not be specific, showing potential in silico amplifications with other fish genera.Since the SNP is in the 3`end of the forward primer in all loci, stringent PCR conditions are required to avoid non-specific amplifications.Due to the forward primer for the locus 12198 is short (12 mer in length), its length may be increased using adaptors attached on its 5`end such as those described by Blacket et al. (2012) or the M13 tail.
Whether the differentially expressed genes influence the fitness and distribution of these species is currently unknown.Given the wider distribution of P. grosskopfii, it is expected that this species has a better ability to face different environments, suggesting that a SNP within a gene related to the development of skeletal and cardiac muscle could be associated with its distribution.Scott, Johnston (2012) showed that embryonic temperature in zebrafish can have dramatic and persistent effects on acclimatization ability at multiple levels of biological organization.Performance differences following temperature acclimatization from cold to hot were partially explained by fiber type composition in swimming muscles.In addition, differences were observed in the expression of genes involved in energy metabolism, angiogenesis, cell stress, contraction and muscle remodelling, and apoptosis.Although the Claudin-4 gene has a pivotal role in salinity adaptation in the flounder species Paralichthys lethostigma Jordan & Gilbert, 1884(Tipsmark et al., 2008), it does not seem to provide differential tolerance to physicochemical changes (atmospheric pressure, dissolved oxygen, temperature, total solutes, and pH) in P. grosskopfii and P. yuma, as synonymous mutations were found in both species.However, since the effect of genetic drift may also explain our results, these hypotheses must be assessed in more individuals and with suitable experimental designs that allow to understand the role of these genetic polymorphisms in the physiology of these species.
In line with the above-mentioned statements, the SNP with non-synonymous mutations 9294 (WNK1) exhibited high genotypic frequencies for the "AA" genotype in P. grosskopfii from upper and middle sections of the Cauca and Magdalena rivers (AA = 0.833; AC = 0.125; CC = 0.042), whereas the "CC" genotypes were found exclusively in the lower section of the Cauca River.Although samples from the lower section of the Magdalena River remains to be examined and more studies are needed, these genotypic frequencies seem to be a marker for the upper/middle (Cauca River: AA; Magdalena River: AA > AC) and lower sector (Cauca River: CC) environments of the Magdalena and Cauca rivers.This outcome differs from that observed in P. yuma, where the genotype "CC" found in all assessed individuals (n = 10) was present in the sectors evaluated of the Cauca and Magdalena rivers.
In conclusion, this study demonstrates that four methodological approaches can discriminate the congeners P. grosskopfii and P. yuma: traditional and geometric morphometrics, phylogenetic analysis of the partial sequence of the cox1 gene and SNPs.These approaches, which can be used separately or in combination by researchers with differing areas of expertise, support the recognition of these species as valid and can help to assign individuals from wild and captive populations to the correct species.That in turn, provides a crucial step for the preservation of these animals and proper management of the fishery that depends upon them.

FIGURE 1 |
FIGURE 1 | Sampling sites of Pimelodus grosskopfii and P. yuma in the Cauca and Magdalena rivers.

FIGURE 6 |
FIGURE 6 | Plot of the two first principal components (PC1: 97.2%, PC2: 1.2%) based on traditional measurements of Pimelodus grosskopfii (Pg) and P. yuma (Py).PC3 (0.6%) does not help to discriminate the species and may represent minor measurement error.

FIGURE 9 |
FIGURE 9 | Bayesian phylogenetic tree based on partial cox1 gene sequences, showing the phylogenetic relationships of individuals analysed in this study with the barcode-confirmed species Pimelodus grosskopfii, P. crypticus, and P. yuma.

FIGURE 10 |
FIGURE 10 | STRUCTURE results based on unlinked SNP analysis for Pimelodus grosskopfii and P. yuma.

TABLE 1 |
Principal components (PC) loading from traditional measurements.*Majorcontribution to the species discrimination.

TABLE 2 |
Linear measurements of Pimelodus grosskopfii and P. yuma expressed as percentages of standard length or of the length of the adipose fin.SD = Standard deviation.