Hydrography rather than lip morphology better explains the evolutionary relationship between Gymnogeophagus labiatus and G. lacustris in Southern Brazil (Cichlidae: Geophagini)

Pedro Ivo C. C. Figueiredo1,2, Luiz R. Malabarba1,3 and Nelson J. R. Fagundes1,2,4

PDF: EN    XML: EN  | Supplementary: S1  | Cite this article

Abstract​


EN

Gymnogeophagus labiatus and G. lacustris have been long recognized as sister species exhibiting different ecological requirements. Gymnogeophagus labiatus occurs in rock bottom rivers in the hydrographic basins of Patos Lagoon (HBP) and Tramandaí River (HBT), while G. lacustris is exclusive from sand bottom coastal lagoons of the HBT. In this study, we used molecular markers, morphological measurements and data from nuptial male coloration to investigate the evolutionary relationship between these species in each hydrographic basin. We found, for all data sets, a closer relationship between G. labiatus and G. lacustris from the HBT than between G. labiatus populations from HBT and HBP. In particular, lip area had a large intraspecific plasticity, being uninformative to diagnose G. lacustris from G. labiatus. Molecular clock-based estimates suggest a recent divergence between species in the HBT (17,000 years ago), but not between G. labiatus from HBP and HBT (3.6 millions of years ago). Finally, we also found a divergent G. labiatus genetic lineage from the Camaquã River, in the HBP. These results show that the current taxonomy of G. labiatus and G. lacustris does not properly represent evolutionary lineages in these species.

Keywords: Ecological divergence, Iterative taxonomy, Mitochondrial DNA, Phenotypic plasticity, Tramandaí-Mampituba ecoregion.

PT

Gymnogeophagus labiatus e G. lacustris vêm sendo consideradas espécies irmãs que possuem diferentes exigências ecológicas. Gymnogeophagus labiatus ocorre em rios de fundo de pedra nas bacias hidrográficas da Laguna dos Patos (HBP) e do rio Tramandaí (HBT), enquanto G. lacustris é exclusivo da HBT, ocorrendo em lagoas costeiras de fundo de arenoso. Nesse estudo, foram usados marcadores moleculares, medidas morfológicas e dados sobre a coloração nupcial em machos para investigar a relação evolutiva entre estas espécies em cada bacia hidrográfica. Para todos os conjuntos de dados foi observada uma relação mais próxima entre G. labiatus e G. lacustris da HBT do que entre as populações de G. labiatus da HBP e HBT. Em particular, a área do lábio teve uma grande plasticidade intraespecífica, não sendo informativa para diagnosticar G. lacustris de G. labiatus. Estimativas baseadas no relógio molecular sugeriram uma divergência recente entre as espécies da HBT (17.000 anos atrás), mas não entre as populações de G. labiatus da HBP e HBT (3,6 milhões de anos atrás). Finalmente, também foi encontrada uma linhagem genética de G. labiatus divergente no rio Camaquã, na HBP. Esses resultados mostram que a taxonomia atual de G. labiatus e G. lacustris não representa adequadamente as linhagens evolutivas nessas espécies.

Palavras-chave: Divergência ecológica, DNA mitocondrial, Ecorregião Tramandaí-Mampituba, Plasticidade fenotípica, Taxonomia iterativa.

Introduction​


Cichlids comprise 1727 valid species showing a wide morphological diversity in the shape of the body and also in its color pattern (Fricke et al., 2021). Many morphological features contributing to such diversity evolved repeatedly in the adaptive radiations of both African and Neotropical cichlids, including characters with known adaptive roles, such as hypertrophic lips, which facilitate foraging in rocky substrate, and body color, which has an important role in sexual selection (Henning, Meyer, 2014; Meier et al., 2017). In the African cichlids, several studies have revealed extensive parallelism at the morphological level with relatively frequent hybridization events, making species delimitation a difficult task (Meier et al., 2017; Salzburger, 2018). In Neotropical cichlids, the role of hybridization has not been studied to the same extent, but genetic data for Cichla Bloch & Schneider, 1801 have also indicated extensive introgression, with important impacts for species delimitation and identification (Willis et al., 2012; Mourão et al., 2017; Diamante et al., 2021).

Cichlinae, the clade representing Neotropical cichlids, contains about one third (566 species) of all valid species for the family (Fricke et al., 2021). Among these, approximately half belong to Geophagini, a well-supported clade within Cichlinae with approximately 18 genera and 250 species distributed in South America and southern Panama (López-Fernández et al., 2010; Ilves et al., 2018). Gymnogeophagus Miranda Ribeiro, 1918 contains 19 valid species (Turcati et al., 2018; Alonso et al., 2019) showing two synapomorphies: a spinous process directed forward on top of the first dorsal pterygiophore (unique among Neotropical cichlids), and the absence of supraneurals (Reis, Malabarba, 1988). Representatives of this genus are distributed among the Paraná, Paraguay and Uruguay basins and in small coastal drainages of Uruguay and southern Brazil, with the exception of G. balzanii (Perugia, 1891), which also occurs in the Guaporé River, in the Amazon drainage (Malabarba et al., 2015; Loureiro et al., 2016; Casciotta et al., 2017).

