Genetic and morphological differentiation of the Lusitanian toadfish ( Halobatrachus didactylus ) between estuarine and coastal areas in Portugal

The Lusitanian toadfish, Halobatrachus didactylus (Bloch and Schneider, 1801) (Batrachoididae) is a sedentary benthic species that lives down to about 50 m, inhabiting the eastern Atlantic from Ghana to the Iberian Peninsula (about 39oN) (Roux, 1986). On the Portuguese coast this species is particularly abundant in estuaries south of Sado (38oN) (Sobral and Gomes, 1997) and has a peculiar distribution: SCIENTIA MARINA 70(4) December 2006, 749-758, Barcelona (Spain) ISSN: 0214-8358


INTRODUCTION
The Lusitanian toadfish, Halobatrachus didactylus (Bloch and Schneider, 1801) (Batrachoididae) is a sedentary benthic species that lives down to about 50 m, inhabiting the eastern Atlantic from Ghana to the Iberian Peninsula (about 39ºN) (Roux, 1986).On the Portuguese coast this species is particularly abundant in estuaries south of Sado (38ºN) (Sobral and Gomes, 1997) and has a peculiar distribution: on the western coast it occurs only in estuaries but on the south coast it occurs in both estuaries and coastal areas (Costa, 1993).Together with its distribution, this species' life history also suggests a low flux of individuals, even between close areas, since males nest under rocks or in crevices defending the clutch (Santos et al., 2000).This is expected to result in a differentiation by genetic drift, inbreeding or, as suggested by Beheregaray and Sunnucks (2001) for species that inhabit coastal and estuarine regions, "divergence-with-gene-flow".
Like H. didactylus, the widely distributed sand goby, Pomatoschistus minutus (Pallas, 1770), inhabits both estuarine and coastal areas and adults have little or no migration.Research on the genetic diversity and differentiation of this species using allozymes suggested that the morphological differentiation observed in the lagoon of Venice was related to population structuring (Stefanni et al., 1996;2003).Studies conducted at four locations in the North Sea using microsatellite markers also revealed differentiation between estuarine, coastal and marine samples of P. minutus (Pampoulie et al., 2004).Within the family Batrachoididae population sub-structuring has been described for Opsanus tau Linnaeus, 1758 and Opasanus beta Goode and Bean, 1879 inhabiting each side of cape Hatteras on the western Atlantic coast (Avise et al., 1987).Considerable genotypic diversity and differentiation was observed within each species due to the constrained gene flux between populations.
In studies considering population differentiation the selection of marine areas where major biogeographic shifts occur are as important as species' life history and ecology (Nielsen et al., 2004).The Portuguese coast is therefore a suitable area for the study of population interactions due to its transitional character as a border area between the warm-temperate and cool-temperate Atlantic zoogeographic areas.Here, Castilho and McAndrew (1998) found a significant population sub-structuring for Dicentrarchus labrax (Linnaeus, 1758) juveniles in five estuaries.However, for two species of soles with a similar life history pattern, Solea solea (Linnaeus, 1758) and Solea senegalensis Kaup, 1858, Cabral et al. (2003a)  2002), the differentiation patterns depend on whether the species only use them as a special habitat during their life cycle or they actually support distinct populations.
In the present study, morphological and enzymatic characters of H. didactylus inhabiting 10 areas along the Portuguese coast were used in order to: 1) identify and quantify the phenotypic and genotypic variability; 2) analyse the morphological and genetic differentiation of H. didactylus; and 3) infer the importance of estuaries in population structuring.

Sampling procedures
A total of 466 individuals of H. didactylus were sampled from 10 areas along the Portuguese coast (Fig. 1).Sampling sites on the western coast (the Tejo, Sado and Mira estuaries) were widely separated, whereas most of those on the south coast comprised estuaries and their adjacent coastal areas.This set of samples allowed the detection of regional and local patterns of differentiation as well as the comparison between estuarine and marine environments.Beam trawl fishing was carried out during the non-matting period, between September and April 2000, and all individuals from each sample were caught in the same trawl.Liver and muscle tissue samples were taken from all specimens and stored at -80ºC until electrophoretic analysis.

