Genetic connectivity and population expansion inferred from multilocus analysis in Lutjanus alexandrei (Lutjanidae–Perciformes), an endemic snapper from Northeastern Brazilian coast

Sampling

Samples from 120 specimens of L. alexandrei were obtained between 2013 and 2018. Some of these samples were kindly provided by colleagues from muscle tissues stored at their facilities. The remaining samples were obtained from dead individuals, procured through direct purchase from specimens locally caught and commercially sold by artisanal fishermen. Therefore, there was no need to apply for a license for collection or approval by the Animal Ethics Committee. All samples were collected/transported with authorization of the Instituto Chico Mendes de Conservação da Biodiversidade (ICMBio) (SISBIO license n. 62250-1).

The analyzed specimens were obtained from four localities in Northeastern Brazilian coast: Salvador (Bahia–BA) (n = 14; all individuals were juveniles and were captured near beaches and provided by artisanal fishermen, in 2018); Aracaju (Sergipe–SE) (n = 28; all individuals were juveniles and sampled in 2018, most were obtained from fish markets; three specimens were captured by artisanal fishermen in the Vaza Barris river estuary); Tamandaré (Pernambuco–PE) (n = 31; all samples were provided from the analysis material of Rodrigo Bastos, a postdoctoral student supervised by professor Beatrice Padovani. The fish were captured in 2017, in the estuary of the Ariquindá River, on the coastal reefs, in Vau de Tamandaré and Sirinhaém beach) and Nísia Floresta (Rio Grande do Norte–RN) (n = 48; most individuals were juveniles and were collected in the estuary of the Pium River or in coastal reefs of the Praia de Pirangi, by artisanal fishermen, in 2013, 2014 and 2018. Some adult specimens were sampled in markets on Ponta Negra Beach in 2014) (Fig. 1). Their taxonomic identification was carried out using the specialized literature (Moura & Lindeman, 2007) and the samples of muscle tissue were tagged and stored in microtube with 100% ethanol at −20 °C up to laboratory procedures. Voucher specimens were deposited in the ichthyological collection from the Laboratory of Applied Genetics at Universidade Federal do Pará, Campus of Bragança, Pará.

Laboratory procedures

Total DNA was isolated using the commercial Wizard Genomic® kit (Promega, Madison, WI, USA), according to the manufacturer’s instructions. The quality of the extract DNA was verified by electrophoresis in 1% agarose gel stained with Gel Red™ (Biotium, Fremont, CA, USA) and exposed to UV light. The amplification of the fragments of interest was performed via PCR (polymerase chain reaction). The analyzed sequences comprised mitochondrial DNA (mtDNA) regions—Cytochrome b (Cyt b), NADH dehydrogenase subunit 4 (ND4) and Control Region (CR)—and nuclear DNA genes (intron 1 of the Myostatin and intron 1 of the S7 ribosomal protein).

The PCRs comprised 2.4 μL of dNTPs (1.25 mM), 1.5 μL of buffer (10x), 0.6 μL of MgCl₂ (50 mM), 0.6 μL of each primer (50 ng/μL), 0.1 μL of Taq DNA polymerase (U/μL), about 50 ng of template DNA and ultrapure water to a final volume of 15 μL. The primers used to amplify each fragment of interest are described in Table 1.

The PCR conditions encompassed a first denaturation step at 94 °C for 3 min; followed by 35 cycles of denaturation at 94 °C for 40 s, annealing at specific temperature for each primer for 40 s (Table 1), and extension at 72 °C for 2 min plus a final extension step at 72 °C for 2 min.

The amplicons were purified in PEG 8000 (polyethylene glycol) according to Paithankar & Prasad (1991) and sequenced by the dideoxy method (Sanger, Nicklen & Coulson, 1977) using the Big Dye kit (ABI Prism™ Dye Terminator Cycle Sequencing Reading Reaction; Thermo Fisher, Waltham, MA, USA). The final products were read in ABI 3500 (Thermo Fisher, Waltham, MA, USA) automatic sequencer.

Dataset