Gymnogeophagus labiatus (Hensel, 1870) and G. lacustris Reis & Malabarba, 1988 have been recovered as sister species in morphological and mitochondrial DNA phylogenies (Reis, Malabarba, 1988; Wimberger, Reis, Thornton, 1998), and have been differentiated by the presence of hypertrophied lips in G. labiatus. This trait possibly represents an adaptation for foraging in rocky environments (Elmer et al., 2010; Burress et al., 2013), as shown by experimental studies in other cichlid species (Baumgarten et al., 2015; Henning et al., 2017). Indeed, G. labiatus has been recorded in rock bottom streams associated with lotic environments with transparent water (Reis, Malabarba, 1988) in the hydrographic basins of Patos Lagoon (HBP) and Tramandaí River (HBT) (Reis, Malabarba, 1988), with additional records from the Mampituba River, an isolated coastal drainage (Malabarba et al., 2013). In contrast, G. lacustris is restricted to lentic environments, in sandy bottom coastal lagoons in the HBT (Reis, Malabarba, 1988), where it shows a preference for environments with little submerged or emergent vegetation or with no vegetation at all (Malabarba et al., 2013). Both species are oral incubators and have marked sexual dimorphism in the reproductive season, when the dominant males develop a hump in the anterior region of the head. This feature may be used for display during the reproductive period or as an energy reserve during the period of parental care (Lowe-McConnell, 1999).

The HBT can be divided into two main ecological subregions (Malabarba, Isaia 1992) that match the habitats of G. labiatus and G. lacustris, respectively. The first is represented by the ancient rocky bottomed rivers and streams that flow in the Serra Geral slopes, where the Maquiné and Três Forquilhas rivers are the major water courses. The second is represented by the sandy lagoons of the Coastal Plain formed in the last 5,000 years following recurrent episodes of Laguna-Barreira marine transgressions that have been occurring since the last 400,000 years (Schwarzbold, Schäfer 1984; Villwock, 1984; Tomazelli, Villwock, 2000). The HBP also possesses rocky bottomed rivers and streams, widely inhabited by G. labiatus, but also sandy lagoons where no populations assignable to either G. labiatus or G. lacustris are found. These ecological differences make these species interesting biological models to study the emergence of ecological specialization. However, understanding such specialization requires a more thorough characterization of the evolutionary scenario involving these lineages.

For example, if the divergence between G. labiatus and G. lacustris represent reciprocally monophyletic lineages, the occurrence of the former species in HBT and HBP would probably represent a recent dispersal event from one drainage to another. Alternatively, the late emergence of HBT lagoons may suggest that G. lacustris originated from a relatively recent speciation event, in which populations related to G. labiatus occurring in HBT slopes colonized and adapted to this lentic environment. In this case, HBT populations assigned to G. labiatus and G. lacustris would be more related, making G. labiatus paraphyletic relative to G. lacustris. Thus, to be recognized as separate species, G. lacustris and G. labiatus populations from HBT must represent recognizable independent lineages (Wiley, 1978; de Queiroz, 2007; Malabarba et al., 2020). In either case, if enough genetic, ecological and morphological differences have accumulated in each of these three lineages (G. labiatus HBP, G. labiatus HBT, G. lacustris), they could be recognized as three different species.

However, because hypertrophic lips, the main feature separating G. labiatus from G. lacustris, have independently evolved several times among cichlids (Henning, Meyer, 2014), and is a polymorphic trait in some species (Machado-Schiaffino et al., 2014), another hypothesis is that the HBT populations could constitute a single cohesive evolutionary lineage, so that the morphological differences between species could be interpreted as ontogenetic and ecological differences between populations. In this case, G. labiatus and G. lacustris populations from the HBT should not be recognized as independent evolutionary lineages. Here, we used genetic markers, morphological measurements and color pattern to discuss these alternatives, which may have important taxonomic implications.

Material and methods


Sampling and laboratory methods. All individuals used in the study are deposited in the scientific collection of the Laboratory of Ichthyology, Departamento de Zoologia, Universidade Federal do Rio Grande do Sul, Porto Alegre (UFRGS). The taxonomic identification of all specimens was based on their external morphology (Reis, Malabarba, 1988). We used 78 individuals (60 G. labiatus and 18 G. lacustris) for the morphological analysis, and tissue samples from 61 individuals (36 G. labiatus and 25 G. lacustris) for the genetic analysis (Fig. 1; Tab. 1).

