Introduction
Gyrodactylidae comprises a group of ectoparasites known mainly from bony fishes. Among the viviparous Gyrodactylidae, Gyrodactylus von Nordmann, 1832 is the most diverse genus, with more than 400 described species worldwide (Boeger et al., Reference Boeger, Kritsky and Pie2003; Bakke et al., Reference Bakke, Cable and Harris2007; Harris et al., Reference Harris, Shinn, Cable, Bakke and Bron2008), of which some are highly pathogenic and represent a serious threat to populations of wild and cultured fishes (McVicar, Reference McVicar1997; Bakke et al., Reference Bakke, Cable and Harris2007). In the Neotropical region, approximately 20 species of Gyrodactylus are currently known (Boeger et al., Reference Boeger, Vianna, Thatcher, Adis, Arias, Rueda-Delgado and Wantzen2006; Cohen and Kohn, Reference Cohen and Kohn2008). In South America, more than 10 species of Gyrodactylus are known from bony fishes (Kohn and Cohen, Reference Kohn and Cohen1998; Cohen and Kohn, Reference Cohen and Kohn2008), of which seven species are known to parasitize callichthyid catfishes in Brazil (Boeger et al., Reference Boeger, Vianna, Thatcher, Adis, Arias, Rueda-Delgado and Wantzen2006; Bueno-Silva and Boeger, Reference Bueno-Silva and Boeger2009, Reference Bueno-Silva and Boeger2014). The high diversification of Gyrodactylus is associated to a combination of biological characteristics (e.g. hyperviviparity and continuous host transfer) that contribute to the establishment of the parasites in a new host species and dispersion to distinct geographic regions (Bakke et al., Reference Bakke, Harris and Cable2002; Cable and Harris, Reference Cable and Harris2002; Morand et al., Reference Morand, Simková, Matejusová, Plaisance, Verneau and Desdevises2002; Boeger et al., Reference Boeger, Kritsky and Pie2003). Also, these characteristics seem to play a role in the diversification of lineages of Gyrodactylidae associated with distinct host species and geographic localities. However, the capacity of a lineage to diversify is dependent upon opportunities provided by the environment, as predicted by the Stockholm Paradigm (Araujo et al., Reference Araujo, Braga, Brooks, Agosta, Hoberg, von Hartenthal and Boeger2015; Hoberg and Brooks, Reference Hoberg and Brooks2015), which postulates that ecological factors (e.g. opportunity to exploit an alternative host species under changing environments) may exert influence on host-parasite relationships. Hence, realized diversification is often dependent on historical environmental events.
Indeed, despite the fact that abiotic and biotic components play an important role in the evolution of species, some authors have argued that coevolutionary interactions between parasites and hosts are dependent not only on ecological and genetic factors but also on geography (Kaltz and Shykoff, Reference Kaltz and Shykoff1998; Thompson, Reference Thompson1999; Gandon and Michalakis, Reference Gandon and Michalakis2002; Dybdahl and Storfer, Reference Dybdahl and Storfer2003; Prugnolle et al., Reference Prugnolle, Theron, Pointier, Jabbour-Zahab, Jarne, Durand and de Meeûs2005). The geographic structuring defines the diversity of ecological patterns and coevolutionary dynamics of species, which promotes the local adaptation of populations (Thompson, Reference Thompson1999, Reference Thompson2005). In addition, geographic structuring is considered one of the most important factors for the diversification, local adaptation, and composition of communities of parasites (Poulin and Morand, Reference Poulin and Morand1999; Thompson, Reference Thompson1999; Gandon and Michalakis, Reference Gandon and Michalakis2002; Dybdahl and Storfer, Reference Dybdahl and Storfer2003; Prugnolle et al., Reference Prugnolle, Theron, Pointier, Jabbour-Zahab, Jarne, Durand and de Meeûs2005; Nuismer, Reference Nuismer2006). Although there is some evidence of the factors that drive diversification of Gyrodactylidae (Morand et al., Reference Morand, Simková, Matejusová, Plaisance, Verneau and Desdevises2002; Ziętara and Lumme, Reference Ziętara and Lumme2002; Boeger et al., Reference Boeger, Kritsky and Pie2003; Meinilä et al., Reference Meinilä, Kuusela, Ziętara and Lumme2004; Bueno-Silva et al., Reference Bueno-Silva, Boeger and Pie2011), little is known about the role of the geographic structuring in this context (Malmberg, Reference Malmberg1970; Dmitrieva and Dimitrov, Reference Dmitrieva and Dimitrov2002; Meinilä et al., Reference Meinilä, Kuusela, Ziętara and Lumme2004; Pettersen et al., Reference Pettersen, Mo, Hansen and Vøllestad2015; Xavier et al., Reference Xavier, Faria, Paladini, van Oosterhout, Johnson and Cable2015). According to Thompson (Reference Thompson1999), the analysis of the host-parasite relationship is an important way to evaluate the local adaptation of species in geographically heterogeneous environments because parasitism (e.g. abundance and virulence) varies between host populations. The spatial variation in host-parasite interactions is one of the major driving forces of biodiversity and contributes to geographic mosaic selection (Thompson, Reference Thompson1999). The patterns of genetic differentiation at a small spatial scale could provide information on the mode of dispersal and evolutionary dynamics of parasites.
In South America, four species of Gyrodactylus are known to parasitize syntopic host species of Corydoras (Siluriformes: Callichthyidae), which represent an important host-parasite model (Boeger et al., Reference Boeger, Kritsky, Pie and Engers2005; Bueno-Silva et al., Reference Bueno-Silva, Boeger and Pie2011) to understand the diversification of species of Gyrodactylus (see Boeger et al., Reference Boeger, Kritsky and Pie2003). Bueno-Silva et al. (Reference Bueno-Silva, Boeger and Pie2011) found out that morphological and genetic differences between suprapopulations of Gyrodactylus corydori Bueno-Silva and Boeger, Reference Bueno-Silva and Boeger2009 corresponded to lineages of parasites from distinct host species and localities, respectively. This could provide a basis for understanding the role of the geographic structuring on the host-parasite interaction and the morphological and genetic diversity of populations. In the present study, an integrative eco-evolutionary (Pelletier et al., Reference Pelletier, Garant and Hendry2009) and comparative phylogeographic approach was adopted in order to understand the role of the historical processes on the geographic isolation, molecular evolution, and morphological differentiation of populations of Gyrodactylus and host species of Scleromystax (Siluriformes: Callichthyidae) in southern coastal Brazil. These Neotropical species were selected for this study based on an analogous host-parasite model analyzed by Bueno-Silva et al. (Reference Bueno-Silva, Boeger and Pie2011), of which one or more species of Gyrodactylus exploit two syntopic congeneric species of callichthyid catfish. Thus, based on previous studies on the host-parasite coevolutionary dynamics and diversification (Thompson, Reference Thompson1999; Desdevises et al., Reference Desdevises, Morand, Jousson and Legendre2002; Boeger et al., Reference Boeger, Kritsky and Pie2003; Huyse et al., Reference Huyse, Audenaert and Volckaert2003; Nieberding et al., Reference Nieberding, Morand, Libois and Michaux2004, Reference Nieberding, Durette-Desset, Vanderpoorten, Casanova, Ribas, Deffontaine, Feliu, Morand, Libois and Michaux2008; Huyse and Volckaert, Reference Huyse and Volckaert2005; Bueno-Silva et al., Reference Bueno-Silva, Boeger and Pie2011), the present study aimed to understand the spatial and temporal dynamics of coevolution of a host-parasite system, exploring eventual differences in processes of genetic and morphological differentiation of each group of associates. This article assessed whether: (i) The host-parasite interaction could be influenced by host transfer events; (ii) The low or absent gene flow between geographically distant subpopulations could contribute to events of evolutionary co-divergence; (iii) The level of morphological and genetic diversification could be related to the geographic distance between subpopulations, in accordance with the isolation-by-distance model (Wright, Reference Wright1943); and (iv) Local environmental and historical processes could have an effect on the past demography of parasites and its hosts. This study was carried out using a combination of molecular markers, such as mitochondrial DNA (mtDNA) and nuclear DNA (nDNA), and morphometrics.
Materials and methods
Study area and sample collection
The study area included two sub-basins of the coastal basin of the state of Paraná, southern Brazil: sub-basin of the Nhundiaquara River and sub-basin of the Paranaguá Bay (Fig. 1). These two areas were selected based on the assumption that the geomorphological characteristics and Quaternary paleoenvironmental changes in the region (Angulo and Lessa, Reference Angulo and Lessa1997; Ab'Saber, Reference Ab'Saber2000; Bigarella, Reference Bigarella2001; Maack, Reference Maack2001; Angulo et al., Reference Angulo, Pessenda and Souza2002) would make possible to test hypotheses about the genetic and morphological diversification of host-parasite populations. The extension of the coastal plain of Paraná ranges from 10 to 20 km wide, reaching the maximum of 50 km in the Paranaguá Bay, and an altitude of up to 10 m above sea level; at the innermost areas, it can reach up to 20 m of altitude (Bigarella, Reference Bigarella2001). Also, the coastal plain of Paraná is bordered on the east by the Atlantic Ocean and to the west by the mountainous border known as Serra do Mar (Bigarella, Reference Bigarella2001), including elevations reaching until more than 1000 m of altitude (Maack, Reference Maack2001).
Fig. 1. Study area at the coastal basin of the state of Paraná, southern Brazil. Sampling localities are indicated by black dots: (1) Marumbi River, municipality of Morretes; (2) Pinto River, municipality of Morretes; (3) Ribeirão River, municipality of Paranaguá; and (4) Fortuna River, municipality of Pontal do Paraná. The dashed line corresponds to the geographical limit between distinct sub-basins.
Two species of Scleromystax are known from coastal rivers in the south and southeast Brazil (Reis, Reference Reis, Reis, Kullander and Ferraris2003; Ferraris, Reference Ferraris2007), Scleromystax barbatus (Quoy and Gaimard, 1824) and Scleromystax macropterus (Regan, 1913). Adult specimens of S. barbatus and S. macropterus were collected from the above-mentioned sub-basins between 2007 and 2010. Scleromystax barbatus was found in the sub-basin of the Paranaguá Bay (n = 63), Ribeirão River (25°36′02″S, 48°37′19″W) at 21 m above sea level, municipality of Paranaguá, and Fortuna River (25°39′54″S, 48°31′05″W) at 8 m above sea level, municipality of Pontal do Paraná; and in the sub-basin of Nhundiaquara River (n = 44), Pinto River (25°30′50″S, 48°50′34″W) at 14 m above sea level, and Marumbi River (25°30′54″S, 48°52′03″W) at 31 m above sea level, municipality of Morretes. Scleromystax macropterus was found in the sub-basin of the Paranaguá Bay (n = 34), Fortuna River (25°39′54″S, 48°31′05″W), municipality of Pontal do Paraná. Method of parasite collection and preparation were the same as described by Bueno-Silva et al. (Reference Bueno-Silva, Boeger and Pie2011). Fish were captured by electrofishing, sacrificed by pithing and kept individually in plastic fish bags with warm water (60 °C) for a few seconds before being fixed in 95% ethanol. Host specimens were brought to the laboratory and screened for Gyrodactylus under a stereomicroscope. Tissue samples were taken from parasites and hosts and kept, respectively, in 95% ethanol and DMSO salt-saturated solution at −20 °C.
Two sets of morphometric and molecular techniques were used as independent sources of evidence. Morphometric and molecular analyses were conducted for parasites and hosts (see subsections below). Ecological terminology follows that described by Bush et al. (Reference Bush, Lafferty, Lotz and Shostak1997). Parasitological indices, such as prevalence (P) and mean intensity (MI ± s.e.), were obtained for each species of Gyrodactylus.
Morphometrics
The haptor of each parasite specimen was excised from the body and mounted on slides with Hoyer's medium to be used in morphometric analyses. The parasite trunk was used for DNA extraction. Haptoral structures of all parasites were photographed with a digital camera (Olympus QColor 5) connected to a phase-contrast microscope (Olympus BX 51), and the images were used to obtain the measurements of the haptoral sclerites using the software SigmaScan Pro 5.0 (SPSS, Inc.). Morphological features of the haptoral sclerites were measured as described in Bueno-Silva and Boeger (Reference Bueno-Silva and Boeger2009). Morphometrics of nine morphological features of the anchors and hooks were analyzed: anchor aperture distance (AAD), anchor proximal shaft width (APSW), anchor point curve angle (APCA), anchor inner curve length (AICL), anchor total length (ATL), hook shank length (HSHL), hooklet length (HL), hooklet aperture (HAD), and hooklet distal width (HDW). All measurements are in micrometres, except the anchor point curve angle (in degrees). For the APCA measurement, the cosine was calculated to transform the values to a linear function as described in Shinn et al. (Reference Shinn, Gibson and Sommerville2001).
Morphometric distances of the host specimens were obtained from the left side of each fish and measured according to Shibatta and Hoffmann (Reference Shibatta and Hoffmann2005) using an analogue calliper with an accuracy of 0.05 mm. All morphometric measurements by Shibatta and Hoffmann (Reference Shibatta and Hoffmann2005) were previously analyzed herein to avoid high misclassification rate, measurement error, and redundant measurements, as described in previous studies on morphometrics (Geets et al., Reference Geets, Appleby and Ollevier1999; Shinn et al., Reference Shinn, Gibson and Sommerville2001; Du Preez and Maritz, Reference Du Preez and Maritz2006; Bueno-Silva and Boeger, Reference Bueno-Silva and Boeger2009). Based on these parameters, only nine out of sixteen morphometric features described in Shibatta and Hoffmann (Reference Shibatta and Hoffmann2005) were measured for specimens of Scleromystax: snout length (SLE), head length (HLE), predorsal distance (PD), depth of body (DB), length of dorsal-fin base (LDB), dorsal to adipose distance (DAD), length of anal-fin base (LAB), depth of caudal peduncle (DCP), and standard length (SL).
All morphometric measurements were log-transformed to correct for increasing variance with increasing mean size of the measured variables (Lleonart et al., Reference Lleonart, Salat and Torres2000; Shinn et al., Reference Shinn, Gibson and Sommerville2001) and were analyzed by principal component analysis (PCA), multivariate analysis of variance (MANOVA) and stepwise discriminant analysis using Statistica 10 (Statsoft, Inc.). These analyses were conducted to assess comparatively the interspecific and interpopulational morphometric differentiation for parasites and hosts. The PCA was performed based on the covariance matrix (for hosts) and correlation matrix (for parasites) of the measured variables.
DNA extraction, PCR amplification, and sequencing
Total DNA was extracted from individual parasites using the DNeasy Tissue kit (Qiagen). The primer pairs cox2F (5′-TACAYAYCGCCCGTCAAYYTCG-3′), cox2R (5′-AATAMWKATWGGCATRWAAGARTG-3′), cox2F2 (5′- TTTCACTGAGATAAGTCGTAAC-3′) and cox2R2 (5′-TTACCGCTTCCYTGAACACG-3′) (Bueno-Silva and Boeger, Reference Bueno-Silva and Boeger2014) were used to amplify and sequence a fragment of approximately 570 bp of the mitochondrial cytochrome oxidase II gene (COII). The polymerase chain reaction (PCR) program for COII followed Bueno-Silva and Boeger (Reference Bueno-Silva and Boeger2014): 5 min at 95 °C, after which 40 cycles of 30 s at 94 °C, 45 s at 42–50 °C (gradient), 45 s at 72 °C, and finally 4 min at 72 °C. The nuclear zinc metalloproteinase gene (ZMP) and the internal transcribed spacer 1 (ITS1) and ITS2 were sequenced to confirm the mitochondrial data. The primers ZMPf1 (5′-CATGARRTTGGMCATAAYTTTGGATC-3′) and ZMPr1 (5′-AKCCACARTCRCATTCYTCDCC-3′) (Jarman et al., Reference Jarman, Ward and Elliott2002) were used to amplify and sequence a fragment of approximately 320 bp of the zinc metalloproteinase gene (ZMP), including intron and flanking exons. The PCR program for ZMP was the following: 5 min at 95 °C, followed by 40 cycles of 30 s at 94 °C, 45 s at 50 °C, 30 s at 72 °C, and finally 1 min at 72 °C. The primers ITS1 (5′-TTTCCGTAGGTGAACCT-3′) and ITS2 (5′-GGTAATCACGCTTGAATC-3′) (Ziętara et al., Reference Ziętara, Arndt, Geets, Hellemans and Volckaert2000) were used to amplify and sequence a fragment of approximately 1400 bp of the internal transcribed spacers ITS1 and ITS2. The PCR program for ITS was the same by Bueno-Silva et al. (Reference Bueno-Silva, Boeger and Pie2011): 5 min at 95 °C, followed by 40 cycles of 1 min at 95 °C, 45 s at 50 °C, 1 min at 72 °C, and finally 5 min at 72 °C. Each amplification reaction (COII, ZMP and ITS) was composed by 3 µL of template DNA (20–30 ng), 0.4 mm dNTP, 3 mm MgCl2, 1 U Platinum Taq polymerase (Invitrogen), 1 × PCR-buffer (Invitrogen), 0.4 pmol each primer and autoclaved water to complete 25 µL final volume.
Total DNA was extracted from individual host specimens using the EZ-DNA kit (Biosystems) or DNeasy Tissue kit (Qiagen). The primer pairs DLA-D (5′-TCCYACCCCTAACTCCCAAAG-3′), DLA-R (5′-AGTCAGGACCAARCCTTTGTGC-3′), DLB-D (5′-AGCRYCGGTCTTGTAATCCG-3′) and DLB-R (5′-GGYYCATCTTGACATCTTCAG-3′) (Montoya-Burgos, Reference Montoya-Burgos2003) were used to amplify and sequence a fragment of approximately 1000 bp of the mtDNA control region (CR). The PCR program for mtDNA CR was as follows: 5 min at 95 °C, after which 35 cycles of 30 s at 94 °C, 30 s at 63 °C, 35 s at 72 °C, and finally 2 min at 72 °C. The nuclear signal recognition particle 54 kDa gene (SRP54) was sequenced to confirm the mitochondrial data. The primers SRP54f1 (5′-ATGGGTGAYATYGAAGGACTGATWGATAAAGTCAA-3′) and SRP54r1 (5′-TTCATGATGTTYTGGAATTGYTCATACATGTC-3′) (Jarman et al., Reference Jarman, Ward and Elliott2002) were used to amplify and sequence a fragment of approximately 260 bp of the SRP54 gene, including intron and flanking exons. The PCR program for SRP54 was the following: 5 min at 95 °C, followed by 35 cycles of 30 s at 95 °C, 30 s at 55 °C, 30 s at 72 °C, and finally 2 min at 72 °C. The PCR mix for amplification of the mtDNA CR and SRP54 contained 1 µL of template DNA (~30 ng), 0.4 mm dNTP, 3 mm MgCl2, 1 U Taq Polymerase (Invitrogen), 1 × PCR-buffer (Invitrogen), 0.4 pmol each primer and autoclaved water to complete 25 µL final volume.
All PCR products were electrophoresed on a 1.5% agarose gel and then purified using MinElute Purification kit (Qiagen). Precipitated PCR products were sequenced using BigDye 3.1 chemistry on an ABI 3130 automatic sequencer (Applied Biosystems). Sequences were edited using the Staden Package 2.0 (Staden, Reference Staden1996) and aligned using MUSCLE (Edgar, Reference Edgar2004) implemented in the European Bioinformatics Institute website (McWilliam et al., Reference McWilliam, Li, Uludag, Squizzato, Park, Buso, Cowley and Lopez2013). All obtained sequences were deposited in GenBank: COII (MK395551-MK395607), ZMP (MK405282-MK405336, MK726393), mtDNA CR (MK395608-MK395653), and SRP54 (MK430079-MK430134).
Phylogenetics and molecular evolution
Phylogenetic relationships based on the COII and ZMP sequences were inferred using Bayesian analysis implemented in the software BEAST 2.5.2 (Bouckaert et al., Reference Bouckaert, Heled, Kühnert, Vaughan, Wu, Xie, Suchard, Rambaut and Drummond2014). The software jModelTest 2.1.4 (Darriba et al., Reference Darriba, Taboada, Doallo and Posada2012) was used to select the most appropriate model of DNA evolution (HKY + G for COII and ZMP; HKY + I for ITS). Posterior probabilities were estimated with non-partitioned data, four independent runs over 10 000 000 generations, sampled at each 1000th tree, a burn-in of 25% of the sampled trees, and estimated base frequencies. Two sequences of COII (KF751723) and ZMP (MK726393) of Aglaiogyrodactylus ctenistus Kritsky, Vianna et Boeger, 2007 and G. corydori were used as an outgroup, respectively. The molecular evolutionary relative rates of COII (2% per My) and ZMP (0.92% per My) were calibrated based on the nucleotide divergence between local suprapopulations of Gyrodactylus scleromystaci Bueno-Silva and Boeger, Reference Bueno-Silva and Boeger2014 from distinct sub-basins and the Quaternary fluvial sediments of the coast of Paraná (see Angulo, Reference Angulo2004). These fluvial sediments were found in present alluvial plains, which originated before and after the last post-glacial marine transgression, namely between the Pleistocene and Holocene Epoch (see Angulo, Reference Angulo2004). A comparative calibration of molecular clocks was performed to test the hypothesis that subpopulations were isolated before and after the last post-glacial maximum transgression. This test revealed unrealistic overestimated rates of molecular evolution for subpopulations hypothetically isolated since the Late Pleistocene up to after the last Holocene maximum transgression (data not shown). Since estimates suggest that S. barbatus is present in the Paranaguá Bay at more than 1 Mya (see Tschá, Reference Tschá2016), molecular clocks were calibrated herein based on the postulate that subpopulations from distinct localities were possibly isolated during the Early Pleistocene, until 2.5 Mya (see Discussion). The estimated evolutionary rate of COII was used to calculate the divergence dates of haplotypes and species of Gyrodactylus from Scleromystax under an uncorrelated lognormal relaxed molecular clock model (Drummond et al., Reference Drummond, Ho, Phillips and Rambaut2006) implemented in BEAST 2.5.2. The past demographical history was inferred based on the ZMP sequences using Bayesian skyline plots produced with the software BEAST 2.5.2 and Tracer 1.7.1 (Rambaut et al., Reference Rambaut, Drummond, Xie, Baele and Suchard2018) with the parameters described above. Bayesian skyline analysis was performed separately for suprapopulations of Gyrodactylus collected from different sub-basins in order to avoid a biased signal over the demographical history of parasites. Trees were simulated under the coalescent Bayesian skyline with a constant population model. Historical demographic changes in population size of parasites were also evaluated by mismatch distribution analysis (Rogers and Harpending, Reference Rogers and Harpending1992) and R 2 test (Ramos-Onsins and Rozas, Reference Ramos-Onsins and Rozas2002) using the software DnaSP 5.1 (Librado and Rozas, Reference Librado and Rozas2009). The R 2 test was performed conducting coalescent simulations (n = 10 000 replicates) on the number of segregating sites of COII sequences. The R 2 test is a powerful statistical test for detecting past population growth (Ramos-Onsins and Rozas, Reference Ramos-Onsins and Rozas2002). Analysis of molecular variance (AMOVA) and F-statistic (F st) were performed to estimate the degree of intra- and interpopulation genetic differentiation using the software Arlequin 3.5.1.3 (Excoffier and Lischer, Reference Excoffier and Lischer2010). The AMOVA and F st values were calculated based on the Kimura 2-parameter (K2P) distance matrix for COII and ZMP haplotypes under 10 000 permutations. Mean nucleotide divergence (d = mean ± s.d.) based on the COII and ZMP sequences was calculated using the software MEGA 6 (Tamura et al., Reference Tamura, Stecher, Peterson, Filipski and Kumar2013) under K2P substitution model with variance estimation by bootstrap method (1000 replicates). The compound DHEW neutrality test (Zeng et al., Reference Zeng, Shi and Wu2007) was used to detect positive selection on the COII and ZMP sequences using the software DH (available at http://zeng-lab.group.shef.ac.uk/wordpress). One sequence of Gyrodactylus from S. barbatus or S. macropterus was used as an outgroup for each parasite species. The significance level of tests was determined through Watterson's theta (θ w) under 50 000 simulations. The DHEW is a compound test of the Ewens-Watterson test (Watterson, Reference Watterson1978), Tajima's D (Tajima, Reference Tajima1989), and Fay and Wu's H (Fay and Wu, Reference Fay and Wu2000). The DHEW test is robust against the presence of recombination and is also relatively insensitive to purifying selection and demography, which allows detecting recent positive selection with high specificity (Zeng et al., Reference Zeng, Shi and Wu2007). Non-synonymous (Ka) and synonymous (Ks) substitution rate ratio (ω) was inferred for COII by an empirical Bayesian method under codon substitution models M8 (Yang et al., Reference Yang, Nielsen, Goldman and Pedersen2000) and M8a (Swanson et al., Reference Swanson, Nielsen and Yang2003) implemented in the web server Selecton 2.4 (Stern et al., Reference Stern, Doron-Faigenboim, Erez, Martz, Bacharach and Pupko2007). This analysis was performed to test the hypothesis for adaptive and neutral selection in each codon of the COII gene between and within species. Bayesian phylogeny inferred from COII sequences of Gyrodactylus was used to test for the existence of a biased distribution of parasite lineages between host species using the program Mesquite 2.7.4 (Maddison and Maddison, Reference Maddison and Maddison2010, available at http://www.mesquiteproject.org.) following the same method described by Bueno-Silva et al. (Reference Bueno-Silva, Boeger and Pie2011).
Phylogenetic relationships based on the mtDNA CR and SRP54 sequences were inferred using Bayesian analysis implemented in the software BEAST 2.5.2. The software jModelTest 2.1.4 was used to select the most appropriate model of DNA evolution (HKY + G for mtDNA CR; TrN for SRP54). Posterior probabilities were estimated with non-partitioned data, four independent runs over 10 000 000 generations, sampled at each 1000th tree, a burn-in of 25% of the sampled trees, and estimated base frequencies. Two sequences of mtDNA CR (GQ178156) and SRP54 (NM_200988) of Trichomycterus areolatus Valenciennes, 1846 and Danio rerio (Hamilton, 1822) were used as an outgroup, respectively. The molecular evolutionary relative rates of the mtDNA CR (0.32% per My) and SRP54 (0.4% per My) were calibrated based on the nucleotide divergence between local subpopulations of S. barbatus from different sub-basins and the geomorphological history of the coast of Paraná as described above. This evolutionary rate of the mtDNA CR was used to estimate the divergence dates of haplotypes and species of Scleromystax under an uncorrelated lognormal relaxed molecular clock model implemented in BEAST 2.5.2. The past demographical history of host populations was inferred based on the SRP54 sequences using Bayesian skyline plots produced with the programs BEAST 2.5.2 and Tracer 1.7.1 with the parameters described above. Bayesian skyline analysis was performed separately for subpopulations of Scleromystax collected from different sub-basins, following the same approach described above for parasites. Trees were simulated under the coalescent Bayesian skyline with a constant population model. Historical demographic changes in population size of host species were also assessed by mismatch distribution analysis and R 2 test using the software DnaSP 5.1 with the parameters described above. The AMOVA and F st values were used to estimate the degree of intra- and interpopulation genetic differentiation of Scleromystax in Arlequin 3.5.1.3. The AMOVA and F st values were calculated based on the K2P and TrN distance matrix for mtDNA CR and SRP54 haplotypes, respectively, under 10 000 permutations. Mean nucleotide divergence (d) based on the mtDNA CR and SRP54 sequences were calculated using the software MEGA 6 with the parameters described above. Differences in the length of mtDNA CR sequences were analyzed using the software TRF 4.04 (Benson, Reference Benson1999) in order to locate patterns of tandem repeats in DNA. The compound DHEW neutrality test was used to detect positive selection on the mtDNA CR and SRP54 sequences using the software DH. One sequence from each of the two species of Scleromystax was used as an outgroup for the other. The significance level of tests was calculated as described above.
On the assumption that morphologic variations between populations of parasites and hosts could result from the genetic divergence between geographically isolated subpopulations, it would be expected to find a positive correlation between the degree of genetic and morphological diversification in populations. This hypothesis was evaluated herein using the partial Mantel correlation (Mantel, Reference Mantel1967; Mantel and Valand, Reference Mantel and Valand1970) based on the morphometrical, genetic and geographic distance matrices obtained from populations of parasites and hosts. Matrices were permuted (n = 1 000 000) using the software ZT 1.1 (Bonnet and Van de Peer, Reference Bonnet and Van de Peer2002). Morphometric distance matrices were calculated based on the Euclidian distance of the morphometric measurements. Genetic distance matrices, inferred as the percentage of nucleotide differences, were based on the sequences of COII (for parasites) and mtDNA CR (for hosts). Geographic distance matrices (distance by land measured in Km) were obtained from the geographic coordinates of the host-parasite populations (see Study area and sample collection subsection for geographic coordinates) using the software GDMG 1.2.3 (Ersts, available at http://biodiversityinformatics.amnh.org/open_source/gdmg).
Given that S. macropterus was found exclusively in a single locality at the study area, the comparative phylogeographic analysis was carried out only for populations of S. barbatus and its respective parasites. Nevertheless, morphological and genetic data from S. macropterus and its respective parasites were taken and analyzed for comparison purposes, including the reconstruction of the historical demography, molecular evolution, and distribution of parasite lineages between host species.
Results
Identification of parasite specimens and occurrence
A total of 313 specimens of Gyrodactylus was recovered from the body surface of S. barbatus and S. macropterus. Four species of Gyrodactylus were recognized: Gyrodactylus bueni Bueno-Silva and Boeger, Reference Bueno-Silva and Boeger2014 (n = 124), Gyrodactylus major Bueno-Silva and Boeger, Reference Bueno-Silva and Boeger2014 (n = 90), G. scleromystaci (n = 92), and Gyrodactylus sp. (n = 7) (Table 1). Due to the limited number of Gyrodactylus sp., which seems to be a rare species (Bueno-Silva and Boeger, Reference Bueno-Silva and Boeger2014), these specimens were used only for estimating parasitological indices and genetic data. Sequences of the COII gene and ITS1-ITS2 of Gyrodactylus from Scleromystax obtained by Bueno-Silva and Boeger (Reference Bueno-Silva and Boeger2014) were analyzed herein (see Table 1).
Table 1. Occurrence, parasitological indices, and DNA sequences of four species of Gyrodactylus collected from Scleromystax barbatus and Scleromystax macropterus, southern Brazil
Prevalence (P) and mean intensity (MI ± s.e.) were estimated for each species of Gyrodactylus. Asterisks indicate host species whose parasites were collected: *Scleromystax barbatus and **Scleromystax macropterus.
Morphometric data
Interspecific and suprapopulational variation in parasites
The PCA performed with anchor and hook measurements of G. bueni, G. major, and G. scleromystaci revealed that the two first components accounted for 70% of the interspecific total variance. In this analysis, five variables showed the highest coefficient (factor loadings greater than 0.7) on the first principal component: anchor aperture distance (AAD), anchor proximal shaft width (APSW), anchor total length (ATL), hooklet length (HL), and hooklet aperture (HAD). On the second principal component, two variables presented the highest coefficient: anchor inner curve length (AICL) and anchor point curve angle (APCA). In the second PCA (performed only with anchor and hook measurements of G. scleromystaci), the two first components accounted for 55% of the suprapopulational total variance. In this analysis, three variables showed the highest coefficient on the first principal component: anchor aperture distance (AAD), anchor inner curve length (AICL), and anchor point curve angle (APCA). On the second principal component, two variables presented the highest coefficient: hooklet length (HL) and hooklet aperture (HAD). These five variables (AAD, AICL, APCA, HL, and HAD) were the same that accounted for the interspecific morphometric variation of Gyrodactylus from Scleromystax.
The MANOVA revealed significant interspecific differences in the morphology of anchor and hook of G. bueni, G. major, and G. scleromystaci (F = 94.6, P < 0.001). According to the discriminant analysis, parasites could be correctly assigned to their respective species in 100% of cases (Wilks’ Lambda = 0.15; F = 106.71; P < 0.001) by differences in eight morphological features of anchors and hooks (AAD, ATL, APSW, APCA, HSHL, HL, HDW, and HAD), while lineages of G. scleromystaci could be distinguished from distinct sub-basins in 71% of cases (Wilks’ Lambda = 0.74; F = 3.7; P < 0.01) by differences in five measurements of haptoral sclerites (AICL, APCA, APSW, HDW, and HAD) (Fig. 2). The univariate analysis of the morphometric features of anchors and hooks indicated that seven measurements are significantly different between species of Gyrodactylus: AAD (F = 231.6, P < 0.001), APSW (F = 3.09, P < 0.05), APCA (F = 3.91, P < 0.05), ATL (F = 242.5, P < 0.001), HSHL (F = 329.7, P < 0.001), HL (F = 6.9, P < 0.01) and HAD (F = 15.6, P < 0.001). In the second MANOVA (performed only with measurements of the anchor and hook of G. scleromystaci), considerable morphologic differences were found between suprapopulations from distinct sub-basins (F = 2.2, P < 0.05) and locations (F = 2.5, P < 0.01). In addition, the discriminant analysis based on these measurements of G. scleromystaci indicated that the suprapopulations could be distinguished from different locations (Wilks’ Lambda = 0.54; F = 6.45; P < 0.001) in 66% of cases. The univariate analysis of the morphometric features of G. scleromystaci revealed significant differences in the morphology of anchor and hook between suprapopulations from distinct sub-basins (AAD: F = 4.88, P < 0.05; AICL: F = 11.56, P < 0.01; APCA: F = 9.59, P < 0.01) and locations (AAD: F = 3.31, P < 0.05; AICL: F = 15.18, P < 0.001; APCA: F = 6.78, P < 0.01; HDW: F = 5.64, P < 0.01). In this analysis, the morphometric variation of suprapopulations of G. scleromystaci was significant between sub-basins and locations for the same measurements of anchor (AAD, AICL and APCA), with exception of the measurement HDW (hooklet distal width). Moreover, five out of eight variables that allowed delimitation of species of Gyrodactylus were the same that accounted for the morphometric variation between suprapopulations of G. scleromystaci from distinct sub-basins (HDW, HAD, APCA, and APSW) and locations (HDW and HSHL). Notwithstanding these morphometric variations, there was little correspondence between the degree of morphological variation and geographic distance for suprapopulations of G. scleromystaci (Mantel test; r 2 = 0.05, P < 0.001).
Fig. 2. Scatter plot of canonical discriminant components based on the morphometrical features of anchors and hooks of species of Gyrodactylus collected from Scleromystax barbatus and Scleromystax macropterus from the sub-basin of the Nhundiaquara River and the sub-basin of the Paranaguá Bay, state of Paraná, Brazil. Symbols correspond to parasite species: Gyrodactylus bueni (open squares), Gyrodactylus major (filled squares), Gyrodactylus scleromystaci from the sub-basin of the Nhundiaquara River (open circles), and Gyrodactylus scleromystaci from the sub-basin of the Paranaguá Bay (filled circles). Each ellipse represents a 95% confidence interval.
Interspecific and interpopulational variation in hosts
The standard length (SL) of S. barbatus from the sub-basin of the Paranaguá Bay (4.9 ± 1.3 cm) was higher than that of the S. barbatus from the sub-basin of Nhundiaquara River (4.5 ± 0.9 cm), whereas S. macropterus presented the lower standard length among host species (3.3 ± 0.6 cm). The PCA based on the morphometric features of S. barbatus and S. macropterus revealed that the two first components accounted for 94% of the interspecific total variance. In this analysis, all measured variables showed factor loadings greater than 0.7 on the first principal component: snout length (SLE), head length (HLE), predorsal distance (PD), depth of body (DB), length of dorsal-fin base (LDB), dorsal to adipose distance (DAD), length of anal-fin base (LAB), depth of caudal peduncle (DCP), and standard length (SL). In the second PCA (performed exclusively with measurements of S. barbatus), similar results were obtained in comparison to the first PCA: the two first components accounted for 94% of the interpopulational total variance and all measured variables showed factor loadings greater than 0.7 on the first principal component.
The MANOVA revealed significant interspecific differences in the morphology of S. macropterus and S. barbatus (F = 34.6, P < 0.001). According to the discriminant analysis, host species could be correctly distinguished in 99% of cases (Wilks’ Lambda = 0.2; F = 53.61; P < 0.001) by differences in six morphological features of the body (SLE, DCP, SL, LDB, DB, and LAB), while subpopulations of S. barbatus could be distinguished from different sub-basins in 92% of cases (Wilks’ Lambda = 0.28; F = 20.48; P < 0.001) by differences in seven measurements of the body (LAB, SL, DCP, DB, LDB, PD, and HLE) (Fig. 3). The univariate analysis of the morphometric features of hosts revealed seven measurements that were significantly different between species of Scleromystax: PD (F = 28.7, P < 0.001), SLE (F = 80.28, P < 0.001), HLE (F = 29.96, P < 0.001), DB (F = 18.53, P < 0.001), LDB (F = 16.65, P < 0.001), SL (F = 38.81, P < 0.001), and DAD (F = 28.83, P < 0.001). The second MANOVA (performed only with measurements of S. barbatus) revealed significant morphologic differences between subpopulations from distinct sub-basins (F = 15.46, P < 0.001) and locations (F = 7.85, P < 0.001). Also, the discriminant analysis based on the measurements of S. barbatus revealed that these subpopulations could be distinguished from different locations (Wilks’ Lambda = 0.18; F = 10.29; P < 0.001) in 83% of cases. In this analysis, some variables provided a significant unique contribution for distinction between subpopulations of S. barbatus from different sub-basins (LAB: F = 7.28, P < 0.01; SL: F = 27.32, P < 0.001; DCP: F = 36.94, P < 0.001; LDB: F = 27.92, P < 0.01) and locations (LAB: F = 4.32, P < 0.05; SL: F = 11.92, P < 0.001; DCP: F = 14.58, P < 0.001; SLE: F = 14.2, P < 0.001; DB: F = 10.29, P < 0.05). Among the morphometric features that allowed differentiation of subpopulations of S. barbatus from distinct sub-basins and locations, six measurements were the same as for both grouping variables (LAB, SL, DCP, DB, LDB, and PD). In addition, five out of six variables that allowed delimitation of species of Scleromystax were the same that accounted for the morphometric variation between subpopulations of S. barbatus from distinct sub-basins and locations simultaneously (DCP, SL, LDB, DB, and LAB). Despite these morphometric differences, there was little correlation between the degree of morphological variation and geographic distance for subpopulations of S. barbatus (Mantel test; r 2 = 0.005, P < 0.05).
Fig. 3. Scatter plot of canonical discriminant components based on the morphometrical features of the body of Scleromystax barbatus and Scleromystax macropterus from the sub-basin of the Nhundiaquara River and the sub-basin of the Paranaguá Bay, state of Paraná, Brazil. Symbols correspond to host species: Scleromystax barbatus from the sub-basin of the Nhundiaquara River (filled squares), Scleromystax barbatus from the sub-basin of the Paranaguá Bay (open squares), and Scleromystax macropterus (open circles). Each ellipse represents a 95% confidence interval.
Molecular data
DNA sequences of Gyrodactylus
A total of 81 sequences of the COII gene were obtained for four species of Gyrodactylus (20 sequences from G. bueni, 12 sequences from G. major, 45 sequences from G. scleromystaci, and four sequences from Gyrodactylus sp. from S. barbatus). Sequences of COII varied in length among species and suprapopulations of parasites: for G. bueni and G. major, 564 bp; for G. scleromystaci, 564 bp and 576 bp (for haplotypes from the sub-basins of the Paranaguá Bay and the Nhundiaquara River simultaneously); for Gyrodactylus sp. from S. barbatus, 564 bp (for a haplotype from the sub-basin of the Paranaguá Bay) and 570 bp (for haplotypes from the sub-basin of the Nhundiaquara River). All COII sequences analyzed herein present TTG start codon and TAA or TAG as stop codons. The TAG stop codon was found only in COII sequences of G. scleromystaci with 576 bp. In the analysis of the interspecific genetic variation, the COII sequences revealed 225 polymorphic sites, 278 mutations (113 synonymous substitutions, 51 non-synonymous substitutions and 90 sites with complex codons), and 16 sites with indels. In the analysis of the intraspecific genetic variation, the COII sequences of G. bueni showed 11 polymorphic sites, 11 mutations (seven synonymous substitutions and four non-synonymous substitutions) and average nucleotide differences between sequences (k) of 2.56; for G. major, 10 polymorphic sites, 10 mutations (three synonymous substitutions and seven non-synonymous substitutions) and an average number of nucleotide differences of 3.78; for G. scleromystaci from the sub-basin of the Paranaguá Bay, seven polymorphic sites, six mutations (four synonymous substitutions and two non-synonymous substitutions) and an average number of nucleotide differences of 2.29; and for G. scleromystaci from the sub-basin of the Nhundiaquara River, 16 polymorphic sites, 16 mutations (seven synonymous substitutions and nine non-synonymous substitutions) and an average number of nucleotide differences of 2.73.
The ZMP gene (intron and flanking exons) was sequenced for 55 specimens of Gyrodactylus (18 sequences from G. bueni, 16 sequences from G. major, and 21 sequences from G. scleromystaci). Nevertheless, no amplicons were obtained for Gyrodactylus sp. from S. barbatus. Sequences of ZMP varied in length among species of parasites: for G. bueni, 319 bp (200 bp in first exon, 35 bp in intron, and 84 bp in second exon); for G. major, 346 bp (200 bp in first exon, 59 bp in intron, and 87 bp in second exon); and for G. scleromystaci, 312 bp (200 bp in first exon, 28 bp in intron, and 84 bp in second exon). All ZMP sequences analyzed herein present an intronic region with typical dinucleotides GT (at the 5′ end of the intron) and AG (at the 3′ end of the intron). In the analysis of the interspecific genetic variation, the ZMP sequences revealed 70 polymorphic sites, 73 mutations (50 synonymous substitutions and 14 non-synonymous substitutions in the exons, and nine mutations in the intron), and 34 sites with indels. In the analysis of the intraspecific genetic variation, the ZMP sequences of G. bueni showed seven polymorphic sites, seven mutations (four synonymous substitutions and one non-synonymous substitution in the exons, and two mutations in the intron), and an average number of nucleotide differences of 1.64; for G. major, 10 polymorphic sites, 10 mutations (three synonymous substitutions and seven non-synonymous substitutions in the exons), and an average number of nucleotide differences of 0.58; for G. scleromystaci from the sub-basin of the Paranaguá Bay, two polymorphic sites, two mutations (two non-synonymous substitutions in the first exon), and an average number of nucleotide differences of 0.22; and for G. scleromystaci from the sub-basin of the Nhundiaquara River, nine polymorphic sites, nine mutations (six synonymous substitutions and three non-synonymous substitutions in the exons), and an average number of nucleotide differences of 1.77.
Phylogenetics and molecular evolution of Gyrodactylus
Phylogenetic relationships based on the COII and ZMP sequences of Gyrodactylus supported the identification of the species as recognized by morphology. In addition, molecular data allowed distinguishing between lineages of G. scleromystaci from the sub-basins of the Paranaguá Bay and the Nhundiaquara River (Figs 4 and 5). A total of 36 haplotypes were detected among the COII sequences of Gyrodactylus (Fig. 4). Eleven haplotypes were detected for G. bueni (H1–H11), of which two were exclusive to S. barbatus, and nine were exclusive to S. macropterus. Eight haplotypes were detected for G. major (H12–H19), of which one was exclusive to S. barbatus, and seven were exclusive to S. macropterus. Thirteen haplotypes were found for G. scleromystaci from S. barbatus (H20–H32), of which three were exclusive to suprapopulations of the sub-basin of the Paranaguá Bay, and 10 were exclusive to suprapopulations of the sub-basin of the Nhundiaquara River. Three haplotypes of G. scleromystaci from Pinto River (H21–H23) were grouped with haplotypes from Marumbi River (H20 and H24–H27), in the sub-basin of the Nhundiaquara River. Four haplotypes were found for Gyrodactylus sp. from S. barbatus (H33–H36), of which one was exclusive to the sub-basin of the Paranaguá Bay, and three were exclusive to the sub-basin of the Nhundiaquara River. The nucleotide divergence (d) of COII sequences between species of Gyrodactylus ranged from 22% to 54%. Mean nucleotide divergence of COII between suprapopulations of G. scleromystaci from distinct sub-basins (d = 5 ± 1%) was higher than the mean divergence values measured for suprapopulations of G. scleromystaci within the sub-basin of the Nhundiaquara River (d = 0.8 ± 0.2%) and the sub-basin of the Paranaguá Bay (d = 0.9 ± 0.3%). The DHEW test suggested that the genetic variation of the COII and ZMP sequences was not a result of positive selection (Table 2). The ratio (ω) of non-synonymous (Ka) and synonymous (Ks) substitution rates suggests that the COII gene is essentially evolving under purifying selection (ω < 1), including some sites under positive selection (ω > 1) (Supplementary Fig. S1). Likelihood ratio test (LRT) between the two codon substitution models (M8 and M8a) showed the interspecific positive selection was significant for species of Gyrodactylus from Scleromystax (P < 0.001). Although the positive selection was detected in some sites of COII for G. bueni and G. major, intraspecific positive selection was significant only for G. scleromystaci, even between haplotypes from the same sub-basin (P < 0.01) or from distinct sub-basins (P < 0.001). Interestingly, the number of sites evolving under positive selection on COII was different between suprapopulations of G. scleromystaci from distinct sub-basins (Supplementary Fig. S1). Among the ZMP sequences of Gyrodactylus, a total of 27 haplotypes were found (Fig. 5). Nine haplotypes were detected for G. bueni (H37–H45), of which one was found simultaneously on S. barbatus and S. macropterus, and eight were exclusive to S. macropterus. Ten haplotypes were detected for G. major (H46–H55), of which one was exclusive to S. barbatus, and nine were exclusive to S. macropterus. Eight haplotypes were found for G. scleromystaci from S. barbatus (H56–H63), of which two were exclusive to suprapopulations of the sub-basin of the Paranaguá Bay, and six were exclusive to suprapopulations of the sub-basin of the Nhundiaquara River. The interspecific genetic divergence of Gyrodactylus based on the ZMP sequences ranged from 2.5% to 17%, of which the lower value of genetic differences (2.5%) was found between G. bueni and G. scleromystaci. The genetic divergence of ZMP between suprapopulations of G. scleromystaci from distinct sub-basins (d = 2.3 ± 0.8%) was higher than the mean divergence values obtained for suprapopulations of G. scleromystaci within the sub-basin of the Nhundiaquara River (d = 0.6 ± 0.3%) and the sub-basin of the Paranaguá Bay (d = 0.7 ± 0.5%).
Fig. 4. Rooted phylogram from Bayesian analysis based on the mitochondrial cytochrome oxidase II (COII) gene from 36 haplotypes of four species of Gyrodactylus collected from Scleromystax barbatus (letters in grey) and Scleromystax macropterus (letters in black) from the sub-basin of the Nhundiaquara River (grey background) and the sub-basin of the Paranaguá Bay (no background), state of Paraná, Brazil. Values above branches correspond to posterior probabilities.
Fig. 5. Rooted phylogram from Bayesian analysis based on the nuclear zinc metalloproteinase (ZMP) gene from 27 haplotypes of three species of Gyrodactylus collected from Scleromystax barbatus (letters in grey) and Scleromystax macropterus (letters in black) from the sub-basin of the Nhundiaquara River (grey background) and the sub-basin of the Paranaguá Bay (no background), state of Paraná, Brazil. The asterisk (*) indicates haplotype of Gyrodactylus found simultaneously on S. barbatus and S. macropterus. Values above branches correspond to posterior probabilities.
Table 2. Estimates of the DHEW test performed for COII and ZMP sequences of Gyrodactylus from host species of Scleromystax
The values indicate Watterson's theta (θ w) and P value.
Reconstruction of divergence dates (the Time to Most Recent Common Ancestor – TMRCA) suggested that the haplogroup of G. scleromystaci is the oldest one, and its origin dates to approximately 2.3 Mya (95% HPD = 5.8–0.2 Mya). The origin dates of the remaining haplogroups were estimated as follows: haplogroup of G. bueni, 533 000 years ago (95% HPD = 1.3–0.07 Mya); haplogroup of G. major, 682 000 years ago (95% HPD = 1.7–0.07 Mya); and haplogroup of Gyrodactylus sp. from S. barbatus, 1.4 Mya (95% HPD = 3.5–0.12 Mya). Estimates of TMRCA for haplogroups of G. scleromystaci suggested that the suprapopulations from the sub-basin of the Nhundiaquara River originated about 624 000 years ago (95% HPD = 1.5–0.07 Mya), whereas suprapopulations of the sub-basin of the Paranaguá Bay originated approximately 530 000 years ago (95% HPD = 1.3–0.04 Mya).
The AMOVA based on the COII sequences revealed significant genetic structuring between suprapopulations of G. scleromystaci from distinct sub-basins (F st = 0.89; P < 0.001). The percentage of variation between the suprapopulations from different sub-basins (89%) was higher than within them (11%). Also, there was significant genetic structuring between suprapopulations of G. scleromystaci from distinct rivers, independently of its localities (F st = 0.89; P < 0.001). The percentage of variation between the suprapopulations from different rivers (89%) was higher than within them (11%). The AMOVA based on the ZMP sequences also revealed significant genetic structuring between suprapopulations of G. scleromystaci from distinct sub-basins (F st = 0.74; P < 0.01). The percentage of variation between the suprapopulations from different sub-basins (74%) was higher than within them (26%). Bayesian skyline plot analysis and mismatch distribution revealed a population expansion signal in the genetic data of species of Gyrodactylus (Fig. 6). All suprapopulations of parasites presented a similar demographic history, independently of host species or geographic locality, which suggests the occurrence of a common event that favoured simultaneously the population growth for distinct species of Gyrodactylus in southern coastal Brazil in the past 250 000 years between the Middle and Late Pleistocene, as well as induced a slight deceleration in population growth between the end of the Pleistocene and Holocene (Fig. 6). The R 2 test also detected a significant population expansion signal for species of Gyrodactylus (Table 3).
Fig. 6. Bayesian skyline plots and mismatch distribution for three species of Gyrodactylus collected from Scleromystax barbatus and Scleromystax macropterus inferred from nuclear ZMP sequences. (A) Gyrodactylus bueni. (B) Gyrodactylus major. (C) Gyrodactylus scleromystaci from the sub-basin of the Nhundiaquara River. (D) Gyrodactylus scleromystaci from the sub-basin of the Paranaguá Bay. Dark lines represent median inferred effective population size (Ne) in Bayesian skyline plots; grey lines indicate the 95% confidence intervals.
Table 3. Estimates of the R 2 test performed for COII sequences of Gyrodactylus from host species of Scleromystax
The values represent the average value of R 2, 95% confidence interval (CI), and P value.
a Sub-basin of the Nhundiaquara River.
b Sub-basin of the Paranaguá Bay.
The Mantel test indicated that there was little relationship between the degree of genetic and morphological variation between suprapopulations of G. scleromystaci (r 2 = 0.05; P < 0.001), suggesting that the differences in morphology are not simply due to neutral variation in genetically different suprapopulations. On the other hand, there was a significant positive correlation between the geographic distance and genetic variation between suprapopulations of G. scleromystaci (r 2 = 0.95; P < 0.001), indicating that geographical isolation has played a role in the genetic differentiation of parasites. Although there was a weak correlation between the degree of morphological variation and geographic distance, it was detected a significant positive correlation between geographic distance, morphological variation, and genetics for parasites (r 2 = 0.94; P < 0.001).
DNA sequences of Scleromystax
The mtDNA CR was sequenced for 46 specimens of Scleromystax, of which 15 were sequenced for S. macropterus and 31 were sequenced for S. barbatus (18 specimens from the sub-basin of the Paranaguá Bay and 13 specimens from the sub-basin of the Nhundiaquara River). The total length of the mtDNA CR sequences varied between host species (947 bp for sequences of S. macropterus and 1068–1289 bp for sequences of S. barbatus). In addition, there was variation in the total length of mtDNA CR sequences between subpopulations of S. barbatus. Sequences of S. barbatus from the sub-basin of the Paranaguá Bay revealed fragments with 1068–1289 bp, whereas sequences of specimens from the sub-basin of the Nhundiaquara River varied between 1086 bp and 1206 bp. These intra- and interpopulational differences in the length of mtDNA CR sequences resulted from tandem repeats of nucleotides in DNA fragments. Two patterns of tandem repeats were found for sequences of S. barbatus, which were referred herein as ‘pattern I’ (5′-ACATATATGTATGTACTAAATACATAATATGTATAATATT-3′), with 40 bp, and ‘pattern II’ (5′-TATGTATAACATCACATATATATA-3′), with 24 bp. Sequences of S. barbatus from the sub-basin of the Paranaguá Bay showed these two patterns of tandem repeats, whereas sequences of those specimens from the sub-basin of the Nhundiaquara River presented only tandem repeats of the pattern I. The number of copies of the pattern I in sequences of S. barbatus from the sub-basin of the Nhundiaquara River (6–9 copies) was different from that obtained from specimens of the sub-basin of the Paranaguá Bay (5–11 copies). On the other hand, there was no variation in the number of copies of the pattern II (two copies) between sequences of S. barbatus from the sub-basin of the Paranaguá Bay. All tandem repeats in DNA sequences were found closely at the 5′ end of the mtDNA CR, between the sixtieth (60th) and five hundredths (500th) sites. In the analysis of the interspecific genetic variation, the mtDNA CR sequences presented 107 polymorphic sites, 112 mutations, and 343 sites with indels. In the analysis of the intraspecific genetic variation, the mtDNA CR sequences of S. macropterus revealed 11 polymorphic sites, 11 mutations and an average number of nucleotide differences of 2.31; for sequences of S. barbatus from the sub-basin of the Paranaguá Bay, 25 polymorphic sites, 25 mutations, 221 sites with indels, and an average number of nucleotide differences of 5.49; and for sequences of S. barbatus from the sub-basin of the Nhundiaquara River, 21 polymorphic sites, 21 mutations, 121 sites with indels, and an average number of nucleotide differences of 5.14.
The SRP54 gene (intron and flanking exons) was sequenced for 56 specimens of Scleromystax, of which 20 were sequenced for S. macropterus and 36 were sequenced for S. barbatus (20 specimens from the sub-basin of the Paranaguá Bay and 16 specimens from the sub-basin of the Nhundiaquara River). All SRP54 sequences presented 253 bp, including two exons (85 bp in the first exon and 48 bp in the second one) and an intron with 120 bp. Moreover, these sequences presented an intronic region with characteristic dinucleotides GT and AG. In the analysis of the interspecific genetic variation, the SRP54 sequences revealed 15 polymorphic sites and 16 mutations (five synonymous substitutions and three non-synonymous substitutions in the first exon, and eight mutations in the intron). In the analysis of the intraspecific genetic variation, the SRP54 sequences of S. macropterus showed two polymorphic sites, two mutations (one synonymous substitution in the first exon and one mutation in the intron), and an average number of nucleotide differences of 0.84; for sequences of S. barbatus from the sub-basin of the Paranaguá Bay, seven polymorphic sites, eight mutations (four synonymous substitutions and two non-synonymous substitutions in the first exon, and two mutations in the intron), and an average number of nucleotide differences of 2.0; and for sequences of S. barbatus from the sub-basin of the Nhundiaquara River, six polymorphic sites, six mutations (three synonymous substitutions and one non-synonymous substitution in the first exon, and two mutations in the intron), and an average number of nucleotide differences of 1.73.
Phylogenetics and molecular evolution of Scleromystax
Phylogenetic reconstructions based on the mtDNA CR and SRP54 sequences of Scleromystax supported the taxonomic identity of the species and the recognition of lineages of S. barbatus from sub-basins of the Paranaguá Bay and the Nhundiaquara River (Figs 8 and 9). A total of 33 haplotypes were recognized in the mtDNA CR sequences of Scleromystax (Fig. 7). Seven haplotypes were found for S. macropterus (H64–H70) and 26 haplotypes for S. barbatus (H71–H96), of which nine haplotypes were exclusive to S. barbatus from the sub-basin of the Nhundiaquara River, 16 haplotypes were exclusive to S. barbatus from the sub-basin of the Paranaguá Bay, and one haplotype was found simultaneously in distinct sub-basins. The mean nucleotide divergence of mtDNA CR sequences between S. barbatus and S. macropterus was 11 ± 2%. Mean nucleotide divergence of mtDNA CR between subpopulations of S. barbatus from different sub-basins (d = 0.8 ± 0.1%) was higher than the mean value found for subpopulations of S. barbatus within the sub-basin of the Nhundiaquara River (d = 0.5 ± 0.1%), but a few smaller than the mean nucleotide divergence obtained within subpopulations of the sub-basin of the Paranaguá Bay (d = 0.9 ± 0.2%). Among the SRP54 sequences (Fig. 8), four haplotypes were detected for S. macropterus (H97–H100) and 11 haplotypes for S. barbatus (H101–H111), of which five haplotypes were exclusive to S. barbatus from the sub-basin of the Paranaguá Bay, one haplotype was exclusive to S. barbatus from the sub-basin of the Nhundiaquara River, and five haplotypes were detected in both the sub-basins. The interspecific mean genetic divergence of SRP54 sequences between S. barbatus and S. macropterus was 3 ± 1%. The intraspecific mean nucleotide divergence of SRP54 between subpopulations of S. barbatus from distinct sub-basins (d = 1 ± 0.3%) was slightly higher than the mean values measured for subpopulations within the sub-basin of the Nhundiaquara River (d = 0.7 ± 0.3%) and the sub-basin of the Paranaguá Bay (d = 0.8 ± 0.4%). The DHEW test suggested that the pattern of genetic variation of the mtDNA CR sequences was consistent with positive selection for S. barbatus but marginally not for S. macropterus, however, no evidence of positive selection was detected in the SRP54 sequences for both host species (Table 4). Additionally, the DHEW test performed separately for subpopulations of S. barbatus showed significant positive selection within subpopulations from the sub-basin of the Paranaguá Bay (θ w = 4.2; P < 0.001) while, conversely, no signal of positive selection was identified within populations from the sub-basin of the Nhundiaquara River (θ w = 4.64; P = 0.6).
Fig. 7. Rooted phylogram from Bayesian analysis based on the mitochondrial control region from 33 haplotypes of Scleromystax barbatus (letters in grey) and Scleromystax macropterus (letters in black) from the sub-basin of the Nhundiaquara River (grey background) and the sub-basin of the Paranaguá Bay (no background), state of Paraná, Brazil. The asterisk (*) indicates haplotype of Scleromystax collected simultaneously from both sub-basins. Values above branches correspond to posterior probabilities.
Fig. 8. Rooted phylogram from Bayesian analysis based on the signal recognition particle 54 kDa gene (SRP54) from 15 haplotypes of Scleromystax barbatus (letters in grey) and Scleromystax macropterus (letters in black) from the sub-basin of the Nhundiaquara River (grey background) and the sub-basin of the Paranaguá Bay (no background), state of Paraná, Brazil. The asterisk (*) indicates haplotype of Scleromystax found simultaneously in both sub-basins. Values above branches correspond to posterior probabilities.
Fig. 9. Bayesian skyline plots and mismatch distribution for host species Scleromystax barbatus and Scleromystax macropterus inferred from nuclear SRP54 sequences. (A) Scleromystax barbatus from the sub-basin of the Nhundiaquara River. (B) Scleromystax barbatus from the sub-basin of the Paranaguá Bay. (C) Scleromystax macropterus. Dark lines represent median inferred effective population size (Ne) in Bayesian skyline plots; grey lines indicate the 95% confidence intervals.
Table 4. Estimates of the DHEW test performed for CR and SRP54 sequences of Scleromystax barbatus and Scleromystax macropterus
The values indicate Watterson's theta (θ w) and P value.
*Significant P value (P < 0.01).
Reconstruction of divergence dates (TMRCA) indicated that the divergence of the two haplogroups of S. barbatus from distinct sub-basins occurred about 1.2 Mya (95% HPD = 3.3–0.08 Mya). The origin dates for these two haplogroups were estimated as 799 000 years ago (95% HPD = 2.2–0.06 Mya) for haplotypes from the sub-basin of the Nhundiaquara River, and 1.1 Mya (95% HPD = 2.9–0.06 Mya) for haplotypes from the sub-basin of the Paranaguá Bay. Estimates of TMRCA for S. macropterus suggested that haplotypes originated approximately 588 000 years ago (95% HPD = 1.6–0.03 Mya).
The AMOVA based on the mtDNA CR sequences indicated low genetic structuring between subpopulations of S. barbatus from distinct sub-basins (F st = 0.09; P < 0.05). The percentage of variation between the subpopulations from different sub-basins (9%) was lower than within them (91%). On the other hand, the AMOVA revealed high genetic structuring (F st = 0.26; P < 0.001) between subpopulations of S. barbatus from different rivers (i.e. when sub-basins were not taken into account). The percentage of variation between the subpopulations from distinct rivers (26%) was lower than within them (74%). The AMOVA based on the SRP54 sequences showed no evidence of genetic structuring between subpopulations of S. barbatus, however when shared haplotypes were excluded from the analysis, a high genetic structuring (F st = 0.63; P < 0.001) was detected between subpopulations of S. barbatus from different sub-basins. In this analysis, the percentage of variation between the subpopulations from distinct sub-basins (63%) was higher than within them (37%). Also, there was high genetic structuring (F st = 0.59; P < 0.001) between subpopulations of S. barbatus from different rivers. The percentage of variation between the subpopulations from distinct rivers (59%) was higher than within them (41%). Bayesian skyline plot analysis and mismatch distribution indicated a population expansion signal in the genetic data of S. barbatus and S. macropterus (Fig. 9). For both host species, a population growth signal was detected over the last 250 000 years between the Middle and Late Pleistocene followed by a slight signal of deceleration in population expansion between the end of the Pleistocene and Holocene (Fig. 9). The R 2 test also presented a significant population expansion signal for species of Scleromystax (Table 5).
Table 5. Estimates of the R 2 test performed for SRP54 sequences of Scleromystax barbatus and Scleromystax macropterus
The values represent the average value of R 2, 95% confidence interval (CI), and P value.
a Sub-basin of the Nhundiaquara River.
b Sub-basin of the Paranaguá Bay.
According to Mantel test, there was no correlation between the degree of genetic and morphological variation between subpopulations of S. barbatus (r 2 = 0.002; P > 0.05), suggesting that the interpopulational morphologic differences could not be explained only by neutral variation in genetically distinct subpopulations. Geographic distance was weak but significantly correlated with genetic variation between subpopulations of S. barbatus (r 2 = 0.10; P < 0.001), suggesting the influence of historic and environmental factors on the genetic diversification of populations. Despite there was little correlation between morphological variation and geographic distance, no correlation was found between morphological variance, geographic distance and genetics for hosts (r 2 = 0.005; P > 0.05).
Distribution of parasites between host species
The optimization of host species (S. barbatus and S. macropterus) on the Bayesian cladogram inferred from COII sequences by the Parsimonious Ancestral State Reconstruction method of Mesquite 2.7.4 resulted in four dispersal events, which were significantly different from the distribution of host transfer events obtained from a casual association (P < 0.001). Randomization of the host-parasite association with subsequent reconstruction of character states onto the Bayesian cladogram resulted in an average of 22 host transfer events with a confidence interval of 95% of 20 to 28 events. Gyrodactylus bueni and G. major were highly associated with S. macropterus, with exception of a few parasite specimens (n = 3) found on S. barbatus. Conversely, G. scleromystaci was almost exclusively associated with S. barbatus, with exception of three parasite specimens collected from S. macropterus. All specimens of Gyrodactylus sp. were found only on S. barbatus.
Discussion
Morphometric analysis on the haptoral sclerites of Gyrodactylus revealed significant differences in anchors and hooks between suprapopulations of G. scleromystaci from distinct rivers and sub-basins, which is consistent with population-specific morphological variation. Nonetheless, not all parasites could be assigned to their respective geographic locality with absolute accuracy based only on morphometrics, as opposed to molecular data, which reflected perfectly the geographical distribution of haplotypes of G. scleromystaci from distinct sub-basins. This relative incongruence between molecular and morphologic data suggests that phenotypic plasticity in anchors and hooks of Gyrodactylus may be related to environmental and biotic factors (Malmberg, Reference Malmberg1970; Mo, Reference Mo1991a, Reference Mo1991b; Harris, Reference Harris1998; Dmitrieva and Dimitrov, Reference Dmitrieva and Dimitrov2002; Dávidová et al., Reference Dávidová, Jarkovský, Matejusová and Gelnar2005; Robertsen et al., Reference Robertsen, Hansen, Bachmann and Bakke2007; Bueno-Silva et al., Reference Bueno-Silva, Boeger and Pie2011). As reported by some authors (Olstad et al., Reference Olstad, Bachmann and Bakke2009; Bueno-Silva et al., Reference Bueno-Silva, Boeger and Pie2011), host species may also exert influence on the morphometrical variance of anchors and hooks in species of Gyrodactylus, suggesting the adaptation of parasites to its hosts. Therefore, the morphometrical variance of haptoral sclerites seems not to be associated simply to geographic distance and genetic divergence between suprapopulations of G. scleromystaci but to a combination of genetic and environmental components, which seems to drive local patterns of development. Anchors were the most variable structure in size and contributed significantly to distinguish between species and suprapopulations of G. scleromystaci from different geographic localities. Five out of eight variables that allowed delimitation of species of Gyrodactylus were the same that accounted for the morphometric variation between suprapopulations of G. scleromystaci, of which three measurements were taken from hooks: hooklet aperture (HAD), hooklet distal width (HDW), and hook shank length (HSHL). This result reinforces the significant role of hooks in morphological differentiation of species and suprapopulations of Gyrodactylus, as previously reported by other authors (Malmberg, Reference Malmberg1970; Dmitrieva and Dimitrov, Reference Dmitrieva and Dimitrov2002; Bueno-Silva and Boeger, Reference Bueno-Silva and Boeger2009; Bueno-Silva et al., Reference Bueno-Silva, Boeger and Pie2011).
Morphometric analysis on the measurements of the host species also revealed significant differences in subpopulations of S. barbatus between distinct rivers and sub-basins. This result suggests that the morphological divergence found between subpopulations of S. barbatus may be due to phenotypic plasticity related to historic and environmental factors, as it is known that geographic isolation and environmental conditions play an important role on the morphological variability of fish populations (Strauss, Reference Strauss1985; Robinson and Wilson, Reference Robinson and Wilson1994; Smith and Skúlason, Reference Smith and Skúlason1996; Langerhans et al., Reference Langerhans, Layman, Langerhans and Dewitt2003; Shibatta and Hoffmann, Reference Shibatta and Hoffmann2005). Based on the morphometric analysis, Shibatta and Hoffmann (Reference Shibatta and Hoffmann2005) were able to distinguish upland subpopulations of Corydoras paleatus (a closely related species to Scleromystax) from different hydrographic basins in the state of Paraná, southern Brazil. On the other hand, Takahashi and Koblmüller (Reference Takahashi and Koblmüller2011) argued that differences in body shape of cichlid fish might be due to random genetic drift, which can evolve phenotypes without adaptation. In the present study, it was postulated that subpopulations of S. barbatus were geographically isolated because of marine transgressions on the coast of Paraná (see further discussion below). Marine transgressions possibly confined freshwater fish to distinct areas, which may have reduced populations of S. barbatus. Under this scenario, one would expect that the morphological variation in S. barbatus could be explained by genetic drift associated with a reduction of subpopulations which were isolated during past marine transgressions. Although there is a genetic component for morphological variance in fish (Robinson and Wilson, Reference Robinson and Wilson1994; Smith and Skúlason, Reference Smith and Skúlason1996), the present study suggests that historic and local environmental factors may have contributed to phenotypic variations between subpopulations of S. barbatus, since morphometric variance was little related to geographic distance and genetic divergence of the subpopulations. Environmental differences detected between collection sites could induce some morphological variation on S. barbatus. According to Guimarães et al. (Reference Guimarães, Menezes and Peret2010), species of Siluriformes, including S. barbatus, exhibit habitat preference, since they mainly occur in waters with a predominance of rocks and pebbles. In the present study, the two sample points from the Marumbi and Pinto River (sub-basin of the Nhundiaquara River) were mainly composed by rocks and pebbles, whereas sampling sites from Fortuna and Ribeirão River (sub-basin of the Paranaguá Bay) were mostly composed by sand. In addition, it is known that the pH and flow of water vary between sampled sites. For example, the Marumbi and Pinto River present alkaline clear water (e.g. pH 7.4–7.7; see Horodesky et al., Reference Horodesky, Abilhoa, Zeni, Neto, Castilho-Westphal and Ostrensky2015) and sampling sites with moderate to fast current, while rivers that flow at lower altitudes in the coastal plain of Paraná present very slow current and acidic dark water (e.g. pH 6.5) due to the high concentration of organic matter (Bigarella, Reference Bigarella2001), as is the case with Fortuna River. Given that variations in pH and flow of water may induce morphological differentiation on fish (see Trippel and Harvey, Reference Trippel and Harvey1987; Imre et al., Reference Imre, McLaughlin and Noakes2002; Reddon and Hurd, Reference Reddon and Hurd2013), it is possible to hypothesize that different environmental conditions may also have led to some variation in morphology of S. barbatus between distinct localities. Some studies have reported significant morphometric variation between subpopulations of catfish from distinct geographical localities (Turan et al., Reference Turan, Yalçin, Turan, Okur and Akyurt2005; Gunawickrama, Reference Gunawickrama2007; Solomon et al., Reference Solomon, Okomoda and Ogbenyikwu2015; Chaklader et al., Reference Chaklader, Siddik, Hanif, Nahar, Mahmud and Piria2016). Remarkably, some authors have noticed that the head measurements (e.g. head length, predorsal distance, and snout length) present significant differences between subpopulations of catfish (Turan et al., Reference Turan, Yalçin, Turan, Okur and Akyurt2005; Gunawickrama, Reference Gunawickrama2007; Chaklader et al., Reference Chaklader, Siddik, Hanif, Nahar, Mahmud and Piria2016). In the present study, head measurements contributed to distinguishing between species and subpopulations of Scleromystax. Also, five out of six variables that allowed delimitation of species of Scleromystax were the same that accounted for the morphometric variation between subpopulations of S. barbatus from different geographic locations, indicating that these measurements are relevant for morphological differentiation of fish populations, i.e. depth of caudal peduncle (DCP), standard length (SL), length of dorsal-fin base (LDB), depth of body (DB), and length of anal-fin base (LAB). Therefore, differences in morphology between subpopulations of S. barbatus seem to reflect population fragmentation and local patterns of development caused by local environmental conditions.
The morphometric variation in G. scleromystaci is comparable to the findings by Agrawal (Reference Agrawal2001) and Thompson (Reference Thompson2005), who stated the morphological variance between different subpopulations because of local adaptation in the host-parasite relationship (e.g. influence of environmental factors and spatial or temporal differences in the host-parasite association). In accordance with Agrawal (Reference Agrawal2001), the evolution of adaptive phenotypic plasticity has led to the success of organisms in novel habitats and potentially contributes to genetic differentiation and speciation. Nevertheless, the morphological diversification seems to be more related to local ecological factors (e.g. local patterns of development) and historical processes (e.g. population fragmentation and demographic changes) rather than the geographic distance between subpopulations.
The Bayesian phylogenetic analysis based on the COII and ZMP sequences allowed recognizing of clades of parasite species, as previously identified by morphology, as well as lineages of G. scleromystaci corresponding to different geographic localities. For parasites, the genetic differentiation was related to the geographic distance between suprapopulations. This result is consistent with the isolation-by-distance model, which suggests the influence of geographic structuring on the genetic variation of local lineages of Gyrodactylus. Unexpectedly, three haplotypes of G. scleromystaci from Pinto River were grouped with haplotypes from Marumbi River, in the sub-basin of the Nhundiaquara River. Two scenarios could explain this result: (1) the alluvial plain of the municipality of Morretes is subjected to periodic floodings (Bigarella, Reference Bigarella2001), which could result in the exchange of haplotypes of G. scleromystaci between different rivers; and (2) these are the oldest haplotypes of G. scleromystaci, which could have expanded and colonized distinct rivers from the same sub-basin. For hosts, the Bayesian phylogenetic reconstruction based on the mtDNA CR and SRP54 sequences confirmed the taxonomic identity of species of Scleromystax and detected haplotypes of S. barbatus derived from distinct sub-basins. In contrast to parasites, the genetic differentiation was poorly related to the geographic distance between subpopulations of S. barbatus. While some authors have found a positive correlation between genetic variation and geographic distance among fish populations (Raeymaekers et al., Reference Raeymaekers, Maes, Audenaert and Volckaert2005; Crispo et al., Reference Crispo, Bentzen, Reznick, Kinnison and Hendry2006; Crookes and Shaw, Reference Crookes and Shaw2016), other detected a lack of relationship between genetics and geography (Hänfling et al., Reference Hänfling, Durka and Brandl2004; Caldera and Bolnick, Reference Caldera and Bolnick2008). The data obtained by Caldera and Bolnick (Reference Caldera and Bolnick2008) suggest that isolation by distance played less a role in genetic structure than physical barriers. These findings give support to the results obtained herein since there is no connectivity between distant subpopulations of Scleromystax due to geographic barriers. In the present study, the genetic variation could be only partially explained by the isolation-by-distance model, which suggests that historic and local environmental factors (e.g. population fragmentation and physical barriers) may have contributed more significantly to genetic differentiation of Scleromystax than the geographic distance between subpopulations.
Unlike the AMOVA calculated for hosts, the genetic data based on mtDNA and nDNA showed a high genetic structuring between suprapopulations of G. scleromystaci from distinct geographic localities. Such divergences between molecular data of parasites and hosts may be related to biological characteristics of Gyrodactylus, which result in high levels of genetic diversification, such as hyperviviparity, a short generation time (about four days), and dispersion as an adult throughout its lifecycle (Cable and Harris, Reference Cable and Harris2002; Bakke et al., Reference Bakke, Harris and Cable2002, Reference Bakke, Cable and Harris2007; Boeger et al., Reference Boeger, Kritsky and Pie2003).
On the other hand, the AMOVA estimated from mtDNA CR and SRP54 sequences of S. barbatus indicates the simultaneous occurrence of exclusive and shared haplotypes between subpopulations from distinct sub-basins, which would suggest partial genetic structuring for hosts. Nonetheless, the hypothesis of active gene flow between distant subpopulations was rejected based on the topographic and geomorphological characteristics of the coastal basin of the state of Paraná (Ab'Saber, Reference Ab'Saber2000; Bigarella, Reference Bigarella2001; Maack, Reference Maack2001), since the formation of geographic barriers would prevent gene flow. Given that the formation of the coastal plain of the state of Paraná is geologically recent, between the Late Pleistocene and Holocene (Maack, Reference Maack2001; Angulo et al., Reference Angulo, Pessenda and Souza2002; Angulo, Reference Angulo2004), one would expect that the exclusive haplotypes of S. barbatus would possibly be the result of recent post-colonization genetic differentiation, as it is known that historic events (e.g. post-glacial colonization) may exert influence on the genetic structure and variability of freshwater fish (Nesbø et al., Reference Nesbø, Fossheim, Vøllestad and Jakobsen1999; Kotlík and Berrebi, Reference Kotlík and Berrebi2001; Gagnon and Angers, Reference Gagnon and Angers2006; Caldera and Bolnick, Reference Caldera and Bolnick2008; Walter et al., Reference Walter, Blum, Snider, Paterson, Bentzen, Lamphere and Gilliam2011). Nevertheless, a comparative calibration of molecular clocks performed in the present study rejects the hypothesis that subpopulations were recently isolated, since results revealed unrealistic overestimated rates of molecular evolution for subpopulations hypothetically isolated since the Late Pleistocene up to after the last Holocene maximum transgression, at approximately 5100 years ago (data not shown). This result suggests that subpopulations of S. barbatus from distinct localities were possibly isolated before the Late Pleistocene. A scenario of Early Pleistocene isolation of subpopulations is proposed herein based on the following: (1) a historical reconstruction on the origin of freshwater fish species from the coast of Paraná suggests that S. barbatus is present in the Paranaguá Bay at more than 1 Mya (Tschá, Reference Tschá2016); (2) the Pliocene Epoch was a time of generalized epeirogenic elevation of the entire eastern region of Brazil, which resulted in an intense erosive process and consequent sedimentation in the coastal regions – during this epoch, a tectonic pulsation affected the Serra do Mar system, and this process was accompanied by an accentuation of the relief of mountains and rifts which at that time assumed their modern aspects (Almeida, Reference Almeida1976); (3) since some authors have registered occurrence of S. barbatus at different altitudes on the coastal basin of Paraná, ranging from coastal plain to elevated terrains at approximately 80 m of altitude (this study; Guimarães et al., Reference Guimarães, Menezes and Peret2010; Wolff and Hahn, Reference Wolff and Hahn2017), it would be expected that subpopulations of S. barbatus would be confined to elevated areas during past marine transgressions (see Suguio and Martin, Reference Suguio and Martin1978; Ab'Saber, Reference Ab'Saber2000; Siddall et al., Reference Siddall, Rohling, Almogi-Labin, Hemleben, Meischner, Schmelzer and Smeed2003); and (4) during past marine transgressions on the coast of Paraná, elevated terrains (10 m above sea level) may have acted as ‘islands’, in some way, analogous to barrier-islands (see Suguio and Martin, Reference Suguio and Martin1978; Martin et al., Reference Martin, Suguio, Flexor, Dominguez and Bittencourt1996), by isolating subpopulations of S. barbatus from distinct localities. At the same time, the occurrence of shared haplotypes of S. barbatus between distant geographic locations indicates that distinct subpopulations were founded by a common ancestral population. The occurrence of shared haplotypes of S. barbatus between different sub-basins was not expected because there is no connectivity between the sub-basin of the Nhundiaquara River and the sub-basin of the Paranaguá Bay. Nonetheless, according to Ab'Saber (Reference Ab'Saber2000), the sea level suffered a regressive decrease in southern coastal Brazil until approximately 100 m below the present elevation between 23 000 and 12 700 years ago. Under this scenario, it would be expected that freshwater fish populations would be able to expand and colonize new habitats which were previously submerged along the coastal plain of Paraná. In the same way, Tschá et al. (Reference Tschá, Baggio, Marteleto, Abilhoa, Bachmann and Boeger2017) suggested that low sea-levels during the Pleistocene may have promoted the expansion of coastal freshwater basins of Paraná, which may have caused an expansion of isolated populations of freshwater fishes, including S. barbatus, and it is most likely that distinct subpopulations were combined into a large panmictic group. Therefore, the presence of shared haplotypes of S. barbatus between distinct sub-basins may be due to past low sea-level periods in the coastal plain of Paraná. On the other hand, the incongruence in the number of shared haplotypes inferred from mtDNA CR and SRP54 sequences could be due to differences in the evolutionary rate between mitochondrial and nuclear markers, since the mtDNA exhibits extensive intraspecific polymorphism and evolves faster than nDNA (Brown et al., Reference Brown, George and Wilson1979; Avise, Reference Avise2000).
The genetic variation of the mtDNA CR appears to be driven by a combination of positive selection (particularly for S. barbatus) and population-specific length variation in the sequences. Length variations in the mtDNA CR sequences could be due to: (1) insertions/deletions in DNA; and (2) occurrence of tandem repeats in DNA. Such variations are known to occur in various species of fish (Árnason and Rand, Reference Árnason and Rand1992; Faber and Stepien, Reference Faber and Stepien1998; Ludwig et al., Reference Ludwig, May, Debus and Jenneckens2000; Hoarau et al., Reference Hoarau, Holla, Lescasse, Stam and Olsen2002; Chen et al., Reference Chen, Ablan, McManus, Bell, Tuan, Cabanban and Shao2004; Ponce et al., Reference Ponce, Infante, Jiménez-Cantizano, Pérez and Manchado2008; Jamandre et al., Reference Jamandre, Durand and Tzeng2014; Satoh et al., Reference Satoh, Miya, Mabuchi and Nishida2016). In addition, Lee et al. (Reference Lee, Conroy, Howell and Kocher1995) discussed that variations in the mtDNA CR may have little effect on its functionality. The pattern of tandem repeats in mtDNA CR sequences of S. barbatus is consistent with the findings of Lee et al. (Reference Lee, Conroy, Howell and Kocher1995), since selective pressure acts to maintain insertions/deletions near at the start (5′) or end (3′) position of the mtDNA CR in order to avoid structural modifications in the DNA-protein interaction sites, which are situated at the central domain of the control region. Comparative analysis on the mtDNA CR sequences of S. barbatus revealed that tandem repeats in DNA may be useful to phylogenetic inference, allowing distinguishing between subpopulations from distinct geographic localities, as reported by Faber and Stepien (Reference Faber and Stepien1998) on fish species of Perciformes. Moreover, the pattern of genetic variation estimated from mtDNA CR sequences of S. barbatus and S. macropterus corroborates the observations by Sbisà et al. (Reference Sbisà, Tanzariello, Reyes, Pesole and Saccone1997) that insertions/deletions of different sizes makes this region subject to mechanisms of DNA expansion, which may be very useful to evaluate the genetic diversification and geographic distribution of fish populations. The DHEW test revealed that the genetic variation in the mtDNA CR was consistent with positive selection. This result was unexpected because it is known that selection mostly acts in the mtDNA coding regions (see Castellana et al., Reference Castellana, Vicario and Saccone2011; Fischer et al., Reference Fischer, Koblmüller, Gülly, Schlötterer, Sturmbauer and Thallinger2013; Consuegra et al., Reference Consuegra, John, Verspoor and Leaniz2015; James et al., Reference James, Piganeau and Eyre-Walker2016). Evidence for positive selection was detected in animal mtDNA genes, such as NADH dehydrogenase genes, ATP8, and Cytb (see Shen et al., Reference Shen, Liang, Zhu, Zhou, Irwin and Zhang2010; Garvin et al., Reference Garvin, Bielawski and Gharrett2011; Melo-Ferreira et al., Reference Melo-Ferreira, Vilela, Fonseca, Fonseca, Boursot and Alves2014; Consuegra et al., Reference Consuegra, John, Verspoor and Leaniz2015). Such positive selection in mtDNA may be associated with local adaptation of fish lineages between distinct geographic localities (Consuegra et al., Reference Consuegra, John, Verspoor and Leaniz2015). Despite the lack of evidence for positive selection in the mtDNA CR (Hardouin and Tautz, Reference Hardouin and Tautz2013), Bazin et al. (Reference Bazin, Glémin and Galtier2006) and James et al. (Reference James, Piganeau and Eyre-Walker2016) found evidence that mtDNA undergoes substantial amounts of adaptive evolution. Although it is not clear here whether genetic variation of the mtDNA CR has resulted only from natural selection or hitchhiking effect, it is possible to hypothesize that inter- and intrapopulational positive selection inferred from mtDNA CR sequences could be related to some adaptation of local lineages of S. barbatus between different geographic localities because of recent post-glacial colonization. However, further studies will be necessary to elucidate the evolution of the mtDNA CR in Scleromystax.
The low genetic structuring found between subpopulations of S. barbatus from different sub-basins may be related to environmental and historical factors which are believed to influence the genetic structure and diversity of fish populations (Montoya-Burgos, Reference Montoya-Burgos2003; Comesaña et al., Reference Comesaña, Martínez-Areal and Sanjuan2008; Rodrigues et al., Reference Rodrigues, Schneider, Santos, Vallinoto, Sain-Paul and Sampaio2008; von der Heyden et al., Reference von der Heyden, Lipinski and Matthee2010). As the mtDNA CR presents a highly variable rate of evolution among fish (see McMillan and Palumbi, Reference McMillan and Palumbi1997; Montoya-Burgos, Reference Montoya-Burgos2003; Liu et al., Reference Liu, Gao, Yokogawa and Zhang2006), it is necessary to take caution when mtDNA CR sequences reveal little or absence of genetic structuring between geographically isolated fish subpopulations (sensu von der Heyden et al., Reference von der Heyden, Lipinski and Matthee2010). In accordance with Bazin et al. (Reference Bazin, Glémin and Galtier2006), mitochondrial lineages with low diversity might sometimes correspond to recently selected, well-adapted haplotypes. Considering that subpopulations of S. barbatus from distinct sub-basins are totally separated by a mean distance of 32 km, it is possible to presume that the low genetic divergence between distant subpopulations of S. barbatus (0.7–0.9%) could be due to historical events, including changes in population size (sensu Zhang et al., Reference Zhang, Cai and Huang2006; Rodrigues et al., Reference Rodrigues, Schneider, Santos, Vallinoto, Sain-Paul and Sampaio2008) or recent diversification of populations (sensu Bazin et al., Reference Bazin, Glémin and Galtier2006). Despite the changes in population size of S. barbatus (see further discussion below) and the Pleistocene diversification of subpopulations, it seems that a past connection between isolated subpopulations of S. barbatus, as suggested by Tschá et al. (Reference Tschá, Baggio, Marteleto, Abilhoa, Bachmann and Boeger2017), may have contributed mostly to formation of a large panmictic population, which could explain the low genetic structuring found between subpopulations of S. barbatus from different sub-basins.
Comparative analysis of the demographic history of populations revealed congruencies between demographic fluctuations of Gyrodactylus and Scleromystax. A population expansion signal was detected for all sampled species of Gyrodactylus and Scleromystax, which occurred approximately 250 000 years ago. These fluctuations in population size of parasites and hosts could be related to environmental effects of glacial-interglacial cycles between the Middle Pleistocene and part of Holocene (5000– 250 000 years ago; sensu Imbrie et al., Reference Imbrie, Boyle, Clemens, Duffy, Howard, Kukla, Kutzbach, Martinson, McIntyre, Mix, Molfino, Morley, Peterson, Pisias, Prell, Raytoo, Shackletons and Toggweile1992; Petit et al., Reference Petit, Jouzel, Raynaud, Barkov, Barnola, Basile, Bender, Chappellaz, Davisk, Delaygue, Delmotte, Kotlyakov, Legrand, Lipenkov, Lorius, Pépin, Ritz, Saltzmank and Stievenard1999; Glennon et al., Reference Glennon, Perkins, Chisholm and Whittington2008; Hickerson et al., Reference Hickerson, Carstens, Cavender-Bares, Crandall, Graham, Johnson, Rissler, Victoriano and Yoder2010). Past sea-level variations seem to be related to such glacial-interglacial cycles (see Siddall et al., Reference Siddall, Rohling, Almogi-Labin, Hemleben, Meischner, Schmelzer and Smeed2003). Some authors discuss that extreme climate changes at the Pleistocene Epoch may have exerted a strong influence on the distribution and demographic fluctuation of fish (Hewitt, Reference Hewitt2000; Lecomte et al., Reference Lecomte, Grant, Dodson, Rodríguez-Sánchez and Bowen2004; Liu et al., Reference Liu, Gao, Yokogawa and Zhang2006). Ribeiro (Reference Ribeiro2006) argued that marine transgressions facilitated the dispersion of several species of fish between different geographic areas of the coast of Brazil. In fact, paleoenvironmental variations were important in defining the morphology of the coast of southern Brazil, including marine transgressions and regressions (Suguio and Martin, Reference Suguio and Martin1978; Angulo and Lessa, Reference Angulo and Lessa1997; Ab'Saber, Reference Ab'Saber2000; Angulo et al., Reference Angulo, Pessenda and Souza2002). According to Ab'Saber (Reference Ab'Saber2000), sea level reached a maximum of approximately 10 m above the present elevation between 130 000 and 23 000 years ago in southern Brazil. In this period, the coast of the state of Paraná was submerged by seawater, including its sub-basins. This sea-level rise was followed by a regressive decreasing until approximately 100 m below the present elevation between 23 000 and 12 700 years ago, which originated coastal plains with streams and medium-sized rivers (Ab'Saber, Reference Ab'Saber2000). According to Angulo and Lessa (Reference Angulo and Lessa1997), a post-glacial marine transgression increased the sea level to approximately 3.5 m above the present elevation in Paranaguá region at 5100 years ago. After this period, the sea level declined until the present elevation. This regressive decreasing of sea level in the last 5000 years is consistent with the stabilization of sea level for Southern Hemisphere at the Holocene (Dawson, Reference Dawson1992; Angulo and Lessa, Reference Angulo and Lessa1997; Angulo et al., Reference Angulo, Pessenda and Souza2002).
Given that S. barbatus and S. macropterus are exclusively found in coastal rivers from the south and southeast coastal Brazil (Reis, Reference Reis, Reis, Kullander and Ferraris2003; Ferraris, Reference Ferraris2007), it may be possible to hypothesize that these fish species were under influence of sea-level fluctuations during the Quaternary. The fact that different populations of Gyrodactylus and Scleromystax presented similar demographic history suggests thereby the common influence of paleoenvironmental variations on population expansions and retractions, mainly due a regressive decreasing of sea level at the last 130 000 years in southern Brazil, which possibly allowed the populations to expand by colonizing new habitats in the coastal plain of the state of Paraná. In the same way, an earlier marine regression that occurred approximately between 150 000 and 200 000 years ago (see Siddall et al., Reference Siddall, Rohling, Almogi-Labin, Hemleben, Meischner, Schmelzer and Smeed2003) could possibly have contributed to population expansions. On the other hand, it was detected a recent slight signal of deceleration in population growth of Gyrodactylus and Scleromystax, which suggests that a sea-level rise between the end of the Pleistocene and Holocene (see Siddall et al., Reference Siddall, Rohling, Almogi-Labin, Hemleben, Meischner, Schmelzer and Smeed2003) may have reduced the size of host-parasite populations. These findings suggest that environmental factors potentially played an important role in the population fluctuations of parasites and hosts simultaneously. Hence, local geomorphological characteristics and past eustatic sea-level changes may have caused effects on the historical demography of parasites and its hosts. Similarly, Tschá et al. (Reference Tschá, Baggio, Marteleto, Abilhoa, Bachmann and Boeger2017) argued that past sea-level variations may have exerted influence in the historical demography of freshwater fish in the coast of Paraná. However, estimates of population expansion of S. barbatus obtained by Tschá et al. (Reference Tschá, Baggio, Marteleto, Abilhoa, Bachmann and Boeger2017) – less than 50 000 years ago – were different from that of the present study (250 000 years ago). This may be due to fact that Tschá et al. (Reference Tschá, Baggio, Marteleto, Abilhoa, Bachmann and Boeger2017) reconstructed historical demography of fishes by using strict molecular clock and a single rate of evolution, simultaneously, for marine and freshwater fish species, while in the present study relaxed molecular clocks were calibrated based on the genetic divergence between local subpopulations of S. barbatus and the geomorphological history of the coast of Paraná.
The expansion of the suprapopulations of Gyrodactylus seems to reflect an increase of synonymous and non-synonymous substitutions in the COII gene, which suggest that the diversification of these parasites may be related to evolutionary radiation, as previously discussed by Boeger et al. (Reference Boeger, Kritsky and Pie2003). According to Kita et al. (Reference Kita, Hirawake and Takamiya1997), parasitic helminths exhibit wide diversity in energy metabolism (via cytochrome oxidase), which plays an important role in the adaptation of parasites to their natural habitats. Furthermore, Kita et al. (Reference Kita, Hirawake and Takamiya1997) commented that some mutations in amino acid sequences of cytochrome oxidase could exert influence on the degree of parasite resistance to certain inhibitors of the cell respiratory chain. Although there are no known studies to date that assess effects of genetic mutations in cytochrome oxidase of Monogenoidea (e.g. fitness and functional consequences), the molecular data by Rawson and Burton (Reference Rawson and Burton2002) provide evidence for evolution and adaptation in intertidal copepods based on COII sequences. These authors revealed that some mutations in amino acid sequences of COII were under purifying and positive selection. Such substitutions could modify the specificity of the respective protein of cytochrome between different populations of copepods, which led to the geographic diversification of these populations. In the present study, the number of sites evolving under positive selection on COII was different between suprapopulations of G. scleromystaci from distinct sub-basins, which is consistent with population-specific selection. Bankers et al. (Reference Bankers, Fields, McElroy, Boore, Logsdon and Neiman2017) found evidence for population-specific genetic responses to parasite infection in the Potamopyrgus-Microphallus interaction between three distinct localities, which represented the first genome-level evidence for the type of population-specific response expected under local adaptation in this snail-trematode system. These findings seem to offer support for the local adaptation of G. scleromystaci. Moreover, according to Storz (Reference Storz2005), when species are distributed across ecologically heterogeneous landscapes, spatially separated populations that occupy distinct environments may be subject to different selection modes. As an adaptive mutation is driven to fixation by positive selection (Storz, Reference Storz2005), it may be possible to hypothesize that the positive selection on COII may reflect local adaptation of suprapopulations of Gyrodactylus between hosts from different geographic localities, but further studies are necessary to address this question.
Additionally, the topologies inferred from COII and mtDNA CR sequences are consistent with rapid divergence event, i.e. when lineages rapidly diverged from a common ancestor during an evolutionary short time period (Hayakawa et al., Reference Hayakawa, Culleton, Otani, Horii and Tanabe2008). Rapid divergence events have been associated with the diversification of fish species due to the influence of geological events, climate changes, and population isolation (see Beheregaray et al., Reference Beheregaray, Sunnucks and Briscoe2002; Allender et al., Reference Allender, Seehausen, Knight, Turner and Maclean2003; Seehausen, Reference Seehausen2006; Jamandre et al., Reference Jamandre, Durand and Tzeng2014). The study carried out by Beheregaray et al. (Reference Beheregaray, Sunnucks and Briscoe2002) sustained that sea-level changes at the Pleistocene and Holocene may be related with rapid radiation of a freshwater fish species in the coastal plain of southern Brazil. These findings seem to corroborate the results obtained herein. Thus, the geographic structuring has played a role in the genetic diversification of host-parasite populations, which may have favoured events of evolutionary co-divergence.
The most parsimonious reconstruction of the host species on the Bayesian phylogeny of parasites indicated a high host specificity of Gyrodactylus spp. from S. barbatus and S. macropterus, however the data suggest that local differences in the host-parasite relationship seem to be mostly explained by host transfer events (e.g. the presence of closely related host species in syntopy). This finding suggests that the host-parasite interaction could differ between subpopulations from different geographic localities if local conditions facilitate host transfer events, as postulated by the Stockholm Paradigm (Araujo et al., Reference Araujo, Braga, Brooks, Agosta, Hoberg, von Hartenthal and Boeger2015; Hoberg and Brooks, Reference Hoberg and Brooks2015). Although host transfer events be common in species of Gyrodactylus (see Ziętara and Lumme, Reference Ziętara and Lumme2002; Boeger et al., Reference Boeger, Kritsky and Pie2003; Huyse et al., Reference Huyse, Audenaert and Volckaert2003; Meinilä et al., Reference Meinilä, Kuusela, Ziętara and Lumme2004; Bueno-Silva et al., Reference Bueno-Silva, Boeger and Pie2011), G. bueni and G. major presented high host specificity to S. macropterus, independently of the occurrence of S. barbatus in syntopy. Based on other studies on Gyrodactylus (see Boeger et al., Reference Boeger, Kritsky and Pie2003, Reference Boeger, Kritsky, Pie and Engers2005; Bueno-Silva et al., Reference Bueno-Silva, Boeger and Pie2011), it could be suggested that the low incidence of G. bueni and G. major in S. barbatus may be related to casual transfer events of parasites between syntopic host species.
Species of Callichthyidae generally exhibit social behaviour and are grouped in shoals (Paxton, Reference Paxton1997; Kaatz and Lobel, Reference Kaatz and Lobel1999), which may facilitate the transmission of parasites between distinct host species, mainly when fish are grouped in mixed shoals (as is the case with S. barbatus and S. macropterus). On the other hand, when S. barbatus and S. macropterus were not found in syntopy, populations of S. barbatus harboured exclusively two parasite species, G. scleromystaci and Gyrodactylus sp., regardless of the geographic location of host populations. These results corroborate the findings by Bakke et al. (Reference Bakke, Harris and Cable2002, Reference Bakke, Cable and Harris2007), which suggest that some species of Gyrodactylus can exhibit high specificity to its host even when the parasites have an opportunity to exploit other close host species in sympatry. On the other hand, the occurrence of closely related host species does not necessarily imply that parasites may exploit them since host-parasite interaction is influenced by ecological factors, such as opportunity and compatibility (see Araujo et al., Reference Araujo, Braga, Brooks, Agosta, Hoberg, von Hartenthal and Boeger2015; Hoberg and Brooks, Reference Hoberg and Brooks2015). According to this view, the high abundance of a parasite species on a single host species does not necessarily mean that this parasite is highly host-specific, as the absence of alternative compatible hosts may reduce the opportunity for parasites to exploit new host species and expand to new geographic areas. Thus, the interaction between species of Gyrodactylus and Scleromystax seems to be defined by a combination of host compatibility and casual host transfer events.
Regarding the results of the present study, it is possible to conclude that: (i) Local variances in host-parasite interaction seem to be mostly influenced by casual host transfer events; (ii) The geographic structuring is an important factor on the genetic diversification of parasites and its hosts, which may have favoured events of evolutionary co-divergence as a response to common climatological and environmental processes; (iii) Evidence of population-specific selection in DNA sequences suggests local adaptation of lineages of parasites and hosts; (iv) The morphological diversification of parasites and hosts may be influenced more by historic and local environmental factors (e.g. population fragmentation and local patterns of development) than isolation by distance; (v) The genetic variation of local lineages of Gyrodactylus is consistent with the isolation-by-distance model while, conversely, the genetic variation of Scleromystax could be partially explained by isolation by distance, which implies that local environmental factors (e.g. geographic barriers and population fragmentation) may have contributed more significantly to genetic differentiation than geographic distance between subpopulations of hosts; (vi) The molecular diversification of Gyrodactylus and Scleromystax seems to reflect a rapid evolutionary divergence of host-parasite populations between distinct geographic localities, which may have occurred in the Early Pleistocene; (vii) The molecular data of Gyrodactylus and Scleromystax revealed a congruent historical demography for parasites and its hosts, which suggest that these populations were influenced by common paleoenvironmental variations, at least, in the last 250 000 years in southern Brazil; (viii) Marine regressions that occurred in the last 200 000 years may have contributed to population expansions, especially the regressive decreasing of sea level at the last 130 000 years in southern Brazil, which possibly allowed the populations to expand by colonizing new habitats in the coastal plain of the state of Paraná; (ix) A recent slight signal of deceleration in population growth of Gyrodactylus and Scleromystax suggests that a sea-level rise between the end of the Pleistocene and Holocene may have reduced the size of host-parasite populations.
Author ORCIDs
Marlus Bueno-Silva, 0000-0002-9215-3797.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182019000556
Acknowledgements
We wish to thank the staff of the Laboratório de Ecologia Molecular e Parasitologia Evolutiva (LEMPE/UFPR), Dr Rogério T. Vianna (FURG) and Dr Jean Vitule (UFPR) for field and laboratory support, and the Programa de Pós-graduação em Zoologia (PPGZOO/UFPR).
Financial support
The Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq, Brazil) provided financial support for this study and granted support to M.B.S. (grant number 140946/2007-1) and W.A.B. (grant number 471165/2009-3).
Conflicts of interest
None.
Ethical standards
The Instituto Brasileiro do Meio Ambiente e dos Recursos Naturais Renováveis (IBAMA, Brazil) granted the licence for the collection of biological material (licence 087/2004, 10811/2007 and 11617/2008-2010).