The sequences were organized into individual datasets for each marker and visualized, edited and corrected in the software BioEdit (Hall, 1999). Each sequence dataset was aligned using CLUSTAL W (Thompson, Higgins & Gibson, 1994), available in Bioedit.

Heterozygous individuals were detected when double peaks were identified on both strands in the electropherograms. The gametic phase was determined using Phase v. 2.1 (Stephens, Smith & Donnelly, 2001), available in DNAsp v. 5.10 (Librado & Rozas, 2009). Three runs with distinct random seeds, 1,000 iterations per chain, thinning interval equal to 1, and a burn-in of 1,000 were carried out. Only the haplotypes with probability values above 0.59 were selected for the further analyses.

The occurrence of intragene recombination was estimated by the PhiW test (Bruen, Philippe & Bryant, 2006) using the software Splits Tree v. 4.6 (Huson & Bryant, 2006). Whenever present, a dataset free of recombination was obtained using the software Imgc (Woerner, Cox & Hammer, 2007).

Population structure

The genetic polymorphism in populations of L. alexandrei was estimated based on the values of haplotype (h) and nucleotide (π) diversity (Nei, 1987) using the software Arlequin v. 3.5 (Excoffier & Lischer, 2010).

To infer about the spatial distribution of haplotypes and their genealogical relationships, haplotype networks were built for each genomic region in the software Haploviewer (Salzburger, Ewing & Von, 2011). For that, we used as input a tree topology inferred from Maximum Likelihood, available in the DNAml package and implemented in Phylip 3.6 (Felsenstein, 2005).

The analysis of molecular variance (AMOVA) (Excoffier & Lischer, 2010) was carried out to test the population subdivision of L. alexandrei along their range based on the frequency and number of mutations among haplotypes (Excoffier & Lischer, 2010). The AMOVA was performed considering the samples from distinct localities as a single group. The significance of fixation indexes was estimated using 10,000 non-parametric permutations. The level of genetic differentiation among the sampled localities was estimated from pairwise FST values as calculated by Weir & Hill (2002) and the significance of this result was based on 10,000 bootstrap permutations assuming a significance level of 0.05.

The population structure was also evaluated using a Bayesian clustering analysis in the software Bayesian analysis of genetic population structure (BAPS) v. 6.0 (Corander et al., 2013). This method allows identifying the genetic groups independently on the geographic distribution of samples. In this study, BAPS was performed assuming an admixture model among specimens from all collection sites with linked loci and up to 20 putative clusters (k = 20).

Neutrality tests and demographic history

The neutrality of the analyzed loci was evaluated using Tajima’s (1989) D and Fu’s (1997) Fs tests in Arlequin v. 3. 5 (Excoffier & Lischer, 2010) based on 10,000 permutations for estimating the p significance values.

The Extended Bayesian Skyline Plot (EBSP) (Heled & Drummond, 2008) available in the software Beast v.1.8 (Drummond et al., 2012) was also used to evaluate the population dynamics through time, allowing integrating multilocus data. Each run comprised 500 million MCMC iterations, sampled at every 10,000 iterations, with a burn-in of 10%.

The best evolutionary models were estimated using the software PartitionFinder 2.1 (Lanfear et al., 2016), as follows: HKY+I+G (CR) and TRN+I+G (Cyt b; ND4; Myostatin and S7). We assumed a relaxed clock with an uncorrelated lognormal distribution and found that the markers’ ucld.stdev values were close to zero, indicating that the uniform replacement rate hypothesis was not rejected. Then strict clock and evolutionary rate of 3.6% per million years for CR (Donaldson & Wilson, 1999; da Silva et al., 2020); 1% for Cyt b and ND4 (Bowen et al., 2001; Lessios, 2008; Gaither et al., 2010; Delrieu-Trottin et al., 2017; Da Cunha et al., 2020); 0.4% for S7 intron I (Janko et al., 2007; Delrieu-Trottin et al., 2017); and for Myostatin we estimate the rate. Was employed in the software Beast v.1.8 using the Cipres Science Gateway v. 3.3 platform (Miller, Pfeiffer & Schwartz, 2010) available in http://www.phylo.org/index.php/portal/cite_us. After finishing the run, ESS >200 values were verified in TRACER 1.7.1 (Rambaut et al., 2018) and the graph was generated using the Python script provided by Heled (2010).

