Iterative taxonomic study of Pareiorhaphis hystrix (Siluriformes, Loricariidae) suggests a single, yet phenotypically variable, species in south Brazil

Pareiorhaphis hystrix is a widely distributed species, occurring in the upper and middle Uruguay River and in the Taquari River basin, Patos Lagoon system, southern Brazil. Morphological variation has been detected throughout the distribution of P. hystrix, and this work seeks to test the conspecific nature of populations in several occurrence areas. Specimens from six areas in the Uruguay River basin and three in the Taquari River basin were compared. Variance analysis (ANOVA) was performed for the meristic data, and Principal Component Analysis (PCA) and Linear Discriminant Analysis (LDA) were conducted for morphometric data. Molecular analyses used coI, cytb, 12S and 16S mitochondrial genes, examining nucleotide diversity, haplotype diversity, genetic distance, and delimitation of possible multiple species through the Generalized Mixed Yule Coalescent (GMYC) method. Phylogenetic relationships of studied populations were also investigated through Bayesian inference. While PCA indicated a tendency of overlap between areas, ANOVA and LDA detected a subtle differentiation between populations from the two hydrographic basins. Yet, both latter analyses recovered the population from Pelotas River, a tributary to Uruguay River, as more similar to populations from Taquari River, which is congruent to morphological observations of anterior abdominal plates. The molecular data indicated a nucleotide diversity lower than the haplotypic diversity, suggestive of recent expansion. The concatenated haplotype network points to slight differentiation between areas, with each locality presenting unique and non-shared haplotypes, although with few mutational steps in general. The species delimitation by coalescence analysis suggested the presence of a variable number of OTUs depending on the inclusion or exclusion of an outgroup. In general, the morphological data suggest a subtle variation by river basin, while the genetic data indicates a weak population structuration by hydrographic areas, especially the Chapecó and Passo Fundo rivers. However, there is still not enough differentiation between the specimens to suggest multiple species. The iterative analyses indicate that Pareiorhaphis hystrix is composed of a single, although variable, species.


Introduction
Among catfish families, Loricariidae is the most diverse with above 1,000 species currently valid [1], a number that continues to raise steadily. Fishes in this family, popularly known as Cascudos in Brazil, are widely distributed in neotropical freshwaters from the Pacific drainages of southern Costa Rica to northeastern Argentina, with six subfamilies currently recognized [1,2].
Pareiorhaphis Miranda Ribeiro, 1918 was removed from the synonymy of Hemipsilichthys Eigenmann, 1889 by Pereira [3] to reflect the phylogenetic relationships of the Delturinae and Neoplecostominae [4]. Pareiorhaphis currently has 26 valid species [5] and is endemic to Brazil, occurring in main coastal river drainages of south, southeast and northeast Brazil, in addition to eastern versants tributaries to the Paraná and São Francisco rivers [5,6]. Species of Pareiorhaphis inhabit streams of strong water current and rocky bottom, usually being abundant where they occur, with greater diversity in the Doce River and coastal rivers of Santa Catarina State [5,7].
The genus was recently demonstrated to be monophyletic by Pereira & Reis [8], and its species show remarkable morphological variation in several aspects like body size (maximum size 34 mm SL in P. nudula (Reis & Pereira, 1999) to 116 mm SL in P. azygolechis (Pereira & Reis, 2002)), color pattern, secondary sexual dimorphism, nature of abdominal cover, head, snout, and lip shapes, and morphometric and meristic features [3,8]. Even the hypertrophied cheek odontodes of adult males, which represent a synapomorphy for the genus [8] are highly variable in size, thickness, density, direction, and position on the sides of the head.
Pareiorhaphis hystrix (Pereira & Reis, 2002) (Figs 1 and 2) was described almost two decades ago from the upper Pelotas River, Uruguay River basin, south Brazil, and later had its distribution expanded, being currently known from the upper and middle Uruguay River and the Taquari River basin, a tributary to the Patos Lagoon system. Its wide and disjoint distribution across two river basins diverges of most Pareiorhaphis species, which are usually restricted to a single drainage basin or to small basins historically connected in recent geological history.
Over the years, specimens of Pareiorhapis hystrix were collected in different environments in southern Brazil and subtle phenotypic variation among individuals has been noticed regarding color pattern, size and position of hypertrophied odontodes, and size of the cheek fleshy lobe. The geographic distribution of such variation has not been investigated, and it is thus not clear whether these represent intraspecific variation or could indicate the presence of multiple species. The combination of phenotypical data with molecular markers in an iterative approach, as opposed to the traditionally used morphology or molecules alone, has been efficient in fish population studies, and their use can contribute to the identification of cryptic species [9,10].
The objective of this study is to perform a comparison between populations of Pareiorhaphis hystrix in the Uruguay and Taquari river basins, southern Brazil, to test the existence of undetected cryptic species. The hypothesis tested is the co-specificity of P. hystrix in the different areas of occurrence, with the null hypothesis being that P. hystrix constitutes a single, widespread species, while the alternative hypothesis predicts that P. hystrix aggregates different cryptic and closely related species lineages.