FIGURE 1 | Map with the distribution of all individuals used in the present study. The hydrographic basins of the Patos Lagoon (HBP) and the Tramandaí River (HBT), as well as the Camaquã Sub-Basin (CSB), are indicated on the map. The gray lines represent the watershed between basins or sub-basins in the study area. The red dots correspond to Gymnogeophagus lacustris, the green to G. labiatus HBT, the yellow to G. labiatus HBP and the orange to G. labiatus CSB. 

TABLE 1 | Individuals and locations of Gymnogeophagus labiatus and G. lacustris included in the present study. 1Sample size for the molecular markers consider individuals characterized for both cytB and D-loop. 2Located in the Camaquã River sub-basin.

Voucher UFRGS

Species

Basin

Sampling size

Coordinates

Locality

 

 

 

Molecular1

Morphology

 

 

3904

G. labiatus

Patos

2

31°37’00”S 53°41’00”W

Bagé, RS

5924

G. labiatus

Patos

1

28°32’24”S 51°33’36”W

Guaporé, RS

6381

G. labiatus

Patos

4

28°57’40”S 51°45’17”W

Cotiporã, RS

6402

G. labiatus

Patos

1

28°53’04”S 51°47’28”W

Guaporé, RS

6409

G. labiatus

Patos

2

28°56’23”S 51°46’47”W

Dois Lajeados, RS

6464

G. labiatus

Patos

1

28°56’01”S 51°28’01”W

Vila Flores, RS

6596

G. labiatus

Patos

1

29°33’54”S 53°17’09”W

Agudo, RS

6606

G. labiatus

Patos

1

31°32’59”S 53°46’17”W

Candiota, RS

6962

G. labiatus

Patos

2

Rio Carreiro, RS

8402

G. labiatus

Patos

2

31°43’10”S 52°53’59”W

Pedro Osório, RS

8800

G. labiatus

Patos

1

29°22’08”S 52°03’30”W

Lageado, RS

8807

G. labiatus

Patos

1

29°19’21”S 52°14’03”W

Lageado, RS

9978

G. labiatus

Patos

3

28°56’15”S 51°27’54”W

Veranópolis, RS

10099

G. labiatus

Patos

1

28°57’47”S 51°45’43”W

Dois Lajeados, RS

10746

G. labiatus

Patos

1

1

30°06’02”S 51°41’40”W

Eldorado do Sul, RS

14174

G. labiatus

Patos

4

29°17’35”S 52°03’44”W

Travesseiro, RS

14317

G. labiatus

Patos

2

29°15’44”S 52°08’50”W

Marques de Souza, RS

19659

G. labiatus

Patos

2

29°20’47”S 50°42’04”W

Canela, RS

20360

G. labiatus

Patos

4

1

29°32’56”S 53°27’50”W

Faxinal do Soturno, RS

20396

G. labiatus

Patos

1

29°22’12”S 52°07’01”W

Forquetinha, RS

20407

G. labiatus

Patos

2

29°43’31”S 53°09’39”W

Paraíso do Sul, RS

22378

G. labiatus

Patos

2

2

30°23’00”S 51°26’00”W

Barra do Ribeiro, RS

22379

G. labiatus

Patos

8

3

29°59’15”S 51°14’24”W

Eldorado do Sul, RS

22380

G. labiatus

Patos

1

30°17’53”S 51°41’03”W

Barão do Triunfo, RS

22524

G. labiatus

Patos

2

28°42’12”S 51°50’57”W

Serafina Corrêa, RS

22719

G. labiatus

Patos

1

29°02’08”S 51°05’16”W

São Marcos, RS

23018

G. labiatus

Patos

4

6

30°17’00”S 51°48’00”W

Barra do Ribeiro, RS

22132

G. labiatus

Patos2

1

30°57’35”S 53°28’52”W

Caçapava do Sul, RS

17753

G. labiatus

Tramandaí

3

2

29°34’13”S 50°16’49”W

Maquiné, RS

17761

G. labiatus

Tramandaí

2

5

29°32’56”S 50°04’13”W

Três Forquilhas, RS

18279

G. labiatus

Tramandaí

1

29°43’42”S 50°08’06”W

Maquiné, RS

18437

G. labiatus

Tramandaí

1

29°34’14”S 50°16’49”W

Maquiné, RS

18457

G. labiatus

Tramandaí

1

29°40’09”S 50°04’59”W

Maquiné, RS

18464

G. labiatus

Tramandaí

2

