Population differentiation of the shore crab Carcinus maenas ( Brachyura : Portunidae ) on the southwest English coast based on genetic and morphometric analyses

1 Centro de Oceanografia, Laboratório Marítimo da Guia, Faculdade de Ciências da Universidade de Lisboa, Avenida Nossa Senhora do Cabo 939, 2750-374 Cascais, Portugal. E-mail: micsilva@fc.ul.pt 2 The Marine Biological Association of the UK, The Laboratory, Citadel Hill, Plymouth PL1 2PB, United Kingdom. 3 Museu Nacional de História Natural, Universidade de Lisboa, Rua da Escola Politécnica 58, 1269-102 Lisboa, Portugal. 4 Centro de Biologia Ambiental, Faculdade de Ciências da Universidade de Lisboa, Campo Grande, 1749-016 Lisboa, Portugal. 5 School of Ocean Sciences, University of Wales Bangor, Menai Bridge, Anglesey LL59 5AB, United Kingdom.


INTRODUCTION
Population genetic analyses can reveal important aspects of evolution and ecology of a species, such as natural selection, gene flow, genetic variation and differentiation.As an example, the amount of genetic variability can determine the adaptation to environmental changes, while gene flow allows connectivity between populations, preventing genetic erosion and inbreeding (Rousset, 2001).Furthermore, from the distribution of molecular genetic markers, important ecological traits, such as mating system and spatial population boundaries, can be inferred (Rousset, 2001).However, understanding these processes in marine species is particularly difficult because barriers to gene flow are far less obvious compared to continental species (Patarnello et al., 2007).Population genetic studies of marine invertebrate species have shown that highdispersal potential due to planktonic larvae is often associated with mild genetic differentiation over large scales, which implies high levels of gene flow within populations (e.g.Reuschel and Schubart, 2006).
A number of recent studies have detected extensive phenotypic variability in shore crabs within relatively restricted geographical areas (e.g.Brian, 2005;Todd et al., 2006).Given that population divergence is prevented by gene flow, patterns of phenotypic variability are likely to reflect differences between the local environmental conditions, as was suggested by Brian et al. (2006).Phenotypic plasticity (i.e.induced changes resulting in different phenotypes in different environments) is today recognized as an important evolutionary mechanism that modulates inherited differences of individuals (Hollander et al., 2006) and can also be an important adaptive strategy in variable or changing environments (Schlichting and Pigliucci, 1998).In fact, if organisms encounter predictable environments, fixed development is expected, whereas in organisms that cannot predict their future environment, phenotypic plasticity would be optimal to increase local adaptation (Hollander et al., 2006).As an example, inducible defences are a ubiquitous form of plasticity that involves the production of chemicals, morphologies, or behaviours by prey species in response to predator cues (Trussell and Smith, 2000).These changes reduce prey vulnerability and examples include the diel vertical migration in marine zooplankton (Bollens et al., 1992) and shell shape and thickness in gastropods and bivalves (Trussell, 1996).Despite improved understanding of the cues inducing these changes and their immediate adaptive value (Schlichting and Pigliucci, 1998), our understanding of how this phenomenon contributes to broader temporal and spatial patterns of phenotypic variation remains poor.
The shore crab Carcinus maenas is common in the intertidal throughout Europe and has been well studied due to the ease with which it can be found, identified, sexed and measured.This highly adaptable crab has recently gained notoriety due to its globally invasive nature associated with drastic ecological and economic effects.Once established, it becomes the dominant intertidal crab in some areas (Yamada, 2001), affecting the abundance, size structure and defence response of native species (Roman and Palumbi, 2004).Thus, investigations of the population structure of this crab within its native range can provide insights into processes of invasions in other regions.
Carcinus maenas has a long planktonic larval phase (up to 50 days in the plankton; Tresher et al., 2003), and its offshore dispersal may account for considerable exchange of individuals between local populations (Peliz et al., 2007).According to Queiroga (1996), the zoeal stages I and II are concentrated in the surface layer, but a gradual ontogenic displacement to deeper waters is observed from then on.Horizontally, there is a clear association of the first zoea with the sites where hatching occurs, while the older zoeal stages are dispersed progressively offshore (Queiroga, 1996).Therefore, in accordance with the population genetics theory that suggests that marine animal species with long planktonic larvae have less genetic structure and higher connectivity than those with direct development (Palumbi, 2003), C. maenas would be expected to have low levels of population differentiation.However, previous studies of the shore crab C. maenas have suggested several different conclusions.Bulnheim and Bahns (1996) reported a slight geographic cline from north to south Europe in one allozyme locus.Bagley and Geller (1999) found no population structure in a microsatellite DNA study of Atlantic European crabs, and a recent study by Pascoal et al. (2009), also with microsatellite data, revealed very weak but significant structuring between populations on the Portuguese coast.A more extensive study by Roman and Palumbi (2004), with the mitochondrial cytochrome oxidase I gene, detected a slight population structure between the central North Sea and populations to the south, and a break between populations from the Faeroe Islands and Iceland and continental populations.
The shore crab Carcinus maenas (Linnaeus, 1758) is a typical inhabitant of the European coastline, which lives in the tidal and in the subtidal zone (Bulnheim and Bahns, 1996).In its native range, C. maenas has a distribution from northern Norway and Iceland to Mauritania, and it has successfully colonized Australia, Tasmania, South Africa, Japan and both coasts of North America (Roman and Palumbi, 2004 and references therein).On the English coast, C. maenas occurs throughout the coastline, the estuarine populations being linked by open coast populations.It has a complex life cycle comprising an exported planktonic larval phase that develops in shelf waters and takes four to six weeks to reach the megalopa stage (Tresher et al., 2003).Megalopae migrate back to the coast during spring tides, and transport is accomplished by selective tidal stream transport (Moksnes et al., 1998).Once they reach a suitable environment, the megalopae of C. maenas settle in a variety of intertidal and subtidal habitats, though they show a preference for those that are structurally complex (Paula et al., 2006).
Our study explored local patterns of population structure of onshore populations of C. maenas on the coastline of southwest England, using eight microsatellite loci as genetic markers.Morphological differences among these populations were assessed to show whether genetic and morphometric tools provide coinciding or contrasting results in population differentiation.

