The environment is an important driver for organismic natural selection. Environmental changes or climatic fluctuations can make organisms evolve rapidly into different morphologic or taxonomic groups or create new functions specialized in different individual living environments [1]. Organismic evolution is the process by which an organism must repetitiously overcome new conditions and create a new set of unique metabolic reactions to a particular environmental stress. This adaptive response to changing environmental conditions results in an acceleration of evolutionary rate of the lineage and the functional evolution of specific stress-response proteins, which favors organismic fitness to a new environment. All evolution of phenotypes results from minor increases in mutation rates of genes related to a particular stress [2]–[4]. The field of modern molecular evolution provides powerful tools for us to study the relationships between the functional changes of proteins and the rates of nucleotide and amino acid substitution [5], [6]. Comparisons of the ratio of non-synonymous (Ka)/synonymous (Ks) substitutions (ω = Ka/Ks) have become a useful means for quantifying the impact of natural selection on molecular evolution [7], [8]. Synonymous mutations are generally neutral in the course of evolution and do not result in changes to the amino acids in a protein, while non-synonymous mutations can occur under strong selective pressure and result in the altering of the amino acids in a protein. A value of ω = 1 denotes a neutral mutation, ω less than 1 purifying selection which describes selection against new variants, while ω greater than 1 denotes positive selection (adaptive molecular evolution) in that non-synonymous mutations offer fitness advantages to the protein [9]. A purifying selection may aid in the detection of regions or residues of functional importance. However, much interest in evolution focuses on positive selection because it is associated with adaptation and the evolution of new forms or functions.

Pikas are small non-hibernating, diurnal lagomorphs (rabbits and relatives; order Lagomorpha) that belong to the family Ochotonidae. Pikas are endemic to the modern Holarctic Region [10], [11]. There are approximately 26 species throughout the world; most are restricted to Asia with only three pika species that presently live outside of Asia. In North America, the only two species (the American pika, O. princes, and the collared pika, O. collaris) are discontinuously distributed across the mountainous areas of western North America, from the southern Sierra Nevada and Rocky Mountains to central British Columbia. In Europe, the range of the steppe pika (O. pusilla) extends west to the Ural Mountains. In Asia, pikas are found throughout central Asia, in the Himalayan massif and associated ranges, and across Western Siberia to Sakhalin Island and onto Hokkaido Island, Japan. Among the species in Asia, approximately 18 species are concentrated in the Qinghai-Tibet Plateau and adjacent areas, accounting for greater than 70% of all pikas in the world [12]–[15]. Obviously, most of the species are confined to those regions with either high altitudes or high latitudes and live in cold climates. The species diversification of pikas on Qinghai-Tibet Plateau implies that pikas are particularly fitted for survival under the environment or climate of the Qinghai-Tibet Plateau. Hypoxia and cold are the two most remarkable climate characteristics of the Qinghai-Tibet Plateau, known to be the highest plateau in the world. Thermoregulation is very important for animals' survival in cold environment. Under these extreme conditions, animals show apparent alteration in energy expenditure in order to meet the varying metabolic requirements imposed by environmental stresses [16]. During evolution, plateau pikas (O. curzoniae), the keystone species in the Qinghai-Tibet Plateau ecosystem [17], have become highly hypoxia- and low temperature-tolerant mammals with markedly high resting metabolic rates (RMR), non-shivering thermogenesis (NST), and a high ratio of oxygen utilization to cope with the cold and hypoxic plateau environment [18]–[20]. Furthermore, plateau pikas also show marked seasonal changes in thermogenic capacities with enhanced nonshivering thermogenesis (NST), cytochrome c oxidase activity and increased production of mitochondrial uncoupling protein 1 (UCP1) in brown adipose tissues (BAT) to deal with the cold in winter [21].

