De novo Assembly, Annotation, and Analysis of Transcriptome Data of the Ladakh Ground Skink Provide Genetic Information on High-Altitude Adaptation

The Himalayan Arc is recognized as a global biodiversity hotspot. Among its numerous cryptic and undiscovered organisms, this composite high-mountain ecosystem harbors many taxa with adaptations to life in high elevations. However, evolutionary patterns and genomic features have been relatively rarely studied in Himalayan vertebrates. Here, we provide the first well-annotated transcriptome of a Greater Himalayan reptile species, the Ladakh Ground skink Asymblepharus ladacensis (Squamata: Scincidae). Based on tissues from the brain, an embryonic disc, and pooled organ material, using pair-end Illumina NextSeq 500 RNAseq, we assembled ~77,000 transcripts, which were annotated using seven functional databases. We tested ~1600 genes, known to be under positive selection in anurans and reptiles adapted to high elevations, and potentially detected positive selection for 114 of these genes in Asymblepharus. Even though the strength of these results is limited due to the single-animal approach, our transcriptome resource may be valuable data for further studies on squamate reptile evolution in the Himalayas as a hotspot of biodiversity.


Introduction
The Himalayan arc represents one of the world's most important faunal and floral hotspots with high species diversity and endemism [1], which result from the Tertiary orogeny of this mountain chain, its complex topography as well as its great climatic heterogeneity and isolation. The genesis of the Tibetan highlands and the Himalayas since the Paleogene, with the Greater Himalayas starting to rise presumably the earliest in the post-Eocene (for a review, see the supplementary in Hofmann et al. [2]), triggered the evolution of unique biodiversity under gradual high-altitude adaptation, as already shown for anurans [3][4][5][6][7]. Besides amphibians, there are also several reptiles that can cope with life at high altitude in those regions, e.g., Thermophis [8], Phrynocephalus [9], and some Laudakia species [10]. Potential constraints to upslope migration of reptiles (and amphibians) to high-elevation environments are the substantial UV-radiation, the thermal extremes, and especially the oxidative stress, referred to as high-altitude hypoxia, which interacts with temperature in a context-dependent manner to influence thermal performance and limits in terrestrial ectotherms [11,12]. Recent advances in high-throughput sequencing technologies have led to a growing number of genomic studies that address the molecular basis of Figure 1. Map of Asymblepharus species based on GBIF (www.gbif.org; accessed on 20 July 2021) records of preserved specimens and georeferenced localities in the taxonomic reptile database (https://reptile-database.reptarium.cz/; accessed on 20 July 2021). The location of our RNA sample of the female A. ladacensis (photo) is indicated by a green circle with a dot in the middle and an arrow. *Note, according to a large-scale phylogeny of squamates, A. sikimmensis is nested within Scincella; however, it remains unclear whether this single "A. sikimmensis" specimen, on which the sequence data are based, had been taxonomically correctly identified. Therefore, we also show the GBIF records of specimens collected as A. sikimmensis. Records of A. eremchenkoi in the databases could not be georeferenced due to insufficient information on the collection site. Photo credit: Sylvia Hofmann.

Sample Collection and Ethics Statement
A single gravid female Asymblepharus ladacensis was collected in Central Nepal, in the Dhaulagiri range, west of the Kali Gandaki valley (28.68° N, 83.59° E; 2714 m a.s.l.; Figure  1). Samples were collected in accordance with regulations for the protection of terrestrial wild animals under the permits of the Nepal expeditions of the Natural History Museum of Erfurt, Germany [27,28]. All treatments were carried out in accordance with approved guidelines and according to the permit as well as the local animal welfare committee's instructions (VNME 17,(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30). Tissues were transferred into RNAlater (Thermo Fisher), kept at ambient temperature during the time of the fieldwork, and later stored at −30 °C.