Sampling sites and procedures
Adult shore crabs were collected using baited traps from the intertidal zone at nine estuarine sites on the southwest English coast during the Autumn of 2004 and 2005 (Fig. 1).The sites were the Axe, Exe, Teign, Tamar, Looe, Fowey, Hayle, Camel and Torridge estuaries.After collection, specimens were transported to the laboratory, where they were frozen at -20ºC for preservation, and the fourth and fifth pereiopods were preserved in absolute ethanol.Specimens were separated by size and sex and labelled, and the carapaces and chelae were processed to eliminate those that were not suitable for the study (broken carapaces and regenerated and/or broken claws).Females of the species C. maenas with a carapace width of less than 1.5 cm and males with one of less than 2.5 cm were also removed from the morphometric analyses.These sizes correspond to the minimum sizes of sexual mature specimens (Neal and Pizzolla, 2008).Thus, this procedure eliminated most allometric growth variation.Morphometric analyses of the carapaces were made following Silva and Paula (2008), using the major chelae in males and right chelae in females.

Genetic data analysis
A total of 270 individuals were collected for genetic analyses.Forty-five specimens from each of the Exe, Teign, Tamar, Looe, Camel and Torridge populations were analysed.The Axe, Fowey and Hayle populations were not considered for the genetic analyses due to logistical constraints.
Total genomic DNA was extracted from muscle tissue of pereiopods using sequential phenol-chloroform extraction steps as described by Hillis et al. (1996).The DNA obtained was resuspended in low TE buffer (Tris-HCl 10 mM pH 8.0, EDTA 0.1 mM pH 8.0) and used as a template in polymerase chain reactions (PCR).Before PCR amplification, DNA concentration was estimated using the Thermo Scientific NanoDrop TM 1000 Spectrophotometer.Genetic diversity was screened at 8 microsatellite loci (Cma10EPA, Cma11EPA, Cma12EPA, Cma14EPA, developed by Tepolt et al. (2006); and SP107, SP229, SP280 and SP495 developed by Pascoal et al., 2009).Microsatellite polymorphism was detected using a fluorescent detection method, the forward primer for each locus being 5'-labelled with one of three fluorophores.The labelled PCR products from the 8 loci were then divided into three sets, taking into account the size of the fragments (Cma10EPA + Cma11EPA + SP107; Cma12EPA + Cma14EPA + SP229; SP280 + SP495).PCR amplifications were carried out in 15-μl reactions with 10 pmol of each primer set, 10-100 ng of genomic DNA, 1x PCR buffer, 1.3-1.8mM MgCl 2 , 0.2 mM of each dNTP (Promega) and 0.5 U Taq DNA polymerase (Promega).Thermal cycling parameters were: initial denaturation at 94ºC for 5 min followed by 26-35 cycles of denaturation at 94ºC for 1 min, primer-specific annealing temperature (Table 1) for 1 min, extension at 72ºC for 1 min and final extension at 72ºC for 10 min.Following PCR amplification, the extension products were resolved on 1.5% agarose gels.PCR products were run on a Beckman Coulter TM sequencer.Quality of PCR products and allelic length were determined using the CEQ TM Genetic Analysis System software.Alleles were designated according to their sizes and the Microsatellites toolkit, a specific tool for Microsoft Excel, was used as a conversion program in order to obtain the genepop, fstat and arlequin input data files.
Raw data were analysed with micro-checker software (van Oosterhout et al., 2004) to check microsatellite data for null alleles and scoring errors.The observed and expected heterozygosity for each locus and population were calculated using genepop 1.2 (Raymond and Rousset, 1995).The number of alleles and the allelic richness for each locus and population were estimated using fstat 2.9.3.2 (Goudet, 2002).To test for deviations from Hardy-Weinberg expectations and linkage disequilibrium a probability test was performed, with unbiased p-values and associated standard errors (SE) estimated by the Markov Chain algorithm.Bonferroni correction of the p-values was used to correct multiple tests.As suggested by Raymond and Rousset (1995), number of batches (B) and number of iterations per batch (C) were increased to SE ≥0.01 (B =400, C =4000).These estimations were achieved with genepop 1.2 (Raymond and Rousset, 1995).The Weir and Cockerman (1984) analogues of Wright's F-statistics, F ST , and their significance per locus were estimated to assess levels of population differentiation over all loci, using fstat 2.9.3.2 (Goudet, 2002).Isolation-by-distance was examined employing Mantel's t-test (Mantel, 1967), computed with mantel 1.18 (Cavalcanti, 2005).Pairwise comparisons were made of geographic distances and F ST .A Bayesian approach was used to assess population structure by model-based clustering methods with the program baps 3.2 (Corander et al., 2005), which identifies the optimum number of partitions among groups of samples.Furthermore, genotypes were plotted in a multidimensional space, using a principal-component analysis (PCA) with the program pca-gen 1.2 (Goudet, 1999) and labelled as 1 (Exe), 2 (Teign), 3 (Tamar), 4 (Looe), 5 (Camel) and 6 (Torridge).