Leptin, the product of the ob gene, is a 16-kDa cytokine-like hormone primarily secreted by adipose tissue that acts on hypothalamic centers and regulates the energy balance by determining changes in both food intake and energy expenditure [22]–[24]. Physiologically, plasma leptin concentrations reflect adipose tissue deposits [25] and appear to be negatively regulated by fasting [22], protein kinase C [26], androgens [27] and β-adrenergic agonists [28], whereas feeding [22], [29], oestrogens [30], glucocorticoids [30], [31] and insulin act by stimulating both ob gene expression and leptin secretion [32]. Leptin is considered as an important regulator of adaptive thermogenesis in response to environmental temperature or diet. Leptin administration can increase an animal's body temperature, basal metabolic rates (BMR), and nonshivering thermogenesis (NST) which enhance its ability to tolerate cold stress [33]. Leptin increases energy expenditure by increasing uncoupling protein 2 (UCP2) expression in white adipose tissue (WAT), increasing sympathetic activation of UCP1 gene expression in BAT, and upregulating LPL gene expression in brown adipose tissue (BAT) [34], [35]. Cold exposure can reduce leptin expression by directly acting on adipocytes or indirectly acting via the sympathetic nervous system [36], [37], while hypoxia can increase leptin expression [38]. Obviously, regulation of leptin plays an important role in meeting the fluctuating requirements of energy expenditure and energy intake under different physiological states. Due to the extraordinary environment in which pikas live, and their apparent changes in energy metabolism under cold conditions, we hypothesized that leptin, acting as a cold stress-response protein, plays an important role in the pikas' ecological adaptation to the harsh plateau environment. However, the sequence characteristics of leptin specialized in extreme environmental stress, which is very important for us to understand animals' ecological adaptation mechanism, remains unknown. In our previous study, only leptin from the plateau pika was cloned and was identified as exhibiting divergent mRNA expression at different altitudes [39]. Together, these data led us to ask: is pika leptin sensitive to the cold and hypoxic plateau environment; does pika leptin itself functionally evolve under this environmental stress?

To examine the extent of leptin variation within the Ochotona family, we carried out both DNA and amino acid sequence analyses in molecular evolution and compared modeled spatial structures. To identify the environmental factor for driving evolution of pika leptin, we compared leptin sequences of pikas from different altitudes. Therefore, our study first explored the relationships between environmental stresses on the Qinghai-Tibet Plateau and molecular evolution of stress-response genes.

A 646-bp fragment in pikas and a 565-bp fragment from both Lepus oiostolus and Oryctolagus cuniculus, which contained the complete coding region, were cloned, respectively. These sequences were submitted to GenBank and were assigned the following accession numbers: DQ983189 (Ochotona curzoniae); EF091861 (Ochtona nubrica); EF091863 (Ochotona cansus cansus 1); EF091864 (Ochotona cansus cansus 2); EF091862 (Ochotona annecten); EF091860 (Ochotona daurica bedfordi); DQ983190 (Lepus oiostolus), and DQ983191 (Oryctolagus cuniculus). The deduced amino acid sequences were composed of 167 amino acids and encoded an apparent signal peptide sequence of 21 amino acids with the signal cleavage site between Ala-21 and Val-22. Thus, the mature secreted protein had a predicted molecular weight of 16.086 kDa and a pI of 6.3. To reveal the evolutionary divergence in leptin sequences among lineages, we assembled and analyzed 20 sequences representing samples from different lineages in vertebrates. The result of the multiple alignments for amino acid sequences is shown in Figure 1.

In the present study, we have compared the entire coding sequence of leptin from different lineages of representative mammals in order to help us identify the variation of functional sites and to understand the mechanism of functional evolution of pika leptin. The phylogenetic tree of leptin yielded similar topology to that of the mitochondrial cytb gene, and thus, was consistent with traditional morphological assortment. However, comparing divergent distance with the inner branch length in the phylogenetic tree of leptin, we found the pika branch (shown as branch B in Figure 2A) to be significantly longer than any of the other lineage branches, implying that the variation of leptin sequence in the pika lineage is great. Next, to identify the evolutionary rate of leptin to be a result of neutral evolution or natural selection, the relative rate test was performed. The results showed that pika branch was significantly different (P<0.05) from all other lineages in nonsynonymous (Ka) and amino acid variation, but not in synonymous (Ks). Interestingly, in comparison of evolutionary rate between pika and rabbit lineages we found that the differences of the evolutionary rate between these two lineages were significant, but the differences between rabbit and all other lineages were not significant. In traditional taxonomy, pika and rabbit have a closer kinship than any other compared lineages and belong to the same family (Lagomorpha). However, for leptin evolution, the taxonomic relationship between these two lineages is changed and the divergence between pika and rabbit lineages appears to be significantly greater than that between rabbit and all other lineages. That is to say, the evolutionary rate of leptin between rabbit and other selected lineages was stable and conformed to the molecular clock model. While for pika, after diverging with rabbit, the evolutionary rate of pika leptin may be accelerated.