RNA Isolation, Library Preparation, and Sequencing
We followed the same procedure as previously described [16]. In brief, total RNA was isolated from the brain, an embryonic disc, and from pooled tissues (including lung, muscle, and heart) using TRIzol Reagent (Thermo Fisher Scientific, Waltham, MA, USA) according to the supplier's recommendation, in combination with the RNeasy Mini Kit (Qiagen, Hilden, Germany) and adjusted to equal concentrations. RNA quality was assessed by RNA concentration, RIN (RNA Integrity Number) value, 28S/18S and fragment length distribution using an Agilent Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA).
Complementary DNA (cDNA) library preparation and paired-end sequencing were carried out by BGI (BGI-Hongkong Co., Ltd., Tai Po District, Hong Kong), using Illumina records of preserved specimens and georeferenced localities in the taxonomic reptile database (https://reptile-database.reptarium.cz/; accessed on 20 July 2021). The location of our RNA sample of the female A. ladacensis (photo) is indicated by a green circle with a dot in the middle and an arrow. * Note, according to a large-scale phylogeny of squamates, A. sikimmensis is nested within Scincella; however, it remains unclear whether this single "A. sikimmensis" specimen, on which the sequence data are based, had been taxonomically correctly identified. Therefore, we also show the GBIF records of specimens collected as A. sikimmensis. Records of A. eremchenkoi in the databases could not be georeferenced due to insufficient information on the collection site. Photo credit: Sylvia Hofmann.

Sample Collection and Ethics Statement
A single gravid female Asymblepharus ladacensis was collected in Central Nepal, in the Dhaulagiri range, west of the Kali Gandaki valley (28.68 • N, 83.59 • E; 2714 m a.s.l.; Figure 1). Samples were collected in accordance with regulations for the protection of terrestrial wild animals under the permits of the Nepal expeditions of the Natural History Museum of Erfurt, Germany [27,28]. All treatments were carried out in accordance with approved guidelines and according to the permit as well as the local animal welfare committee's instructions (VNME 17,(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30). Tissues were transferred into RNAlater (Thermo Fisher), kept at ambient temperature during the time of the fieldwork, and later stored at −30 • C.

RNA Isolation, Library Preparation, and Sequencing
We followed the same procedure as previously described [16]. In brief, total RNA was isolated from the brain, an embryonic disc, and from pooled tissues (including lung, muscle, and heart) using TRIzol Reagent (Thermo Fisher Scientific, Waltham, MA, USA) according to the supplier's recommendation, in combination with the RNeasy Mini Kit (Qiagen, Hilden, Germany) and adjusted to equal concentrations. RNA quality was assessed by RNA concentration, RIN (RNA Integrity Number) value, 28S/18S and fragment length distribution using an Agilent Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA).

Assembly and Assessment of Transcriptome Quality and Completeness
First, we controlled our data for rRNA quantity using SortMeRNA 4 [30]. The three de novo assemblies were then created following the Oyster River Protocol (ORP; Docker image 2.2.8) best practices [31]. This protocol implements both pre-assembly procedures and a number of different kmer lengths and assemblers, finally merging these assemblies into a single, comprehensive transcriptome. The rationale behind it is that assembling RNAseq reads with different assembly tools increases assembly quality and mapping rate and, in turn, the ability to draw conclusions from that fraction of the sample [31]. Thus, merging the contigs resulting from several assemblers and parameter configurations to combine the advantages of certain assembly mechanisms and to overcome their different disadvantages seems to be the best way to obtain a comprehensive de novo transcriptome assembly [32,33].
Illumina sequencing adapters and nucleotides with quality Phred ≤ 2 were removed using Trimmomatic v0.36 [34], then the reads were error corrected by Rcorrector version 1.0.2 [35]. These reads were then assembled using Trinity release 2.8.4 [36] with default settings (k = 25), two independent runs of SPAdes assembler version 3.11 with kmer lengths of 55 and 75 [37], and the assembler Shannon version 0.0.2 with a kmer length of 75 [38]. The resulting four distinct transcriptome assemblies were then merged to a single, comprehensive transcriptome using Ortho-Fuser [31]. This final transcriptome was evaluated with TransRate version 1.0.3 [39], which is modified for and packaged with the ORP, and with BUSCO (Benchmarking Universal Single-Copy Orthologs) version 4.1.4 [40,41]. Searching the assembly for conserved single-copy orthologs found in orthologous sets of genes constructed from genomes representing eukaryotes (70 species: 255 BUSCOs), vertebrates (67 species: 3354 BUSCOs), and tetrapods (38 species: 5310 BUSCOs).
A detailed quality assessment of the assembly with respect to known genes was further obtained with rnaQUAST version 2.2.0 [42] using the reference genomes of Anolis carolinensis, Gekko japonicus, and Python bivittatus.

Functional Annotation
Transcripts were functionally annotated as previously described [16]. Briefly, sequence homology searches were conducted against seven databases (Gene Onthology, GO; Kyoto Encyclopedia of Genes and Genomes, KEGG; EuKaryotic Orthologous Groups, KOG; InterPro; the non-redundant nucleotide database, NT; the non-redundant protein database, NR; SwissProt). To align our data to KEGG, KOG, NR, and NT and SwissProt, we used Diamond v0.8.31 [43] or the BLASTx [44] algorithm; matched transcripts were filtered by using a cut-off e-value of 1 × 10 −25 . Transcripts that aligned to the NR database were transferred to the GO database with Blast2GO v2.5.0 [45] and assigned into the following three groups: biological process, cellular components, and molecular functions. InterProScan5 v5.11-51.0 tool [46] was used to annotate against the InterPro databases. Blast v2.2.23 was used to search in SwissProt and hmmscan v3.0 [47] for search against Pfam database (for each sample individually). Candidate coding areas within the transcript sequences were predicted by TransDecoder v.3.0.1 [36]. For a coding sequence with multiple open reading frames (ORF), the longest one was selected. We also used the getorf program of the EMBOSS v6.5.70 package [48] to find the ORF of each transcript and mapped them to the Animal Transcription Factor DataBase (AnimalTFDB2.0). The threshold of transcript lengths used for annotation and downstream analyses was ≥200 bp.

Positively Selected Genes Related to Mechanisms of High-Altitude Adaptation
We selected transcripts of genes reported to exhibit molecular adaptation to high elevations in lizards and anurans as provided in the literature [7,15]. This comprised a total of 143 genes identified to be under positive selection (PSG) in high-elevational lineages of lizards (Phrynocephalus vlangalii) [15]. We included further 1481 PSGs of these lizards (P. erythrurus, P. putjatia, P. vlangalii) and of dicroglossid frogs (Nanorana, Quasipaa) [7], genes that were identified across an elevational gradient (~1000 m to 4500 m). These additional, individual genes were grouped according to the phylogenetic tree branches across the elevational gradient as presented in Sun et al. (2018) [7]: PSGs attributed to branches that represent (i) lowland species (~1000 m), (ii) species distributed in colline zones up to about 2000 m, (iii) submontane and montane species (2000 and 3500 m), and (iv) subalpine and alpine species (>3500 and 4500 m). Given the vertical distribution of Asymblepharus ladacensis in the Himalayas between~2500 m and 4500 m [49], we expected to find primarily genes under positive selection reported as PSGs for the submontane and montane, as well as the subalpine and alpine species.
To identify scincid orthologs to the genes under positive selection (high-altitude adaptations) in reptiles and amphibians, we used the corresponding coding sequences from Anolis carolinensis (AnoCar2.0, gene build from ensemble 104.2; we used Anolis gene numbers that were mentioned in the publications cited in the last paragraph) as a reference for blast searches. We performed blast searches with the newly sequenced and assembled transcriptome of A. ladacensis, as well as with publicly available transcriptome data from the following scincid lizards: Scincella lateralis (NCBI-SRA SRR629642), Lampropholis guichenoti (NCBI-SRA SRR4293354), Lepidothyris fernandi (NCBI-SRA SRR10360868), and Pseudemoia entrecasteauxii (NCBI-SRA SRR3099521). Only the best reciprocal hit between Anolis and each scincid species (and only if the e-value was below 10 −25 ) were used for subsequent analyses.
Alignments of orthologs to the Anolis reference were done with mafft v.7.455 [50] making use of the "adjustdirection" and "keeplength" options to get alignments that are in the same direction and keep reading frames intact. Alignments were inspected for good representation of all species under study. When one or more species had substantially shorter, incomplete contigs as best reciprocal hits, we omitted those from the alignment. A phylogenetic tree for each alignment was produced using FastTree v. 2.1.10 [51] with default settings, except for the nucleotide option.
Alignments and trees were analyzed with the HyPhy (Hypothesis Testing using Phylogenies) package v. 2.5 [52,53] using the following methods (for details see [54]): (i) The BUSTED (Branch-Site Unrestricted Statistical Test for Episodic Diversification) [55] model, to test whether a given gene has been subject to positive, diversifying selection at any site, at any time (we tested all lineages for positive selection); (ii) FUBAR (Fast, Unconstrained Bayesian AppRoximation) [56], a Bayesian approach to infer which site(s) in a gene are subject to pervasive, i.e., consistently across the entire phylogeny, diversifying selection (we considered a posterior probability of at least 0.90 as significant); and (iii) the branch-site model aBSREL (adaptive Branch-Site Random Effects Likelihood) [57,58], to test whether codon sites and individual branches are subject to positive selection across the phylogeny. The threshold for significance in BUSTED and aBSREL was set at a p-value lower than 0.05.
GO term and metabolic pathway enrichment analysis was done using pantherdb (pantherdb.org; accessed on 2 August 2021, [59]), using the genome of Anolis carolinensis as reference.

Sequencing and Transcriptome Assembly
All RNA was reasonably high quality; A260/280 ratios ranged between 1.81 (brain), 1.76 (embryonic disc), and 2.09 (pooled tissues); RIN values were between 7.90, 8.20, and 9.10, respectively. RNA-seq libraries of the two tissues yielded a total of 74.78, 73.56, and 64.20 million raw sequence reads (Table 1). Pre-processing of reads via read trimming and read error correction removed approximately 2-4% of the raw data, resulting in 73.14, 71.89, and 61.40 million clean reads for the brain, embryonic disc, and pooled tissues (Table 1). GC content of these clean reads was 48%. Table 1. Summary of sequencing data used to obtain the de novo transcriptome assemblies of Asymblepharus ladacensis based on paired-end Illumina sequencing. Final assemblies based on four unique assemblies per sample generated by ORP using different assemblers and k-mers.  Figure S1). A total of 48,884 (32.22%; brain), 44,358 (42.19%; embryonic disc), and 25,772 (38.64%; pooled tissues) transcripts were longer than 500 bp (Supplementary Table S2).