Morphometric data analysis
Fifty males and fifty females from each population were used in the morphometric analysis.To quantify possible shape variation in the carapaces and in the chelae, landmark-based geometric morphometrics was used.Images of each specimen were taken using a Canon Power Shot A85 digital camera with a resolution of 4.0 megapixels and consistent zoom and distance, in order to avoid false perspectives and allow accurate comparisons.Fifteen carapace landmarks and nine claw landmarks were digitized using the program tpsDig 2.10 (Rohlf, 2006) (Fig. 2).Only one side of the carapace (the left) was used to avoid duplication of equivalent landmarks and computation problems (Rufino et al., 2004).
The digitized landmark configurations were then subjected to a generalized Procrustes analysis (GPA; Rohlf and Slice, 1990) to remove the effects of size, position and orientation in the digital images.The aligned landmark configurations produced by this analysis were used to generate the shape variables, the relative warp scores and the uniform component, from a thin-plate spline analysis (TPS).These variables were obtained with tpsRelw 1.45 (Rohlf, 2007a).The software tpsRegr 1.34 (Rohlf, 2007b) was used to detect the presence of allometry, and when detected, it was removed by the regression of each relative warp against a measure of body size (centroid size), thus estimating residual shape variation.
To understand patterns of morphometric differentiation among populations, a multivariate analysis of variance (MANOVA) was performed on the relative warps scores.After assessment of the degree of variation among populations, pairwise comparisons among these populations were done with the post-hoc Tukey HSD test.In addition, discriminant function analysis was used to determine Mahalanobis distances between each pairs of sites, for a more accurate differentiation between populations.Finally, a cluster analysis was generated from the Squared Mahalanobis distance matrix to assess phenotypic relationships among the populations.These procedures were performed in statistica 6.0.

