extfieldsofstudy: Biology Medicine

s2fieldsofstudy:
Biology
Environmental Science

text:
Tracing glacial refugia of Triturus newts based on mitochondrial DNA phylogeography and species distribution modeling

Introduction The major climatic oscillations during the Quaternary Ice Age heavily influenced the distribution of species and left their mark on intraspecific genetic diversity. Past range shifts can be reconstructed with the aid of species distribution modeling and phylogeographical analyses. We test the responses of the different members of the genus Triturus (i.e. the marbled and crested newts) as the climate shifted from the previous glacial period (the Last Glacial Maximum, ~21 Ka) to the current interglacial. Results We present the results of a dense mitochondrial DNA phylogeography (visualizing genetic diversity within and divergence among populations) and species distribution modeling (using two different climate simulations) for the nine Triturus species on composite maps. Conclusions The combined use of species distribution modeling and mitochondrial phylogeography provides insight in the glacial contraction and postglacial expansion of Triturus. The combined use of the two independent techniques yields a more complete understanding of the historical biogeography of Triturus than both approaches would on their own. Triturus newts generally conform to the ‘southern richness and northern purity’ paradigm, but we also find more intricate patterns, such as the absence of genetic variation and suitable area at the Last Glacial Maximum (T. dobrogicus), an ‘extra-Mediterranean’ refugium in the Carpathian Basin (T. cristatus), and areas where species displaced one another postglacially (e.g. T. macedonicus and western T. karelinii). We provide a biogeographical scenario for Triturus, showing the positions of glacial refugia, the regions that were postglacially colonized and the areas where species displaced one another as they shifted their ranges.

Introduction
The Quaternary Ice Age (~2.59 Ma-present) is associated with large climatic oscillations [1][2][3][4]. Long, cold and dry glacial cycles are alternated by relatively short, warm and wet interglacials. The transition between the two takes place in a geological blink of an eye [5,6]. The climatic oscillations have a major impact on species distribution patterns [1]. Generally, impacts of glacial cycles are more extreme further away from the equator (and at higher elevations): areas at a higher latitude (and elevation) become inhospitable, whereas climate at a lower latitude (and elevation) remains habitable [1,7]. Populations at higher latitude may cope with climate change in situ, through adaptation or phenotypic plasticity, or alternatively they may track suitable habitat [6,8]. However, a more likely outcome is that such populations go extinct [8][9][10]. At lower latitudes, populations can endure glacial periods relatively unimpaired, in so called glacial refugia [2,10,11]. During subsequent interglacials, species can reclaim their former distribution, by rapidly recolonizing the large tracts of still uninhabited, but freshly habitable land from glacial refugia [1].
The fluctuating climate during the Quaternary has left its mark on patterns of genetic diversity [1][2][3]5,12]. Populations persisting in glacial refugia have a relatively long and stable demographical history compared to those in areas claimed postglacially. As a result, populations in glacial refugia are characterized by high levels of genetic diversity, whereas populations established after the most recent glacial cycle typically show little genetic variation. This is the conceptdevised from a northern hemisphere perspectiveof 'southern richness and northern purity' [1]. Furthermore, species displacement after postglacial secondary contact regularly coincides with introgression of genetic material (especially of cytoplasmic DNA) [12,13]. By uncovering the spatial structuring of genetic lineages between glacial refugia and along recolonization routes, and by detecting mismatches between genetic markers and species boundaries, phylogeographical surveys provide insights in past distribution rearrangement [9].
Independent from genetics, species distribution modeling can be applied to answer historical biogeographical questions [14]. Species distribution modeling involves the approximation of the ecological requirements of a species, based on the range of environmental conditions experienced at known localities [15]. The constructed model can then be extrapolated on current climate layers, to determine the species' potential distribution. Similarly, the model can be projected on climatological reconstructions of the past. Niches evolve over time, questioning the validity of predicting past distributions based on present day models. However, considering the relatively short time period spanning the shift from the last glacial cycle to the present day, niche conservatism is a realistic assumption [16,17]. Comparison of present and past potential distribution provides information on range shifts.
The distribution of amphibians is tightly linked to environmental conditions and thus has to follow suit in the face of rapid climate change [18][19][20][21]. The marbled and crested newts (genus Triturus) are a group of nine closely related species [22,23] (Table 1). Triturus newts are distributed across most of Europe and adjacent Asia ( Figure 1). They are found in regions generally regarded as important glacial refugia, such as the Iberian, Italian and Balkan Peninsulas, Anatolia, Caucasia and the southern Caspian basis [3,9,10]. On the other hand, Triturus newts occupy large tracts of land which would have been uninhabitable during glacial periods, particularly temperate Eurasia [24]. Thus, within this single model system, we expect to observe varying responses to glaciation.
We first conduct a dense phylogeographical survey: employing mitochondrial DNA, we determine geographical genetic structuring (i.e. diversity within and divergence among populations) for each of the different Triturus species. We also delineate areas where mitochondrial DNA has been asymmetrically introgressed. Subsequently, we approximate the distributions of the different Triturus species at the Last Glacial Maximum (~21Ka) using species distribution modeling. The outcome of the two independent approaches is visualized on composite maps, which show the nine different species side by side. Our hypothesis is that those areas predicted suitable at the Last Glacial Maximum based on species distribution models are also the ones showing the highest genetic diversity. Vice versa, we expect areas suggested to have been unsuitable at the time are genetically depleted. We summarize the signatures of past distribution dynamics as inferred from mitochondrial DNA phylogeography and species distribution modeling on a map.

Genetic approach
Generally the identification of Triturus newts is straightforward based on geography ( Figure 1). For individuals occurring near the contact zone more care should be taken. The Triturus individuals included were in the first place identified based on morphology [26,27] and for a subset nuclear genetic data confirming their identity [23] (Arntzen et al., submitted; Wielstra et al., unpublished data).
We included genetic data (658 bp of subunit 4 of the NADH dehydrogenase gene complex; ND4) for 2470 Triturus newts, representing 493 populations ( Figure 2 and Additional file 1). A large proportion of these individuals (n = 1795) was taken from previous studies (published [25,[28][29][30][31][32] or submitted [Arntzen et al.]; see Additional file 1 for details). The remainder (n = 675) *The taxon traditionally referred to as T. karelinii actually comprises three distinct species for which the taxonomy still has to be worked out fully [22,23]. For now these three species are referred to as the western, central and eastern T. karelinii species.
was newly sequenced for the current paper, using the primers in Table 2 and following the protocol outlined in [33]. Sequences were manually aligned and identical ones merged into haplotypes using MacClade 4.08 [34]. Our purpose was 1) to test for mitochondrial DNA introgression between species, and 2) to infer the spatial genetic variation within each species.
To detect which individuals retain their 'original' mitochondrial DNA and which possess introgressed mitochondrial DNA (derived from one of the other species), we first assigned mitochondrial DNA haplotypes to species by conducting a Neighbor Joining analysis in MEGA 5.05 [35]. Calotriton asper was used as outgroup (sequence taken from [28]; GenBank accession number: GU982378).
To infer interspecific geographical structuring, we first excluded introgressed mitochondrial DNA (as it does not properly reflect the evolutionary history of either 'host' or 'donor'). Subsequently, we determinedfor each species -1) the genetic diversity within populations and 2) the genetic distance among populations. The mean number of pairwise differences among haplotypes (π), as determined with Arlequin 3.5 [36], was used as a measure of genetic diversity within populations. To determine the genetic distance among populations we used Alleles in Space 1.0 [37]. Alleles in Space connects the populations in a network, based on Delaunay triangulation. Subsequently, the program produces values of average genetic distance among the populations that are connected by the network (the proportion of mismatched nucleotide sites; Z i ). These values are positioned at the midpoints of the connections in the network. For the measure of genetic diversity within populations we only included populations for which more than one sequence was available (π will always be zero for populations with only one sequence).
We interpolated the values for π and Z i across geographical space using inverse distance weighting in the spatial analyst extension of ArcGIS (www.esri.com) [38]. The output for each species was cropped according to its distribution range ( Figure 1). For both π and Z i , we compiled a single composite map from the nine crops (i.e. one for each Triturus species). We used a color scheme running from blue to red to reflect low to high values of π and Z i . Using a single scale for all Triturus species facilitates comparison among them, but runs the risk that variation in genetically poor species is overshadowed by that of genetically rich species. Hence, we additionally apply a species specific scale to better express intraspecific structuring. This way Figure 1 The distribution of the genus Triturus. Species mostly meet at parapatric contact zones, but note the area of sympatry of T. marmoratus and T. cristatus in western and central France. This map is adapted from [25].
we managed to visualize 1) regions of relatively high and low intraspecific genetic diversity and 2) the relative genetic divergence among populations within species.