29°40’08”S 50°12’24”W

Maquiné, RS

19594

G. labiatus

Tramandaí

1

29°32’13”S 50°14’45”W

Barra do Ouro, RS

21101

G. labiatus

Tramandaí

3

2

29°39’07”S 50°12’34”W

Maquiné, RS

21912

G. labiatus

Tramandaí

2

29°13’44”S 50°01’18”W

Praia Grande, SC

3885

G. lacustris

Tramandaí

1

29°42’00”S 50°05’59”W

Capão da Canoa, RS

3894

G. lacustris

Tramandaí

2

29°22’60”S 49°49’59”W

Torres, RS

10751

G. lacustris

Tramandaí

6

30°32’26”S 50°25’12”W

Mostardas, RS

16751

G. lacustris

Tramandaí

6

30°32’22”S 50°25’18”W

Mostardas, RS

16915

G. lacustris

Tramandaí

1

29°35’56”S 49°58’44”W

Terra de Areia, RS

17246

G. lacustris

Tramandaí

6

1

30°32’22”S 50°25’18”W

Mostardas, RS

17345

G. lacustris

Tramandaí

2

29°36’01”S 49°58’53”W

Terra de Areia, RS

17481

G. lacustris

Tramandaí

1

1

30°09’23”S 50°14’04”W

Cidreira, RS

18273

G. lacustris

Tramandaí

1

29°39’07”S 50°12’33”W

Maquiné, RS

18405

G. lacustris

Tramandaí

1

29°45’44”S 50°05’02”W

Capão da Canoa, RS

18406

G. lacustris

Tramandaí

2

30°32’26”S 50°25’16”W

Mostardas, RS

19081

G. lacustris

Tramandaí

4

29°57’56”S 50°13’46”W

Osório, RS

19511

G. lacustris

Tramandaí

4

29°47’03”S 50°11’02”W

Osório, RS

19570

G. lacustris

Tramandaí

3

29°58’16”S 50°13’07”W

Osório, RS

19586

G. lacustris

Tramandaí

1

29°57’04”S 50°13’02”W

Osório, RS

 

Specimens were obtained from the Patos Lagoon and Tramandaí-Mampituba ecoregions (Abell et al., 2008). The first includes all interconnected rivers and lagoons draining to the Patos Lagoon estuary. The second includes all interconnected lagoons draining to the Tramandaí Lagoon estuary, and can be further divided into two subsystems of interconnected lagoons (Schwarzbold, Schäfer 1984): one to the north of the Tramandaí Lagoon (formed by the Tramandaí River itself, Itapeva Lake, Quadros Lake, and a set of small lagoons in the municipality of Osório), and one to the south (including a string of interconnected lagoons reaching the Porteira Lake). We also included individuals sampled from isolated coastal lagoons (Bacopari Lake and Corvina Lake) that are presently mapped to the Patos Lagoon ecoregion (Abell et al., 2008) in order to understand the origin of its fish fauna.