Characterization of the dataset and genetic diversity

Five different fragments from mtDNA and from nuclear DNA were sequenced in the present study, totaling 2,997 base pairs (bp), being 2,084 bp related to mtDNA and 913 pb from nuclear DNA.

The verification of the occurrence of intragenic recombination showed that, in contrast to the myostatin intron 1, recombination was detected for the fragment of the intron 1 of S-7 ribosomal protein.

As for the mitochondrial fragments, for CR (mtDNA), a bank of sequences with a length of 806 bp was obtained for 111 specimens (accession numbers: MZ327145–MZ327255). From this total, we recovered 110 haplotypes and high levels of genetic diversity across the sampled area which can be observed by the values obtained for haplotype diversity (h = 1.0 for the entire range of samples) and for nucleotide diversity (π ranging from 0.033 in Bahia to π = 0.041 in Rio Grande do Norte). For the CytB, sequences of 808 bp were obtained for 110 specimens (accession numbers: MZ343890–MZ343999), comprising 25 haplotypes and haplotype diversity (h) ranging from 0.1818 in Bahia to h = 0.8376 in Pernambuco while the nucleotide diversity (π) varied from 0.001 and 0.002 for the same localities. The ND4 fragment, woth470 bp for 117 individuals (accession numbers: MZ356036–MZ356152), encompassed 23 haplotypes and h values from 0.6527 for the samples from Pernambuco to 0.9103 in Bahia, locations in which nucleotide diversity for this fragment ranged from π = 0.002 to π = 0.003, respectively (Table 2).

With regard to nuclear genome fragments, for intron I S7 a total of 36 haplotypes were recovered from the 506 bp sequenced for 70 specimens (accession numbers: MZ355926–MZ356035), in addition to allelic diversity values above 0.9 for all populations and nucleotide diversity ranging from π = 0.006 in Bahia to π = 0.008 in Pernambuco. At last, the myostatin intron 1 comprised sequences of 407 bp for 112 specimens (accession numbers: MZ344000–MZ344114), distributed into nine haplotypes. The h values for this fragment varied from 0.5404 to 0.7101, while the π values varied from 0.002 to 0.003 for the samples from Rio Grande do Norte and Bahia, respectively (Table 2).

Population genetic structure

The haplotype network inferred from both mitochondrial and nuclear markers revealed a remarkable admixture pattern with homogeneous distribution across the sampled area in Northeastern Brazil. For the CR, a single haplotype was shared by two specimens while the others represented unique haplotypes. A total of 25 haplotypes were recovered using Cyt B with the three most frequent haplotypes being shared by 50, 18 and 10 specimens, respectively, while the remaining haplotypes showed low frequency or were unique, similarly to the pattern observed in ND4. Three haplotypes were also the most frequent ones for the Myostatin intron 1, being present in 131, 56 and 20 specimens, respectively, from all populations sampled in this study. The haplotype network for the intron 1 of S7 ribosomal protein revealed a more patched scenario, with two frequent haplotypes shared among entre 16 and 12 specimens, while the others were observed at low frequencies (Fig. 2).

This lack of population subdivision among the sampled populations was confirmed by AMOVA (Table 3), since most of genetic variation was observed within populations, as also reinforced by F/θST values. In addition, these results were supported by BAPS inasmuch as a single cluster (K = 1) presented the maximum probability value.

Finally, the pairwise FST values among the localities corroborated the homogeneity among populations of L. alexandrei, as revealed by low and non-significant values for all studied fragments (Table 4).

Neutrality and demographic history