Morphological analysis
Characters easily repeatable from individual to individual, including distances between clearly recognisable landmarks of the species anatomy, were selected for the morphological analysis of the 10 samples.This resulted in a total of 20 morphometric-including horizontal as well as vertical dimensions of the body-and 16 meristic characters (Table 1).Analyses of these variables were carried out separately because they are different in their statistical nature and might therefore be responding differently to environmental and genetic factors (Ihssen et al., 1981).
All morphometric characters, measured to the nearest 0.1 mm, were transformed to logarithms to approximate multivariate normality.Given that body measurements are generally strongly correlat-ed with the individuals' total length (TL), morphometric characters were standardised to the overall mean total length using Hurlbut and Clay's (1998) expression: in which TL is the total length, M x is the original measurement, TL -is the overall mean total length and b is the slope, within areas, of the geometric mean regression (Ricker, 1973) on the logarithms of M x and TL.This regression model was chosen because any of the variables could be considered as depending on another.The resulting variables were then considered shape discriminators (Humphries et al., 1981).After size effect removal, variables remaining highly correlated were considered redundant and eliminated from the analysis.Meristic  characters have been reported as independent of fish size (e.g.Tudela, 1999;Murta, 2000) and were therefore not standardised.
For each variable, differences between males and females were evaluated by t or Mann-Whitney tests (morphometric and meristic variables, respectively).Individuals were considered as multivariate observations in order to account for the joint effect of variables (Misra and Carscadden, 1987).A stepwise multivariate discriminant analysis was performed separately for morphometric and meristic data in order to identify the combinations of variables that best separate the samples of H. didactylus (Hair et al., 1996).The effectiveness of the discriminant analysis was evaluated by means of the percentage of individuals correctly classified in the original sample and by a one-way ANOVA followed by a posteriori Bonferroni test procedures.
Morphological distances between samples were computed as Euclidean distances and their correlation with geographic distances was evaluated through a Mantel test.
All calculations were performed using SPSS (SPSS Inc.) and 0.05 as the level of significance.

Genetic analysis
Samples of muscle and liver tissue were hydrated and then disrupted using ultrasound prior to horizontal starch gel electrophoresis.The technical procedures for allozyme analysis followed Murphy et al. (1996) and the histochemical staining methods were adapted from Harris and Hopkinson (1976).Alleles were designated by their electrophoretic mobility relative to the anodal mobility of the most common allele, which was designated as 100.From the wide range of enzyme stains tested, only 10 (corresponding to 11 loci) showed variation between samples and were therefore assayed.This electrophoretic variation that was observed conformed to the simple Mendelian model of co-dominance and to the known quaternary structure of each enzyme (Pasteur et al., 1985;Alayse, 1987).
The proportion of polymorphic loci, mean number of alleles per locus (MNA) and observed (H O ) and expected (H E ) heterozygizities were calculated in GENETIX 4.02 (Belkhir et al., 1996(Belkhir et al., -2001)).A locus was considered polymorphic when the frequency of its most common allele did not exceed 0.95.
Wright's F-statistics and Nei's genetic distance (Nei, 1978) were used to analyse genetic differenti-ation and computed in GENETIX 4.02.Wright's Fstatistics were estimated according to Weir and Cockerham (1984) and their significance was tested using permutations (1000 replicates).F IT , F IS and F ST measure the departure from Hardy-Weinberg proportions at the level of the whole sample, individual samples and differentiation between samples, respectively (Wright, 1965).Confidence intervals (95%) for multilocus F ST estimates were calculated by bootstrapping over loci.Deviations from Hardy-Weinberg equilibrium (HWE) were also tested according to the Marckov-chain method in GENEPOP 3.1 (Raymond and Rousset, 1995).
The model of isolation by distance (IBD) was evaluated by the correlation between genetic and geographical distances using the Mantel test (Mantel, 1967) as implemented in GENETIX 4.02.Geographic distances were computed as the shortest coastal distances between sites and genetic distances as F ST /(1-F ST ).To compare morphological and molecular differentiation, the Mantel test was also performed in GENETIX 4.02 using morphological distances, defined as Euclidean distances, instead of geographic distances.

Morphological analysis
No statistical differences were found between males and females (t>0.282 or U>20854, P>0.05, for all tests), so sexes were pooled.After size adjustment calculations, all correlations between variables were lower than 0.6, and, as so, none of the variables were discarded prior to the discriminant analysis.The stepwise discriminant analyses performed on morphome- tric and meristic data revealed that the first two discriminant functions explained 60.1% and 80.0% of total variance, respectively (Table 2).For the morphometric data, the first discriminant function was mainly correlated with DN2 (distance between the second pair of nostrils) and HW (head width), whereas the second function was mainly associated with CPH (caudal peduncle height) (negative correlation) and IOD (interorbital distance).For the meristic data, the variables presenting the highest correlations with the discriminant functions were LP (number of rays on the left pectoral fin) and SB (number of sublabial barbells) (function 1) and RP (number of rays on the right pectoral fin) (function 2) (Table 2).
The ordination diagrams obtained from the discriminant analysis of the morphometric data (Fig. 2a) revealed some differentiation between estuarine (dark symbols, left side) and coastal (light symbols, right side) samples along the first discriminant function.These differences were mainly related to the higher values of DN2 and HW presented by the coastal samples as indicated by the high and positive correlation coefficients of these variables with the first discriminant function (Table 2).The one-way ANOVA revealed significant differences between samples (F>29.1,P<0.05, in all tests) and the a posteriori Bonferroni test results showed that estuarine and coastal samples were always significantly different (P<0.05).For the meristic analysis it was not possible to identify a distinctive pattern of differentiation between estuarine and coastal samples in the discriminant analysis diagrams (Fig. 2b).However, estuarine samples showed a more homogeneous pattern, being distributed on the right side of the diagram, whereas coastal samples were distributed throughout the diagram.This pattern was probably related to the higher counts of LP, SB and RP in estuarine individuals, as indicated by the high and positive correlation coefficients of these variables with the first discriminant function (Table 2).Oneway ANOVAs and Bonferroni post-hoc tests showed that these variables had statistically different values between the 10 samples (F>1.3,P<0.05, in all tests) and between coastal and estuarine samples (P<0.05),respectively.