To determine the nature of variation sites occurring in pika leptin, a set of evolutionary analysis was performed. A comparison of the one-ratio vs. the two-ratio in branch-specific models revealed the ω ratio along the pika lineage was significantly different from all other lineages. In site-specific models, both the M8 (beta & ω) and M3 (discrete) models demonstrated unconsentaneous ω ratios among sites, yielding a ω ratio of 1.22961 and 1.25333, respectively, and predicted one common site under positive selection, 92F. Under branch-site models, model B provided a ω ratio of 1.31007 and identified the following sites to be under positive selection: 2S, 4W, 7R, 28H, 29A, 36I, 44A, 59V, 60L, 62K, 63H, 92A, 94Q, 95G, 98P, 103D, 106S, 108N, 111E, 113I. The ancestral sequence reconstructed by the models of Goldman and Yang [43] and of Yang, Kumar, and Nei [44] suggested the following amino acid changes along the pika branch: 2SP, 4WQ, 7RQ, 28HQ, 29AS, 36IV, 44AG, 59VA, 60LV, 62KQ, 63HQ, 94QK, 95GS, 98PL, 103DG, 108NE, 113IV, and 134HQ. Obviously, sites under positive selection, inferred by Bayes prediction, were highly consistent with those from the reconstruction of the ancestral sequence. Therefore, our evolutionary analysis confirmed the previous hypothesis that adaptive evolution occurred in pika leptin.

The establishment of new or modified function of a protein under specific stress is derived from the adaptive evolution in this protein. We speculate the possible effect of positive selection sites on the functional evolution of pika leptin interpreted from the analysis of literature on the evolutionary, functional-structural, and biochemical information concerning the leptin protein. Previous investigations indicated that segment 85-119 in leptin protein was of special functional significance and underlied the functional differences between human and other non-hominoid leptins [45]–[48]; the key binding sites between leptin and its receptor were 9D, 12T, 15K, 16T, 82N, 85D, and 86L; the important sites for leptin signaling activation were 20R 29S, 30V, 31S, 34Q, 35R 41F, 75Q 115E, 117S, 120S, 121T, 122E 138Q, 139Q, and 142V. Mutations at these sites showed strongly decreased binding affinity or signaling activation for its receptor [49]–[52]. In pika leptin, positive selection (PS) mutations occur at the following 20 sites: P2S, Q4W, Q7R, Q28H, S29A, V36I, G44A, A59V, V60L, Q62K, Q63H, F92A, K94Q, S95G, L98P, G103D, T106S, D108N, G111E, and V113I. Nine PS sites locate within the functionally significant segment 85-119 of leptin: F92A, K94Q, S95G, L98P, G103D, T106S, D108N, G111E, and V113I, implying an important conclusion that the function of pika leptin may be divergent from the other leptins in all compared lineages. Additionally, the G44A mutation results in a unique motif appearing only in the pika lineage, the ATP synthase α and β subunit signature site, composing of the functional sites (nucleotide-binding site for ATP and ADP in the α subunit and catalytic activity in the β subunit) of the ATP synthase complex that take part in energy transduction in living cells [53]. This motif predicted in pika leptin appears to be consistent with its general role in energy regulation and also seems to be associated with the requirement of varying energy expenditure under extreme environmental stress. The exact experimental evidence for the existence and function of the ATP synthase α and β subunit signature site in pika leptin needs to be further studied. Substitution by Ala at site 29 occurs in a key site for signaling, changing from a polar (Ser) to a non-polar (Ala) residue. What effect this substitution has on receptor signal activation requires further study. Taken together, these analyses mentioned above led us to suppose that for pika leptin, nine sites under positive selection were within the functionally significant segment 85-119 and one unique motif (ATP synthase α and β subunit signature site) appeared only in the pika lineage, indicating the functional variation of pika leptin.