Material and methods
The concept of iterative taxonomy is used in the sense of Yeates et al. [11], where species boundaries are treated as hypotheses to be tested by the comparison of multiple lines of evidence analyzed iteratively, in opposition to integrative taxonomy, where different sources of evidence should be integrated and analyzed together.
Biodiversidade (ICMBio) of the Ministry of Environment. Fishes were euthanized through an over-exposition to a solution of clove oil, according to the protocol of Lucena et al. [12].
Morphological observations were made searching for variations on cheek fleshy lobe, hypertrophied male odontodes on cheek and pectoral-fin spine, color pattern, and abdominal plates with the specimens submerged in alcohol. Further morphological analyses were performed on adult specimens (above 60 mm) preserved in 70% ethanol. Counts and measurements were obtained according to Pereira et al. [7] from 211 specimens (S1 Table) representing each of the nine populations/areas, with body measurements being expressed as percent of the standard length, while those of the cephalic region expressed as percent of the head length. The comparisons between the different areas were performed through statistical analyses of the morphometric and meristic values, presenting the minimum, maximum, mean, and standard deviation by area. The distribution map of Fig 3 was prepared using Quantum-GIS (v3.8), with shape and raster files available at IBGE (Instituto Brasileiro de Geografia e Estatística: http://mapas.ibge.gov.br/bases-e-referenciais), and ANA (Agência Nacional de Á guas: http://www.snirh.gov.br/hidroweb), following the map preparation tutorial of Calegari et al. [13].
Power regressions (y = a.SL b ), where SL stands for standard length, were applied to morphometric data to convert measurements into residues and eliminate bias of individual size and allometric growth. Conversion was performed by nonlinear regression in SPSS v22.0 software using the Levenberg-Marquardt algorithm [14], seeking to eliminate the log-transformation bias. The same residues were compared morphometrically by multivariate analysis, Principal Component Analysis (PCA), which finds the main axes of variation of a dataset in any dimension, not discriminating whether these data belong to the same class) and Linear  Discriminant analysis (LDA), which considers the existence of classes in the data, highlighting a linear separation in case it exists), using PAST v3.12 software [15]. Analysis of variance (ANOVA), box plots, and Tukey's tests were used to compare meristic data and were performed using PAST v3.12.