Species distribution modeling approach
We composed a database of 4532 Triturus localities (incorporating the 493 populations included in the genetic analysis; see Figure 2 and Additional file 2). This database is based on [39] and updated in the light of the taxonomic developments [23] (see Table 1). For most species the majority of localities concern accurate locations, but in particularly for the wide ranging T. cristatus we mostly only had the name of the nearest village of a locality or digitized atlas data to our disposal. Species distribution models are still expected to perform relatively accurately if spatial error in localities is introduced and the method we use (Maxent, see below) is robust against such errors [40]. We used the bioclimatic variables available at a 2.5 arcminute resolution (c. 5 × 5 km) from the WorldClim database version 1.4 [41] to calibrate our species distribution models. In order to prevent model overfitting, which would negatively influence transferability [16,42], we minimized multicolinearity Figure 2 Maps showing the Triturus localities used for the mitochondrial DNA phylogeography and species distribution modeling. The inset shows part of the localities used in the genetic analyses (white circles) and the additional ones used for species distribution modeling (black circles). Two cut-outs (A and B) show more details on introgressed mitochondrial DNA: populations containing foreign mitochondrial DNA are labeled with a colored star (with the color denoting the 'donor' species) and populations containing two mitochondrial DNA types (of which one mostly, but not exclusively, is of the original species) as a black star; populations containing original mitochondrial DNA are again labeled with white circles. The colors used correspond to Figure 1. Details on localities can be found in Dataset S1 and S2. among data layers by selecting a subset that showed a Pearson's correlation of r < 0.7. Furthermore, focusing on climate layers that are deemed biologically meaningful based on life history knowledge of the model system yields the most appropriate species distribution models [43]. Variation in the availability of standing water during the breeding season appears to have driven the ecological radiation among Triturus newts [25]. Taking these two points into consideration, we included a set of layers that encompasses seasonal variation in evaporation and precipitation, in casu bio10 = mean temperature of warmest quarter, bio11 = mean temperature of coldest quarter, bio15 = precipitation seasonality, bio16 = precipitation of wettest quarter, and bio17 = precipitation of driest quarter. Bioclimatic variables are also available from the WorldClim database for the Last Glacial Maximum (~21 Ka). These data are derived from the Paleoclimate Modelling Intercomparison Project phase 2 [44]; [http://pmip2. lsce.ipsl.fr/] and based on two climate simulations: the Model for Interdisciplinary Research on Climate version 3.2 (MIROC) and the Community Climate System Model version 3 (CCSM) [45].
Species distribution models were created with Maxent 3.3.3 k [46]. We restricted the feature type to hinge features as this produces smoother model fits, so forcing models to be focused on key trends rather than potential idiosyncrasy in the data. This approach facilitates extrapolation to a different time (or place) [47]. The environmental range covered by the pseudo-absence data, used to discriminate presence data from background, should neither be too narrow nor too broad [48][49][50]. Too narrow a range results in complex models that do not generalize well, whereas too broad a range results in too simple models that focus on coarse-scale and neglect fine-scale variation. A practical solution is to focus on area that is potentially accessible to the species of interest if it were not for abiotic factors, i.e. an area where spread would not be hampered by major physical barriers. However, competition could lead to further exclusion of the target species, even though the area would be suited in the absence of such biotic interactions [51]. Taking these considerations into account, we restricted the area from which pseudo-absence was drawn to the distribution of the entire genus Triturus. This area was broadly defined as a 200 km buffer zone [48] around known Triturus localities (Additional file 2).
The species distribution models were tested for statistical significance against a null model derived from random localities [52]. For each tested species distribution model, we created a null distribution of 99 AUC values. These AUC values were derived from species distribution models, based on as many random localities as used for the tested species distribution model. The AUC value of the tested species distribution model was treated as a 100 th value and deemed statistically significant if it ranked higher than the 95 th value (i.e. above the 95% confidence interval). Random point data were created with ENMTools 1.3 [53]; [http://enmtools.blogspot.com/]. The null model approach prevents interpreting model quality based on an arbitrary AUC threshold and precludes the requirement to set aside part of the localities for model testing [52].
The species distribution models were projected on the current and Last Glacial Maximum climate layers. Composite maps were created for each set of climate layers (in the same way as explained for the genetic approach). Maxent provides predicted probability values between zero and one and we used a color scheme running from blue to red to reflect these values.

Results
The 2470 Triturus sequenced newts comprise 315 haplotypes (see Additional file 1 for details and Additional file 3 for GenBank accession numbers). The Neighbor Joining phylogeny shows that haplotypes group in nine reciprocally monophyletic lineages, corresponding to species (Figure 3). A large number of individuals (n = 408; c. 16.5%) that belong to one species, possess mitochondrial DNA characteristic of another ( Figure 2 and Additional file 1). This mitochondrial DNA introgression is mostly restricted to near the contact zones; only T. pygmaeus, T. macedonicus, T. cristatus and central T. karelinii have extended ranges in which foreign mitochondrial DNA is present (derived from T. marmoratus for the first one and from western T. karelinii for the other three).
The number of populations included to determine genetic diversity within (π) and genetic divergence among (Z i ) populations was 44 and 48 out of 52 populations for T. carnifex, 60 and 88 out of 104 for T. cristatus,   Figure 5 (e.g. in the case of two allopatric populations). Uncropped maps for each species, showing a more comprehensive picture per species but at the cost of conciseness, can be found in Additional file 5. Triturus species generally show considerable spatial variation in their genetic composition (reflected by 'warm' and 'cold' areas in Figures 4 and 5). Species distribution models perform significantly better than random (Additional file 6). Composite maps depicting the predicted suitability of each species' distribution range at the Last Glacial Maximum, based on the two different climate simulations (MIROC and CCSM), are provided in Figure 6. As species were not necessarily bound to their current ranges through time, projections for each species on a wider area (roughly the distribution of the entire genus Triturus) are provided in Additional file 7. The ranges of all Triturus species are predicted to have been restricted at the Last Glacial Maximum. A composite map showing species distribution models projected on present day climate layers is provided in Figure 6. Predicted suitability of the distribution range of each species under current conditions shows a considerable overlap with the sketched outlines of their ranges (Figure 1).

Discussion
Before we discuss our results, we first consider where care should be taken when interpreting genetic or species distribution modeling results in isolation and explain how interpretation of the results of the two independent techniques together at least partially negates these problems. Then we identify the probable glacial refugia and postglacially colonized area for each Triturus species. Finally, we summarize our scenario regarding the biogeographical response of Triturus newts to the climate shift from the previous glacial period to the current interglacial on a map.

Considerations and combination of the genetic and species distribution modeling results
A high value for genetic diversity within populations ( Figure 4) can result from two distinct processes: postglacial secondary contact from distinct source populations or long term presence in situ throughout glacial-interglacial cycles (potentially enforced by a low connectivity between populations due to e.g. xeric environment) [9]. To separate these two processes, phylogenetic relationships among haplotypes provide some insight. For example, the obvious 'hotspot' found in T. cristatus (SE Poland) likely reflects the comingling of haplotypes belonging to the two distinct clades present in this species ( Figure 3; Additional file 1). On the other hand, the hotspot found in eastern T. karelinii (SE Azerbaijan) likely represents standing genetic variation, given the lack of phylogeographic structure among the haplotypes involved. Species distribution modeling can further help to distinguish both processes by establishing whether an area became suitable only recently or whether it has continuously been so on the long term. Figure 4 The geographical distribution of genetic variation for the different Triturus species. These are composite maps for all nine Triturus species. For each species, the genetic variation within each population (π) was determined and subsequently interpolated across its distribution range. In Figure 4a we use a single scale for all Triturus species (allowing direct comparison among species) whereas in Figure 4b we use a species specific scale (better expressing genetic structure in genetically relatively poor species). Warmer colors refer to a higher genetic diversity. The insets show the situation for T. cristatus (left) and T. marmoratus (right) in their area of sympatry (shown in dark gray on the main maps).
Maps showing genetic divergence among populations ( Figure 5) reflect a different aspect of geographical genetic structuring. On the one hand, areas with a high genetic overturn among populations (e.g. the southern part of the range of T. macedonicus) show many hotspots. Such areas would be expected to have been habitable on the long term. On the other hand, genetically homogenous areas, due to high levels of gene flow and/or demographic expansion stand out as cold areas (e.g. both of these processes likely determine the absence of genetic structuring among T. dobrogicus populations). These two processes are more difficult to disentangle: whereas demographic expansion would particularly be expected in areas that became suitable postglacially, gene flow could reduce genetic structuring both inside and outside of glacial refugia. In some cases geographical barriers to gene flow stand out, such as the Greater Caucasus mountain range separating eastern T. karelinii on its northern and southern side. Such barriers reflect climatically unsuitable regions and can be identified as such in the species distribution models.
The reliability of genetic diversity values for populations depends on sample size, and some populations are relatively poorly sampled in our study. The genetic turnover between populations provides additional insight that counters the adverse effects of low sample sizes. If populations within an area are genetically diverse, but sample sizes per population are low, genetic diversity per population might be underestimated due to chance effects in sampling. However, these effects are random per population. Therefore, the more populations that are being compared within that area, the less likely it is that actual genetic diversity remains hidden by chance. For Figure 5 The genetic (dis)similarity among populations within the different Triturus species. These are composite maps for all nine Triturus species. For each species, the genetic divergence among populations (Z i ) was determined and subsequently interpolated across its distribution range. In Figure 5a we use a single scale for all Triturus species (allowing direct comparison among species) whereas in Figure 5b we use a species specific scale (better expressing genetic structure in genetically relatively poor species). Warmer colors refer to a higher genetic divergence. The insets show the situation for T. cristatus (left) and T. marmoratus (right) in their area of sympatry (shown in dark gray on the main maps).
this same reason determining the genetic turnover between populations is less susceptible to sample size to begin with: low sample size per population adds up to a big sample sizes within a region. So, in short, if multiple populations within a region are found to be genetically similar, even though samples sizes per population are low, together they provide a good indication that genetic diversity in the region really is low. We thus recommend exploring genetic diversity within populations and genetic divergence among populations together to improve Figure 6 The predicted suitability of each Triturus species' range at the Last Glacial Maximum (MIROC and CCSM model) and at the present. This is a composite map for all nine Triturus species. For each species, its species distribution model was projected on MIROC (a) and CCSM (b) Last Glacial Maximum climate layers and on current climate layers (c), cut according to its current distribution range. Warmer colors refer to a higher predicted suitability. The insets show the situation for T. cristatus (left) and T. marmoratus (right) in their area of sympatry (shown in dark gray on the main map). the estimate of spatial genetic variation, especially when sample sizes are low.
We present the results of a single (mitochondrial) genetic marker. Although we consider our results indicative (and use the species distribution modeling results as an independent check), it should be taken into account that a single marker may not accurately reflect species history [54]. In fact, we explicitly exploit the mismatch between species history and mitochondrial DNA in the case of asymmetrically introgressed mitochondrial DNA. Asymmetrically introgressed mitochondrial DNA may suggest that species displaced one another (e.g. [29,55]). For Triturus, the similarity of the introgressed mitochondrial DNA to that present in the 'donor' species suggests a recent (postglacial) transference across the species boundary. Given the climatic oscillations during the Quaternary Ice Age, the Triturus newt species are presumed to have been in periodic contact through time. However, we did not identify ancient mitochondrial DNA introgression events (as in e.g. [56][57][58]). A potential explanation is that introgressed mitochondrial DNA was limited to areas of postglacial expansion and never penetrated those areas that repeatedly acted as glacial refugia. In effect, any mitochondrial DNA that became introgressed during an interglacial would be erased by the next glacial period. Species distribution modeling allows us to test this hypothesis.
A disadvantage of species distribution models is that they cannot distinguish between potential and realized distribution: areas that had suitable conditions at the Last Glacial Maximum may not have actually been inhabited. Genetic information can narrow down the position of glacial refugia (e.g. both France and the Carpathian Basin are suggested to have had suitable conditions during the Last Glacial Maximum for T. cristatus but only the latter shows the high genetic variation expected to be present in a former glacial refugium). Another issue is that species distribution models are based on the climatic conditions currently experienced by a species, but the actual conditions under which they could occur might be broader (e.g. no suitable conditions were predicted to have been present at the Last Glacial Maximum for T. dobrogicus, but yet it did survive to the present day).

Glacial refugia and area postglacially colonized for each Triturus species
The two marbled newt species illustrate the different effect of glaciations on species with a relatively southern and relatively northern distribution. The northern T. marmoratus is predicted to have withdrawn its range at the Last Glacial Maximum more extensively than the southern T. pygmaeus (especially according to the CCSM climate simulation; Figure 6). In T. marmoratus genetic structuring is highest in the south of its range, but for T. pygmaeus it is not (Figures 4 and 5). The presence of T. marmoratus mitochondrial DNA in the northern part of the range of T. pygmaeus suggests that in this area the former was replaced by the latter [31,59,60]. This is confirmed by the distribution models projected on Last Glacial Maximum climate data outside of the current species ranges (Additional file 7): at the time, area predicted suitable for T. marmoratus (but not for T. pygmaeus) was present in the northern part of the area currently occupied by T. pygmaeus. Striking is that both T. marmoratus and T. pygmaeus mitochondrial DNA occurs throughout this area (Figure 2). Although a species will often simply show foreign mitochondrial DNA where it replaced another [13] see below for examples in Triturus, exceptions similar to the T. marmoratus-T. pygmaeus case are known (e.g. [61]).
The species distribution models based on the two climate simulations agree that most of the range of T. cristatus was unsuitable at the Last Glacial Maximum ( Figure 6). Suitable area was suggested to be present at the Last Glacial Maximum in both the Carpathian Basin and in France. The genetic data suggest that only the Carpathian Basin acted as a glacial refugium at the Last Glacial Maximum as genetic diversity here is high (Figures 3, 4, 5). This patterns corresponds to the geographical distribution of genetic variation based on nuclear DNA markers found in previous studies (MHC genes, allozymes and microsatellites [62,63]). This also suggests that the remainder of T. cristatus' current range (including France) was colonized post-glacially, as the species is genetically depleted outside of the Carpathian Basin. Although T. cristatus shows asymmetrically introgressed western T. karelinii mitochondrial DNA south of the Danube river, relatively distinct T. cristatus haplotypes also occur here (Figures 2 and 3, Additional file 1), suggesting a complex pattern of both range expansion at the cost of western T. karelinii but also long term presence south of the Danube.
For T. carnifex there is a discrepancy between the two climate simulations, with MIROC showing a more extensive predicted distribution at the Last Glacial Maximum than CCSM. Particularly CCSM suggests a more considerable range reduction and no suitable area in the current northern Balkan range ( Figure 6). The genetic data better correspond to the modeled distribution based on MIROC, as it shows distinct genetic clusters distributed on both sides of the Adriatic Sea and high genetic diversity in Italy (Figures 3, 4, 5), supporting long term presence in the Balkans and a stable population in Italy (see [32] for more detail). Both climate simulations agree that T. carnifex colonized the part of its range east and north of the Alps postglacially. The T. carnifex mitochondrial DNA introgressed into T. cristatus where the two species meet suggests that Balkan T. carnifex were involved in this range expansion.
Compared to the other Triturus species, intraspecific genetic diversity (i.e. phylogenetic depth) is low in T. dobrogicus (Figure 3). Triturus dobrogicus shows a 'starburst' phylogenya large number of similar haplotypessuggesting a demographic expansion after a bottleneck. At the Last Glacial Maximum, the range of T. dobrogicus is predicted to have been unsuitable ( Figure 6). This is in line with a population bottleneck. However, the species distribution models do not make it possible to pinpoint a potential refugium: the entire currently occupied range was predicted unsuitable at the Last Glacial Maximum. Suitable area was also not available outside the current range (Additional file 7). Evidently, T. dobrogicus must have had a refugium somewhere and given its specialization on river floodplains (e.g. [64]) this was likely positioned within its current range. Furthermore, compared to the other Triturus species T. dobrogicus does not show reduced genetic diversity of nuclear DNA (allozymes [65]). This suggests that perhaps the actual ecological niche of T. dobrogicus is broader than is currently realized (and used to construct the species distribution models). Another option would be that the niche of T. dobrogicus has evolved in a brief time interval to adapt to the environmental shift [66]. We consider this less likely as T. dobrogicus must have survived multiple glacial-interglacial cycles in situ.
The western part of the current range of T. macedonicus is predicted to have been suitable during the Last Glacial Maximum ( Figure 6). This does not fully correspond to mitochondrial DNA data, which shows structuring in the southern part of the range, including the southeast, but not in the northwest (Figures 4 and 5). Genetic variation in T. macedonicus is the highest of all Triturus newts (Figure 3). In a large part of its range, T. macedonicus possesses western T. karelinii mitochondrial DNA (Figure 2). This can be explained by T. macedonicus displacing western T. karelinii there after the Last Glacial Maximum as the area only became unsuitable for T. macedonicus after the Last Glacial Maximum (Figure 6; see [29] for a detailed analysis).
In western T. karelinii, a basal split separates a clade occurring in European Turkey and extreme southeastern Bulgaria (Figure 3). The species distribution modeling approach suggests suitable area was present here during the Last Glacial Maximum (Figure 6). The genetic distinction of this clade is expressed by the 'barrier' in Figure 5 surrounding the populations in which it is represented. The other clade shows extensive structuring in Asiatic Turkey, suggesting the region acted as a refugium. Suitability at the Last Glacial Maximum is confirmed by the species distribution modeling approach. A single sub-clade derived from (i.e. nested within) Asiatic Turkey has recently colonized the Balkan part of the range; it shows a starburst phylogeny in line with demographic growth during a range expansion [12]. The species distribution modeling approach agrees that the area only became hospitable after the Last Glacial Maximum.
Central T. karelinii shows a fragmented distribution at the Last Glacial Maximum (Figures 6 and 7). This is reflected by an east west divide (the 'barrier' in Figure 5) between two distinct clades separated by a basal split (Figure 3), which are currently in secondary contact (the 'warm' spot in Figure 4). Central T. karelinii contains western T. karelinii mitochondrial DNA in northwestern Asiatic Turkey (Figure 2), suggesting postglacial species displacement (see [23] for a detailed analysis).
For eastern T. karelinii, the southern Caspian basin and western Georgia (Colchis) are suggested to have harbored suitable area at the Last Glacial Maximum according to species distribution models ( Figure 6). Genetic diversity is particularly high along the southern Caspian Sea shore (Figures 3, 4, 5). Relatively recently a single nested clade colonized the Caucasus and Crimea from the Caspian part of the range [28]. However, this clade is still distinct (reflected by a relatively deep coalescence; see Figure 3), suggesting that colonization happened during a previous interglacial and thus that the clade persisted here in a glacial refugia during one or several glacial cycles.

Historical biogeographical scenario
We provide a historical biogeographical scenario for the genus Triturus in Figure 7, summarizing glacial refugia, areas of postglacial expansion and regions where species displaced one another. The genetic and species distribution modeling results provide insight into the distribution dynamics of Triturus in response to climate change during the Quaternary Ice Age. Triturus conforms to the general pattern of the Iberian, Italian and Balkan Peninsulas functioning as glacial refugia [1,67]. Furthermore, several regions that are increasingly appreciated as glacial refugia are also identified as such for Triturus: Anatolia, the southern Caspian Basin and the Caucasus region (e.g. [10,68]). The novel idea of extra-Mediterranean refugia, positioned more to the north [69], also applies to Triturus: the glacial refugium of T. cristatus was situated in the Carpathian Basin.
Similarly, areas typically identified as having been postglacially colonized from glacial refugia are predicted as such for Triturus: T. marmoratus, T. carnifex and (to a lesser extent) 'eastern T. karelinii' extended their ranges into temperate Europe from their southerly positioned refugia. However, postglacial expansion is best exemplified by T. cristatus. The contemporary range of T. cristatus outside its glacial refugium in the Carpathian Basin encompasses a huge (c. 4.75 million km 2 ) postglacially acquired area, stretching all the way from western Europe to Scandinavia and central Russia. This gives an indication of the speed with which postglacial colonization can be accomplished.

Conclusions
The combination of phylogeography and species distribution modeling aids locating of glacial refugia [70]. The two approaches to visualize intraspecific geographical genetic structuring together illustrate which areas are genetically rich and thus suspected to reflect long term inhabited area and which areas are genetically poor and thus probably only recently became habitable (Figures 4  and 5). Additionally, asymmetrically introgressed mitochondrial DNA suggests where species displaced one another upon post-glacial secondary contact ( Figure 2). Independently, the two reconstructions of potential distribution at the Last Glacial Maximum (based on the MIROC and CCSM climate simulations) provide an indication of which areas were habitable at the time and which were unsuitable for Triturus ( Figure 6). We provide a qualitative comparison of the results of the two independent techniques (summarized in Figure 7). We anticipate that future developments in the synthesis of phylogeography and species distribution modeling will make it possible to integrate the two approaches into a single analysis even further.


extfieldsofstudy: Chemistry Medicine

s2fieldsofstudy:
Biology
Materials Science
Medicine

text:
Activation of the ERK1/2 Signaling Pathway during the Osteogenic Differentiation of Mesenchymal Stem Cells Cultured on Substrates Modified with Various Chemical Groups

The current study examined the influence of culture substrates modified with the functional groups –OH, –COOH, –NH2, and –CH3 using SAMs technology, in conjunction with TAAB control, on the osteogenic differentiation of rabbit BMSCs. The CCK-8 assay revealed that BMSCs exhibited substrate-dependent cell viability. The cells plated on –NH2- and –OH-modified substrates were well spread and homogeneous, but those on the –COOH- and –CH3-modified substrates showed more rounded phenotype. The mRNA expression of BMSCs revealed that –NH2-modified substrate promoted the mRNA expression and osteogenic differentiation of the BMSCs. The contribution of ERK1/2 signaling pathway to the osteogenic differentiation of BMSCs cultured on the –NH2-modified substrate was investigated in vitro. The –NH2-modified substrate promoted the expression of integrins; the activation of FAK and ERK1/2. Inhibition of ERK1/2 activation by PD98059, a specific inhibitor of the ERK signaling pathway, blocked ERK1/2 activation in a dose-dependent manner, as revealed for expression of Cbfα-1 and ALP. Blockade of ERK1/2 phosphorylation in BMSCs by PD98059 suppressed osteogenic differentiation on chemical surfaces. These findings indicate a potential role for ERK in the osteogenic differentiation of BMSCs on surfaces modified by specific chemical functional groups, indicating that the microenvironment affects the differentiation of BMSCs. This observation has important implications for bone tissue engineering.

Introduction
Bone tissue engineering is an attractive approach with great potential for repairing bone defects. The basic strategy of tissue engineering for the structural and functional restoration of damaged bone involves the use of artificial biomaterials as scaffolds for the culture of specific types of cells, which are then stimulated with growth factors [1]. On the basis of their promising osteogenic differentiation ability, bone marrow stromal cells (BMSCs) have been widely used as seed cells in bone tissue engineering, either in vitro or in vivo [2,3]. The hypothesis that stem-cell differentiation can be induced by parameters, including structure, chemistry, mechanic, and molecule delivery of biomaterials, has already been substantiated [4][5][6][7][8][9][10][11]. Biomaterial technologies provide the exciting possibility of deconstructing and then reconstructing niches, allowing quantitative analysis of stem-cell behavior in a manner not previously possible [4].
Current research has confirmed that the surface chemistry of the substrates influences the adsorption and the conformation of extracellular matrix protein, such as fibronectin (FN) [12]. Binding of specific integrin adhesion receptors to this protein modulates the activity of focal adhesion kinase (FAK) and the intracellular signaling cascades of osteoblast-and myoblast-like cells [13,14]. As one of the primary transducers of integrin signals to the cell nucleus, the mitogen-activated protein kinase (MAPK) pathway provides a plausible link between cell surface integrin activation and subsequent stimulation of core binding factor alpha 1 (Cbf -1) in osteoblast-like cells [15,16]. Extracellular signal-related kinase (ERK1/2) is an important member of the MAPK family of serine/threonine protein kinases [17]. Appropriate soluble and insoluble signals can activate ERK1/2, which then activates important cellular processes, including key transcriptional and phenotypic differentiation programs via sequential activation of protein kinase cascades [18]. The activation of the ERK1/2 pathway has been reported to be involved in the osteogenic differentiation of BMSCs and differentiation of the 3T3L1 preadipocyte cell line into mature adipocytes [19][20][21]. The regulation of ERK1/2 also contributes to the balance between the osteogenesis and adipogenesis of human adipose-derived stem cells (hASCs) [18,22]. However, the mechanism by which the surface chemistry of the substrate influences osteogenesis of BMSCs is still unknown, and the effects of ERK1/2 pathway on the stem-cell fate are also unclear. To investigate these processes, we used self-assembled monolayers (SAMs) of alkanethiolates on gold surfaces to evaluate the effects of surface chemistry on the stem cells. Previous studies of biomaterials have routinely used protein-coated dishes to mimic the cellular microenvironment and to elucidate the functions of the matrix in regulating cellular processes. However, the adsorption of proteins is a complex process and often proceeds with a lack of control over the orientation and conformation of proteins at the surface. It is therefore difficult to control the biological activities of proteins adsorbed to man-made materials, and it is impractical to use these substrates as models of the extracellular matrix (ECM) [23]. SAMs, a recently developed technology, offer a set of monolayers which surface chemistry characteristics are well suited to the studies of cell-ECM interactions. Using SAMs technology, previous studies have demonstrated that surface chemistry regulates the specific differentiation of stem cells [24][25][26][27][28].
In the current study, we used a novel SAMs technology to precisely control the density of specific functional groups on the substrate. Multiple phenotypic features of rabbit BMSCs, including initial adhesion, morphology, long-term growth, gene expression, and functional differentiation, were examined. We also analyzed the relationship between the osteogenic differentiation of BMSCs and activation of the ERK1/2 pathway as a function of culture on SAMs substrates presenting different surface chemistries. These results provide support for future studies to identify surface chemistry parameters to achieve the optimal control over stem-cell osteogenic differentiation.
Ambient air-water-substrate contact angle measurements (4 L of ultrapure water) were taken with the OCA20 contact angle system (Dataphysics, German) fitted with a digital camera and analyzed using in-house image analysis software. Photoelectron Spectroscopy (XPS) analysis was carried out on Escalab 250 Xi (VG, UK) with Carbon (284.8 eV) as a marker, and the examined area was 300 m × 700 m. The surfaces with different terminal functional groups were observed using an atomic force microscope (MFP-3D-S, Asylum Research, USA) in contact mode under an air atmosphere. The photographs were taken under the openloop condition, and an Olympus AC240TS probe with Si3N4 cantilevers was used.

Culture and Harvest of BMSCs.
Rabbit bone marrow stromal cells (BMSCs) were isolated from the tibias of 4month-old New Zealand white rabbits. Briefly, the tibia of an anesthetized rabbit was punctured with a 16-gauge needle, and 4-5 mL of bone marrow was aspirated into a 10 mL syringe containing 5000 units of heparin. The bone marrow was filtered through a cell strainer (40 m), and the cells were resuspended and cultured in Dulbecco's modified Eagle's medium (DMEM, Gibco, Rockville, MD, USA) containing 10% fetal bovine serum (FBS, Hyclone, Logan, UT, USA), 100 IU/mL penicillin G, and 100 g/mL streptomycin. The medium was changed every 2 days. When 70-80% of the bottom of the dish was covered by the cell monolayer, the cells were enzymatically harvested and centrifuged at 800 g for 5 min. All experiments were conducted using cells below passage 4. The cells were then reseeded in standard growth medium at 5 × 10 4 cells/cm 2 onto substrates modified with different chemical functional groups and cultured to permit attachment and spreading.

Cell Viability Assay.
Cell viability was determined with the CCK-8 cell viability assay (Dojindo Laboratories, Kumamoto, Japan) according to the manufacturer's protocol. Briefly, after incubation at 37 ∘ C/5% CO 2 for 1, 3, 5, and 7 days, the culture medium was replaced with 300 L of medium containing 30 L of CCK-8 per well, and the cells were incubated for 3 h in the incubator. Finally, the cell viability was measured as the absorbance at 450 nm using a microplate reader (Bio-Rad, Model 680, Hercules, CA, USA).

2.5.
Real-Time RT-PCR. After 1, 3, 5, 7, 10, and 14 days of culture, the cells were harvested, and total RNA was extracted using the Trizol reagent (Invitrogen). The RNA concentrations were determined spectrophoto-metrically at 260 nm prior to reverse transcription. Real-time RT-PCR was performed on a quantitative real-time amplification system (MxPro-Mx3000, Stratagene, La Jolla, CA) using a SYBR PrimeScript RT-PCR Two-Step Kit (TaKaRa, Otsu, Shiga, Japan) according to the manufacturer's instructions. Equal amounts of total RNA (1 g) from each sample were converted into cDNA in a reverse transcription reaction. PCR for each gene was carried out in a 20 L reaction mixture containing 2 L cDNA template. The PCR conditions were 95 ∘ C for 1 min followed by 40 cycles of 95 ∘ C for 10 s and 60 ∘ C for 40 s. The 2 −ΔΔCT method was used to calculate the relative levels of gene expression using Strata gene MxPro QPCR software v3.00 (Stratagene). The housekeeping gene GAPDH expression was used as the internal control. A standard melting-curve cycle was used to confirm the quality of amplification. The reactions were performed in triplicate for each sample. The gene-specific primers were designed with the Primer premier 5 software (Premier Biosoft International, Palo Alto, CA, USA), as presented in Table 1.

Statistical Analysis.
All experiments were performed in triplicate, and representative results are presented as the means ± standard deviations. The statistical analyses for quantitative assays were performed using the SPSS 11.0 software (SPSS, Inc., Chicago, IL, USA).
The ratio of S/Au atoms analyzed by XPS and the XPS spectra demonstrated that the four types of functional groups successfully self-assembled. The four surfaces had similar S/Au ratios with values of approximately 0.22, indicating the similar surface density of functional groups. AFM images displayed the same morphology of these functional groups surfaces. The distance between neighboring functional group was approximately 0.5 nm, corresponding with the results of Widrig et al. [33]. The AFM images indicated that the four functional groups consisted of the well-known ( √ 3× √ 3) R30 ∘ structure unit, and the densities of four functional groups were 5 × 10 18 /m 2 [31].

Cell Viability on the Different Chemical Functional Groups.
The effects of modifying the substrates with the different chemical functional groups on the proliferative activity of BMSCs were assessed by CCK-8 analysis. Figure 1 showed the CCK-8 conversion by BMSCs cultured on the various types of SAMs on days 1, 3, 5, and 7, in conjunction with a clean glass (TAAB) control. On day 1, the levels of cell viability on the -NH 2 , -OH, and -COOH substrates were significantly higher than that on TAAB; the levels of cell viability on the -CH 3 substrates were significantly lower than that on TAAB. On days 3, 5, and 7, the cell viability significantly upregulated compared to their respective values on day 1 for cells grown on all of the chemically modified substrates. The greatest increase in cell viability was observed on the -NH 2 substrate, for which the value on day 7 was approximately 4fold higher than the value observed on day 1. In contrast, the smallest upregulation in cell viability was observed on the -CH 3 substrate, for which the value after 7 days of cultures was approximately twice that on day 1. The BMSCs exhibited substrate-dependent cell viability of -

Cytoskeleton and Focal
Adhesions on the Different Chemical Functional Groups. Figures 2 and 3 showed the typical morphologies of BMSCs on the different chemical substrates. The cytoskeletal structure was examined using a greenfluorescent F-actin stain, and FAK structure was examined using anti-vinculin antibody stained using a red fluorophore.
The cells plated on -NH 2 -and -OH-modified substrates were well spread and homogeneous with the actin organized into stress fibers and demonstrated long microtubules and distinct focal adhesions. In contrast, the cells on the -COOHand -CH 3 -modified substrates showed a more rounded phenotype. After 1 day, the cells cultured on the -NH 2 -and -OHmodified substrates demonstrated good surface coverage and clear evidence of stress fiber and focal contact formation. The images clearly showed that the focal contacts were located at the ends of the stress fibers (Figures 2(h) and 2(p)). The cells cultured in contact with the -NH 2 -modified substrates demonstrated a well-defined cytoskeleton at both day 1 and day 3 (Figures 2(a) and 2(e)), comparable to the cells cultured on the -CH 3 -modified substrate (Figures 3(a) and 3(e)). The morphology of the BMSCs cultured on the -OH-modified surface was similar to that observed on the -NH 2 -modified substrate during the test period (Figures 2(l) and 2(p)).
Cells cultured in contact with the -CH 3 -modified surface demonstrated evidence of a spindled morphology at day 1 (Figures 3(a) and 3(d)), whereas there were increases in cell coverage, stress fiber, and focal contact formation by day 3 on this surface (Figures 3(e) and 3(h)). The morphology of the cells on the -CH 3 -modified surface was less well spread than these on -NH 2 -and -OH-modified surfaces (Figure 3(h)).
The individual cells cultured in contact with the -COOHsubstrate displayed a smaller spread area at day 1 (Figures 3(i) and 3(l)) than those on -NH 2 -and -OH-modified surfaces, with strong evidence of focal contact formation at both days 1 and 3 (Figures 3(l) and 3(p)).

Chemical Functional Groups Affected the Osteogenic
Differentiation of BMSCs. The osteogenic differentiation of BMSCs is manifested by numerous specific markers. The most commonly used markers are the synthesis of calciumbinding glycoproteins, osteocalcin and osteopontin; collagen I; the activity of ALP; and the bone tissue mineralization [34]. Because Cbf -1, a potent transcription factor involved in the regulation of several osteoinductive genes, is essential for the differentiation of osteoblasts from mesenchymal precursors and for bone formation [35,36], we used real-time RT-PCR to examine the expression of Cbf -1 and ALP in BMSCs cultured on the different chemical substrates (Figure 4). All results were normalized to the appropriate GAPDH values of the test samples and are represented as an increase/decrease in cellular expression.
Significant upregulation of both Cbf -1 and ALP was observed for the cells cultured on the -NH 2 -modified substrate relative to those cultured on the -CH 3 -modified surface. Increased Cbf -1 expression was observed on -NH 2modified surface on days 7, 10, and 14. In addition to the effect of the -NH 2 -modified surface, the -OH-and -COOHmodified surfaces caused increases in the expression levels of Cbf -1 and ALP on certain days. For example, increases in both markers were observed in the cells cultured on the -OH-modified surface on day 10 and on the -COOH-modified surfaces on days 10 and 14. However, these trends were not as obvious as those observed for the -NH 2 -modified surface (Figure 4(a)). The expression of ALP was increased for cells cultured on the -NH 2 -modified surface on days 7, 10, and 14.
The ALP expression was also increased for the cells cultured on the -OH-modified surface on day 10 and on the -COOHmodified surface on days 10 and 14 (Figure 4(b)). chemical modifications of the substrate can affect the adhesion, spreading, and osteogenic differentiation of BMSCs. The -NH 2 -modified surfaces promoted osteogenesis of the BMSCs in the absence of biological stimuli as previous described [24][25][26]. To investigate the downstream signaling pathways in BMSCs that may be involved in the regulation of osteogenesis during culture on the -NH 2 -and -CH 3modified substrates (the modifications with the strongest and weakest ability to promote the osteogenic differentiation of BMSCs, resp.), we focused on the key proteins of the ERK1/2 signaling pathway, as suggested by earlier studies implicating these pathways [37]. Cells were cultured on the -NH 2 -and -CH 3 -modified surfaces for up to 14 days. The analysis of the mRNA and lysates of cells revealed that chemical groups influenced the ERK1/2 pathway. Integrins are a family of transmembrane proteins that bind the cell to extracellular matrix proteins, attach to the cytoskeleton at focal adhesion sites, and induce intracellular signals [38][39][40][41]. The results of the evaluation of the mRNA expression of integrins 1, 5, and in cells cultured on -NH 2 -and -CH 3 -modified substrates are shown in Figure 5. The expression level of integrin 1 in cells cultured on the -NH 2 -modified surface was significantly greater than that for cells cultured on the -CH 3 -modified surface at 15 min ( Figure 5(a)); integrin 5 on the -NH 2 -modified surface was higher than on the -CH 3 -modified surface at 5 min and significantly higher at 15 min and 30 min (Figure 5(b)); the expression of integrin was significantly higher on the -NH 2 -modified surface than on the -CH 3 -modified surface at 15 min and 30 min ( Figure 5(c)). The peak of the expression of integrin 1 occurred at 15 min ( Figure 5(a)), while the peaks of integrin 5 and integrin expression occurred at 30 min ( Figure 5(b)).

Activation of the ERK1/2 Pathway in BMSCs
Next, we examined integrin activities at 5, 15, 30, and 45 min. As shown in Figure 6, the peak level of integrin 1 activity in cells cultured on substrates modified with the -NH 2 group was observed at 15 min and was followed by a decline in the levels between 15 and 45 min; the peak level of integrin 1 activity in cells cultured on substrates modified with the -CH 3 group was observed at 30 min (Figure 6(a)); the activity of integrin 5 on the -NH 2 -modified surface was significantly higher than on the -CH 3 -modified surface at 15 min and 30 min (Figure 6(b)); the activity of integrin was significantly higher on the -NH 2 -modified surface than on the -CH 3 -modified surface at 30 min (Figure 6(c)). The and subunits of the integrins do not have any intrinsic enzymatic activity. However, these proteins interact with cytoplasmic proteins that nucleate the formation of large protein complexes. These complexes in turn activate intracellular signaling cascades, leading to alterations in gene expression. Focal adhesion kinase (FAK) is one of the numerous signaling molecules present in focal adhesion sites [42]. We examined FAK activity at 5, 15, 30, and 45 min. As shown in Figure 7, the peak level of FAK activity in cells cultured on substrates modified with the -NH 2 group was observed at 15 min and was followed by a decline in the levels between 15 and 45 min. The peak level of FAK activity in cells cultured on substrates modified with the -CH 3 group 8 BioMed Research International was observed at 45 min, but almost no detectable activation of FAK could be measured at 5 and 15 mins. Taken together, these results indicated that the -NH 2 -modified substrate promoted the activation of FAK.
To investigate whether the osteogenic differentiation of BMSCs on the -NH 2 -modified surface paralleled the activation of ERK1/2, we examined ERK1/2 activity on days 1, 3, 5, 7, 10, and 14. As shown in Figure 8, the chemically modified substrates did not affect ERK1/2 phosphorylation until 7 days. A high level of ERK activity was observed on day 10, which was followed by a decline. The time period during which phosphorylation of ERK1/2 was the highest was associated with high expression of Cbf -1 and ALP. Meanwhile, the activation of ERK on the -CH 3 -modified substrate at each BioMed Research International Cells were cultured on the -NH 2 -and -CH 3 -modified substrates, and lysates were prepared at the indicated times. Lysates were subjected to immunoblot analysis using integrin antibodies. An aliquot of each lysate was subjected to a kinase detection assay using SDS-PAGE. GAPDH was used as a control. , and lysates were prepared at the indicated times. Lysates were subjected to immunoblot analysis using phosphospecific and nonactivated FAK antibodies. An aliquot of each lysate was subjected to a kinase detection assay using SDS-PAGE. GAPDH was used as a control.
time point was less than on the -NH 2 -modified substrate. Taken together, these results indicated that the osteogenic differentiation of BMSCs on the -NH 2 -modified surface correlated well with the activation of ERK1/2. To clarify whether ERK activation is necessary for osteogenic differentiation of the BMSCs cultured on the -NH 2modified substrate, we blocked the ERK1/2 signaling pathway using different amounts of PD98059, a potent and selective inhibitor of MEK1 (a kinase upstream of ERK1/2) [43]. To determine whether PD98059 suppressed the increases in ERK1/2 activity observed during culture of the cells on the -NH 2 -and -CH 3 -modified substrates, we measured the levels of ERK1/2 and p-ERK1/2 in the presence of PD98059 at concentrations of 0, 10, 25, and 50 M on day 10. As shown in Figures 9(a) and 9(b), treatment of BMSCs with PD98059 (0, 10, 25, and 50 M) resulted in concentration-dependent inhibition of ERK activation on day 10.
Osteogenic differentiation of BMSCs was assessed on the basis of Cbf -1 and ALP expression using real-time RT-PCR. As shown in Figures 10(a) and 10(b), treatment of BMSCs with PD98059 at day 10 resulted in concentration-dependent inhibition of Cbf -1 and ALP activities. However, on day 10, the Cbf -1 and ALP activities in each set of samples cultured on the -NH 2 -modified substrate in the presence of the same concentrations of PD98059 were still significantly higher than that in the corresponding groups cultured on the -CH 3 -modified substrate ( < 0.05). On the basis of the western blot analysis, the inhibitory effect of PD98059 on the phosphorylation of ERK1/2 paralleled the downregulation of Cbf -1 and ALP activities in a similar dose-dependent manner (Figures 9 and 10).

Discussion
Cell-material interactions are critical for many biomedical and biotechnological applications. What makes these cells different when they are adhered to different materials? To answer that question, we must first fully understand how the cells translate the cues they receive from the extracellular microenvironment to the intracellular signaling, which eventually leads to activation of transcription factors and altered gene expression [44].
A series of reviews recently summarized numerous studies that clearly show that the biochemical characteristics of the substrate, as well as its rigidity and spatial organization, are recognized by cells through differential signaling from  Figure 8: The -NH 2 -modified substrate induces the activation of ERK1/2 in BMSCs. Cells were cultured on the -NH 2 -modified substrate (a) and the -CH 3 -modified substrate (b), and lysates were prepared at the indicated times. Lysates were subjected to immunoblot analysis using phosphospecific and nonactivated ERK1/2 kinase antibodies. An aliquot of each lysate was subjected to a kinase detection assay using SDS-PAGE. GAPDH was used as a control. integrin-based molecular complexes [45][46][47]. Cells interact with solid substrates such as extracellular matrix (ECM) proteins through integrins, which, upon phosphorylation, initiate signaling cascades, many of which alter the activity of transcription factors that control gene expression. Integrin-ECM binding induces a wide range of intracellular signaling events such as the activation of Ras, FAK, Src, Rac/ Rho/CDC42 GTPases, MAPK, PKC, and PI3K in different cell systems [48][49][50]. The pathway from integrin binding to FAK to MAPK is a particularly well-studied pathway, but there are a large number of pathways by which microenvironmental cues are transmitted into the cell and eventually alter gene expression. The chemical characteristics of the substrate surface have been shown to affect the biological behavior in different cell systems in general, and the osteogenic differentiation potential of BMSCs in vitro in particular, in a number of studies [24][25][26]. However, the molecular mechanism underlying the substrate-dependent osteogenic commitment of these cells remains unclear. Elucidation of these mechanisms could improve the ability to promote the osteogenic differentiation of BMSCs for the purposes of cell therapy and tissue engineering. As one of the major components of the MAPK family, ERK1/2 has been associated with cellular survival, proliferation, and differentiation and especially the osteoblastic differentiation of BMSCs [51][52][53]. Thus, the objective of this study was to investigate whether activation of the ERK1/2 pathway was involved in the differentiation of BMSCs cultured on differentially modified surfaces into an osteogenic lineage.
It is important to recognize that the same intracellular signaling molecules that are activated by integrins are also activated by the stimulation of various cytokines and other cues [54]. Cell behavior is often dependent on synergy between growth factor/cytokine-mediated and integrinmediated signaling pathways [55][56][57]. The well-defined regulatory effect of dexamethasone (Dex), present in osteogenic medium (OM), on cellular growth and differentiation is usually assumed to contribute to the mitogenesis and osteogenic differentiation of BMSCs. This mediator has previously been reported to activate the ERK1/2 pathway [19]. To eliminate the osteogenic differentiation effect of Dex, we cultured BMSCs in the absence of Dex and evaluated the osteogenesis of BMSCs. The biological behavior of BMSCs has been previously described [58].
SAMs are currently the best available class of model organic surfaces with which to control the interfacial structure and properties. Previous studies have demonstrated that terminal chemical groups control cell behavior including adhesion, migration, and differentiation in various cell types including human osteoblast-like cells [59], human fibroblasts [60], mesenchymal stem cells [24][25][26], and neural stem cells [27]. We have developed molecules with specific functional groups by a novel SAM technology and confirmed the structures produced by these functional groups by AFM and XPS. XPS spectra demonstrated that the four types of functional groups successfully self-assembled, and AFM images displayed the same morphology of these functional group surfaces [30,31]. Because the different chemical functional groups had same structures and same elements, the only reason of the different behaviors of cells seeded on different chemical functional groups was their chemistry.
It is well established that cell adhesion and morphology influence subsequent activities [61]. In accordance with previous publications, our study confirms that cell viability depends on surface properties and varies with individual functional groups rather than general surface properties such as hydrophilicity [62][63][64]. For example, the -OH-modified substrate was the most hydrophilic surface (Figure 1), but similar cell contact areas were observed on this substrate and on the -NH 2 -modified surface (Figure 2). Figure 2 showed that the cells plated on -NH 2 -and -OH-modified substrates were well spread and homogeneous with the actin organized into stress fibers and demonstrated long microtubules and distinct focal adhesions. On the other hand, the cell viability results demonstrated that all of the substrates maintained levels of viable cell adhesion throughout the test periods in all of the mediums tested. However, BMSCs cultured on the -NH 2 substrate maintained increased viable cell adhesion on day 7.
Consistent with the studies of JM Curran, which indicated that the surface properties of the culture substrate influence the phenotype of BMSCs [24,25], our gene expression analysis using real-time RT-PCR suggests that BMSCs cultured on -NH 2 -modified substrates in the absence of osteoinductive medium upregulate the expression of Cbf -1 and ALP during the first 14 days of culture ( Figure 4). The function of Cbf -1 dominates the osteoblast differentiation, which is characterized by regulation of the expression of major osteogenesis-related genes, including COLI, ALP, and OCN. Cbf -1 also maintains the functions of early-stage differentiation of the osteoblasts [21]. As one of the important markers of osteogenesis assessment, ALP was shown to regulate organic or inorganic phosphate metabolism via the hydrolysis of phosphate esters and to function as a plasma membrane transporter for inorganic phosphates [18]. Based on the results of cell viability, cell adhesion, and the mRNA expression, we affirmed that -NH 2 chemical functional group could contribute to the osteogenic differentiation of BMSCs. So we selected the -NH 2 -modified substrate to examine the role of ERK1/2 signal pathway in osteogenic differentiation of BMSCs as a function of the chemical groups used to modify the substrate.
Integrin-mediated multiprotein adhesion complexes link the extracellular matrix to the actin cytoskeleton. Molecular analyses of these adhesion sites indicate that the integrin adhesome consists of ∼160 distinct components. Examination of the molecular interactions that takes place among the different constituents of the adhesome points to an extraordinary connectivity. The biological activities of the adhesome components are diverse and include several actin regulatory elements that affect the organization of the attached cytoskeleton. Many of the adaptor proteins directly or indirectly link actin to integrins and a wide range of signaling molecules, such as kinases, phosphatases, and G proteins and their regulators. It seems likely that the tight association between the structural and signaling elements of the adhesome provides its unique properties as a sensitive environmental sensing system to the adhesion machinery [65]. Therefore, optimizing the osteogenic differentiation potential may not simply involve providing a uniform coating to engage a distinct receptor, but instead the identification of the surface properties that control the presentation of integrin-specific epitopes associated with specific chemical groups may be required [37].
Based on our findings, 1, 5, and integrins-FAK appear to mediate the adhesion of BMSCs mediated by specific chemical groups. These associations may perhaps even be required for osteogenic differentiation (Figures 5,  6, and 7). Figures 5 and 6 showed the results of the expression of integrins 1, 5, and in cells cultured on -NH 2 -and -CH 3 -modified substrates. The expression level of integrins 1, 5, and in cells cultured on the -NH 2 -modified surface was significantly greater than that for cells cultured on the -CH 3 -modified surface at different times; Figure 7 showed the FAK activity in cells which were cultured on substrates modified with the -NH 2 group and -CH 3 group. The peak level of FAK activity in cells cultured on substrates modified with the -NH 2 group was observed at 15 min and was followed by a decline in the levels between 15 and 45 min, but the peak level of FAK activity in cells cultured on substrates modified with the -CH 3 group was observed at 45 min, but almost no detectable activation of FAK could be measured at 5 and 15 mins. Taken together, these results indicated that the -NH 2modified substrate promoted the activation of integrin and FAK.
Our results show that the activation of ERK1/2 in BMSCs cultured on the -NH 2 -modified substrate was greater than that on the -CH 3 -modified substrate during the period between 7 and 14 days (Figure 8). The chemically modified substrates did not affect ERK1/2 phosphorylation until 7 days. The expression level of the phosphorylation of ERK1/2 in cells cultured on the -NH 2 -modified surface was significantly greater than that on the -CH 3 -modified surface at 10 and 14 days. In this study, we found that ERK1/2 was activated in a time-dependent manner consistent with the expression of Cbf -1 and ALP in BMSCs cultured on the -NH 2 -modified substrate. Hence, it can be concluded that ERK1/2 activation parallels the osteogenic differentiation of BMSCs on the -NH 2 -modified substrate.
To further investigate the role of the ERK1/2 pathway in BMSC osteogenic differentiation on the -NH 2 -modified substrate, the activation of the ERK1/2 pathway was blocked with PD98059, and the corresponding osteogenic differentiation of BMSCs was examined on day 10. PD98059 is a specific inhibitor of MEK1, a kinase upstream of ERK1/2 and is able to bind to the inactive forms of MEK1 and thereby block the ERK1/2 signaling pathway [66]. PD98059 inhibited expression of the phosphorylation of ERK1/2 in a dose-dependent manner in both substrates. At a concentration 50 M of PD98059, expression of the phosphorylation of ERK1/2 was almost completely inhibited ( Figure 9). Figure 10 showed that PD98059 inhibited expressions of Cbf -1 and ALP in a dosedependent manner also. But at the concentration 50 M of PD98059, the expressions of Cbf -1 and ALP were only partly inhibited ( Figure 10).
Thus, these results demonstrated that blocking the activation of the ERK1/2 pathway leads to the partial inhibition of the osteogenic differentiation of BMSCs on -NH 2 -modified substrates and further suggests that activation of the ERK1/2 pathway is primarily responsible for the osteogenesis of BMSCs on -NH 2 -modified substrates.

Conclusion
The data presented here suggest that -NH 2 chemical functional group could contribute to the osteogenic differentiation of BMSCs, and the influences depend on chemical functional group rather than general surface properties such as hydrophilicity.
The BMSCs which cultured on -NH 2 -modified substrates alter their expression of integrins and FAK and activate the downstream ERK1/2 signaling pathway that is required for osteogenic differentiation. It is believed that the balance of molecular events responding to various extracellular stimuli determines the process of cell type-specific differentiation of multipotent mesenchymal stem cells. Therefore, other signaling pathways, such as c-Jun N-terminal kinase, p38 and PI3K, as well as cross-talk between these pathways, need to be further characterized to explore the detailed mechanisms responsible for the osteogenic differentiation of BMSCs cultured on various chemical groups.