It is known that cold and hypoxia are the two most remarkable climatic characteristics of the Qinghai-Tibet Plateau. To identify the environmental factor driving the adaptive functional variation of pika leptin, we collected pikas from different altitudes, five species from the Qinghai-Tibet Plateau (average altitude >3000 m) and the other from the Inner Mongolia steppe (altitude of 1300 m). We also collected plateau pikas from three altitudes (3200 m, 3900 m and 4790 m). Because most of the literature concerns cold survival environment in pika species, it is difficult to find a pika species living in a warmer climate. We therefore used trial animals, rabbit and rat as controls living under the environment of warmer temperature and lower altitude. Mean rates of synonymous substitution (Ks), non-synonymous substitution (Ka) and amino acid substitution (Aa) relative to the outgroup-human were investigated and stepwise multiple regression was used to determine how mean January actual temperature (Tjanu, °C) and altitude (Al, in m) affected these substitution rates. The results of stepwise multiple regression showed that pikas with relatively lower mean January actual temperature reached relatively higher substitution rates of Ka and Aa, while altitude was not included in the model, and thus, did not significantly affect the substitution rates of Ks, Ka and Aa. Multiple alignment of sequences of plateau pika leptin from different altitudes also showed that there was only one synonymous mutation in the nucleotide sequences, but no changes in amino acid sequences. The relation between altitude and barometric pressure or inspired oxygen pressure is negatively correlated [54], [55]. Namely, as the altitude increases, the oxygen content of the air decreases, and thus hypoxia aggravates. Therefore, hypoxia differs at the different altitudes. In addition, the evolutionary history of the pika family is also indicative of the cold adaptation in the Ochotona family. Fossil records show that this family had an Asiatic origin with the earliest emergence of pikas hypothesized to be during the late Miocene in central Asia and migrated to North America through the Bering Strait during the Mid-Pliocene [56], [57]. Time of divergence within the Ochotona family was consistent within the historical episodes of both geologic and climatic changes. The frequent uplift of the Qinghai-Tibet Plateau from 3.4–1.6 mya [58], [59] resulted in strong environmental changes (climatic fluctuations), potentially producing strong selective pressures that resulted in varied clades of Ochotona of the Qinghai-Tibet Plateau. With gradual global warming, there were apparent distribution trends toward increased elevations and contracted ranges in pika populations, eventually resulting in the global extinction of some species [13]. Obviously, pikas have been living in a harshly cold environment both in ancient and present time. Together, all these data indicate that cold, and probably not hypoxia, may be the primary environmental factor for driving adaptive evolution of pika leptin. Furthermore, similar changes of pika leptin in the two regions (Qinghai-Tibet Plateau and Inner Mongolia steppe) lead us to suppose that the adaptive evolutionary variation of leptin may have occurred within the entire Ochotona family throughout the world due to their common evolutionary history of Asiatic origin and their similar cold environment for survival in all ages.

Adaptive thermogenesis is the main way of heat production for small mammals in response to cold environmental stress and is produced mainly by means of nonshivering thermogenesis (NST) associated with an increase in BAT weights, mitochondria protein concentrations, and uncoupling protein 1 (UCP1) mRNA expression [60], [61]. Leptin administration increases an animal's body temperature, basal metabolic rates (BMR), NST, and UCP1 mRNA expression in BAT [33]–[35]. Compared with other small mammals that live in warmer environments, pikas generally show markedly high levels of resting metabolic rates (RMR) and non-shivering thermogenesis (NST) to cope with the harsh plateau environment [20]. Therefore, we speculate that adaptive functional evolution in pika leptin may play an important role in contributing to fitness enhancement for pikas' survival in a stressful environment. Leptin also regulates lipid and glucose metabolism and insulin action [62], [63]. Leptin stimulates fatty acid oxidation and glucose uptake in skeletal muscle [64]–[67], and inhibits glucose output, lipogenesis in liver [68], [69], and insulin secretion [70], which increase insulin sensitivity. Leptin can regulate osteoblastic osteocalcin production via insulin response. Bromage speculates that natural selection may be acting on insulin production, achieving this by altering leptin in such a way as to enhance hypothalamic sensitivity to leptin and the resultant osteocalcin production (Bromage, pers comm.). These facts mentioned above make it plausible that adaptive evolution in pika leptin may be targets of study in order to help us to clarify the adaptive mechanism for small mammals living under extremely stressful environments, and to identify potential therapeutic strategies for human's disease associated with metabolic disorders.

In summary, our study confirmed the previous hypothesis that leptin is a cold stress-response protein and that cold probably is the primary environmental factor for driving the adaptive functional evolution of leptin within the native cold-adapted Ochotona family, contributing to fitness enhancement for the pikas' survival in a stressful environment. We also put forward an important hypothesis: adaptive functional evolution of pika leptin may be a common characteristic of the entire Ochotona family throughout the world.