Assembly Completeness
TransRate's optimal assembly score (min 0.0, max 1.0) is considered to be a good parameter of the quality of an assembly [32]; it captures the confidence and completeness of the assembly. The TransRate scores of our final assemblies were high, ranging between 0.44 (optimized score 0.49) for the brain sample, 0.45 (optimized score 0.57) for the embryonic disc sample, and 0.44 (optimized score 0.52) for the pooled tissues. A transRate score >0.22 is generally thought to be acceptable [31,39]. More than 90% of the reads were used to assemble the transcriptomes, and 87% (brain, pooled tissues), as well as 90% (embryonic disc) of the fragments, were considered as good mappings, while only 1.6% (brain), 3.7% (embryonic disc), and 2.7% (pooled tissues) of assembled contigs had no coverage (Supplementary Table S1).
The assessment of completeness of our assemblies by the BUSCO pipeline resulted in a moderate to high percentage of complete eukaryotic orthologues (from 69.8 to 98.0% of 255 BUSCOs) but also a significantly higher percentage of putative paralogues; in the vertebrate (3354 BUSCOs) and tetrapod (5310 BUSCOs) databases, more than half of the markers were recovered completely in the brain and embryonic tissue ( Table 2). The fraction of missing BUSCOs ranged between 0.8% (Eukaryota; embryonic disc) and 54.9% (Tetrapoda; pooled tissues). These BUSCO values are comparable to recent de novo transcriptome studies in many vertebrates [60][61][62]. BUSCO recovery rate tends to be highest when full organism and/or multiple developmental stages were used to generate the transcriptomes, compared to those assembled from specific organs or tissues [63]. Coverage of specific gene databases (Anolis carolinensis, Gekko japonicus, Python bivittatus) ranged between 7% and 29%, being highest for the Gekko reference genome (Supplementary Table S2). Up to 40,000 transcripts (22-39%) could be aligned to one of the three databases. The duplication ratio varied between 1.3 and 1.5 and was in the range reported for vertebrate transcriptomes [32,64]. Although the proportion of misassembled contigs was low (<2%), the assemblies showed only a small number of 95%-assembled genes and isoforms (Supplementary Table S2). We assume that this might reflect biological novelty in the study species rather than fragmentation of the assemblies [41]. For example, the low proportion of completeness against the three reptile gene databases contrasts with a >40% completeness score for the vertebrate gene set and could be the result of an overrepresentation of gene sets from more intensively studied lineages [65].
Alternatively, suboptimal sample quality and a resulting higher proportion of fragmented genes due to potential RNA degradation could be a reason for the relatively low number of assembled genes and lower Eukaryota BUSCO scores, especially in the pooled tissues sample.