DNA extraction, amplification, sequencing and alignment
Total DNA extraction was performed with the DNeasy Blood & Tissue Extraction Kit (Qiagen, Valencia, CA, USA), for 47 specimens of Pareiorhaphis hystrix plus five outgroup specimens (S1 Table) following the manufacturer's protocol. DNA amplification was performed by Polymerase Chain Reaction (PCR) for the following mitochondrial genes: Cytochrome c oxidase subunit I (coI), Cytochrome b (cytb), 16S, and 12S. These fragments were amplified using previously published and available primers ( Table 2). The PCR protocol for coI included: initial denaturation at 94˚C for 2 min, followed by 40 cycles of denaturation at 94˚C for 20s (or 30s when using primers LCOI490 and HCO2198), annealing at 50˚C, 48˚C, 46˚C, 44˚C, 42˚C and 40˚C with 20s each end and 5s for intermediates, extension at 72˚C for 2 min and final extension at 72˚C for 10 min. For cytb: initial denaturation at 94˚C for 2 min, followed by 35 cycles of denaturation at 94˚C for 30s, annealing at 57˚C, 55˚C and 53˚C with 20s each, extension at 72˚C for 1.5 min and final extension at 72˚C for 10 min. For 16S: initial denaturation at 94˚C for 2 min, followed by 35 cycles of denaturation at 94˚C for 30s, annealing at 48 or 40˚C for 20s, extension at 72˚C for 60s and final extension at 72˚C for 10 min. Finally, for 12S: initial denaturation at 94˚C for 2 min, followed by 35 cycles of denaturation at 95˚C for 30s, annealing at 52˚C, 50˚C and 48˚C with 20s each, extension at 72˚C for 80s and final extension at 72˚C for 10 min. Negative controls were used in all procedures. The amplified DNA fragments were stained with Gelred and visualized on agarose gel. All plates with PCR products were shipped to Functional Bioscience Inc., United States, for purification and sequencing. Geneious R8 6.0.5 [22] was used for editing and obtaining consensus sequences and automatic sequence alignment (coI, cytb), using the MUSCLE algorithm [23]. SATE v2 was employed for aligning 12S and 16S genes, also using MUSCLE. The sequences of the four mitochondrial genes were subsequently visually inspected and the alignments concatenated in Geneious.
Delimitation of possible multiple species was also tested using the General Mixed Yule Coalescent (GMYC) model, available on the GMYC server (https://species.h-its.org/gmyc; [30]). The coalescence trees were produced by Bayesian inference for the gene coI. The HKY model was used for all positions, considering relaxed molecular clock using a lognormal time distribution and birth-death prior through the softwares BEAUTi and BEAST. BEAST was programmed to run 100 million generations, sampling every 1,000 trees. TreeAnnotator [31] was used to build a consensus tree, with 10% of the trees discarded. To perform the analysis on the online server, the tree generated in TreeAnnotator was converted to Newick format using the software R v3 [32].

Morphological observations
Hypertrophied odontodes. Several development patterns of cheek hyperthrophied odontodes and associated lateral fleshy lobe were observed in males. These odontodes vary in size, thickness, density, direction, and position on the sides of the head, while the fleshy lobe varies in thickness and shape (Fig 4). Such variation, however, occurs independently of the geographic area and is found concurrently within each area, not being useful to distinguish populations. Conversely, the odontodes on the pectoral-fin spine were more delicate and homogeneous than the pattern detected for cheek odontodes, without significant variation. However, it is possible to observe syntopic adult males with small and barely visible odontodes on the pectoral-fin spine (Fig 5A), and others with slightly larger odontodes ( Fig 5B).
Coloration. The comparative analysis evidenced subtle differences in color pattern among individuals between and within areas, without a detectable geographical pattern. In general, the coloration varied from gray to light brown shades dorsally, with many scattered spots, and the ventral region usually with shades of pale yellow. Some specimens possess dark small or medium-sized dots clearly visible along the head and dorsal surface of the trunk ( Fig  6A and 6B). Other specimens show a darker dorsal coloration, with a more mottled pattern of fine to coarse dark vermiculations (Fig 6C and 6D).
Abdominal plates. Small, granular plates were observed in specimens from some of the study areas. Specimens from the Uruguay River basin showed a predominance of anterior abdominal plates, with the exception of those from Pelotas. Specimens from that area and those from areas in the Taquari River basin lack visible plates in the abdomen. The anterior abdominal plates are small and bear minute odontodes. They are located in the middle or at the edges of the pectoral girdle, usually just posterior to gills slits. In general, anterior abdominal plates occur only in adult specimens, both male and female (Fig 7A to 7E). One individual from Passo Fundo ( Fig 7A) had a large number of abdominal plates on the pectoral girdle.

PLOS ONE
The principal component analysis (PCA) did not differentiate the nine areas, which were generally overlapping. However, it evidences a slight separation by river basin, except for Prata and Pelotas, which were closer to their respective adjacent basins (Fig 9). Principal Component 1 (PC1) accounted for 28.9% of the total variance, PC2 for 20.1% and PC3 for 7.8%.
The Linear discriminant analysis (LDA) (Fig 10) revealed better differentiation between basins in comparison to PCA, with Pelotas clustering with areas of the Taquari basin. The percent separation achieved by each discriminant function was 51.8% for LDA1, 20.2% for LDA2, 10.3% for LDA3 and 6.5% for LDA4. The loads for discriminant function LDA1 suggest that the most significant measures were: (1) caudal-peduncle length (0.22), for which specimens from Middle Antas, Middle Uruguay, Ijuí, and Canoas presented higher length, with Middle Antas also grouping the shortest lengths along with Pelotas.     Plates at base of dorsal-fin 5 7 6.2 0.6 6 7 6.6 0.

Haplotype networks
The concatenated haplotype network (Fig 11) presented 22 haplotypes and one subtle genetic structure by geographic area, with unique and non-shared haplotypes. Pelotas presented proximity with specimens from areas in the Taquari River basin, while Passo Fundo and Chapecó presented greater distance from the other areas. For each individual haplotype network (Fig 12), in general, haplotype sharing was observed between the different areas, although exclusive haplotypes are also observed. There is low differentiation between populations, with few mutations separating them.  For coI (Fig 12) 12 haplotypes were observed in the network. There is a central haplotype shared between the different river basins, from which several others depart, with 1-4 mutational steps. The cytb (Fig 12) presented 23 haplotypes in the network, and a greater genetic structure per area, with only haplotype 17 shared between the separate basins: Pelotas (Uruguay River basin) and Middle Antas (Taquari River basin). The 16S (Fig 12) showed little variation among specimens, consisting of only six haplotypes with few mutations between them (1-3 mutational steps), with a central haplotype shared by all areas except Chapecó. Finally, for the 12S (Fig 12), 15 haplotypes were obtained, the central haplotype being shared between the two watersheds, originating several others with few mutations (1-5 steps).

Genetic distance
The genetic distance of the concatenated data (Table 6) points to a species-level separation with a percent above 3%, as observed for Pareiorhaphis azygolechis and P. stendachneri, or approximately 1% for P. vestigipinnis. In general, each area of P. hystrix, when compared to the others, presented genetic distance below 1%, with Passo Fundo (0.60 to 0.95%) and Chapecó (0.66 to 0.95%) presenting highest genetic distance when compared to other areas. For coI alone (Table 7), the genetic distance indicated a species-level separation ranging from 4-6%, except for P. vestigipinnis, which also showed lower differentiation (less than 2%). The comparison among areas of P. hystrix never exceeded 1.37%, and Prata presented the greater genetic distance relative to the other areas (some values higher than 1%).

Phylogeny
The phylogenetic analysis of the concatenated genes revealed most basal branches with high posterior probability (PP) values, and found Pareiorhaphis hystrix as monophyletic at 99% PP (Fig 13). Chapecó grouped the specimens in a clade sister to all other areas (100% PP), while Passo Fundo had a longer branch length when compared to the others, indicating a larger number of transformations. Despite some individual clades depict high support (96-100% PP), most internal nodes are weakly supported and highly polytomic, producing no separation between areas or even between watersheds.

Coalescence analysis
The coalescence analysis-GMYC with the entire data matrix suggested the presence of 14 operational taxonomic units (OTUs), 12 of which grouping mixed populations from different areas and both river basins, one of which including specimens of P. hystrix and P. vestigipinnis. The maximum likelihood for the null model was 344.29, and 349.43 for the GMYC model, with significant difference (0.005 � ). When the genetically most distant outgroup, P. steindachneri, was removed (Fig 14), the GMYC model suggested three separate OTUs, with low support values. The first OTU was represented by specimens from different areas and basins, the second by specimens from Middle Uruguay, and the third by P. vestigipinnis and P. hystrix from Canoas and Middle Uruguay (Uruguay River basin), and Prata (Taquari River basin). The maximum likelihood for the null model was 337.37, and 345.52 for the GMYC model, with significant difference (0.0002 �� ).

Discussion
Specimens of Pareiorhaphis hystrix with both small and delicate or large and strong odontodes in the cheek fleshy lobe are common and were compared in this study. The cheek fleshy lobe also shows a large variation in size and shape. When the fleshy lobe is well developed, odontodes tend to be poorly visible, since they are structures attached to the bone and emerging through the fleshy lobe. Both the variation in odontodes and the associated fleshy lobe did not show any association with area or river drainage. Similarly, observed variation in coloration does not seem to be related with area or river basin. The color pattern of the species as described by Pereira & Reis [33], with the dorsal surface of body and head dark gray, sometimes brownish gray, covered by dark brown to black blotches, and the ventral region characterized by a light, pale yellow, presented variation between the specimens studied. This variation is expected at intraspecific level and is not sufficient to distinguish specimens at interspecific level. The original description of the species [33] did not mention the presence of plates in the anterior region of the abdomen as a characteristic of Pareiorhaphis hystrix, and the presence of such plates is a novelty herein described. Some species such as P. nasuta Pereira, Vieira & Reis, 2007 and P. scutula Pereira, Vieira & Reis, 2010, described from the Doce River basin, P. ruschii Pereira, Lehmann & Reis, 2012, from tributaries of the Piraquê-Açu and Reis Magos rivers, and P. parmula Pereira, 2005, from the Iguaçu River basin, possess minute plates embedded in the abdominal skin. These species have minute plates distributed on each side of the pectoral girdle, right posterior to the gill opening, and P. scutula also has plates scattered throughout the abdominal region, from the pectoral girdle to the insertion of the pelvic fin. The function of such plates is still unknown, but specimens possessing a high concentration of plates in the abdomen may have an advantage when fixing themselves to stones in a strong water current, a common environment for these loricariids. In P. hystrix abdominal plates occur in individuals from the Uruguay River basin except the population from Pelotas, which is devoid of plates like the populations from the Taquari River. This pattern is corroborated by some results of the ANOVA and by the LDA which also indicated a closer similarity of the Pelotas and the areas in the Taquari River basin.
The significant variation ANOVA found for some counts, such as the number of premaxillary and dentary teeth, may be due to the fact that specimens often lose and replace their teeth and may be in different stages of tooth replacement [34]. Although significant values were observed, the differences between the counts were not discrete, fitting within the expected variation for a species.
The star-shaped individual haplotype networks with a central haplotype and additional haplotype starting with a few mutations (genes coI, 16S, 12S), may suggest a recent expansion of the parental haplotype [35][36][37][38]. Moreover, a high number of haplotypes was found considering the relatively small sample in some of the areas, followed by the small nucleotide diversity-few mutations among haplotypes, also corroborating the idea of a recent expansion [39,40]. The cytb gene is the most differentiated between the areas, with each area carrying its own haplotypes. The concatenated network, on the other hand, lost the star pattern observed individually, with the areas having their own haplotypes, although the distance between them still being small. In general, the most different areas are Chapecó and Passo Fundo, in agreement with the genetic distance and the concatenated phylogeny obtained.
Recently, Lima et al. [41] conducted a population analysis of Pareiorhaphis garbei in the Macaé and Macacu rivers, southeastern Brazil, where they presented haplotype networks generated for coI and cytb. No haplotype was shared among populations with a large number of mutational steps between populations. Contrarily, in the present study the haplotypes of the separated populations are very close, and there is still a large internal variability within the areas. Although there is not yet enough morphological evidence indicating a different species restricted to the upper-middle portion of the Uruguay basin, a slight differentiation by geographic area is already observed, especially in the Chapecó and Passo Fundo areas.
Although the coalescence species delimitation analysis suggests the existence of separate OTUs, this separation varies largely with the removal of the most genetically distant outgroup, suggesting the revealed genetic structure represents an artifact, a conclusion also supported by the groups not being observed iteratively in this study. The genetic distance for coI fails to indicate the presence of multiple species, considering the 2% threshold suggested to identify separate species lineages [42][43][44][45]. Higher values of genetic differentiation observed for Prata and Canoas specimens may have influenced the identification of OTUs in the coalescence analysis.
A member of the outgroup, Pareiorhaphis vestigipinnis, is poorly differentiated from P. hystrix. The main distinguishing morphological features of the former species are the absence of an adipose fin, which is always present and well-developed in P. hystrix, and the club-shaped pectoral-fin spine [46], which is homogeneously wide in the latter. While the phylogeny recovered P. vestigipinnis as sister to P. hystrix, corroborating the results of Pereira & Reis [8], it was found more closely related, despite weakly supported, to P. hystrix from Canoas, Prata, and Middle Uruguay in the coalescence analysis. Moreover, the concatenated haplotype network corroborates the phylogeny, where a greater separation is observed between P. vestigipinnis and all populations of P. hystrix. The 2% interspecific differentiation threshold commonly used for coI does not appear to apply to P. vestigipinnis and P. hystrix, which have genetic divergence with values even below 1% when compared to some P. hystrix populations, while P. azygolechis and P. steindachneri show distances above 3%. Despite the low genetic distance, two morphological features consistently distinguish P. vestigipinnis from P. hystrix, and they are thus maintained as a separate species.
While morphology indicates a slight trend of watershed separation, this was not corroborated by the molecular data. In general, haplotype sharing occurs between the two river basins with few mutations separating the haplotypes. However, some specimens from Pelotas presented greater proximity to those from Taquari, similarly to the morphology. A closer proximity of the individuals from Pelotas, which belongs to the Uruguay River basin, to specimens from the adjacent Taquari River basin, could be possibly explained by headwater capture events [47][48][49][50][51]. For example, the trichomycterid catfish Cambeva balios (Ferrer & Malabarba, 2013) occurs in the upper portions of the tributaries of the Antas and Caí rivers of the Taquari basin, and in the headwaters of the coastal Mampituba River. Such occurrence at the headwaters of two unconnected drainages could be possibly explained by headwater capture events [48]. The same phenomenon is reported for Astyanax brachypterygium Bertaco & Malabarba, 2001, which occurs in the upper portion of the Uruguay, Jacui, and Taquari basins, and Cambeva poikilos (Ferrer & Malabarba, 2013), with specimens at the upper portions of the Taquari and Jacuí tributaries. Headwater capture is known for moving a creek from one river basin to another, and thus mixing different species groups, considerably increasing regional levels of interspecific richness [52].
Although Pareiorhaphis hystrix is widely distributed in the Uruguay and upper Taquari river basins, this pattern does not seem to be common to other loricariids. The opposite situation occurs in other loricariid genus sharing most of the same geographic distribution of P. hystrix, the hypoptopomatine Eurycheilichthys. The type species, E. pantherinus (Reis & Schaefer, 1992) is widely distributed in the upper reaches of the Uruguay River, E. limulus Reis & Schaefer, 1998 is endemic to the upper portions of the Jacuí River, and seven species occur in tributaries of the Taquari River [53]. Contrary to Eurycheilichthys, the present analyses indicate that P. hystrix retains old haplotypes common to both watersheds, showing little differentiation between the specimens in these geographical areas. If these individuals were separated long enough, there would probably become different species distributed along the two watersheds. Natural selection seems to be acting differently on the morphology and the genetics of P. hystrix, as the iterative analyses suggest a single, although phenotypically variable species, thus not refuting the null hypothesis.
Supporting information S1 Table. Information of specimens used in morphological and molecular analyses of Pareiorhaphis. Specimens used for morphological analyses are marked MORPH in the Gene column. (DOCX)