Genetic data
A total of 270 Carcinus maenas collected from southwest England were typed for eight microsatel- lite loci.All loci were highly polymorphic (number of alleles: 30 for Cma10EPA, 47 for Cma11EPA, 36 for Cma12EPA, 12 for Cma14EPA, 11 for SP107, 24 for SP280 and 10 for SP495), with the exception of locus SP229, for which only one allele was detected.Therefore, locus SP229 was not used in further analyses.All the populations were characterized by a considerable mean number of alleles per locus, ranging from 13.71 to 16.29, and by a number of private alleles ranging from 3 to 8. The expected heterozygosity (H Exp ) per locus ranged from 0.356 to 0.957 and the observed heterozygosity (H Obs ) from 0.347 to 0.795 (Table 1).Three loci showed a significant heterozygosity defect (Table 1, p-value <0.05).Considering the Hardy-Weinberg Equilibrium (HWE) results by locus and by population, loci Cma10EPA, Cma11EPA and SP280 deviated from HWE in all populations.In the other four loci, no evidence of heterozygote deficiencies was found, so it was possible to assume that all populations were in HWE.The observation of many cases of excess homozygosity at a single locus may suggest the presence of null alleles.This was confirmed by the results of micro-checker 2.2.3.After correction for null alleles with the Brookfield null allele estimator 1 (Brookfield, 1996), subsequent analyses were performed on the corrected data set.Genotyping errors due to stuttering and large allele dropout were not found (confidence interval = 95% and 10000 iterations).Fisher exact test for linkage disequilibrium among loci across populations revealed no significant pvalues among 21 possible combinations.
Table 2 shows F ST estimates across all loci for all population pairs.F ST values ranged from 0.0015 to 0.0458, illustrating low levels of genetic structuring.Following Bonferroni correction, only three of the 15 pairwise F ST values were found to be significant (between Exe-Tamar; Tamar-Camel; Tamar-Torridge).No significant isolation-by-distance was shown when pairwise F ST values and geographic distance were compared (Mantel test, 10000 randomizations, r =-0.224, p =0.205).The hierarchical analysis of molecular variance (amova) showed that genetic variation within and among populations amounted to 97.98% and 1.4% of total variation, respectively.The genetic variation among pre-defined groups was close to zero (0.62%).Overall F ST and the F SC values were low (0.0202 and 0.01404 respectively), and not significant (p >0.05).This indicates no genetic differentiation among all populations sampled and among populations within groups, respectively.The F CT value found was also very small (0.00622) and non-significant (p >0.05), indicating low genetic differentiation among the proposed groups.The baps analyses, conducted assuming different numbers of groups of population, indicated that the best partition included all six sampled populations (ln likelihood =-8753.2).The result of PCA on genotypes is  shown in Figure 3. Scores were plotted onto the two principal axes (PC-I and PC-II), which cumulatively explain 74.7% of the total genetic diversity.The pvalue for the overall F ST was 0.01 and the global F ST was 0.02616.Per axis, the F ST values achieved were, respectively, 0.01287, 0.00667 0.00270, 0.00208 and 0.00184.This analysis showed a separation of all populations, so no clear clustering was formed.However, the Exe (1) and Torridge ( 6) populations seemed to be closer together than to the remaining populations, and the Tamar (3) population was further apart from all the other populations.

