Inter- and intraspecific genetic and morphological variation in a sibling pair of carabid species

Background Pogonus littoralis and Pogonus chalceus are very close related species with quite different ecological preferences within salt marshes. We study the evolutionary processes in and between these presumably young species. Therefore, we compare the variation in ecologically relevant characters and the genetic variation within one of the species (intraspecific differentiation) with the variation of the two types of characters between the two species (interspecific variation). Data are compared between two independent sets of populations, one set at a small geographical scale (the ecologically diverse Guérande area in France) and the other set at a Atlantic-Mediterranean scale. Results Body and relative wing size and IDH1 allozyme data show that the intraspecific variation in P. chalceus is high and in the same range as the interspecific variation (P. chalceus versus P. littoralis). Based on neutral markers (other allozymes and mitochondrial DNA) on the other hand, the intraspecific variation in P. chalceus is much lower in comparison to the interspecific variation. Conclusion The different ecotypes in the highly polytypic species P. chalceus are as highly differentiated in ecological characters as true species, but are not recognised as such by screening neutral DNA polymorphisms. This can be interpreted as a case of ongoing speciation driven by natural selection adapting each ecotype to its respective ecological niche. The same ecological process can be recognised in the differentiation between the two sister species, where en plus reproductive isolation between the two gene pools occurred, allowing independent drift and mutation accumulation in neutral genetic characters.


