Abstract
Cyclotus taivanus consists of five subspecies, with clear morphological diversity. The molecular phylogenetic relationships of this group have never been discussed before.
In order to investigate the relationships between C. taivanus ssp., I sequenced part of the mitochondrial COI and the 16S rRNA gene from 27 sampling sites. I also measured 9 shell traits for morphology analysis. Even though the morphology PCA revealed a more or less continuous distribution of individuals in morph-space, the two highly divergent haplotype clades in COI and 16S rRNA analysis indicated the presence of two independently evolving lineages. The sequence divergence between two clades was almost as high as between other Cyclophoridae species(COI: 0.102 in C. taivanus ssp., 0.028–0.172 in other Cyclophoridae species; 16S rRNA: 0.121 in C.
taivanus ssp., 0.041–0.166 in other Cyclophoridae species, see chapter 2, table 2.5).
Therefore C. adamsi should be a valid species. For the environmental analysis, temperature may be a limited element of the distribution of C. adamsi and C. taivanus group (C. t. dilatus, C. t. diminutus, C. t. peraffinis, and C. t. taivanus). The ecological divergence probably appears rule of speciation in my case. The PLS analysis results indicate, phenotypic plasticity may be a key element of variable shell in C. taivanus group. The speciation process is not complete among C. taivanus group, and the adaptation of climatic pressure continuing be a rule of speciation process.
Introduction
Cyclotus taivanus is widely distributed in Taiwan and Okinawa and five subspecies were described. In traditional classification, the north Taiwan subspecies C. t. adamsi were tall spire with uniform yellow or tawny yellow color, few ornamented with reddish brown zigzag pattern. The other four are flat spire. However, C. t. dilatus is wide extended lip, C. t. diminutus is very narrow lip, C. t. taivanus is intermediate, C.
t. peraffinis is similar to C. t. diminutus but with very polish periostracum (Lee & Wu
2001). The evolutionary processes shaping the phenotypic variation among and within relative landsnails have become a focus of research interest recently (Pfenninger &
Magnin 2001, Teshima et al. 2003, Pfenninger et al. 2003).
Two basic approaches have been taken to the study of variation in shell shape. One approach, involves the analysis of the distribution of shapes among taxa and hypothesizing the reasons for low or high frequencies of certain shapes. The second approach involves the study of correlations between variation in shape and environmental variation (Goodfriend 1986). Land snails with low vagility, therefore they can not avoid unfavorable conditions by migration but are doomed to extinction (Akcakaya & Baur 1996). It makes land snails are an ideal model to understand the roles of speciation history. Cyclotus and Cyclophorus are large cyclophorids occurring in East Asia. They occupy almost the same niche. However, the biogeographic pattern of Cyclotus is quite different from Cyclophorus. It may cause by different speciation process. Ring speciation had been described in Taiwan Cyclophorus (see chapter 3).
The phylogeny and speciation model of C. taivanus ssp. have not been investigated.
The objectives of this study were using the genetic information from COI and 16S rRNA haplotype analysis as a background to study quantitative morphologic differences of all described C. taivanus ssp. In particular, I focus on five questions: (1) How do currently biogeographical pattern of C. taivanus ssp. cause by? (2) What is the speciation process due to? (3) Does the observed phenotypic variation correspond to different evolutionary lineages? (4) How is the phenotypic variation distributed within and among populations in relation to the COI and 16S variation? And (5) Can I identify environmental variables that co-vary with population differences in morphology?
Materials and methods Populations sampled
In total, 27 populations were examined in Taiwan (Fig. 4.1) (Table 4.1). I do not get any living specimen (sampling site 2) and adult individual (sampling site 5, 18) in some sampling sites. It makes morphological analysis was performed on adult individuals of a subset from 25 locations, molecular analysis using 89 individuals from 26 locations.
through a Canon D1 digital camera from the aperture view. The paper prints were digitalized with a resolution of 300 pixel/inch. Eight shell morphology characters were measured by digital image analysis: H (height), W (width), AW (aperture width), AH (aperture height), BW (body whorl width), BH (body whorl height), LW (outer lip width), UW (umbilicus width) and HN (whorl number) (Fig. 4.2). H×W, H/W, AW/W, AW/AH, UW/W, LW/AW, BH/BW and (W－UW－LW)/HN were calculated as indices of size, tallness/flatness, aperture proportion, aperture shape, umbilicus proportion, outer lip proportion, flat degree of body whorl and the degree of whorl increasing, respectively. These shell morphology variables were used to perform a principal component analysis (PCA) using PCA option of the package XLSTAT (2007).
Partial least square (PLS) analysis was performed using the PLS option of the package XLSTAT (2007) to assess the correlations between the quantitative shell traits and environmental variables for each population.
DNA preparing and sequencing
The shells and soft parts were separated at library. Shells were well cleaned for identification and measuring its shell characters, soft parts were stored at -80℃ until DNA extraction. DNA was extracted from columellar muscle. I extracted DNA from separate individuals using TEK-based protocol (Jiang et al. 1997) with minor modifications. Tissue was placed in TEK buffer (12.5mM Tris-Cl pH 7.3, 2.5mM EDTA, 0.4% KCl), then ground with glass pestle and incubated at 57℃ with 20µl of proteinase K (20mg/ml) more than 2 hrs. The tissue extract was extracted at least twice with phenol and chloroform. 400µl DNA extract was precipitated by adding 1000µl pure ice-cold ethanol, and was then placed in -20℃ for 20 min. DNA was pelleted by centrifugation for 30 min. After 70% ethanol rinse, DNA was resuspended in distilled water and stored at -80℃ for DNA amplification. An approximately 860bp mitochondrial 16S rRNA gene was amplified by PCR using primers 16SRT (5’–ACA TAT CGC CCG TCA CTC TC–3') and 16SL900 (5’–AAA TGA TTA TGC TAC CTT TGC–3'), exactly 531bp COI using primer LCO1490 (5'–GGT CAA CAA ATC ATA AAG ATA TTG G–3') and HCO2198 (5'–TAA ACT TCA GGG TGA CCA AAA AAT CA–3') (Williams et al. 2003). Some individuals got pool PCR products using primer LCO1490 and HCO2198. For those I use primer CoL90 (5’-TAG TGT TAA AAT TAC GAT CAG T-3') and CoR600 (5’-AAG TCT TCT AAT TCG TGC AGA A-3').
PCR reactions contained template DNA 10–50ng/µl, 10pmol of each primer, 5µl 10×
reaction buffer (10mM Tris-HCl, pH9.0, 50Mm KCl, 1.5Mm MgCl2, 0.1% gelatin, 1% Trinton X-100), 0.4µl 25mM/µl dNTP, 0.2µl 50mM Mg2+ and 0.4µl Taq polymerase (5unit/µl) in a total volume of 50µl. Thermal cycling for 16S rRNA was
performed with an initial denaturation for 5min at 95℃, followed by 30 cycles of 30 sec. at 95℃, 45 sec. at 57℃, 50 sec. at 72℃ and ultimate extension at 72℃ for 10min, final hold in 4℃. Thermal cycling for COI gene was performed similarly, changing only its annealing temperature to 47℃. PCR products were purified using a purification kit (AMP PCR purification, Beckman) and then sequenced using an ABi 3700 autosequencer.
Phylogenetic analyses
Both COI and 16S rRNA sequences were combined with data of the out-group Ptychopoma wilsoni (Pfeiffer, 1865). Sequences were assembled and edited using Bioedit 5.0.9 (Hall 1999). All alignments employed Clustal X (Thompson et al. 1997) and were manually proofread. Codon positions within COI were tested using the incongruence length difference (ILD) test (Farris et al. 1995), as implemented by the partition homogeneity test in PAUP 4.0b10 (Swofford 1998) (100 replicates). COI sequence data were divided into two partitions, first and second codon positions in one and third codon positions in the other. Two parts of COI gene were congruent and all codon positions were combined and used in the following analysis.
All data sets were subject to Neighbor-joining (NJ) using PAUP 4.0b10, to the maximum likelihood (ML) analyses using PHYML 3.0 (Guindon & Gascuel 2003), to the Maximum parsimony (MP) analysis using MEGA 4.0 (Tamura, Dudley, Nei, and Kumar 2007). The substitution model used for COI data set corresponded to the Hasegawa, Kishino, Yano 85 model, and included invariable sites, and rate variation among sites (HKY+I+G); for 16S rRNA data set corresponded to the Transversional model, and included invariable sites, and rate variation among sites (TVM+I+G).
These were the best models found using Modeltest 3.06 (Posada & Crandall 1998).
Before model fitting, the full-length sequences were tested to confirm that there was no significant heterogeneity in base frequencies across taxa (in COI: X2 = 9.47, df = 96, P = 1; in 16S rRNA: X2 = 34.03, df = 117, P = 1). NJ and MP bootstraps consisted of 1000 iterations. Reliability of ML trees were estimated by the approximate likelihood ratio test (aLRT) (using custom define model, base frequencies: A = 0.2800, C = 0.1865, G = 0.1617, T = 0.3719 in COI and A = 0.3425, C = 0.0962 G = 0.1859, T = 0.3754 in 16S) using PHYML 3.0 (Guindon & Gascuel 2003). Gaps (from insertions/deletions) were treated as missing data.
populations were explored through Mantel-tests in XLSTAT (2007) on the pairwise genetic p-distance, pairwise shell characters similarity (using Euclidean distance calculated by Primer 5.1.2) and pairwise geographic distance.
Results
Morphology analysis
The first three principal components account for 49.59%, 22.38% and 9.68 % of the total variation in the matrix. The loadings of characters on the principal components reflect distinct groups of shell traits: umbilicus proportion (PC1), outer lip proportion (PC2), and aperture proportion (PC3). The C. adamsi group (sampling sites 1–8) are fully separated from the other four taxa (sampling sites 9–27) along the two axes, but C. peraffinis, C. dilatus, C. taivanus and C. diminutus show partial overlap (Fig 4.3).
Phylogenetic analysis
The aligned 531 bp COI gene data matrix, included 171 variable sites of which 147 (85.96%) were parsimony informative. No length difference from the out-group was detected among members of Cyclotus. The average p-distance was =0.087. Sequence divergence among the haplotypes ranged from 0.002 to 0.128. The alignment of the 16S rRNA gene fragment data yielded 816 characters of which 353 are variable sites.
Of these 353 positions, 307 (86.97%) were parsimony informative. The average p-distance among haplotypes was 0.094. Sequence divergence among the haplotypes ranged from 0.001 to 0.129.
The inferred phylogenetic trees between the haplotypes of COI gene are shown in Figs. 4.4, 4.5, 4.6 & 4.7 of 16S rRNA gene are shown in Figs. 4.8, 4.9, 4.10 & 4.11.
Both COI and 16S rRNA cladograms construct using different tree building method present similar topology. The topology of the COI and 16S rRNA gene trees strongly suggested the presence of two fundamental haplotype clades named A and B.
Analysing each clade separately, the average sequence divergence of COI among haplotype clades A and B was 0.102 changes per site compared to 0.009 and 0.057 changes within clade A and clade B, respectively. In contrast, the sequence divergence of 16S rRNA among haplotype clades A and B was 0.121 changes per site compared to 0.011 and 0.049 changes within clade A and clade B, respectively. All C. adamsi cluster in clade A. Clade B comprises C. peraffinis, C. dilatus, C. taivanus and C.
diminutus.
Environmental characterization of sites
PCA of the sites shows that the sampled populations primarily comprise (axis 1,
44.03% of total variation) an environmental gradient from north Taiwan, with low temperature, to south Taiwan, with high temperature (Fig. 4.12). The second most important gradient (axis 2, 32.42%) is a gradient with more changeable temperature (annual range of monthly mean temperature), lower warm season precipitation and higher cold season precipitation from north to south Taiwan.
When the results of the molecular species characterization are plotted on the PCA, it becomes clear that the sites inhabited by the two clades species in question are well separated by this analysis, mainly by the second axis. Clade A populations are preferentially situated on changeable temperature, lower warm season precipitation and higher cold season precipitation climate region of the plot. The ecological amplitude of clade A is much larger than that of clade B.
Mantel test
There are four subspecies including variable shell morphology in clade B. Therefore, I perform Mantel test to investigate the correlation between molecular data, morphology and geographic distribution, respectively.
The Mantel test for association between the geographic distance separating populations and pairwise p-distance was significant (COI p-distance / geographic distance Mantel test, r = 0.347, P = 0.0001; 16S rRNA p-distance / geographic distance Mantel test, r = 0.489, P = 0.0001). By contrast, there were no correlation between molecular pairwise p-distance and overall morphological similarity in both COI and 16S rRNA cases (COI p-distance / morphological similarity Mantel test, r = 0.029, P = 0.330; 16S rRNA p-distance / morphological similarity Mantel test, r = 0.019, P = 0.499). Mantel test was also performed between pairwise geographic distance and overall morphological similarity matrix, and got significant correlation between them.
PLS analysis
PLS analysis detected significant correlation between shell traits and environmental variables of clade B for the respective sampling site (Table 4.2). Populations with a few cold season amount of precipitation tend to be composed of individuals showing a broad outer lip, more rapid increasing shell whorls and more whorl number.
though the PCA revealed a more or less continuous distribution of individuals in morph-space, the two highly divergent haplotype clades in COI and 16S rRNA analysis indicated the presence of two independently evolving lineages. The sequence divergence between two clades was almost as high as between other Cyclophoridae species (COI: 0.102 in C. taivanus ssp., 0.028–0.172 in other Cyclophoridae species;
16S rRNA: 0.121 in C. taivanus ssp., 0.041–0.166 in other Cyclophoridae species, see chapter 2, table 2.5).
Despite of an almost continuous morph-space when looking at the overall phenotypic similarity, PCA yielded a highly discrimination model and the genetic clade as classification variable (Fig. 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 4.10, 4.11). This shows that the genetic clade of an individual is a good predictor of its morphology. The presence of two highly divergent genetic clades or evolutionary lineages is therefore in concordance with morphological differences. This led me to the conclusion that there are two separately evolving lineages in C. taivanus ssp. This view is strengthened by the observation that clade A, found only in sampling sites 1–8 from north Taiwan, occurs in parapatry to populations of clade B (Fig. 4.1). The absence of gene-flow between clade A and clade B suggested reproductive isolation of the clades (Table 4.3). The average sequence divergence in COI and 16S rRNA among clade A and B equalled 0.102 and 0.121 changes per site, respectively. Because of lack of calibration point, I have to assume that two clades diverged approximately 4.5–6.1 million years ago, if the Cyclotus species of investigation with similar mutation rate as Taiwan Cyclophorus (see chapter 3, 2.28% and 2% per million years in COI and 16S rRNA, respecitively). Most species concepts (reviewed in Hull 1997) would recognize clade A as separate species in the light of the presented evidence. However, additional evidence, like crossbreeding experiments and analyses of habitat differences, are needed to confirm the species status of the inferred clades. In order to easy describe, and by comparing genetic clade and shell traits PCA data, I will refer to all individuals with haplotype of clade A as C. adamsi, with haplotype of clade B as C. taivanus group (C. t. dilatus, C. t. diminutus, C. t. peraffinis, and C. t. taivanus) hereafter.
Portions of the molecular and phenotypic variation could be shown to be due to an obviously long separation of different evolutionary lineages that I can consider as separate evolutionary entities – C. adamsi and C. taivanus group. As this was a rather unexpected result, I could not adjust the sampling scheme towards the inclusion of more populations of C. adamsi into this study. This lack of data prevented me from exploring further the spatial distribution of phenotypic variation in this taxon and to compare it to the organization in C. taivanus ssp. An initial survey of the geographical distribution of taxon, that there are no distinct differences in habitat preference between the two lineages. However, C. adamsi and C. taivanus group are always
parapatric but never synpatric. In the scatter plot of environmental variables PC2 and shell traits PC1, C. adamsi seems to have a good adaptation to unstable, arid warm season, moist cold season and cold temperature, but not in C. taivanus group (Fig.
4.13). Temperature may be a limited element of the distribution of C. adamsi and C.
taivanus group. Opposite to C. adamsi, C. taivanus group with flatter shell inhabit in warmer and stable temperature environment. The warm and stable climate tends to have luxuriant vegetation. The flat shell would be at an advantage when roam over the ground in the luxuriant vegetation as keeled shell Cyclophorus species case (see chapter 3).
Hypothetic of C. taivanus ssp. speciation model
It appears that a climatic gradient is responsible for the distribution pattern of species (Fig. 4.12). Along this cline, C. adamsi occupies only the sites with the moist winter and arid summer climate. This type of climate seems to exclude C. taivanus group.
The actual level of water stress at a particular site depends strongly on the local microclimate, which may account for the observed intermingled pattern in the contact zone. Even though, the sampling population is few and limited, there is probably a cline in umbilicus proportion between populations near contact zone (Fig. 4.14).
Likewise, there are more or less clines in the other eight shell traits between populations near contact zone (Fig. 4.15). Even though premating or postmating reproductive isolation may have evolved as a byproduct of ecological divergence, the ecological divergence appears sufficient to prevent immediate contact and therefore acts as an effective barrier to mating. The pattern of distribution along a climatic gradient suggests an ecotone divergence. Clines or ecotones are thought to promote the evolution of specialized adaptations in populations on either side of the border (Pfenninger et al. 2003). The speciation between C. adamsi and taivanus group may due to ecological speciation. The alternative hypothesis is due to vicariance events and secondary contact. The alternative model of speciation offer much less likely scenarios. There was no evidence indicate presenting an obstacle for dispersal around the contact zone in historical geology.
Associations between molecular and morphological variability within C. taivanus group
p-distance and overall morphological similarity in shell traits could not be rejected within C. taivanus group. This suggested that the observed shell traits differences have not evolved in concordance with the mitochondrial genome and resulting phenotypic population structure does not reflect the phylogenetic history of the populations. The significant phenotypic divergence among populations must therefore have other explanations.
Associations between the environment and morphological variability within C.
taivanus group
The results from PLS analysis showed that differences between populations in some traits co-varied significantly with long-term climatic conditions of the sampling site. I found precipitation is worth to affect some shell traits of C. taivanus group. Cold season precipitation and out lip proportion are significantly negative correlation. By contrast, warm season precipitation and out lip proportion are significantly positive correlation. It is the same in whorl number, sampling sites with more annual amount precipitation were associated with fewer whorl numbers. This was astonishing at first sight because a humid climate allows more feeding activity and snails should therefore attain more whorl numbers. A similar tendency of smaller snail occurs in more rainfall habitat was reviewed by Goodfriend (1986), and possible explained as a potential adaptation of sexual maturity. There is no further shell growth when Cyclotus reaches sexual maturity. To reproduce earlier in wetter cold season would offer the next generation the possibility to reach a bigger size before summer. Early reproduction could thus constitute a selective advantage because small size juveniles of most snails are high mortality in dry summer of Taiwan (experience of field work).
The ecological divergence probably appears rule of speciation in C. taivanus ssp. case.
The speciation process is not complete among C. t. dilatus, C. t. diminutus, C. t.
peraffinis, and C. t. taivanus, and the adaptation of climatic pressure continuing is a rule of speciation process. However, further experiments are needed to determine whether the size/precipitation correlation in C. taivanus group is due to phenotypic plasticity in response to the prevailing conditions or whether it has an adaptational significance.
In conclusion, I consider the phenotypic and phylogenetic divergence between the two identified lineages as sufficiently large to propose the presence of two distinct evolutionary entities as valid species of C. adamsi and C. taivanus group. The speciation process between C. adamsi and C. taivanus group is complete, as evidenced by a lack of gene-flow (Table 4.3). The reproductive isolation of two taxa may due to divergent selection on traits in different environments. In another word, the speciation model of this case is ecological speciation. However, more sampling
sites around the contact zone should be contained in the further study. For the subordinate members of C. taivanus, there is no resolution in phylogeny analysis of the present investigation. To resolve the relationship of the C. taivanus group, more rapid evolved molecular marker is needed in further study. Besides, the hypotheses about correlation among C. taivanus group divergence and environmental variables
sites around the contact zone should be contained in the further study. For the subordinate members of C. taivanus, there is no resolution in phylogeny analysis of the present investigation. To resolve the relationship of the C. taivanus group, more rapid evolved molecular marker is needed in further study. Besides, the hypotheses about correlation among C. taivanus group divergence and environmental variables