Both neutrality tests (Tajima’s D and Fu’s Fs) showed negative and significant values for most of populations and considering all fragments included in this study. Regarding the mitochondrial DNA fragments included in this study, for the Control Region, Tajima’s D values ranged from −0.88204 in Rio Grande do Norte to 1.00193 in Bahia, while for Fu’s Fs values ranged from −22.48052 in Rio Grande do Norte and −1.71508 for the same locations. For Cytochrome b, the variation was from −1.59303 in Sergipe to −1.42961 in Bahia, and from −8.04383 in Rio Grande do Norte to 0.50645 in Bahia, for Tajima’s D and Fu’s Fs, respectively. The ND4 fragment analyzed reported a variation from −1.60572 in Rio Grande do Norte to −1.07663 in Bahia in Tajima’s D values and in relation to the Fu’s Fs test, the variation was in the range of −6.86688 in Rio Grande do Norte at −3.20917 in Sergipe (Table 2).

As for nuclear fragments, Intron 1 of the S7 ribosomal protein showed Tajima’s D values between −0.67670 in Sergipe and 0.54201 in Pernambuco, while the variation in the Fu’s Fs test was between −17.39503 and −3.65941, for Rio Grande do Norte and Bahia. Finally, Tajima’s D values found for Intron 1 of the Myostatin comprised the range of −0.34041 to −0.05216 for Sergipe and Bahia, respectively; for Fu’s Fs the variation ranged from −2.16023 in Sergipe to −0.50836 in Rio Grande do Norte (Table 2).

The demographic history of L. alexandrei throughout the Northeastern coast of Brazil was investigated using the EBSP approach. The multilocus analysis resulted in a skyline plot indicating increases of the effective population size during Pleistocene, around 800 thousand years ago (Fig. 3).

Genetic diversity

Nuclear introns are expected to be more suitable for phylogenetic analysis than population-genetics, due to the slow rate of evolution, among other issues. However, to obtain a more complete and realistic scenario, the inclusion of polymorphic sequences is necessary, considering that the use of mitochondrial markers and microsatellite sequences, widely used and recognized as useful, have limitations that indicate the need to use other genomics regions to complement their results, leading fragments of nuclear DNA to be applied to studies of population genetics, with introns being the best candidates (Zhang & Hewitt, 2003).

In general, the analyses of the genomic regions in this study revealed high levels of genetic variation in L. alexandrei, particularly in relation to CR. This pattern has been commonly reported in snappers from the Brazilian coast, such as L. purpureus (h = 0.997; π = 0.029) (da Silva et al., 2016), O. chrysurus (h = 0.963; π = 0.0179) (da Silva et al., 2015), L. analis (h = 0.9962; π = 0.0315) and L. jocu (h = 0.9962; π = 0.,0362) (Souza et al., 2019). Nonetheless, the diversity herein reported in L. alexandrei is higher than that observed in the abovementioned species.

In the case of snappers, highly targeted by commercial fisheries, such high values of genetic diversity are unexpected since several overexploited marine species have shown a reduction of genetic variation (Pinsky & Palumbi, 2014). On the other hand, the fisheries pressure on Lutjanidae is a recent process, dating back about 50 years ago (Fonteles-Filho, 1972). Therefore, the negative effects of this activity might remain undetected by genetic analyses even if populations of snappers have been impacted by overexploitation. In addition, some fragments such as the CR are recognized by indicating diversity losses related to past events, such as historic reduction of effective population size, as reported for other marine fishes like Merluccius paradoxus (Heyden, Lipinski & Matthee, 2010) and other snappers from Northeastern Brazil and sympatric to L. alexandrei, like L. synagris (Silva et al., 2018). Consequently, these data are not suitable to diagnosis of pressure by fisheries since this is a recent activity as mentioned earlier (Evangelista-Gomes et al., 2020).

Genetic connectivity

The present results indicated that L. alexandrei is composed of a single population, a trait putatively associated with the species-specific biological features as proposed for other snappers from Northeastern coast of Brazil, such as L. purpureus. In the latter, the genetic homogeneity was hypothetically explained by the presence of pelagic eggs and long period of larval development, as also reported in other Lutjanidae species from Atlantic Ocean (da Silva et al., 2016). Likewise, O. chrysurus, a closely related species to Lutjanus, distributed in sympatry with several species of the genus, is also composed of a panmictic population along the Brazilian coast, what has been associated with ecological traits such as the duration of larval stages (da Silva et al., 2015).