Positive Selection
Among the 143 PSGs previously reported for the high-elevation lineage of the toadheaded agama Phrynocephalus vlangalii [15], we found ten genes to be likewise under positive selection in Asymblepharus based on all three tests that we performed. Two of these genes (IL1RAP, GRK6) were specific to the Asymblepharus branch of the gene tree (Table 4  and Supplementary Table S3). Table 4. Summary of the positive selection analysis for high-altitude candidate genes of a toad-headed agama (Phrynocephalus vlangalii) [15] tested likewise positive in Asymblepharus (A) using BUSTED (B); p-value > 0.05) [55], FUBAR (FB; number of sites under positive selection) [56] and aBSREL (aB) [57,58] methods. PSGs are represented by the last six digits of the anole lizard's (Anolis carolinensis) ENSEMBL gene and transcript identifiers (starting with ENSACAG00000, or ENSACAT00000, respectively).

GeneID
Gene Out of the 410 transcripts reported to be under positive selection in lowland frog and lizard species [7], 321 could be identified in Asymblepharus. Of these, 23 (7.2%) transcripts were found to be under positive selection in the Ladakh Ground Skink in all three tests (Table 5). Similarly, a total of 32 (7.1% of 449 transcripts) PSGs of colline, 24 (7.6% of 314 transcripts) of submontane and montane, and 24 (6.9% of 350 transcripts) of subalpine and alpine frog and lizard species were tested to be under positive selection in Asymblepharus. Moreover, out of the 32 parallel PSGs that were identified in both high-elevation frog and lizard species [7], one gene (PGS1) showed positive selection in the Ladakh Ground Skink (Table 5). Analyzing these genes under positive selection for GO terms and pathways reveals an overrepresentation of genes involved in the following processes: low-density lipoprotein receptors and catabolic processes; mitochondrial citrate transmembrane transport; glycolysis and fructose/galactose metabolism; nucleoside phosphate binding; p53 pathway feedback loop (involved in DNA repair), platelet-derived growth (PDGF) factor binding (involved in blood-vessel formation). Table 5. Summary of the positive selection analysis for candidate genes of lineages of dicroglossid frogs and toad-headed agamas [7] identified across an elevational gradient, tested likewise positive in Asymblepharus (A) using BUSTED (B); p-value > 0.05) [55], FUBAR (FB; number of sites under positive selection) [56] and aBSREL (aB, number of branches with positive selection) [57,58] methods. Ensembl gene and transcript identifier (ENSACAG00000, ENSACAT00000) refers to Anolis carolinensis.