Genetic analysis
Of the 11 loci assayed (Table 4), five (AAT*, GPI-1*, GPI-2*, ME-1* and PGM*) were monomorphic in all samples.A further two loci (IDH* and MDH*) were weakly polymorphic, with common allele frequencies higher than or equal to 0.95 in all the samples, whereas the remaining loci showed common allele frequencies lower than 0.95 in at least one sample (Table 5).Differences in the pattern of enzymatic systems between samples were also found, the most important being those found for IDH*, MDH* and MPI*.IDH* was polymorphic only for Tavira, Quarteira and Sado estuary samples, whereas MDH* was polymorphic only in the Guadiana estuary, Quarteira and the Tejo estuary; MPI* was polymorphic only in the Monte Gordo and Sado estuary samples but presented different alleles in each (MPI 90 and MPI 110 , respectively).The mean number of alleles per locus (MNA) varied between 1.00 and 1.55 and the mean expected heterozygosities (H E ) ranged from a minimum of <0.001 for several loci in different samples to a maximum of 0.62 for G6PDH* in the Sado estuary (Table 6).
The extremely high and positive values of F IS found for most loci in all samples indicated heterozygote deficiency.According to the permutations test only IDH* showed non-significant F IS values in all populations (   et al., 1985); B -tris-citrate pH 8.0 (Pasteur et al., 1985) 3 Tissue: L, liver; M, muscle Overall exact tests of HWE using the Markov-chain method, under the null hypothesis of heterozygote deficiencies, were significant for all populations and loci (P<0.05).Using the same procedure, but only testing for HWE on populations, Olhão and Sagres were the only samples showing non-significant departures from equilibrium.F ST values were generally low (<0.100) but significant, the highest value registered being between Ria Formosa and its coastal area, Olhão (1.000) (Table 7, upper diagonal).Geographically close estuaries showed lower and significant F ST values when compared to those obtained for distant ones: Tejo and Sado (0.001) vs. Tejo and Guadiana (0.007); Ria Formosa and Guadiana (0.381) vs. Ria Formosa and Tejo (0.455).Nei's genetic distances were also generally low (most values <0.010), the highest being registered between Ria Formosa and Monte Gordo (0.033) and Ria Formosa and the Sado estuary (0.031) (Table 7, lower diagonal).Although no consistent pattern was found, Nei's genetic distances were also generally lower between geographically close estuaries than between distant ones.
The Mantel test performed under 10000 permutations showed that genetic and geographic distances were independent, thus refuting the hypothesis of isolation by distance (Z=491.4,P>0.05).The Mantel test performed using genetic and morphological distances also showed that these were not correlated (Z=9030.0,P>0.05).
According to the morphometric analysis, estuarine and coastal samples of H. didactylus showed some degree of differentiation along the Portuguese coast, the distance between the second pair of nostrils, head width, caudal peduncle height and interor-SCI.MAR., 70( 4  bital distance being the most important characters in this distinction.Although showing a reduced discriminating power when compared to the morphometric traits, meristic counts of the number of rays on the pectoral and caudal fins and the number of sublabial barbells were also important factors in sample differentiation. Allozymes revealed that the genetic variability of H. didactylus along the Portuguese coast was similar to that determined by Smith and Fujio (1982) over a large number of marine teleosts (mean heterozygosities of 0.006).The mean F ST value found in the analysis of all loci (0.042) was lower than that reported by Ward et al. (1994) for fishes (0.14) but was in accordance with the low spatial genetic heterogeneity suggested by Gyllensten (1985) for marine teleosts.The Mantel test showed that genetic and geographic distances were not correlated, thus refuting the hypothesis of isolation by distance, and that morphological distances were independent of genetic ones.This pattern of low genetic divergence has been reported for other marine fishes (e.g.Gyllensten, 1985;Kinsey et al., 1994;Chikhi et al., 1998;Cabral et al., 2003a, Pinheiro et al., 2005).
A fine-scale genetic structure was detected between estuarine and coastal areas, as described for other species inhabiting both environments.Using allozyme markers, Ayvazian et al. (1994) detected genetic subdivision in the cobbler, Cnidoglanis macrocephalus (Valenciennes, 1840), a demersal sedentary fish inhabiting estuaries and nearshore marine areas in western and southern Australia.This species' reproductive biology is very similar to that observed in the Lusitanian toadfish, with eggs deposited in nests and males exhibiting parental care, therefore limiting larval dispersion.The authors detected a high genetic divergence between estuaries on a local and regional scale, suggesting low levels of gene flow between estuarine and marine populations.Chaplin et al. (1998) also reported genetically distinct assemblages of the black bream, Acanthopagrus butcheri Gomon et al., 1994, in nine estuaries and one land-locked lake in western Australia, despite the low levels of allozymic variation.However, whereas in these species the greatest divergence was found between estuaries, the estuarine samples of H. didactylus showed low levels of genetic divergence.The high F ST values between the Ria Formosa estuary and its adjacent coastal areas (Olhão and Quarteira) revealed an estuarine vs. coastal sample differentia-tion.A similar pattern was reported for P. minutus along its distribution range in Europe, with a clear distinction of the Venetian lagoons' samples towards the Mediterranean Sea and Atlantic Ocean (Stefanni et al., 1996;2003).The differentiation found in H. didactylus may have arisen from the operation of the barrier islands of Ria Formosa, which physically constrain the flux of individuals.Another plausible explanation is the "divergence-with-gene-flow" (Beheregaray and Sunnucks, 2001) induced by the sedentary behaviour and brood care of toadfish adults, which may contribute to a decrease in local genetic variability and subsequently increase differentiation among populations.
Many of the species found in a given estuary are likely to have resulted from distant immigration and undergone a recent founder effect or bottleneck (Bilton et al., 2002).The genetic difference found between the Tejo and Sado estuaries, which are only a few kilometres apart, might reflect the fluctuation of the toadfish assemblage in the Tejo estuary, characterised by a great abundance period followed by an apparent complete disappearance of this species and a recent re-colonisation (Costa, 1993).Assuming that the heterozygote deficiency and high F IS observed in this location are not the result of misinterpretations, the Whalund effect, inbreeding and directional selection of loci to which these individuals were subjected might explain such values.It would be expected that individuals from the nearest location, the Sado estuary, were the founders of the present day population in the Tejo estuary.The smallest F ST value obtained between the Tejo estuary and the Sado estuary and their small Nei's genetic distance corroborate this hypothesis.However, the low values of these parameters obtained between the Guadiana estuary and the Tejo estuary, the furthest location, are not conclusive, so the process of such colonisation remains intriguing.
Studies performed with marine species whose populations inhabit separate coastal areas suggest an environmental basis for the high morphological differentiation.Tudela (1999) reported a high level of morphological heterogeneity in the European anchovy (Engraulis encrasicolus Linnaeus, 1758) from three sampling sites in the Mediterranean but, according to the electrophoretic analysis, they were genetically homogeneous.Similar results were found by Cabral et al. (2003b) for the Portuguese sole, Synaptura lusitanica Capello, 1868 inhabiting two locations about 250 km apart on the Portuguese coast.Although phenotypic traits are highly variable and often have low heritabilities, reflecting substantial environmental influence (Lindsey, 1988), the overall low genetic differentiation found for H. didactylus should be investigated using more sensitive markers, such as mitochondrial DNA genes or microsatellites.Pampoulie et al. (2004) studied the genetic structure of P. minutus from the North Sea and, although no consistent differentiation was found using allozymes, microsatellite markers revealed a clear differentiation between estuarine and coastal areas samples.Due to their typical high allelic diversity, microsatellites increase the probability of detecting genetic differences among populations (Goudet et al., 1996), and may therefore, detect a genetic structure concordant with the high morphological differentiation of the Lusitanian toadfish on the Portuguese coast.

TABLE 1 .
-Morphometric and meristic variables considered in the morphological analysis of H. didactylus.

TABLE 2 .
-Correlations between the morphometric and meristic variables included in the models and the discriminant functions (see Table1for variables' acronyms).

TABLE 3
. -Percentage of individuals of H. didactylus reallocated in each group in the validation of the discriminant analysis for the morphometric and meristic data.Rows are the original sample group and columns the reallocation group.Meristic data are in italics.GE, Guadiana estuary; MG, Monte Gordo; Tv, Tavira; Ol, Olhão; RF, Ria Formosa estuary; Qt, Quarteira; Sg, Sagres; ME, Mira estuary; SE, Sado estuary; TE, Tejo estuary.

TABLE 4 .
-Enzyme systems, loci studied, separation technique, buffer system and tissue used in the genetic analysis of H. didactylus samples.