Morphometric data
Analyses performed in the male carapaces revealed that the first relative warp (RW1) explained 47.6%, the second (RW2) 12.5% and the third (RW3) 8.1%, summing 68.2% of variance explained.The first axis revealed shape variation in the rostrum and on the overall carapace shape, with the negative end of this axis corresponding to more round shapes.RW2 explained shape variation in the lateral-posterior spines and the distances among them, and RW3 accounted for differences due to the variation in the carapace width (Fig. 4).A similar variation pattern was observed in females, but the rounder carapaces corresponded to the positive values of relative warp 1.In this case, the first three axes explaining 43.3%, 11.5% and 9%, respectively (Fig. 4).Shape variation in the chelae was more marked.For the males, RW1 explained 53.7% of the variation explained by all the relative warps, RW2 explained 12.8% and RW3 explained 9.4%.The first warp accounted for the variation due to the size of the pollex, the height of the manus and the flattening of the posterior part of the claw (Fig. 5).The variation explained by RW2 corresponded to the orientation of the pollex and the total length of the chela.Relative warp 3 accounted for the shape variation due to the size of the gape.In accordance with the shape variation of the carapaces, a similar shape variation pattern was found in the females, with the first three axes explaining 53.3%, 11.3% and 8.9%, respectively (Fig. 5).
The multivariate analysis of variance revealed significant morphometric differences among the nine populations of Carcinus maenas (males: Wilks' λ =0.189, p <0.0001; females: Wilks' λ =0.226, p <0.0001).The comparisons made among populations, obtained with the post-hoc Tukey HSD test, showed that the Torridge population was significantly different from almost all the other populations.Additionally, the post-hoc test revealed the absence of significant differences between the sites that were closer together.In contrast, the sites further apart tended to show significant differences between each other (Table 3).The discriminant function analysis and cluster analysis revealed that there was extensive inter-population variability in the morphology of both male and female C. maenas, but differences were more pronounced in males, as was detected by the squared Mahalanobis distance.Nevertheless, the patterns exhibited by each sex were similar.Distances between the morphology of male crabs showed two well-defined clusters (Exe, Axe and Teign; Hayle, Looe and Tamar).Males from Fowey, Camel and Torridge exhibited the greatest differences from these groups (Fig. 6).The cluster analysis performed in the females revealed three clusters on the north coast (Hayle, Camel and Torridge) and two overlapping clusters on the south coast (Exe, Teign, Looe; and Axe, Fowey and Tamar).However, Hayle, Camel and Torridge females were most distinct in terms of their morphology.