Lowland
Gene

Discussion
This study presents the first transcriptome sequences from different tissues of the Ladakh Ground Skink Asymblepharus ladacensis, a high-altitude reptile species endemic to the Greater Himalayas. We provide high-quality de novo transcript assemblies and wellannotated results, enabling comparisons with transcriptomes of related scincid or higher lizards available at public databases. Although squamate reptiles (lizards and snakes) represent one of the most diverse vertebrate groups with over 10,000 species spanning more than 200 million years of evolution [66], genomic data of squamates are limited and still poorly studied [67]. To our knowledge, no annotated transcriptome has been published for the genus Asymblepharus so far. Although our study is mostly descriptive, it has yielded discoveries with respect to genes known to play roles in the adaptation of vertebrates to high elevations and adds data resources for genomic studies in Himalayan herpetofauna.
Yang et al. [15] used comparative transcriptomic analyses of two closely related lizards, Phrynocephalus przewalskii from low elevations (500-1500 m a.s.l.) and P. vlangalii from high elevations (2000-4600 m a.s.l.), to identify candidate genes that are potentially linked to adaptation to high elevation environments. In addition, Sun et al. [7] tested amphibian and reptile populations at various altitudes in Tibet, which show parallel evolution. These studies demonstrated convergent and continuous adaptation to high elevations in Anura (Ranidae) and Sauropsida (Agamidae). Genes with related functions, especially DNA repair and energy metabolism, exhibit featured rapid changes and are positively correlated to elevation. These data let us assume that a similar genomic high-elevation selection syndrome might be detectable in Asymblepharus, sampled in 2714 m a.s.l. (Methods), and with a vertical distribution between~2500 and 5500 m ( [26] and references therein).
Indeed, we identified a total of 10 out of 143 and 104 out of~1500 key genes [7,15] under positive selection in the Asymblepharus transcriptome. Interestingly, several of the 10 genes (Table 4) have been reported to be under positive selection or significantly enriched or differentially methylated for pathways consistent with physiological compensation for limited oxygen in high elevation dwellers, e.g., IL1RAP [68] (human), MIA3 [69] (pika), and MICU1 [70] (Ladakhi cow). Several genes we identified to be under positive selection have GO terms that suggest their involvement in, e.g., energy metabolism and DNA repair. It is well-known that high UV radiation and hypoxia are major challenges for organisms in high-altitude habitats. The extreme environments in the Greater Himalayas necessitate high energy metabolism, strong resistance to UV by an efficient DNA repair, and adaptation to hypoxia in species endemic to these mountains. However, the proportion of genes we found in Asymblepharus per 'altitude-specific' gene group does not appear to reflect a convergently evolved gene set as previously reported [7]. In other words, we found a comparable number of genes under positive selection from the group of genes identified for lowland species and those identified for colline, montane, or alpine species. Given the idea of convergent evolutionary changes and, thus, a gradual accumulation of high-elevation genetic adaptations, a higher number of candidate genes reported for montane and alpine species would have been expected in Asymblepharus compared to those genes reported in species distributed in lower elevations. The potential reasons for this supposed lack of confirmation of the suggested pattern are complex. A major limitation is that our analyses were restricted so far to a single Asymblepharus transcriptome. Due to limitations of sampling, data were unobtainable from radiation of species or an altitudinal gradient as available for frogs and other lizards [7,15], preventing us from intraspecific comparisons. Moreover, data of a single specimen cannot reflect the breadth of allelic diversity in the selected genes, putatively associated with adaptations to high altitude, especially in the view of wide vertical distribution. Another deficiency results from the fact that only one female but no male could be sequenced. Indeed, many genes that might be sex-specifically expressed might not have been sequenced or characterized with our approach. Therefore, future research with multiple high-elevation species and populations across a larger scale of altitudinal variation should validate genes known to contribute to high elevation adaptation in scincid reptiles and, thus, yield additional evidence for potential convergent evolution.
A promising future genomic approach might be to include populations of Asymblepharus between the species' lower and upper distributional periphery, sampling three populations at each of four elevation levels (e.g., <2000 m, 3000 m, 4000 m, and 5000 m a.s.l.), to investigate the expression of genes presumably related to adaptions to high altitude. It would also be desirable to reveal potential tissue-specific expression patterns across altitude-associated genes, samples of organs sensitive to UV radiation, and for oxygen, e.g., heart, lung, skin, and embryonic structures, which are particularly of interest. Moreover, to address adaptive convergence, additional comparative transcriptomic analyses on Japalura or Laudakia species might be promising since these taxa have a similarly broad vertical distribution as Asymblepharus and often co-occur with sympatric ground skinks. Asymblepharus might even become an excellent model to study local adaptation by reciprocal transplantation experiments between high and low altitude populations. Such studies could enhance our understanding of how organisms might cope with rapid environmental changes in fluctuating demographic contexts. However, such intensive field studies require adequate access to suitable habitats in the Himalayas and, thus, a much higher logistic and financial effort than available in our pilot study.