For the molecular analysis, genomic DNA was obtained using the CTAB method adapted from Doyle, Doyle (1987). Two mitochondrial DNA (mtDNA) markers: Cytochrome B (cytB) and Control Region (D-loop) were amplified using the polymerase chain reaction (PCR) technique and specific primers (Palumbi et al., 2002; Sivasundar et al., 2001). The amplification reactions were prepared with 0.4 mM dNTP, 1.5 mM MgCl2, 0.5 μM of each primer, 1U Taq Polymerase and 40ng of genomic DNA. Amplification conditions for cytB consisted of 94ºC for 5’ and 10 cycles of 94ºC for 1’, 55ºC (−0.5ºC/cycle) for 1’ and 72ºC for 1’30”, followed by 30 cycles of 94ºC for 1’, 50ºC for 1’ and 72ºC for 1’30”, with a final extension of 72ºC for 5’. For D-loop, amplification conditions were 94ºC for 5’ and 10 cycles of 94ºC for 1’, 70ºC (−0.5ºC/cycle) for 1’ and 72ºC for 1’30”, followed by 30 cycles of 94ºC for 1’, 65ºC for 1’ and 72ºC for 1’30”, with a final extension of 72ºC for 5’. The success of the amplification was checked on a 1% agarose gel stained with GelRedTM (Biotium). The PCR products were enzymatically purified with Exonuclease I and Alkaline Phosphatase (ExoSAP) and sequenced by the Sanger method (Ludwig Biotec, Porto Alegre, Brazil). DNA sequencing was carried out in both directions (forward and reverse) for cytB, and only in the forward direction for the D-loop due to the presence of a repetitive region close to the annealing site of the reverse primer that caused poor reading quality. DNA sequencing was repeated from independent PCR amplifications whenever necessary to resolve ambiguities. We used Geneious 11.1.2 program (http://www.geneious.com/) to check the quality of the chromatograms and to assemble the consensus sequence for all individuals. The sequences were automatically aligned using ClustalW (Thompson et al., 1994) and edited in BioEdit (Hall, 1999). Sequences produced in the present study are available from the GenBank (cytB: MZ667483 – MZ667543; D-loop: MZ667544 – MZ667609).

Molecular data analysis. All analyses were performed with the concatenated cytB and D-loop mitochondrial markers. We retrieved from the GenBank cytB and D-loop sequences (GU736952 and MG581478, respectively) for G. setequedas Reis, Malabarba & Pavanelli, 1992, which was used as an outgroup. The list of distinct haplotypes was generated in the DnaSP 5.10.01 software (Librado, Rozas, 2009). The evolutionary relationship between haplotypes was estimated by the median-joining method (Bandelt et al., 1999) in the Network program (http://www.fluxus-engineering.com/) and by the Bayesian phylogenetic framework implemented in the package BEAST2 (Bouckaert et al., 2014). The best evolutionary model for each marker was estimated in Partition Finder (Lanfear et al., 2012) based on the Bayesian information criterion (BIC) (Sullivan, Joyce, 2005), which indicated four different substitution models: one for each cytB codon position (K80, HKY; TrN+G) and another for D-loop (HKY+I+G).

We estimated the age of the major clades by assuming, for both markers, a strict molecular clock model, which is generally well justified for analysis between closely related species (Li, Drummond, 2012; dos Reis et al., 2016). We assumed an evolutionary rate of 0.0024/site/million years (0.0019 – 0.0029/s/My) for cytB, which was inferred for Neotropical cichlids using a dataset with cytB and COI markers and fossil data to calibrate the phylogeny (Tougard et al., 2017). Because there was no prior information for the evolutionary rate of the D-loop, we assigned it to an unlinked clock partition whose rate was estimated relative to the cytB rate (e.g., Lambert et al., 2015; Ramos-Fregonezi et al., 2017). We assumed a lognormal prior distribution (M = 0.0, S = 2.0) for the D-loop rate can be estimated relative to the cytB rate. We used 100,000,000 steps in the Markov Chain Monte Carlo (MCMC) sampling every 1,000 steps and discarding the first 10,000,000 samples as burnin. We checked for convergence and sampling sufficiency in Tracer 1.7 (Rambaut et al., 2018) ensuring an effective sample size (ESS) >200 for all parameters.

The evolutionary relationship among populations was estimated using the STACEY module (Jones, 2017) implemented in BEAST2 (Bouckaert et al., 2014). In this strategy, individuals are classified into N “minimal clusters” a priori, and the number of potential species (ranging from one to N) is estimated based on the depth of coalescence of the genetic lineages within and between potential species, which is controlled by parameters collapse height (e) and collapse weight (w). In other words, the minimal clusters may be merged but not split to form potential species (Jones et al., 2015). Smaller values for e ​​are more sensitive to recent divergences but may inflate the number of potential species. The parameter ω controls the number of potential species and can be used as a proxy for prior taxonomic knowledge (Matos-Maraví et al., 2019). We assumed minimal clusters of individuals: Gymnogeophagus setequedas, G. lacustris, G. labiatus HBT, G. labiatus HBP, and G. labiatus CSB (from the Camaquã River drainage), which showed a very divergent haplotype (see Results). We used different values for e between 0.001 and 0.00001 to evaluate its impact over species delimitation. The prior for w was set using a beta distribution between [0,1], with an initial value of 0.5. Priors for the Yule birth-death model were set to the default distributions. Priorsfor population size parameters were changed according to the program manual. For popPriorScale we used a lognormal prior distribution (M = -4.0, S = 2.0, respectively), while for the parameter popPriorInvGamma we used a mixture of four gamma distributions (a=1.0 b=1.0). The clock model was set as explained previously. We used 100,000,000 steps in the Markov Chain Monte Carlo (MCMC) sampling every 1,000 steps and discarding the first 10,000,000 steps as burnin. As before, we checked for convergence and sampling sufficiency in Tracer 1.7 (Rambaut et al., 2018) ensuring an effective sample size (ESS) >200 for all parameters.

We used the Analysis of Molecular Variance (AMOVA) (Excoffier et al., 1992) to measure the degree of genetic structure among groups based on Ф-statistics. For this analysis, we used the “minimal clusters” (with the exception of G. setequedas, see above) as populations: G. labiatus HBP, G. labiatus CSB G. labiatus HBT and G. lacustris (from HBT). We also estimated pairwise population differentiation based on Ф-statistics. All calculations were performed in Arlequin 3.5 (Excoffier, Lischer, 2010) using 10,000 permutations to assess statistical significance.

Morphological data analysis. For the morphological analyzes we use the linear measurements and counts following Malabarba et al. (2015), which are represented in Fig. 2A. The measurements were taken using a caliper with a precision of 0.05 mm, and the counts were made under the stereomicroscope. For measuring the lip area we took standardized photographs of all individuals in lateral view and measured the lip area using the software ImageJ (Abràmoff et al., 2004) (Fig. 2B). All linear measurements and photographs were taken from the specimens’ left side. To correct for allometry effects, all linear measurements were divided by the body standard length, while the lip area was divided by the total body area. All measurements were also transformed using logarithms to normalize the scale of variation among variables. We did not include sexually dimorphic characters.

FIGURE 2 | Measurements and counts applied to specimens used in the present study, according to Malabarba et al. (2015). 1: Standard length; 2: Body depth; 3: Head length; 4: Dorsal-fin base length; 5: Pectoral-fin length; 6: Caudal peduncle depth; 7: Caudal peduncle length; 8: Eye diameter; 9: Interorbital width; 10: Upper jaw length; 11: Pre-orbital length; 12: Snout length; E1: Scales in a longitudinal series; DL: Scales between the origin of the dorsal fin and the upper part of the lateral line; UL: Scales in the upper portion of the lateral line; LL: Scales in the lower portion of the lateral line; AL: Scales between the origin of the anal fin and the top of the lateral line. Measurements of the lip area applied to the specimens used in the present study, modified by Buckup, Reis, (1985). SL: upper lip area; IL: lower lip area; TL: total lip area; BA: body area. 

For the morphological analyses, three populations were considered: Gymnogeophagus labiatus HBP, G. labiatus HBT and G. lacustris. Gymnogeophagus labiatus CSB was not considered because only one individual was sampled. To check the role of lip area in the morphological differentiation among populations, we performed all analyses considering the linear measurements alone or combined to lip area. A principal component analysis (PCA) was performed to estimate the degree of morphological variation present in the sample, and population differentiation was tested using a non-parametric multivariate analysis of variance (PERMANOVA) from the Mahalanobis distance using 9,999 replications for the post-hoc pairwise test and considering Bonferroni’s correction for multiple tests. The degree of differentiation between populations was also estimated using a Discriminant Analysis or Canonical Variable Analysis (CVA). We used a Kruskal-Wallis test to check for interpopulation differences considering the three most informative morphological variables in CVA. All statistical tests were performed using the Past3.2 program (Hammer et al., 2001). Photographs of live specimens were analyzed to assess the color pattern of these groups.

Results​


Considering both cytB and D-loop markers, we obtained a total alignment of 1367bp (773bp, 594bp, respectively), with 92 variable sites distributed among 32 haplotypes for a global haplotype diversity of 0.884 +/- 0.030. We found 36 sites separating HBP and HBT populations, irrespective of the taxonomic affiliation of the specimens. Samples from Bacopari and Corvina isolated lagoons grouped with HBT populations. As shown in the haplotype network (Fig. 3), G. labiatus and G. lacustris had 10 closely related haplotypes in the HBT, two of which are shared between species, while the remaining haplotypes were observed in a single individual (seven G. lacustris and one G. labiatus). Concerning the HBP, one individual identified as G. labiatus sampled in the Camaquã Sub-Basin (CSB) showed a very divergent haplotype, sister to the clade formed by HBT and the remaining HBP haplotypes, which was considered as a different genetic population in the remaining analyses (Fig. 4). Other individuals from CSB sequenced for CytB showed haplotypes closely related to this divergent lineage (Fig. S1), but were not included in the main analysis because we were not able to get data for their D-loop sequence. These divergent lineages came from good quality sequencing reactions, and do not show any obvious features of nuclear mitochondrial DNA insertions (NUMTs) such as “heterozygous” sites, frameshift and nonsense mutations. In some specimens from the Camaquã River basin, for which sequences of only one of the two markers were obtained, we also found haplotypes from the canonical HBP clade, indicating that both mtDNA lineages occur in this drainage (Fig. S1). The AMOVA corroborated the large genetic structure among populations, which accounted for nearly ~90% of the total variation (ФST = 0.909, P < 0.0001). Pairwise ФST values ranged between -0.034 and 0.997 (Tab. 2), with non-significant values ​between the two species in HBT and for the comparisons involving CSB due to the small sample size for this population.

TABLE 2 | Genetic structure (FST) between Gymnogeophagus populations. Values in bold values were statistically significant (P < 0.05)

Population

G. labiatus CSB

G. labiatus HBP

G. labiatus HBT

G. lacustris

G. labiatus CSB

 

 

 

G. labiatus HBP

0.857

 

 

G. labiatus HBT

0.997

0.895

 

G. lacustris

0.988

0.915

-0.034

 

FIGURE 3 | Network of mitochondrial haplotypes of Gymnogeophagus considering the concatenated D-loop and CytB regions. The size of the circles is proportional to the number of individuals present in each haplotype. The colors represent the analyzed populations, according to the legend. Small circles represent medium vectors, and the number of traits or values associated with each branch corresponds to the mutational distance between haplotypes.

FIGURE 4 | MtDNA genealogy for the individuals analyzed. The values above the nodes correspond to the posterior probability of each clade. Only the values for the clades discussed in the text were presented. The colors represent the analyzed populations, according to the legend.

The population tree (Fig. 5) evidenced a very recent relationship between species in the HBT (~17,000 years, 95% CI 1,000 – 53,000 years). However, the divergence between the HBP and HBT clades dated ~ 3.6 Ma (95% CI 1.1 Ma – 6.0 Ma, while the divergence between the clade formed by these populations and the CSB lineage dated ~ 4.7 Ma (95% CI 2.8 Ma – 7.4 Ma). Under a coalescent criterion, the estimated number of putative species suggested, with high support, that G. labiatus HBT and G. lacustris belong to the same evolutionary lineage, while the other genetic populations of G. labiatus would constitute independent evolutionary lineages. These results were consistent regardless of the specific value of e (Tab. 3).

FIGURE 5 | Evolutionary relationship between populations of Gymnogeophagus based on mitochondrial haplotypes. The x-axis scale corresponds to millions of years. The values above the nodes correspond to the posterior probability of each clade. The bar associated with each node represents the 95% credibility interval for the date of the common ancestor. The colors represent the analyzed populations, according to the legend.

TABLE 3 | Delimitation of candidate species under a coalescent criterion. N species = number of candidate species; PP = Posterior probability for a specific arrangement, considering the number of species and which populations make up each group. The letters assigned to each population correspond to all candidate species suggested in the analysis. For a collapse height value of 0.001, three alternative arrangements were contained within the 95% credibility range. Note that in all arrangements Gymnogeophagus labiatus HBT and G. lacustris were assigned as the same candidate species.

Collapse Height

N species

PP

Populations

G. labiatus CSB

G. labiatus HBP

G. labiatus HBT

G. lacustris

G. setequedas

0.00001

4

0.995

A

B

C

C

D

 

 

 

 

 

 

 

 

0.0001

4

0.962

A

B

C

C

D

 

 

 

 

 

 

 

 

0.001

2

0.387

A

A

A

A

B

 

3

0.276

A

B

B

B

C

 

4

0.199

A

B

C

C

D

 

Because the individual showing the CSB lineage had a small body size (69.7 mm), and to avoid having a group with a single individual in the morphological analysis, only individuals sampled in HBT and HBP (except CSB) were included and assigned to one out of three populations: Gymnogeophagus labiatus HBP, G. labiatus HBT and G. lacustris. Scale counts showed low variation among populations and, for this reason, were excluded from further analyses. All measures are shown in Tab. 4. The PCA did not show a strong separation among populations, either or not considering lip area in addition to linear measurements (Fig. 6). Considering linear measurements only, PC1 accounted for 39.75% while PC2 accounted for 19.83% of the total variance. Considering lip area, PC1 accounted for 81.37% while PC2 accounted for 6.69% of the total variance. In addition to lip area, other variables with heavy loadings in the PCA were the length of the pectoral fin, the size of the snout and the height of the body. PERMANOVA showed significant differences among populations, regardless of whether or not the lip area was included in the analysis (F = 3.359, P = 0.0001; F = 3.643, P = 0.0001, respectively), with all pairs being statistically different from each other in both analyses (P < 0.01 for all pairs).

FIGURE 6 | Principal Component Analysis (PCA) for morphological data (A) disregarding the lip area; (B) considering the lip area. Gymnogeophagus lacustris is represented in red, G. labiatus HBT is represented in green and G. labiatus HBP is represented in yellow. 

TABLE 4 | Morphological measurements for the analyzed Gymnogeophagus specimens. Min = minimum observed value; Max = maximum observed value; SD = standard deviation; *Values relative to standard length; **Values relative to body area; E1 = Scales in a longitudinal series; DL = Scales between the origin of the dorsal fin and the upper part of the lateral line; AL = Scales between the origin of the anal fin and the top of the lateral line.

 

G. labiatus HBP (N = 48) 

G. labiatus HBT (N = 12) 

G. lacustris (N = 18) 

Min 

Max 

Mean 

SD 

Min 

Max 

Mean 

SD 

Min 

Max 

Mean 

SD 

Linear measurements (mm) 

 

 

 

 

 

 

 

 

 

 

 

 

Standard length (SL) 

80.9 

180.2 

111.4 

– 

80.3 

133.5 

100.6 

–  

80.43 

146.2 

100.8 

–  

Body depth* 

31.7 

39.3 

36.1 

1.57 

34.8 

39.6 

37.1 

1.62 

37.0 

41.0 

39.0 

1.18 

Head length* 

28.6 

38.1 

33.9 

1.78 

32.5 

38.0 

34.6 

1.65 

32.0 

35.8 

33.9 

0.92 

Dorsal-fin base length* 

47.7 

56.5 

52.3 

2.02 

48.1 

54.9 

51.5 

2.01 

51.0 

57.1 

52.9 

1.49 

Pectoral-fin length* 

24.0 

34.4 

27.7 

2.31 

27.7 

36.4 

31.0 

2.32 

25.4 

35.2 

32.0 

2.51 

Caudal peduncle depth* 

11.6 

13.8 

12.7 

0.50 

12.1 

13.6 

12.8 

0.42 

12.7 

14.2 

13.4 

0.39 

Caudal peduncle length* 

16.3 

21.2 

18.7 

1.20 

17.4 

20.4 

18.9 

1.08 

16.6 

21.1 

18.7 

1.24 

Eye diameter* 

5.7 

9.7 

7.4 

0.82 

6.8 

9.6 

8.2 

1.17 

6.9 

10.2 

8.6 

0.92 

Interorbital width* 

8.1 

12.5 

9.9 

1.06 

9.0 

12.4 

10.7 

1.11 

8.9 

12.8 

10.1 

1.06 

Upper jaw length* 

6.2 

11.1 

8.3 

1.11 

5.9 

8.9 

7.8 

0.89 

4.8 

8.6 

7.0 

0.94 

Pre-orbital length* 

9.2 

14.5 

12.0 

1.24 

10.6 

15.0 

12.4 

1.33 

10.8 

15.6 

13.1 

1.45 

Snout length * 

11.6 

19.6 

14.9 

1.85 

10.9 

19.0 

14.0 

2.17 

11.3 

15.6 

13.5 

1.22 

Area (cm²) 

 

 

 

 

 

 

 

 

 

 

 

 

Body area 

8.09 

56.77 

20.07 

–  

10.78 

34.80 

18.94 

–  

12.31 

41.90 

19.77 

–  

Upper lip area** 

0.3 

2.5 

1.2 

0.58 

0.5 

1.3 

0.9 

0.24 

0.2 

0.7 

0.4 

0.12 

Lower lip area** 

0.5 

3.2 

1.6 

0.70 

0.7 

1.3 

1.1 

0.20 

0.4 

0.8 

0.6 

0.12 

Total lip area** 

1.0 

5.6 

2.7 

1.25 

1.3 

2.6 

1.9 

0.42 

0.6 

1.4 

1.0 

0.21 

Counts 

 

 

 

 

 

 

 

 

 

 

 

 

Upper lateral line 

25 

30 

27.4 

1.03 

26 

28 

27.0 

0.60 

26 

28 

26.9 

0.73 

Lower lateral line 

11 

15 

12.4 

0.98 

10 

14 

11.8 

1.22 

10 

13 

11.7 

0.84 

E1 

11 

9.3 

0.66 

11 

9.3 

0.62 

11 

9.0 

0.77 

DL 

16 

20 

18.2 

0.93 

16 

20 

17.8 

1.11 

17 

19 

18.0 

0.69 

AL 

4.6 

0.54 

4.6 

0.51 

4.7 

0.46 

Spines 

12 

15 

13.7 

0.69 

12 

14 

12.9 

0.51 

12 

14 

12.8 

0.55 

Rays 

12 

10.1 

0.92 

10 

12 

11.0 

0.60 

10 

12 

10.8 

0.71 

 

The CVA showed high discrimination among groups, with 89.74% and 87.18% of correct classifications, either considering or not lip area, respectively (Tab. 5). Considering each measure individually, no morphological character was able to differentiate among the three populations simultaneously. However, while the size of the pectoral fin and the size of the snout discriminated populations from HBP or HBT (Fig. 7), lip area and body height differentiated between G. lacustris and G. labiatus (Fig. 8).

TABLE 5 | Confusion matrix for the CVA of morphological data of Gymnogeophagus species. 

Original population

Inferred population

G. labiatus HBP

G. labiatus HBT

G. lacustris

Total

Linear measurements + lip area

G. labiatus HBP

42

5

1

48

G. labiatus HBT

0

11

1

12

G. lacustris

1

0

17

18

Total