Asian pika populations track local glaciation events through the Pleistocene

Asian pika populations track local glaciation events through the Pleistocene | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Asian pika populations track local glaciation events through the Pleistocene Harsha Kumar , Nagarjun Vijay , R Nandini doi: https://doi.org/10.1101/2025.06.25.661461 Harsha Kumar 1 Department of Biology, Indian Institute of Science Education and Research Tirupati , Andhra Pradesh 517507, India 3 Tata Institute of Fundamental Research Hyderabad , Telangana 500046, India Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: nandini{at}labs.iisertirupati.ac.in harshakkumar{at}students.iisertirupati.ac.in Nagarjun Vijay 2 Department of Biology, Indian Institute of Science Education and Research Bhopal , Madhya Pradesh 462066, India Find this author on Google Scholar Find this author on PubMed Search for this author on this site R Nandini Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: nandini{at}labs.iisertirupati.ac.in harshakkumar{at}students.iisertirupati.ac.in Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract Background The Pleistocene glaciation cycles (2.6 mya - 11 kya) were major climatic events that shaped diversity and community assemblages on Earth. Cold-adapted, high-elevation specialists are expected to have responded negatively to warm interglacial periods, with populations contracting and being pushed up the elevational gradient, leading to isolated populations on mountain tops. During the cool glacial maximums, they are expected to have been distributed at both low and high elevations, with populations expanding and facilitating gene flow across mountains. Pikas ( Ochotonidae ) are poorly studied high-elevation specialist Lagomorphs being extirpated at alarming rates due to climate change. Insights into their historical demography and current effective population size are crucial in determining their current population status to inform conservation planning. Results We use de novo assembly to construct partial genomes of three Asian pika species ( Ochotona ladacensis, O. nubrica, and O. macrotis ) that differ in their life history, lifestyle, and social behaviors. We use a combination of species distribution modeling (SDM) projected back in time, in addition to inferring historical demographic patterns using a pairwise sequentially Markovian coalescent (PSMC) approach. Our SDMs predicted the largest bio-climatic niches for all three species during the last glacial maximum (LGM ∼20 kya), and the smallest niches at the last interglacial (LIG∼120 kya), agreeing well with the elevational shift hypothesis. Our PSMC models for the populations in the Changthang Biotic province (CBP), revealed largely synchronous oscillations of populations sampled in geological timescales but displayed contrasting temporal patterns of population spikes when compared with SDMs. The largest effective population sizes of all species examined on the CBP were inferred to pre-date both the LGM and LIG and were placed around the MIS-6 (∼ 191 kya). After the MIS-6 local glaciation event, populations were inferred to have been through a deep bottleneck, correlating well with local glaciation patterns on this plateau. Conclusion Our study illustrates the power of using complementary approaches to infer historical demography in mountainous landscapes with complex glaciation histories. While SDMs appropriately describe species responses to historical climate across the geographic range of species, PSMCs are more informative about local population dynamics when subject to large degrees of population differentiation. All three species of pikas responded similarly to historical glaciation events on the CBP that differ from the rest of the world and other parts of the Himalayas and likely have very low effective population sizes at current timescales. This is one of the first few studies to examine the population demography of small mammals on the Indo-Tibetan Plateau and has the potential to inform species and population-specific conservation planning in the Himalayas. It is very likely that small mammals such as pikas, marmots and voles on the CBP have low effective population sizes due to historical inbreeding during the LIG and the lack of a glaciation rescue during the LGM. This highlights the need for conservation measures for small mammals on the landscape. Background Climate change-associated warming and erratic weather patterns are expected to cause rapid declines in biodiversity worldwide [ 1 , 2 ]. Alpine specialists on mountain tops are at differentially higher risk of extinction, given their limited potential to shift ranges [ 3 – 5 ]. Since most species living in remote landscapes have not been assessed for extinction risk, identifying traits that make animals more susceptible to climate change is useful for species management [ 6 – 9 ]. Species at risk today have survived large-scale climatic events during the Pleistocene glacial cycles, and gaining insights into their historical responses might help us predict how they respond to future climate change [ 10 ]. Pikas of the family Ochotonidae are a unique lineage of small mammals adapted to life on high-elevation cold mountains [ 11 , 12 ], and have been recognized as ecologically important sentinels that are under threat due to climate change [ 13 , 14 ]. Species distribution modeling (SDM) approaches offer insights into how species might have responded to historical climate change [ 15 , 16 ]. By understanding ecological factors that govern species distributions today, we can project preferences of species into the past or the future, assuming niche conservatism over evolutionary time scales [ 17 – 21 ]. SDMs in their simplest form don’t account for critical biotic interactions that govern community structuring [ 22 – 24 ] and work best at shallow timescales, given that species can adapt to a changing climate over large evolutionary timescales [ 25 ] . Recent advancements in population genetics and coalescent theory have allowed researchers to capture genomic signatures of historical demographic fluctuations [ 26 – 30 ]. The pairwise sequentially Markovian coalescent (PSMC) analysis can work on a single genome by using nucleotide variation between two haploids to reconstruct demography [ 27 , 30 ]. It formalizes the relationship between the rate of coalescence of alleles at different loci and effective population size under a model of neutral evolution [ 27 , 30 ]. Although PSMCs are sensitive to selection effects, population structure, and the inclusion of mitochondrial/sex chromosomes; they can be reliably used to understand deep timescale (20kya-1000kya) demography when interpreted with caution [ 27 , 30 – 34 ]. Given that the two methods work optimally at different temporal and spatial scales (historic vs current and populations vs species), complementing PSMCs with SDMs allows researchers to understand how species have responded to climate change historically and how they might meet climate change in the Anthropocene [ 35 – 38 ]. The Indo-Tibetan plateau is a natural laboratory for the study of adaptations and species responses to complex glaciation patterns [ 39 ]. Its high elevation, extreme climate and dynamic geological history have challenged species to adapt or migrate to more favourable climates. While many researchers have studied demography on the plateau using species distribution models across geological timescales, it is focused heavily on plants [ 21 , 40 – 44 ]. Inferences of historical demography of species using population genomics from the region however, has been greatly limited to large mammals [ 45 – 49 ]. While there are a handful of studies on other taxa such on small mammals such as pikas[ 50 ], birds [ 36 , 51 ], plants [ 35 ] and, butterflies [ 52 ]; there is scope for a lot more work on the landscape given that glaciation histories are vastly different across this region [ 43 , 53 – 57 ]. Our understanding of pika biogeography has large implications for the conservation of these specialists [ 12 ]. Pika populations are expected to track cool habitats, move up the elevational gradients at warm time points in the Earth’s history, and expand to lower elevations at cool time points (Elevational shift hypothesis) [ 58 ] ( Figure 1 ). This hypothesis is supported by the: a) Eocene pika fossil record: Specimens found exclusively at high elevations (Qinghai-Tibetan Plateau, the Alps, and the Rocky Mountains) corresponding to a warming event [ 12 , 59 ]; b) Miocene - Pliocene fossil record: Specimens found at both low and high elevations corresponding to a cooling event [ 59 ]; c) Pleistocene fossil records: Specimens of O. princeps found at low elevations corresponding to glacial maximums [ 58 ]. In this study, we use a two-pronged approach (SDMs and PSMCs) to examine the historical demographic responses of three sympatric species of cold-adapted mountain specialist pikas from the Changthang biotic province of the Tibetan plateau. Ochotona ladacensis is a social, burrowing, high elevation specialist (>4200 m) occupying alpine meadows , O. nubrica is a social, burrowing species that utilizes a large elevational gradient (3300-6000 m) and is an alpine scrub specialist [ 60 ] . Ochotona macrotis is a rock-dwelling, asocial species that also occupies a large elevational gradient [ 60 ]. Given that pikas are cold adapted specialists, we expect that: a) population sizes and fundamental niches of all three species would be smallest at the last interglacial (LIG; 120kya) and largest at the last glacial maximum (LGM; 20kya) (Hypothesis 1); b) burrowing species would be more resilient to historical climate change when compared to rock-dwelling species since burrows provide better insulation against extreme weather when compared to rock refuges [ 98 ] (Hypothesis 2); c) O. ladacensis that exclusively occupies high elevations (>4500m) at current timescales, is likely to have small effective population sizes (Hypothesis 3). This is due to warming events causing severe bottlenecks, driving populations up the elevational gradient much more than they would for species that have wider elevational distribution. Understanding how small mammals from the Tibetan plateau have responded to historical climate change will likely be useful for managing alpine specialists. Download figure Open in new tab Figure 1: The elevational shift hypothesis describing how pika populations track cool climates during the Pleistocene glacial cycles (adapted from [ 58 ]). Results Draft genome assemblies for pikas Our genome completeness were assessed by both metrics of Megahit N5o and BUSCO. N50 measures were largest for O. ladacensis (25741 bp), followed by O. nubrica (15016 bp) and finally O. macrotis (N50 - 8050 bp). BUSCO scores had a large variance with 53% of complete BUSCO groups being sequenced for O. nubrica , 41% for O. ladacensis (N50 - 25741 bp) and 19% for O. macrotis; with 13%-18% being sequenced partially across species ( Figure 2A ). Our Orthovenn3 analyses indicated a large proportion of shared genes clusters across species examined (830 shared between species sequenced in this study; 445 shared between species sequenced in this study and O. princeps ). We also found equally large proportions of gene clusters being uniquely associated with the two partial genomes of O. nubrica (462) and O. ladacensis (387) ( Figure 2B ). Download figure Open in new tab Figure 2: A - Proportion of BUSCO groups sequenced from the mammalian database (N=9226 BUSCO groups) across the three species. B - Orthologous gene clusters found between species using OrthoVenn3 using the protein prediction tool Augustus. Projecting bio-climatic niches into the past: pikas have the largest niches at LGM and the smallest niches at LIG Ladakh pika ( Ochotona ladacensis ) - The best models constructed for the species had an AUC score of 0.967 with precipitation seasonality (Bio15), mean temperature of warmest quarter (Bio10), mean temperature of the wettest quarter (Bio8) and precipitation of the driest month (Bio 14) being the most critical variables driving the current distribution of the species ( Table 1 , Figure S1). Niche suitability (suitability > 0.4) for this species was inferred to be the smallest at LIG (∼4800 km 2 ) and expanded significantly at LGM (a 160-fold increase). Niches for this species were inferred to have shrunk post the LGM into the mid-Holocene warming (∼ 5-fold decrease) and expanded slightly at current timescales (2 % increase since mid-Holocene) ( Figure 3 ). Download figure Open in new tab Figure 3: Habitat suitability maps for different species (columns) at different time slices (rows): Current, Mid-Holocene (5kya), Last Glacial Maximum (20kya), and Last Interglacial (120kya). Warmer colors indicate higher suitability of habitats for species. View this table: View inline View popup Download powerpoint Table 1: Top variables (Percent contribution > 5%) identified in MAXENT models of the three species of pikas. Nubra pika ( Ochotona nubrica ) - Best models constructed for the species had an AUC score of 0.972, with the mean temperature of the driest quarter (Bio9), precipitation seasonality (Bio15) and temperature annual range (Bio7) being the most critical variables driving current distribution of the species ( Table 1 , Figure S1). Models projected back into the past indicated that a small area was suitable for the species at LIG (∼1600 km 2 ), with a significant expansion of bio-climatic niches at LGM (a 400-fold increase in suitable habitat). Subsequently, there was a significant reduction in niches for the species during the mid-Holocene warming (82% reduction since LGM). Suitable niches have expanded post the mid-Holocene into the current timescales (a 2.5-fold increase) ( Figure 3 ). Large-eared pika ( Ochotona macrotis ) - Best models constructed for the species had an AUC score of 0.918 with annual mean temperature (Bio1), precipitation of driest month (Bio14), mean temperature of the driest quarter (Bio9) and precipitation seasonality (Bio15) being the most critical variables driving current distribution of the species ( Table 1 , Figure S1). Niche suitability for this species was inferred to be the smallest at LIG (∼ 240000 km 2 ), expanded significantly at LGM (a 10-fold increase from LIG), shrunk post the LGM into the mid-Holocene (∼ 20% decrease from LGM) and expanded by a large extent at current timescales (67 % increase since mid-Holocene) ( Figure 3 ). Taken together, this yields support for our hypothesis that pikas have larger niches at LGM when compared to other time points (Hypothesis 1). Univariate response curves (Maxent models run only with one variable) showed that the rock-dwelling species O. macrotis had the largest bioclimatic envelope, which largely explains its extensive distribution today (Figure S1). In particular, O. macrotis presence was associated with a wide range of annual temperatures (Bio7), mean temperature in the wettest quarter (Bio8), precipitation in the driest month (Bio14), and precipitation seasonality (Bio15) when compared to O. nubrica and O. ladacensis . Pair wise sequential markovian coalescent (PSMC) simulations: pikas of the Indo-Tibetan plateau have the largest population sizes around the MIS-6 Estimates of effective population sizes of different species of pikas from the CBP followed consistent trends ( Figure 4 ). We estimated the largest effective population sizes for species during a) the warm LIG (∼120kya) for O. nubrica , and b) the cool MIS-6 (∼191kya) for O. macrotis and O. ladacensis ( Figure 4 ). Effective population sizes of all species were inferred to have continued to drop post the LIG, through the LGM (∼20kya), into the Mid-Holocene, with only O. nubrica showing a small recovery post the LGM ( Figure 4 ). We find that the species that are currently restricted to high elevations; namely O. ladacensis has the smallest effective population sizes (Hypothesis 2). Estimates of effective population sizes were not affected by lifestyle; with the rock-dwelling species examined ( O. macrotis ) being sandwiched between the two burrowing species ( O.ladacensis , O.nubrica ) ( Figure 4 ). Download figure Open in new tab Figure 4: Top left - PSMC plots of three sympatric species of pika from the Changthang Biotic province showing estimates of effective population sizes between 2kya-2mya; Bottom left - Corresponding average atmospheric surface temperature relative to present constructed using GCMs (data from [ 97 ]. Legend: a - Mid-Holocene; b-Global last glacial maxima; c-d- Local last glacial maxima/Changthang Biotic province glaciation event (MIS3 - MIS4); e - Global last inter-glacial period; f - MIS-6; g - MIS-13; h - MIS-21; i - MIS-22. Discussion Studying the historical demography of species in the context of climate change over evolutionary timescales provides a window into past species’ responses. In this study, we aimed to recover the demographic histories of cold-adapted high-elevational specialist pikas from the Indo-Tibetan Plateau using a combination of SDMs and PSMCs. We detected low effective population sizes for pika populations at LGM using PSMCs. This contrasts our species distribution models which suggest large geographic distributions of species at this time point. Accounting for local geology and historical climate experienced by different populations, in addition to consideration of global distribution patterns will be useful for the conservation planning of these species. Pikas populations move up the elevational gradient during LIG and expand to lower elevations during LGM The SDMs predicted and supported our hypotheses that the last glacial maximum (∼20kya) was likely a time point when pika populations were large (Hypothesis 1). We recovered large niches at the LGM and small niches at LIG for ⅔ species ( O. nubrica, O. ladacensis ) ( Figure 2 ). Habitats of moderate-high suitability (>0.4) were typically at lower elevations during the LGM when compared to current time scales following the elevational shift hypothesis [ 58 ] (Table S3). However, suitable habitats for the species Ochotona macrotis alone were at higher elevations during the LIG compared to the LGM (Table S3). This is potentially due to the complex glaciation history and historical climate of the Himalayas detailed below. SDMs for other cold-adapted high-elevation plants and frogs on the Tibetan plateau through the Pleistocene have revealed complex patterns involving both expansions ( Kobresia spp, Toxicodendron vernicifluum, Rhodiola spp, Allium mongolicum ) [ 40 , 61 – 63 ] and contractions ( Kobresia spp, Larix sibirica , Nanorana pleskei, Primula spp, Mirabilis himalaica, Anisodus tanguticus ) [ 41 , 44 , 62 , 64 – 67 ] during the LGM, linked to species ecology and thermal physiology. The historical climate of the Indo-Tibetan plateau and its impact on native flora A biogeographical and geological perspective is needed to understand the population-level responses of pikas to historical climate change. The complex geology of the ongoing upliftment of the Himalayas and the expansion of the Tibetan Plateau that began in the Pliocene has caused large-scale changes in climate, landscape, and vegetation at local scales, driving the evolution of new groups of mammals [ 68 , 69 ]. Our PSMCs placed the largest population sizes of species at MIS6 (∼190kya), pre-dating the LGM and LIG. This can be explained by the large disparity in the climatic records between the ‘global last glacial maxima (GLGM: 18kya to 22kya)’ experienced by the Northern Hemisphere and the ‘local last glacial maxima (LLGM: 40kya - 78kya)’ experienced by the Tibetan Plateau and the Himalayas during the ‘Batal stage’ as evidenced by geological records [ 43 , 53 – 57 ]. The largest glacial maximum events on the Indo-Tibetan plateau, however, are known to have occurred during the late Pliocene (MIS-12 - MIS-6 strata; 0.175mya - 0.5mya) and differ from glacial dynamics from the rest of the world [ 70 – 72 ]. More interestingly, glaciation patterns experienced by different parts of the Indo-Tibetan plateau such as the Changthang biotic province (CBP) and the Lauhul Himalayas are also significantly different [ 43 , 53 , 54 , 56 , 68 ]. The local glacial maxima experienced by the CBP from where the genomes were sampled, has been dated to span much of the MIS-3 (57 kya) - MIS-4 (71 kya) and was characterized by extremely low precipitation (snow) and small disconnected glacial valleys (<10km from current ranges) [ 43 , 53 , 54 , 56 , 68 ] as opposed to large interconnected glacial valleys predicted through the LGM in the Northern Hemisphere [ 73 , 74 ]. This pushes the timeline of the actual ‘tiny’ LLGM at our sites back to 78kya, likely lasting till 40kya [ 43 , 53 , 54 , 56 , 68 ], situating it way closer to the global last interglacial (120kya). Additionally, pollen analysis from the Tsokar Lake Basin on the CBP suggested periodic expansion of Juniperus communities and marshy vegetation (indicators of warm climate) and a decline of Alpine steppes (indicators of cold climate) at different time points including the GLGM [ 75 , 76 ] suggesting that vegetation communities did not support the distribution of pikas at this time points given the mismatch in habitats. Population responses of pikas on the Indo-Tibetan plateau Our PSMCs largely track the glacation history of the CBP with the largest population sizes of pikas either being at LIG (O. nubrica) or MIS-6 ( O. ladacensis and O. macrotis ) before crashing into the Holocene, much like other range restricted pikas such as O. sikkimaria [ 50 ] . Such disparity in the timing of population peaks between species might have arisen either due to O. nubrica showing delayed effective population size responses to changing climate [ 77 , 78 ] or the less likely but incorrect choice of mutation rates and generation times while scaling plots (Figure S2). Species on the Tibetan plateau have typically experienced population size increases dating back to the MIS-6 period, and includes domestic and wild mammals [ 47 , 49 ], birds [ 68 , 74 ] and reptiles [ 79 ]. Of the limited studies that have examined demographic history through PSMCs on the Tibetan plateau, patterns recorded for the Iberian Ibex [ 45 ] and the Tibetan Wolf [ 46 ] show contrasting patterns. Given the absence of any glaciation events post MIS-3 and MIS-4 (which was also very mild on the CBP), pikas are likely to have remained isolated on cool mountain tops causing them to inbreed for long periods, leading to a decline in population size in the absence of out-breeding rescue [ 50 , 58 , 59 ]. Recent studies of the subgenus Ochotona have suggested that outbreeding and adaptive introgressions between O. curzoniae and other members of the sub-clade have rescued populations over historic timescales [ 50 ]. While recent inbreeding is known to affect the estimation of coalescence times and Ne [ 80 , 81 ], inbreeding at older time scales doesn’t necessarily affect these estimates [ 82 , 83 ]. We, therefore, believe that the population demographic signatures obtained by us across species are driven by the lack of severe glaciation events in the CBP in the last 1mya and inbreeding depression post the LIG. This requires further validation through other analyses of formal runs of homozygosity analysis (ROH) [ 80 , 82 ]. Conclusion While our PSMCs are likely to be true depictions of population responses of pikas from the Indo-Tibetan plateau (Changthang Biotic Province), our SDMs are capturing the general patterns of species responses (elevational shift) across their geographic range. We find that lifestyle (burrowing and rock-dwelling), did not interact with population trajectories recovered. We hypothesize that pikas in the Changthang Biotic Province had the largest population sizes during the cool period of MIS-6, were forced up the elevational gradient during the warm period of LIG and remained isolated and inbred in the absence of a glaciation rescue causing the drop in effective population sizes. Sampling demographics of cold-adapted specialists on the Indo-Tibetan plateau in a spatially explicit manner with consideration of local geology should shed more insights for conservation management of these species under threat. Methods Species distribution modeling (SDMs) Occurrence records: We collected primary points of occurrence of all three species of pikas between 2017 and 2023 in the Tibetan Marginal Mountains, Trans-Himalayas, Ladakh, India. We supplemented this data with secondary points from published literature, GBIF (Global Biodiversity Information Facility), iNaturalist and museum records. Only confirmed observations on citizen science platforms with a photo where species identification could be verified were used. Our dataset included 63 occurrence points for O. ladacensis, 203 points for O. macrotis, and 93 points for O. nubrica . Raster data processing Raster layers of bioclimatic variables (Bio1 - Bio19) were sourced from Bioclim [ 84 ] at 1km resolution for four different time slices: LIG (120 kya), LGM (20 kya), Mid-Holocene (5000 ya) and current [ 85 ]. While the LIG, Mid-Holocene, and current data were available at a 30-second resolution, the LGM data was available only at a coarser, 2.5-minute resolution. The LGM data was interpolated and downscaled to a 30-second resolution to match other raster datasets. A conservatively large area (coordinates - 62.2749, 25.5666: 104.9499, 46.8416) encompassing the distribution of the most widely distributed species was chosen as the study area. ( Figure 3 ). The advanced digitizing toolbox was used to create this shape file, and all rasters were clipped to this mask. While doing so, all rasters were converted to the same projection system (EPSG 4326: WGS 84), with no data pixels being set to -9999 and converted to an ‘ASCII’ file format as required by MAXENT. Processed rasters were imported into QGIS, and the spatial extent, pixel sizes, and range of pixel values were inspected before further downstream processing. The ENMevaluate function was run for occurrence data at current timescales to determine the best model settings for MAXENT. Model performance under different parameter settings was evaluated for each species using AIC scores, and the best model settings for each species are presented (Table S1). Occurrence records were spatially thinned (1 point/km) to avoid data over-fitting. Post spatial thinning, the dataset had 33 points for O. ladacensis , 37 points for O. nubrica , and 121 points for O. macrotis . A bias file was also created to account for sampling biases by modifying existing R scripts (Josh Banta Lab, https://sites.google.com/site/thebantalab/tutorials#h.zsfviftgr2 ). Genome sequencing and De-novo assembly One tissue sample was collected for each of the three species using standard trapping methods with baited Tomahawks or from opportunistic roadkills. These tissues were stored in ethanol at -20°C for downstream processing. DNA was extracted from tissues using QIAGEN Blood and Tissue Kits with a modified protocol to include an additional digestion step with AL buffer for 20 minutes at 75°C and a digestion with RNase. Extracted DNA was quantified using a Qubit fluorometer and sequenced using an Illumina platform sequencer at 30x coverage. The reads obtained for each sample were subject to quality checks using FASTQC [ 86 ] and trimmed with a quality score cutoff of 30 using TRIM GALORE [ 87 ]. Trimmed reads were repaired and quality checked before de-novo assembly using the Megahit [ 88 ] pipeline under default settings. Genome quality assessment Quality and completeness of genomes assembled using Megahit were assessed using BUSCO [ 89 ] against the mammalian database using Augustus [ 90 ]. Proteins predicted were merged together for each species, prefixed with species identifiers and prepped for comparative genomic analyses. Orthovenn3 [ 91 , 92 ] was used to identify orthologous gene clusters across the draft genomes assembled. It does so by employing BLASTP in conjunction with constructing a similarity graph followed by clustering this network [ 91 , 92 ]. Pairwise sequential Markovian Coalescent (PSMC) analysis While demographic inferences are robust to genome assembly quality [ 78 ], the presence of repeats in mammalian genomes, mitochondrial DNA, and sex chromosomes are known to affect demographic inference [ 32 , 82 ]. Assemblies were therefore, hard masked for repeats using RepeatMasker [ 93 ], following which, contigs were mapped to reference mitochondrial genomes of O. princeps (NCBI Accession number: NC_005358.1 ) and sex chromosomes of O. cuniculus (NCBI Accession number: NC_067395.1 ) using BLAST to identify and remove such contigs in assemblies. Trimmed reads from the Illumina sequencing experiment were marked for PCR duplicates and mapped back to this assembly to generate a binary alignment file (.bam) using ‘samtools mpileup’ [ 94 ]. A ‘diploid.fq’ file was generated following standard PSMC protocols for downstream processing using default depth cutoffs when variant calling. The fq2psmca tool was used to generate psmcfa files, and PSMCs were performed with 100 bootstraps under settings that detected more than 1000 recombinations at each of the timesteps (ON: -N25 -t5 -r5 -p "26*2+4+7+1"; OL: -N25 -t15 -r5 -p "4+25*2+4+6"; OM: -N25 -t9 -r5 -b -p "26*2+4+7+1") after post ensuring that the average depth of sequencing was >30x [ 27 , 30 ] (Supplementary Figures 1). The psmc_plot.pl function was used to scale PSMCs and generated text output for plotting bootstraps on R. While scaling plots, the mutation rate for all species was set at 2.96*10 -9 following published literature on pikas and rodents (mean mutation rate for nuclear protein-coding genes) [ 50 , 95 , 96 ]. Different generation times were used for different species of pikas based on both field guides and natural history data collected while on the field between 2019-2023 ( O. nubrica - 0.5 years, O. ladacensis - 1 year, O. macrotis - 1 year) after scaling them by 2 [ 11 , 32 ]. We additionally modeled the sensitivity of the PSMC results obtained using various parameters of mutation rates and generation times (Figure S2). Author contributions HK - conceptualization, analysis, writing, sample collection; NV - conceptualization, inputs on analysis, reviewing the manuscript; NR - conceptualization, reviewing the analysis, reviewing the manuscript, sample collection. Acknowledgments We thank the Jammu and Kashmir Wildlife Department and Department of Wildlife Protection, Ladakh for permits to collect small tissue samples from caught animals/ roadkills (471-75/WLWL/E.P/2022-23) for whole genome sequencing, Senan D’Souza & Dhanesh Ponnu for help with trapping animals and mapping distributions of species in Ladakh, Dr. Meghana Natesh, Dr. Nivetha Murugesan and Vinay KL for comments on the work, and IISER Tirupati for funding and logistics. Funder Information Declared Indian Institute of Science Education and Research, Tirupati, https://ror.org/032d0e990 References 1. ↵ Dueñas A , Jiménez-Uzcátegui G , Bosker T . The effects of climate change on wildlife biodiversity of the galapagos islands . Climate Change Ecology . 2021 ; 2 : 100026 . 2. ↵ Bellard C , Bertelsmeier C , Leadley P , Thuiller W , Courchamp F . Impacts of climate change on the future of biodiversity: Biodiversity and climate change . Ecol Lett . 2012 ; 15 : 365 – 77 . OpenUrl CrossRef PubMed 3. ↵ Engler R , Randin CF , Thuiller W , Dullinger S , Zimmermann NE , Araújo MB , et al. 21st century climate change threatens mountain flora unequally across Europe: Climate change impacts on mountain florae . Glob Chang Biol . 2011 ; 17 : 2330 – 41 . OpenUrl CrossRef 4. Rubidge EM , Patton JL , Lim M , Burton AC , Brashares JS , Moritz C . Climate-induced range contraction drives genetic erosion in an alpine mammal . Nat Clim Chang . 2012 ; 2 : 285 – 8 . OpenUrl CrossRef 5. ↵ DeSaix MG , George TL , Seglund AE , Spellman GM , Zavaleta ES , Ruegg KC . Forecasting climate change response in an alpine specialist songbird reveals the importance of considering novel climate . Divers Distrib . 2022 ; 28 : 2239 – 54 . OpenUrl CrossRef 6. ↵ Olden JD , Poff NL , Bestgen KR . Trait synergisms and the rarity, extirpation, and extinction risk of desert fishes . Ecology . 2008 ; 89 : 847 – 56 . OpenUrl CrossRef PubMed Web of Science 7. Kopf RK , Shaw C , Humphries P . Trait-based prediction of extinction risk of small-bodied freshwater fishes . Conserv Biol . 2017 ; 31 : 581 – 91 . OpenUrl CrossRef PubMed 8. Chichorro F , Juslén A , Cardoso P . A review of the relation between species traits and extinction risk . Biol Conserv . 2019 ; 237 : 220 – 9 . OpenUrl CrossRef 9. ↵ Chichorro F , Urbano F , Teixeira D , Väre H , Pinto T , Brummitt N , et al. Trait-based prediction of extinction risk across terrestrial taxa . Biol Conserv . 2022 ; 274 : 109738 . 10. ↵ Oberle B , Schaal BA . Responses to historical climate change identify contemporary threats to diversity in Dodecatheon . Proc Natl Acad Sci U S A . 2011 ; 108 : 5655 – 60 . OpenUrl Abstract / FREE Full Text 11. ↵ Smith AT , Johnston CH , Alves PC , Hacklander K , editors. Lagomorphs . Baltimore, MD : Johns Hopkins University Press ; 2018 . 12. ↵ Wang X , Liang D , Jin W , Tang M , Liu S , Zhang P . Out of Tibet: Genomic Perspectives on the Evolutionary History of Extant Pikas . Mol Biol Evol . 2020 ; 37 : 1577 – 92 . OpenUrl CrossRef PubMed 13. ↵ Beever EA , Brussard PF , Berger J . Patterns of apparent extirpation among isolated populations of pikas (Ochotona princeps) in the great basin . J Mammal . 2003 ; 84 : 37 – 54 . OpenUrl CrossRef Web of Science 14. ↵ Beever EA , Ray C , Wilkening JL , Brussard PF , Mote PW . Contemporary climate change alters the pace and drivers of extinction: CLIMATE CHANGE ALTERS PACE, DRIVERS OF LOSSES . Glob Chang Biol . 2011 ; 17 : 2054 – 70 . OpenUrl CrossRef 15. ↵ Park I-K , Borzée A , Park J , Min S-H , Zhang Y-P , Li S-R , et al. Past, present, and future predictions on the suitable habitat of the Slender racer (Orientocoluber spinalis) using species distribution models . Ecol Evol . 2022 ; 12 : e9169 . OpenUrl CrossRef PubMed 16. ↵ Franklin J . Species distribution modelling supports the study of past, present and future biogeographies . J Biogeogr . 2023 ; 50 : 1533 – 45 . OpenUrl CrossRef 17. ↵ Phillips SJ , Dudík M , Schapire RE. A maximum entropy approach to species distribution modeling. In: Proceedings of the twenty-first international conference on Machine learning . New York, NY, USA : Association for Computing Machinery ; 2004 . p. 83. 18. Phillips SJ , Anderson RP , Schapire RE . Maximum entropy modeling of species geographic distributions . Ecol Modell . 2006 . 19. Hope AG , Waltari E , Malaney JL , Payer DC , Cook JA , Talbot SL . Arctic biodiversity: increasing richness accompanies shrinking refugia for a cold-associated tundra fauna . Ecosphere . 2015 ; 6 :art159. 20. Zhang Y-Z , Zhu R-W , Zhong D-L , Zhang J-Q . Nunataks or massif de refuge? A phylogeographic study of Rhodiola crenulata (Crassulaceae) on the world’s highest sky islands . BMC Evol Biol . 2018 ; 18 : 154 . 21. ↵ Zhang L , Zhu L , Li Y , Zhu W , Chen Y . Maxent Modelling Predicts a Shift in Suitable Habitats of a Subtropical Evergreen Tree (Cyclobalanopsis glauca (Thunberg) Oersted) under Climate Change Scenarios in China . For Trees Livelihoods . 2022 ; 13 : 126 . 22. ↵ Carneiro LR de A , Lima AP , Machado RB , Magnusson WE . Limitations to the Use of Species-Distribution Models for Environmental-Impact Assessments in the Amazon . PLoS One . 2016 ; 11 : e0146543 . OpenUrl CrossRef PubMed 23. Dormann CF , Bobrowski M , Dehling DM , Harris DJ , Hartig F , Lischke H , et al. Biotic interactions in species distribution modelling: 10 questions to guide interpretation and avoid false conclusions . Glob Ecol Biogeogr . 2018 ; 27 : 1004 – 16 . OpenUrl CrossRef 24. ↵ Cosentino F , Seamark ECJ , Van Cakenberghe V , Maiorano L . Not only climate: The importance of biotic interactions in shaping species distributions at macro scales . Ecol Evol . 2023 ; 13 : e9855 . OpenUrl CrossRef PubMed 25. ↵ Leonardi M , Boschin F , Boscato P , Manica A . Following the niche: the differential impact of the last glacial maximum on four European ungulates . Commun Biol . 2022 ; 5 : 1038 . OpenUrl CrossRef PubMed 26. ↵ Drummond AJ , Rambaut A , Shapiro B , Pybus OG . Bayesian coalescent inference of past population dynamics from molecular sequences . Mol Biol Evol . 2005 ; 22 : 1185 – 92 . OpenUrl CrossRef PubMed Web of Science 27. ↵ Li H , Durbin R . Inference of human population history from individual whole-genome sequences . Nature . 2011 ; 475 : 493 – 6 . OpenUrl CrossRef PubMed Web of Science 28. Cornuet J-M , Pudlo P , Veyssier J , Dehne-Garcia A , Gautier M , Leblois R , et al. DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data . Bioinformatics . 2014 ; 30 : 1187 – 9 . OpenUrl CrossRef PubMed Web of Science 29. Schiffels S , Wang K . MSMC and MSMC2: The Multiple Sequentially Markovian Coalescent . Methods Mol Biol . 2020 ; 2090 : 147 – 66 . OpenUrl CrossRef PubMed 30. ↵ Mather N , Traves SM , Ho SYW . A practical introduction to sequentially Markovian coalescent methods for estimating demographic history from genomic data . Ecol Evol . 2020 ; 10 : 579 – 89 . OpenUrl CrossRef PubMed 31. Kozma R , Melsted P , Magnússon KP , Höglund J . Looking into the past - the reaction of three grouse species to climate change over the last million years using whole genome sequences . Mol Ecol . 2016 ; 25 : 570 – 80 . OpenUrl CrossRef 32. ↵ Nadachowska-Brzyska K , Burri R , Smeds L , Ellegren H . PSMC analysis of effective population sizes in molecular ecology and its application to black-and-white Ficedula flycatchers . Mol Ecol . 2016 ; 25 : 1058 – 72 . OpenUrl CrossRef PubMed 33. Natesh M , Vinay KL , Ghosh S , Jayapal R , Mukherjee S , Vijay N , et al. Contrasting trends of population size change for two Eurasian Owlet species—Athene brama and Glaucidium radiatum from south Asia over the late Quaternary . Front Ecol Evol . 2020 ; 8 . 34. ↵ Foote AD , Hooper R , Alexander A , Baird RW , Baker CS , Ballance L , et al. Runs of homozygosity in killer whale genomes provide a global record of demographic histories . Mol Ecol . 2021 ; 30 : 6162 – 77 . OpenUrl CrossRef 35. ↵ Zhao Y-P , Fan G , Yin P-P , Sun S , Li N , Hong X , et al. Resequencing 545 ginkgo genomes across the world reveals the evolutionary history of the living fossil . Nat Commun . 2019 ; 10 : 4201 . OpenUrl CrossRef PubMed 36. ↵ Song K , Gao B , Halvarsson P , Fang Y , Jiang Y-X , Sun Y-H , et al. Genomic analysis of demographic history and ecological niche modeling in the endangered Chinese Grouse Tetrastes sewerzowi . BMC Genomics . 2020 ; 21 : 581 . 37. Patil AB , Vajja SS , Raghavendra S , Satish BN , Kushalappa CG , Vijay N . Jack of all trades: Genome assembly of Wild Jack and comparative genomics of Artocarpus . Front Plant Sci . 2022 ; 13 : 1029540 . 38. ↵ Huynh S , Cloutier A , Chen G , Chan DTC , Lam DK , Huyvaert KP , et al. Whole-genome Analyses Reveal Past Population Fluctuations and Low Genetic Diversities of the North Pacific Albatrosses . Mol Biol Evol . 2023 ; 40 . 39. ↵ Zhou Z , Deng T . The Tibetan Plateau is a natural laboratory for studying organic evolution and environmental change . Sci China Earth Sci . 2020 ; 63 : 169 – 71 . OpenUrl CrossRef 40. ↵ Li Y-C , Zhong D-L , Rao G-Y , Wen J , Ren Y , Zhang J-Q . Gone with the trees: Phylogeography of Rhodiola sect. Trifida (Crassulaceae) reveals multiple refugia on the Qinghai-Tibetan Plateau . Mol Phylogenet Evol . 2018 ; 121 : 110 – 20 . OpenUrl CrossRef PubMed 41. ↵ Zhang H-X , Zhang M-L , Williams DM . Genetic evidence and species distribution modelling reveal the response of Larix sibirica and its related species to Quaternary climatic and ancient historical events . Biochem Syst Ecol . 2014 ; 54 : 316 – 25 . OpenUrl CrossRef 42. Hu Z , Guo K , Jin S , Pan H . The influence of climatic changes on distribution pattern of six typical Kobresia species in Tibetan Plateau based on MaxEnt model and geographic information system . Theor Appl Climatol . 2019 ; 135 : 375 – 90 . OpenUrl CrossRef 43. ↵ Yang W , Feng L , Jiao P , Xiang L , Yang L , Olonova MV , et al. Out of the Qinghai-Tibet plateau: Genomic biogeography of the alpine monospecific genus Megadenia (Biscutelleae , Brassicaceae). Mol Ecol . 2023 ; 32 : 492 – 503 . OpenUrl CrossRef PubMed 44. ↵ Zhang H , Zhao B , Huang T , Chen H , Yue J , Tian Y . Responses of the Distribution Pattern of the Suitable Habitat of Juniperus tibetica Komarov to Climate Change on the Qinghai-Tibet Plateau . For Trees Livelihoods . 2023 ; 14 : 434 . 45. ↵ Robin M , Ferrari G , Akgül G , Münger X , von Seth J , Schuenemann VJ , et al. Ancient mitochondrial and modern whole genomes unravel massive genetic diversity loss during near extinction of Alpine ibex . Mol Ecol . 2022 ; 31 : 3548 – 65 . OpenUrl CrossRef 46. ↵ Zhang W , Fan Z , Han E , Hou R , Zhang L , Galaverni M , et al. Hypoxia adaptations in the grey wolf (Canis lupus chanco) from Qinghai-Tibet Plateau . PLoS Genet . 2014 ; 10 : e1004466 . OpenUrl CrossRef PubMed 47. ↵ Zhang Y , Xue X , Liu Y , Abied A , Ding Y , Zhao S , et al. Genome-wide comparative analyses reveal selection signatures underlying adaptation and production in Tibetan and Poll Dorset sheep . Sci Rep . 2021 ; 11 : 2466 . OpenUrl CrossRef PubMed 48. Martchenko D , Prewer E , Latch EK , Kyle CJ , Shafer ABA. Population Genomics of Ungulates. In: Population Genomics . Cham : Springer International Publishing ; 2018 . p. 185 – 209 . 49. ↵ Wang M-S , Wang S , Li Y , Jhala Y , Thakur M , Otecko NO , et al. Ancient Hybridization with an Unknown Population Facilitated High-Altitude Adaptation of Canids . Mol Biol Evol . 2020 ; 37 : 2616 – 29 . OpenUrl CrossRef PubMed 50. ↵ Ge D , Wen Z , Feijó A , Lissovsky A , Zhang W , Cheng J , et al. Genomic Consequences of and Demographic Response to Pervasive Hybridization Over Time in Climate-Sensitive Pikas . Mol Biol Evol . 2023 ; 40 . 51. ↵ Song K , Gao B , Halvarsson P , Fang Y , Klaus S , Jiang Y-X , et al. Demographic history and divergence of sibling grouse species inferred from whole genome sequencing reveal past effects of climate change . BMC Ecol Evol . 2021 ; 21 : 194 . 52. ↵ Zhao Y , Su C , He B , Nie R , Wang Y , Ma J , et al. Dispersal from the Qinghai-Tibet plateau by a high-altitude butterfly is associated with rapid expansion and reorganization of its genome . Nat Commun . 2023 ; 14 : 8190 . OpenUrl CrossRef PubMed 53. ↵ Owen LA , Finkel RC , Caffee MW . A note on the extent of glaciation throughout the Himalaya during the global Last Glacial Maximum . Quat Sci Rev . 2002 ; 21 : 147 – 57 . OpenUrl CrossRef 54. ↵ Owen LA , Caffee MW , Bovard KR , Finkel RC , Sharma MC . Terrestrial cosmogenic nuclide surface exposure dating of the oldest glacial successions in the Himalayan orogen: Ladakh Range, northern India . GSA Bulletin . 2006 ; 118 : 383 – 92 . OpenUrl Abstract / FREE Full Text 55. Zhang W , Niu Y , Yan L , Cui Z , Li C , Mu K . Late Pleistocene glaciation of the Changbai Mountains in northeastern China . Chin Sci Bull . 2008 ; 53 : 2672 – 84 . OpenUrl CrossRef 56. ↵ Hedrick KA , Seong YB , Owen LA , Caffee MW , Dietsch C . Towards defining the transition in style and timing of Quaternary glaciation between the monsoon-influenced Greater Himalaya and the semi-arid Transhimalaya of Northern India . Quat Int . 2011 ; 236 : 21 – 33 . OpenUrl CrossRef 57. ↵ Das S , Murari MK , Sharma MC , Saini R , Jaiswal MK , Kumar PV , et al. Reconstraining the Quaternary glacial history of Lahaul Himalaya, northern India . Quat Sci Rev . 2023 ; 316 : 108258 . 58. ↵ Galbreath KE , Hafner DJ , Zamudio KR . When cold is better: climate-driven elevation shifts yield complex patterns of diversification and demography in an alpine specialist (American pika , Ochotona princeps). Evolution . 2009 ; 63 : 2848 – 63 . OpenUrl PubMed 59. ↵ Ge D , Wen Z , Xia L , Zhang Z , Erbajeva M , Huang C , et al. Evolutionary history of lagomorphs in response to global environmental change . PLoS One . 2013 ; 8 : e59668 . OpenUrl CrossRef PubMed 60. ↵ Kumar H , Ponnu D , D’souza S , Das D , Thundup S , Rajamani N. Scale-dependent niche segregation along dimensions of rocks, vegetation and elevation in sympatric pikas of Ladakh . Preprints. 2023 . 61. ↵ Zhang Y , Yu Q , Zhang Q , Hu X , Hu J , Fan B . Regional-scale differentiation and phylogeography of a desert plant Allium mongolicum (Liliaceae) inferred from chloroplast DNA sequence variation . Plant Syst Evol . 2017 ; 303 : 451 – 66 . OpenUrl CrossRef 62. ↵ Hu Z , Guo K , Jin S , Pan H . The influence of climatic changes on distribution pattern of six typical Kobresia species in Tibetan Plateau based on MaxEnt model and geographic information system . Theor Appl Climatol . 2019 ; 135 : 375 – 90 . OpenUrl CrossRef 63. ↵ Wang L , Li Y , Noshiro S , Suzuki M , Arai T , Kobayashi K , et al. Stepped Geomorphology Shaped the Phylogeographic Structure of a Widespread Tree Species (Toxicodendron vernicifluum , Anacardiaceae) in East Asia. Front Plant Sci . 2022 ; 13 : 920054 . 64. ↵ Wang B , Xie F , Li J , Wang G , Li C , Jiang J . Phylogeographic investigation and ecological niche modelling of the endemic frog species Nanorana pleskei revealed multiple refugia in the eastern Tibetan Plateau . PeerJ . 2017 ; 5 : e3770 . OpenUrl CrossRef PubMed 65. Ren G , Mateo RG , Liu J , Suchan T , Alvarez N , Guisan A , et al. Genetic consequences of Quaternary climatic oscillations in the Himalayas: Primula tibetica as a case study based on restriction site-associated DNA sequencing . New Phytol . 2017 ; 213 : 1500 – 12 . OpenUrl CrossRef PubMed 66. Ren G , Mateo RG , Guisan A , Conti E , Salamin N . Species divergence and maintenance of species cohesion of three closely related Primula species in the Qinghai-Tibet Plateau . J Biogeogr . 2018 ; 45 : 2495 – 507 . OpenUrl CrossRef 67. ↵ Rana HK , Luo D , Rana SK , Sun H . Geological and Climatic Factors Affect the Population Genetic Connectivity in Mirabilis himalaica (Nyctaginaceae): Insight From Phylogeography and Dispersal Corridors in the Himalaya-Hengduan Biodiversity Hotspot . Front Plant Sci . 2019 ; 10 : 1721 . OpenUrl PubMed 68. ↵ Lei F , Qu Y , Song G . Species diversification and phylogeographical patterns of birds in response to the uplift of the Qinghai-Tibet Plateau and Quaternary glaciations . Curr Zool . 2014 ; 60 : 149 – 61 . OpenUrl CrossRef 69. ↵ Harris SA , Jin H , He R , Yang S . Tessellons, topography, and glaciations on the Qinghai-Tibet Plateau . Sci Cold Arid Reg . 2018 . doi: 10.3724/SP.J.1226.2018.00187 . OpenUrl CrossRef 70. ↵ Wei Z , Zhijiu C , Yonghua L . Review of the timing and extent of glaciers during the last glacial cycle in the bordering mountains of Tibet and in East Asia . Quat Int . 2006 ; 154 – 155 :32–43. 71. Zhou S , Wang X , Wang J , Xu L . A preliminary study on timing of the oldest Pleistocene glaciation in Qinghai–Tibetan Plateau . Quat Int . 2006 ; 154 – 155 :44–51. 72. ↵ Chevalier M-L , Wang S , Replumaz A , Li H. Marine oxygen Isotope Stage (MIS)-6 glacial advances on the Tibetan Plateau more extensive than during MIS-2 due to more abundant precipitation . Acta Geol Sin - Engl Ed . 2022 ; 96 : 1484 – 94 . 73. ↵ Du Y-R , Guo S-C , Wang Z-F , Ci H-X , Cai Z-Y , Zhang Q , et al. Demographic history of the Tibetan antelope Pantholops hodgsoni (chiru) . J Syst Evol . 2010 ; 48 : 490 – 6 . OpenUrl CrossRef 74. ↵ Qu Y , Lei F , Zhang R , Lu X . Comparative phylogeography of five avian species: implications for Pleistocene evolutionary history in the Qinghai-Tibetan plateau . Mol Ecol . 2010 ; 19 : 338 – 51 . OpenUrl CrossRef PubMed 75. ↵ Bhattacharyya A . Vegetation and climate during the last 30,000 years in Ladakh . Palaeogeogr Palaeoclimatol Palaeoecol . 1989 ; 73 : 25 – 38 . OpenUrl CrossRef GeoRef Web of Science 76. ↵ Ranhotra PS , Bhattacharyya A , Kotlia BS . Vegetation and climatic changes around Lamayuru, Trans-Himalaya during the last 35 kyr B.P . Journal of Palaeosciences . 2007 ; 56 : 117 – 26 . OpenUrl CrossRef 77. ↵ Dussex N , von Seth J , Knapp M , Kardailsky O , Robertson BC , Dalén L . Complete genomes of two extinct New Zealand passerines show responses to climate fluctuations but no evidence for genomic erosion prior to extinction . Biol Lett . 2019 ; 15 : 20190491 . 78. ↵ Patton AH , Margres MJ , Stahlke AR , Hendricks S , Lewallen K , Hamede RK , et al. Contemporary Demographic Reconstruction Methods Are Robust to Genome Assembly Quality: A Case Study in Tasmanian Devils . Mol Biol Evol . 2019 ; 36 : 2906 – 21 . OpenUrl CrossRef PubMed 79. ↵ Gao W , Yu C-X , Zhou W-W , Zhang B-L , Chambers EA , Dahn HA , et al. Species Persistence with Hybridization in Toad-Headed Lizards Driven by Divergent Selection and Low Recombination . Mol Biol Evol . 2022 ; 39 :msac064. 80. ↵ Lu C-W , Yao C-T , Hung C-M . Domestication obscures genomic estimates of population history . Mol Ecol . 2022 ; 31 : 752 – 66 . OpenUrl CrossRef 81. ↵ Jaspers C , Ehrlich M , Pujolar JM , Künzel S , Bayer T , Limborg MT , et al. Invasion genomics uncover contrasting scenarios of genetic diversity in a widespread marine invader . Proc Natl Acad Sci U S A . 2021 ; 118 . 82. ↵ Nadachowska-Brzyska K , Li C , Smeds L , Zhang G , Ellegren H . Temporal Dynamics of Avian Populations during Pleistocene Revealed by Whole-Genome Sequences . Curr Biol . 2015 ; 25 : 1375 – 80 . OpenUrl CrossRef PubMed 83. ↵ von Seth J , van der Valk T , Lord E , Sigeman H , Olsen R-A , Knapp M , et al. Genomic trajectories of a near-extinction event in the Chatham Island black robin . BMC Genomics . 2022 ; 23 : 747 . 84. ↵ Fick SE , Hijmans RJ . WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas . Int J Climatol . 2017 ; 37 : 4302 – 15 . OpenUrl CrossRef PubMed 85. ↵ Karger DN , Conrad O , Böhner J , Kawohl T , Kreft H , Soria-Auza RW , et al. Climatologies at high resolution for the earth’s land surface areas . Sci Data . 2017 ; 4 : 170122 . 86. ↵ Andrews S. FastQC: a quality control tool for high throughput sequence data . 2010. 87. ↵ Martin M . Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal . 2011 ; 17 : 10 – 2 . 88. ↵ Li D , Liu C-M , Luo R , Sadakane K , Lam T-W . MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph . Bioinformatics . 2015 ; 31 : 1674 – 6 . OpenUrl CrossRef PubMed 89. ↵ Manni M , Berkeley MR , Seppey M , Zdobnov EM . BUSCO: Assessing genomic data quality and beyond . Curr Protoc . 2021 ; 1 : e323 . OpenUrl CrossRef 90. ↵ Stanke M , Steinkamp R , Waack S , Morgenstern B . AUGUSTUS: a web server for gene finding in eukaryotes . Nucleic Acids Res . 2004 ; 32 Web Server issue:W309–12. 91. ↵ Wang Y , Coleman-Derr D , Chen G , Gu YQ . OrthoVenn: a web server for genome wide comparison and annotation of orthologous clusters across multiple species . Nucleic Acids Res . 2015 ; 43 : W78 – 84 . OpenUrl CrossRef PubMed 92. ↵ Sun J , Lu F , Luo Y , Bie L , Xu L , Wang Y . OrthoVenn3: an integrated platform for exploring and visualizing orthologous data across genomes . Nucleic Acids Res . 2023 ; 51 : W397 – 403 . OpenUrl CrossRef PubMed 93. ↵ Tarailo-Graovac M , Chen N . Using RepeatMasker to identify repetitive elements in genomic sequences . Curr Protoc Bioinformatics . 2009 ;Chapter 4:4.10.1–4.10.14. 94. ↵ Danecek P , Bonfield JK , Liddle J , Marshall J , Ohan V , Pollard MO , et al. Twelve years of SAMtools and BCFtools . Gigascience . 2021 ; 10 . 95. ↵ Teng H , Zhang Y , Shi C , Mao F , Cai W , Lu L , et al. Population Genomics Reveals Speciation and Introgression between Brown Norway Rats and Their Sibling Species . Mol Biol Evol . 2017 ; 34 : 2214 – 28 . OpenUrl CrossRef PubMed 96. ↵ Zeng L , Ming C , Li Y , Su L-Y , Su Y-H , Otecko NO , et al. Out of Southern East Asia of the Brown Rat Revealed by Large-Scale Genome Sequencing . Mol Biol Evol . 2018 ; 35 : 149 – 58 . OpenUrl CrossRef PubMed 97. ↵ Bintanja R , van de Wal RSW . North American ice-sheet dynamics and the onset of 100,000-year glacial cycles . Nature . 2008 ; 454 : 869 – 72 . OpenUrl CrossRef GeoRef PubMed Web of Science 98. ↵ Kumar H. Adaptations of pikas to alpine stressors on the Tibetan Marginal Mountains. PhD Thesis. IISER Tirupati . 2024 . View the discussion thread. Back to top Previous Next Posted June 28, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Asian pika populations track local glaciation events through the Pleistocene Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Asian pika populations track local glaciation events through the Pleistocene Harsha Kumar , Nagarjun Vijay , R Nandini bioRxiv 2025.06.25.661461; doi: https://doi.org/10.1101/2025.06.25.661461 Share This Article: Copy Citation Tools Asian pika populations track local glaciation events through the Pleistocene Harsha Kumar , Nagarjun Vijay , R Nandini bioRxiv 2025.06.25.661461; doi: https://doi.org/10.1101/2025.06.25.661461 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Genomics Subject Areas All Articles Animal Behavior and Cognition (7635) Biochemistry (17697) Bioengineering (13894) Bioinformatics (41951) Biophysics (21455) Cancer Biology (18592) Cell Biology (25507) Clinical Trials (138) Developmental Biology (13380) Ecology (19903) Epidemiology (2067) Evolutionary Biology (24321) Genetics (15610) Genomics (22509) Immunology (17737) Microbiology (40398) Molecular Biology (17182) Neuroscience (88618) Paleontology (667) Pathology (2833) Pharmacology and Toxicology (4825) Physiology (7641) Plant Biology (15158) Scientific Communication and Education (2046) Synthetic Biology (4296) Systems Biology (9825) Zoology (2271)