DISCUSSION
Carcinus maenas exhibited no clear evidence of genetic differentiation among the sites sampled on the southwest English coast.Different analyses, including AMOVA and the Mantel test, revealed no pattern of geographic separation among populations.Therefore, genetic variation was not associated with geographic subdivision.Furthermore, the Bayesian approach used to assess genetic population structure suggested that all sampled individuals comprised a single genetic population.Pairwise F ST low values suggested no restriction of movement among sites and, consequently, no impediment to gene flow between populations.
Migrations during larval development to the marine environment promote transport over large distances, allowing gene flow between populations (Bilton et al., 2002).Carcinus maenas has a long planktonic larval phase (up to 50 days) that can be dispersed extensively in the water column before settling to the substratum and undergoing metamorphosis.Peliz et al. (2007) conducted a study with C. maenas along the west coast of the Iberian Peninsula, and suggested a larval dispersal distance of about 60 km (s ~ 20 km), which is mostly representative of alongshore larval dispersion.Because sampled populations along the studied coast are separated by distances of 15 to 432 km, such a dispersal radius may account for considerable exchange of individuals between local populations.Thus, the long planktonic development mode of C. maenas can allow larval flow between distant populations and, consequently, might promote population homogenization.These results further suggest that there are high degrees of connectivity within the study area with little evidence of reduction in gene flow, either resulting from barriers that block the dispersion potential or failures in the reproduction of migrant individuals within the new population to which they are dispersed.
A recent study by Pascoal et al. (2009) with microsatellite data detected very weak but significant structuring between C. maenas populations on the Portuguese coast, which, according to the authors, might result from hydrodynamic or topographic barriers.On the Portuguese coast, C. maenas forms large populations in estuaries and on their adjacent rocky shores.The geographic distance that separates estuaries may be characterized as a barrier, thus limiting larval flux between populations.On the other hand, C. maenas has a continuous distribution on the English coast, occurring on all types of shores (Neal and Pizzola, 2008), and no evidence of genetic differentiation among the sites sampled was found in our study.These differences detected in the population structure of C. maenas could therefore be due to the contrasting distribution patterns exhibited on both coasts.
The population from the River Tamar Estuary revealed statistically significant levels of population differentiation when compared with the Exe, Camel and Torridge populations.The River Tamar is a large river, and its drainage to the Atlantic Ocean may affect nearshore currents and salinity, and influence larval transport.Newly hatched larvae migrate vertically in synchrony with the tidal cycle, attaining their highest position in the water column during ebb (Queiroga, 1996).This migration, under behavioural control (Zeng and Naylor, 1996), ensures that larvae are exported to the sea shortly after spawning.The influence of estuarine currents and brackish water can also limit the pattern of offshore larval dispersal and onshore megalopal transport, favouring genetic differentiation from the rest of the populations.A similar pattern of population differentiation was found in the Portuguese Tagus Estuary and on its adjacent rocky shores for the shore  crab species Pachygrapsus marmoratus (Silva et al., 2009).The Tagus River drains into the Atlantic, affecting nearshore currents and salinity, and may prevent larval transport across the river mouth.The influence of estuarine currents and brackish water can also limit the pattern of offshore larval dispersal and onshore megalopal transport, favouring genetic differentiation from the rest of the populations.Another hypothesis that could explain this genetic differentiation involves natural selection.Population differentiation resulting from differential mortality of larval recruits was documented in the mussel Mytilus edulis (see Koehn et al., 1980).Kordos and Burton (1993) also observed in the crab Callinectes sapidus that selection against megalopae of specific allozyme genotypes varied from bay to bay, resulting in both differences between larvae and adults in a given bay and differentiation of adult populations among bays.Thus, due to the influence of estuarine brackish water, post settlement selection might have acted on the megalopae, restricting suitable genotypes to this particular environment.
In contrast with the low levels of genetic variability found between populations, extensive morphometric variability was detected for both males and females.This morphometric variability and population differentiation is in accordance with studies already carried out in the UK with this species (e.g.Brian, 2005;Brian et al., 2006).Brian et al. (2006) also showed that the relationship between the phenotypic characteristics and patterns of genetic similarity was not strong.The results achieved in the current study, in accordance with Brian et al. (2006), suggest that genetic and phenotypic characters are in some way independent, and that environmental factors, such as habitat, wave action, predation, food supply and salinity can play an important role in the patterns of shape variability.However, the extent of plasticity in shore crabs is unclear.Two aspects of the life history of the shore crab may have led to the evolution of plastic responses (Todd et al., 2006): due to its pelagic larval stage, high levels of genetic homogeneity are expected, so there is limited opportunity for local adaptation (Roman and Palumbi, 2004); and post-settlement movements can occur, particularly in response to predator pressures, even though settling megalopae have swimming and habitat-selection ability (Moksnes et al., 1998).Furthermore, the morphology of C. maenas may respond plastically to site-specific selection pressures, causing changes in the expression of a particular morphological trait (Brian et al., 2006).
Carcinus maenas is a heterochelous species, and Silva and Paula (2008) have shown that the crusher claw is more suitable for population differentiation studies.In the current study, clearer patterns of shape variation were seen in the crusher claw than in the carapace.Previous studies have shown that claw size is capable of responding to environmental pressures during development and that differential use influences claw size (Smith and Palmer, 1994).Moreover, patterns in claw size, form and performance observed by Smith (2004) and Baldridge and Smith (2008) suggested diet and water temperature responses to geographic differences in prey armour.These studies, as well as the results achieved in this study, have thus demonstrated that morphometric characters of C. maenas can respond in a plastic way to environmental pressures.
Given their generalist genotype and the environmental heterogeneity associated with their native habitat, the ability to respond to challenges in this plastic way is likely to be a highly advantageous strategy for C. maenas (Brian et al., 2006).These factors, combined with the capacity to survive under extreme conditions (Roman and Palumbi, 2004), high fecundity and a long planktonic larval phase, can explain its success as an invasive species.

Fig. 4 .
Fig. 4. -Scatterplot of individual scores from the relative warp analysis.Grids show deformation from consensus configuration of the carapace of Carcinus maenas males (A) and females (B).

Fig. 5 .
Fig. 5. -Scatterplot of individual scores from the relative warp analysis.Grids show deformation from consensus configuration of the major chela of Carcinus maenas males (A) and females (B).

Fig. 6 .
Fig. 6. -UPGMA dendrogram derived from cluster analysis of squared Mahalanobis distances for relative warp scores, for the major chela of Carcinus maenas males.

Table 2 .
-Genetic differentiation among Carcinus maenas populations.F ST values are below diagonal, and geographic distances (km) among populations are above diagonal.* p<0.05: significance of F ST after sequential Bonferroni correction.

Table 3 .
-Pairwise comparisons among populations achieved with the post-hoc Tukey HSD test.Female values are above diagonal and male values are below diagonal.** p<0.001; n.s., non significant value.