Conclusions
In summary, our present study provides the first transcriptomic data for a Himalayan reptile of the genus Asymblepharus and evidence for genes under positive selection for highaltitude adaptation of the Ladakh Ground Skink. Further research is encouraged to validate the key genes confirmed in this study by population genetic and functional genomic approaches. Comparative sequencing analyses for other Asymblepharus species may enable further insights into the adaptive basis of reptiles to different altitude environments in the Himalayas. Our data are available for future investigations on the evolution and environmental adaptation in Himalayan high-altitude vertebrates.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/genes12091423/s1, Figure S1: Basic statistics for the data derived from brain, embryonic disc, and pooled tissues based on the rnaQUAST report, Table S1: TransRate results for Asymblepharus ladacensis ORP assemblies, Table S2: Results of the quality assessment of the transcriptome assemblies from brain tissue, an embryonic disc, and pooled tissues of Asymblepharus ladacensis using rnaQUAST and the reference database of (a) Anolis carolinensis, (b) Gekko japonicus, and (c) Python bivittatus, Table S3: List of positively selected genes in Phrynocephalus species (according to [2]) and Asymblepharus ladacensis, their functional categories and summary of tests for positive selection.  Institutional Review Board Statement: All treatments were carried out in accordance with approved guidelines and according to the permit as well as the local animal welfare committee's in-structions (VNME 17, 15-30).

Informed Consent Statement: Not applicable.
Data Availability Statement: Data generated in this study are publicly available from the NCBI Gen-Bank database under the Bioproject ID PRJNA750278, BioSamples SAMN20458631, SAMN20458632, and SAMN20458631. All sequence data were deposited in the NCBI Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/Traces/sra/; accessed on 2 August 2021) under the accession numbers SRR15283177, SRR15283178, and SRR15283179; assembled sequences were transmitted to NCBI Transcriptome Shotgun Assembly Sequence Database (TSA, http://www.ncbi.nlm.nih.gov/genbank/tsa (accessed on 2 August 2021).