Background
Pogonus chalceus is a wing polymorphic beetle with extremely variable wing size from short to completely developed wings, with all possible intermediates [1]. A recent study presented population genetic results on P. chalceus (Marsham, 1802) beetles from the Guérande saltfields on the French Atlantic coast, based on allozymes and microsatellites, as well as results on wing and body size [2]. In the unique man made Guérande salt-fields, two contrasting habitat types are found mixed on a microscale in hundreds of replicates (sea canal versus salt extraction ponds). Body, relative wing size and IDH1 allozyme alleles are strongly divergent between these two contrasting microhabitats; divergent selection led to two clearly distinguishable ecotypes, respectively adapted to canal and pond habitat. Comparisons between the Guérande region (microscale) and populations along the Atlantic coast (macroscale) confirmed the generality of the hypothesis regarding ecological processes responsible for this differentiation: habitat stability [2]. The Guérande ecotypes are also slightly differentiated based on neutral molecular markers (microsatellites and allozymes), suggesting that partial barriers to gene flow between the two ecotypes are present. Previous work on a wide range of taxa has demonstrated that strong natural selection can lead to divergence in spite of gene flow [3][4][5][6][7]. Our Guérande results can therefore be interpreted as a case of ongoing speciation driven by natural selection adapting each ecotype to its respective ecological niche, i.e. species in status nascenti (see also [8,9]).
In the same Guérande region and along the European Atlantic and Mediterranean coast, another Pogonus species, P. littoralis (Duftschmid, 1812) lives in a third kind of microhabitat: unvegetated, temporary dry salt marsh ponds or creeks, where it lives between cracks in humid sea clay. This species is, in contrast to P. chalceus, constantly macropterous, always with maximally developed wings and functional flight musculature [10]. The beetle is highly mobile because it regularly has to move between temporarily dry salt marsh ponds and creeks during its life cycle. Both species can be hardly distinguished by external morphology (for example large individuals of P. chalceus versus small P. littoralis) but have clearly distinguishable genitalia.
The data in this article are to some extent compiled from previous works [2,11,12]. Nevertheless, the novelty of this study lies in the fact for the first time the two carabid sister species are analyzed jointly allowing for valuable comparisons to be made. In this study, we will first compare the two ecotypes of Pogonus chalceus with the closely related species, Pogonus littoralis at a microscale (Guérande region). Therefore, we will use population data on wing and body size, IDH1 allozyme polymorphism as well as apparently neutral markers (other allozymes and mtDNA). We will also test if the microscale results are valid at a larger scale across Europe by means of an independent data set of different populations of both Pogonus species. In all of these cases, we will evaluate the contribution of intra population, inter population, inter ecotype and interspecific variation to the total variance. Fig. 1  In the Guérande region and considering the two species (nested design ANOVA; six P. chalceus populations (canals and ponds pooled) versus three P. littoralis populations, the major part of variance (based on body size) is found among species (Table 1; 74.24% for males and 51.96% for females). If we consider three groups (three canal populations (P. chalceus), three pond populations (P. chalceus) and three P. littoralis populations, the major part of variance is even more pronouncedly found among groups (84.96% for males and 72.08% for females). Variance among populations within groups considering three groups instead of two drops from 17.5 to 2.35% for males and from 29.57 to 4.28% for females. This indicates that this variance was almost completely due to the differences in body size between populations of P. chalceus from different microhabitats. All variance components are statistically significant.

Body size
On a European scale and considering the two species (25 P. chalceus populations versus six P. littoralis), the major part of variance (based on body size) is found among species (  In the Guérande region and considering the two species (nested design ANOVA; six P. chalceus populations (canals and ponds pooled) and three P. littoralis populations), the major part of variance (based on relative wing size) is found among species (Table 2; 78.84% for males and 80.37% for females). If we consider three groups (three canal populations (P. chalceus), three pond populations (P. chalceus) and three P. littoralis populations), the major part of variance is even more clearly found among groups (95.48% for males and 93.92% for females). Variance among populations within groups considering three groups instead of two drops from 18.64 to 0.64% for males and from 16.21 to 0.8% for females. This indicates that this variance is almost completely due to the differences in relative wing size between populations of P. chalceus from different microhabitats. All variance components are statistically significant.
On a European scale and considering the two species (25 P. chalceus populations versus six P. littoralis populations), the major part of variance (based on relative wing size) is found among species (Table 2; 57.61% for males and 60.16% for females). If we consider four ecological groups (14 temporary (P. chalceus), five intermediate (P. chalceus), six stable (P. chalceus) and six P. littoralis populations), the major part of variance is again even more pronounced among groups (82.22% for males and 82.58% for females). Variance among populations within groups considering four groups instead of two drops from 36.11 to 6.34% for males and from 32.8 to 6.89% for females. This indicates that this variance is again almost completely due to the differences in relative wing size between populations of P. chalceus from different ecological or salt marsh area stability groups. All variance components are statistically significant.

IDH1 allozyme marker
In Guérande, both Idh1-2 and Idh1-4 alleles are frequent in ponds, whereas canals are nearly fixed at Idh1-4 (  (Table 3; 61.93%). If we consider three groups (three canal populations (P. chalceus), three pond populations (P. chalceus) and three P. littoralis populations, the major part of variance is still found among groups (64.25%). Variance among populations within groups considering three groups instead of two drops from 11.49 to 0.1%. This indicates that this variance is almost completely due to differences in IDH1 between populations of P. chalceus from different microhabitats. All variance components are statistically significant.
At a European scale, both Idh1-2 and Idh1-4 alleles are frequent in the intermediate stability populations, whereas the temporary populations are nearly fixed at the Idh1-2 allele and the stable populations at the Idh1-4 allele (Fig. 3B) [see also additional file 2]. P. littoralis populations are fixed at the Idh1-6 allele. Considering two groups (25 P. chalceus populations versus six P. littoralis populations, the major part of variance (based on IDH1) is found among groups (Table 3; 62.29%). If we consider four ecological groups (14 temporary (P. chalceus), five intermediately stable (P. chalceus), six highly stable (P. chalceus) and six P. littoralis populations, the major part of variance is somewhat lower but is still found among groups (53.26%). Variance among populations within groups, considering four groups instead of two, drops from 12.26 to 2.07%. This indicates that this variance is almost completely due to differences in IDH1 between populations of P. chalceus from different ecological groups, as P. littoralis is fixed in a different allele. All variance components are statistically significant.

Other allozymes
The number of studied individuals and allozyme allele frequencies for each population in the Guérande region is given in additional file 3. In the Guérande and considering two groups (AMOVA; six P. chalceus populations (canals and ponds) versus three P. littoralis populations), the major part of variance (based on four neutral allozymes) is found between species (

Mitochondrial DNA
The 459-bp COI mitochondrial sequences revealed two haplotypes in the Guérande region. Haplotype one was shared by individuals of both canal and pond ecotype.
Haplotype two was exclusive to the canal ecotype ( Table  4). The 497-bp 16S sequences revealed only one haplotype in the Guérande region ( Table 5).
The 459-bp COI sequences included 32 variable sites on a European scale (29 parsimony informative) and revealed nine unique haplotypes (four for P. chalceus and five for P. littoralis, with no haplotype shared by the two species). Most haplotypes were exclusive to a particular sampling site, with the exception of haplotype one and three which appeared in eight different localities, and haplotype five, which was found in three localities (Table 4). Selective neutrality was confirmed for this gene (P > 0.1 for all test statistics with D* and F*).
The neighbour joining tree in Figure 5 shows that both species form clearly separated entities (high bootstrap values). Differences are found in 28 positions between haplotypes 5,6,8,9 (P. littoralis) and haplotypes 1 and 2 (P. chalceus; Table 4). 27 base differences are found between haplotypes 5,6,8,9 and haplotypes 3 and 4 (P. chalceus). 27 base differences are found between haplotype 7 (P. littoralis) and haplotypes 1,2 (P. chalceus) and 26 base differences with haplotypes 3,4 (P. chalceus). Intrapopulation differences in both P. chalceus and P. littoralis are very small (only a limited number of individuals studied) and between each haplotype there are only one to at most two bases different.
The 497-bp 16S sequences included 11 variable sites (9 parsimony informative) and revealed three unique haplotypes (Table 5; two for P. chalceus and one for P. littoralis,  with no haplotype shared by the two species). Only one haplotype was exclusive to a particular sampling site (haplotype two). Haplotype one appeared in all 11 P. chalceus localities, and haplotype three appeared in all three P. littoralis sites. Selective neutrality was confirmed for this gene (P > 0.1 for all test statistics with D* and F*). Both species form, as in the case for COI, clearly separated entities (Table 5). 10 base differences are found between haplotype 3 (P. littoralis) and haplotypes 1 and 2 (P. chalceus).
Interpopulation differences in both P. chalceus are very small (between the two haplotypes there is only one base different). There were no interpopulation differences found in P. littoralis.

Discussion
Pogonus littoralis and Pogonus chalceus are closely related species, sometimes relatively hard to identify without dissection of the genitalia. We are interested to study the evolutionary processes in and between these presumably young species. We therefore compare the degree of intraspecific variation (in ecological groups of P. chalceus) and the degree of interspecific variation (P. chalceus versus Allele frequencies for IDH1 in the Guérande region (part A; cnals (P. chalceus) and ponds (P. chalceus) and P. littoralis) Figure 3 Allele frequencies for IDH1 in the Guérande region (part A; cnals (P. chalceus) and ponds (P. chalceus) and P. littoralis). Allele frequencies for IDH1 on a European scale (part B; stable (P. chalceus) and intermediate (P. chalceus), temporary and P. littoralis). Idh1-2: black, Idh1-4: light gray; Idh1-6: dark gray.

P. littoralis) between a variety of morphological characteristics and molecular markers.
In all of these cases, we did this with an ANOVA splitting up the total variance among groups, among populations within groups and within populations (Table 6).
At both geographical scales and considering two groups (P. chalceus populations versus P. littoralis populations), a very large part of the total variance (based on body size, relative wing size, IDH1 and four neutral allozymes) is found between species (summary in Table 6). The study of the two mitochondrial genes also shows that both species form clearly separated entities. It is clear that relative wing size differences as well as genetic differences between the sister species P. chalceus and P. littoralis (interspecific) in this study are very marked and allow an easy species recognition.
Body size, relative wing size and IDH1 allozyme data in the beetle P. chalceus are also strongly divergent between contrasting microhabitats (intraspecific: two ecotypes in Guérande) as well as between three ecological groups at macroscale (highly stable versus intermediately stable and temporary populations; [2] and this study). If we consider four groups on a macroscale (3 groups in P. chalceus + 1 group of P. littoralis) or three groups on a microscale (2 ecotypes in P. chalceus + 1 group of P. littoralis), the variance among populations within groups drops drastically as compared to the analysis of two groups (all P. chalceus populations versus P. littoralis; based again on body size, relative wing size and IDH1; summary in Table 6). This study clearly shows that the intraspecific variation based on those three characteristics in P. chalceus is very high and in the same order of magnitude as the degree of interspecific variation (P. chalceus versus P. littoralis). We have suggested earlier that this huge phenotypic and IDH1 divergence in P. chalceus has been driven by divergent natural selection [2]. As relative wing size is to a large extent genetically determined [1], this indeed suggests divergent selection between populations. And the observation that the IDH1 locus screened within our samples shows alelic differences between habitats strongly suggests a locus undergoing evolution through natural selection. Moreover, the canal and pond microhabitats differ from each other with respect to temperature, salinity and water level fluctuations [2]. Numerous studies based on allozymes have revealed patterns of allelic distribution associated with environmental factors, such as temperature and salinity [13,14]. Regarding the function of IDH1, the enzyme catalyses the rate-limiting step of the tricarboxylate cycle. Possible links with growth, however, are not direct and could be associated with the energy that is produced from the reaction. Divergent selection can lead to reproductive isolation and assortative mating and ultimately to speciation [8,15].
On the other hand, in a previously study was shown that P. chalceus ecotypes in the Guérande region were only slightly differentiated (based on allozyme and microsatellite markers) compared to the results based on adaptive characteristics [2]. The smaller degree of intraspecific divergence is also reflected in the mitochondrial data from this study. Moreover, allozyme and mtDNA data from this study show that the populations of P. chalceus are much more related to each other than to their sister species P. littoralis both on a micro-and macroscale. Often, little or no genetic divergence is found in neutral markers between ecologically and morphologically differentiated populations [3][4][5]7,[16][17][18]. Our results can be interpreted as a case of ongoing speciation in P. chalceus where divergence reflects a balance between selection and gene flow (see also [2]). Several studies suggest that tital marshes may be an appropriate ecotone in which to search for instances of ecological speciation. The studied species show, as is the case in our study, distinct morphological differences despite little divergence in molecular markers [7,[19][20][21].
In view of the above shown analogy between intra-and interspecific variation, it seems reasonable to assume that the same ecological adaptive bifurcation was also the first step in the speciation process of P. chalceus and P. littoralis. The speciation process was here fully accomplished by the reproductive isolation between the two groups, allowing independent drift and mutation accumulation in neutral genetic characters.

Methods
Sampling P. chalceus populations from three different sites in the Guérande region are analysed (microscale; Fig. 5; see also [2]). Each site consists of two drastically differing micro-Neighbour joining tree based on COI haplotypes described in Table 2  Figure 4 Neighbour joining tree based on COI haplotypes described in Table 2. Haplotype 1-4 are found in P. chalceus and haplotype 5-9 in P. littoralis. Bootstrap values are indicated (1000 replicates).   T  A  T  T  T  T  C  C  C  A  T  T  A  G  T  G  T  G  T  C  T  T  T  T  A  A  T  T  T   6  A  A  A  T  A  T  T  T  T  C  C  C  A  T  T  A  G  T  G  T  G  T  C  T  C  T  T  A  A  T  T  T   7  A  A  A  T  A  T  T  T  T  C  C  C  A  T  T  A  A  T  G  T  G  T  C  T  T  T  T  A  A  T  T  T   8  G  A  A  T  A  T  T  T  T  C  C  C  A  T  T  A  G  T  G  T  G  T  C  T  T  T  T  A  A  T  T  T   9  A  A  A  T  A  T  T  T  T  C  C  C  A  T  T  A  G  T  G  T  T  T  C  T  T  T  T  A  A  T  T  T   FREQUENCY    habitats, situated only 10-20 metres from each other and separated by one or two dikes. We compare P. chalceus populations from three canals (CANAL1; CANAL2; CANAL3: Fig. 5) to three adjacent pond populations (POND1; POND2, POND3; see also [2]). Furthermore, we analyse Guérande P. littoralis populations from three different sites (GUE1, GUE2, GUE3; Fig. 5; see also [12]) nearby the aforementioned P. chalceus population couples.
The two related species are also studied on a macroscale with completely independent population samples (Guérande populations not included). Data on P. chalceus  5; see also [11,22]. For P. littoralis, six populations are analysed here, five of them from France (AUT, MSM, TOU, CAM, ROU; [12]). From these sites in France, we also sampled P. chalceus populations (see above). In Belgium, P. littoralis is critically endangered and the previous record went back to 1956 and was from Ostend [23]. Recently, a supposed new P. littoralis population has been discovered in Belgium and is also included here (   [2,11]. Temporary populations of P. chalceus are situated in the Mediterranean part of Europe or occur in small (<4 ha) and young (<400 years) Atlantic salt marshes. Stable and intermediate populations live in larger marshes situated along the Atlantic coast. The only difference between both salt marsh areas is their estimated age (Stable: >1000 years; Intermediate: between 400-1000 years). The age of salt marshes was estimated using historical information [24][25][26][27].

Morphological analysis
Body size (elytral length) and wing size were measured by means of a calibrated ocular under a binocular microscope. Generally, carabid wing size follows an allometric relationship with body size. [28] developed an index that corrects for this allometry, i.e. percentage of maximal realisable relative wing size. Relative wing size is wing length × width divided by elytral length × width. Relative wing size is then expressed as a percentage of the maximal wing size for a beetle of a given size. This index was shown to be an unbiased estimator for comparing different individuals, populations and species of carabid beetles [28]. In ground beetles, females are generally larger than males. Therefore, we analyse male and female body sizes separately. To be complete, we analyse female and male relative wing size also separately. Body size and relative wing size are compared between species and populations with ANOVA's. Total variance is partitioned among groups (species or species and ecotypes), among populations within groups, and within populations by carrying out a nested design ANOVA using STATISTICA (version 7.1; StatSoft Inc., Tulsa, UK) on both a micro-and macroscale.  [29]. Earlier work showed that one locus (IDH1) was non-neutral and we will always analyse it separately [11].

Genetic divergence
Departures from Hardy-Weinberg equilibrium expectation were tested with an exact test using the GENEPOP software (Version 3.2; [30]). Significance levels were adjusted by using sequential Bonferroni correction. Similarly as in the analyses for body and wing size, total genetic variance is partitioned among groups (species, ecotypes), among populations within groups, and within populations by carrying out a hierarchical analysis of molecular variance (AMOVA) using ARLEQUIN (version 3.000; [31]) on both a micro-and macroscale.
PCRs for nucleotide sequencing of COI utilized primers C1-J-1718 and C1-N-2191 and for 16S we utilized primers LR-J-1307 and LR-N-13398 [32]. DNA amplification reactions were performed in a 25 µL final volume. Each reaction mix contained 5 µL of extract, 1× buffer (Sigma), 1.5 mM MgCl 2 , 200 µM of each dNTP, 0.4 µM of each primer and 0.6 U RedTaq DNA polymerase (Sigma). Initial denaturation was for 2 min at 95°C, followed by 35 cycles of 1 min at 95°C, 1 min 30 s at 48°C (and 46°C for 16S), and 2 min at 72°C; 9 min at 72°C completed the program. The reaction was purified with columns following manufacturer's recommendations. Sequencing was done by BigDye Terminator version 3.1 kits on an ABI 3130 sequencer (Applied Biosystems). Sequences were aligned using BioEdit version 5.0.6 [33]. We tested for neutrality of mutations following Fu & Li's method with D* and F* test statistics using DNASP 4.0 [34,35]. A phylogeny of unique haplotypes was constructed from the calculated Kimura two-parameter distances using the neighbour-joining approach within MEGA ( [36]; 1000 bootstrap replicates).