The type of eggs and larval periods, favored by the synergetic action of current tides, are features that could account for the panmixia in both L. purpureus and O. chrysurus from the Brazilian coast (da Silva et al., 2015, 2016). Most likely, the currents along the Brazilian coast could also potentialize the dispersal of larvae of L. alexandrei, a scenario particularly favored by the apparent lack of barriers to population subdivision of these snappers, as also suggested for other marine fishes (Palumbi, 1994).

Besides the features related to type of eggs and larvae and the role of Brazilian currents, the genetic connectivity found in L. alexandrei might be also influenced by the close relationship of this species with estuaries. Juveniles of Brazilian snapper usually inhabit coastal areas such as mangroves while adults are frequently observed in deeper reef environments (Aschenbrenner & Ferreira, 2015). This migration to both types of habitats might increase the admixture pattern reported in L. alexandrei, as also inferred for L. jocu from Northeastern Brazil and sympatric to L. alexandrei (Martins, 2018).

The hypotheses about the causes of the genetic connectivity within L. alexandrei could be more precise if other studies about the life traits of this species were available inasmuch as reliable phylogeographic inferences invariably rely on biological and ecological aspects of the analyzed organisms (Bowen et al., 2014). Unfortunately, there is no information whether the duration of the larval stages in L. alexandrei is extended or not, as previously reported in other sympatric snapper species like L. analis, L. synagris, L. jocu and Rhomboplites aurorubens (Lessa, Nóbrega & Bezerra, 2004; Freitas et al., 2011), or else if this species presents great dispersal abilities that could explain the high admixture among haplotypes.

Demographic history

The negative values obtained for most populations by the neutrality tests based on the fragments used to compare the samples of L. alexandrei indicate deviations from neutrality, what could be related to population expansion events. This suggestion was corroborated by the results from EBSP inasmuch as a 10-fold increase in Ne (effective population size) of L. alexandrei has taken place about 800 thousand years ago.

According to Rabassa, Coronato & Salemme (2005) and Rabassa, Coronato & Martínez (2011), this period corresponds to the Middle Pleistocene, culminating in a glaciation. These authors highlighted that climate changes in this period influence the development of ecosystems in South America.

A study carried out by Costa et al. (2019) about the demographic changes in Neotropical fishes from Pleistocene revealed that the sea level along the Brazilian coast during the last glacial maximum was 120 m below the present level, thus exposing large areas of the continental shelf and putatively increasing the availability of brackish and freshwater habitats. Therefore, freshwater and estuarine species would benefit from new colonization sites and increased resources, thus leading to increased effective population sizes.

Partially similar results in relation to the present data were reported in other snappers that co-occur with L. alexandrei. Events of population expansion in L. jocu and L. analis were estimated between 36 to 55 thousand years ago and 13 to 61 thousand years ago, respectively (Souza et al., 2019). In the yellowtail snapper (O. chrysurus) from the Brazilian coast, a sudden increase of population size was also inferred during Pleistocene (da Silva et al., 2015), similarly to the pattern reported for L. purpureus (Gomes, Sampaio & Schneider, 2012; da Silva et al., 2016) and L. synagris (Silva et al., 2018).

Although other snappers distributed sympatrically with L. alexandrei also experienced population growth during the Pleistocene, these showed a much more recent population expansion than L. alexandrei. This is probably due to the fact that most studies reporting this increase in the effective number were based only on mitochondrial DNA (Souza et al., 2019; Gomes, Sampaio & Schneider, 2012) or used a mutation rate of 10% per million years (da Silva et al., 2015, 2016; Silva et al., 2018).

EBSP has its limitations, mainly in relation to the use of mutation rates between different molecular fragments, which can sometimes make it difficult to calibrate (Heled & Drummond, 2008), however this method has been widely used to estimate the demographic history of fish, including snappers (da Silva et al., 2015, 2016; Silva et al., 2018).