布威海丽且姆·阿巴拜科日，单文娟，马合木提·哈力克，钟林强

（1.新疆维吾尔医学专科学校，和田，848000；2.新疆大学生命科学与技术学院，乌鲁木齐，830017）

物种适应性进化的遗传基础可以通过基因序列或有助于适应性表型的基因表达的变化来表现［1］。物种基因组序列的变异可能通过顺式或反式作用引起基因表达谱的变化，或改变蛋白质编码序列从而改变相应的适应性表型［2］。相对于基因组的变异，基因表达的变化更短暂，能够反映相应适应性表型的变化，可以增强表型的灵活性，因此基因的表达调控是一种特别有用的环境适应机制［3］。在物种环境适应性研究中，转录组学（RNA-Seq）可作为比较基因组学和群体遗传学的补充，能够有效地识别那些作为强烈自然选择目标的潜在特征的适应性基因［2］。利用转录组学对同一物种的不同组织或不同处理条件下的同一个物种或组织进行特异性差异表达分析，可以揭示物种环境适应性机制［3］。此外，转录组学只需对少数个体的血液或组织取样就能够识别与干旱—沙漠适应表型相关的基因［3］。因此，在物种环境适应性机制的研究中转录组学被广泛应用，如基于RNA-Seq 技术分析干旱环境适应性旗尾更格芦鼠（Dipodomys spectabilis）的肾脏和脾脏组织表达谱，揭示该物种肾脏组织中高表达的基因与肾脏滤液中的水重吸收作用有关［4］；对双峰驼（Camelus bactrianus）沙漠环境适应性研究发现，在禁水条件下饲养的双峰驼肾皮质和髓质中上调表达的基因参与调节钠和水分的重吸收、葡萄糖代谢和能量产生、渗透保护和糖酵解等重要的过程，而下调表达的基因主要参与肾髓质的渗透调节，进一步揭示双峰驼在高血糖水平下独特的渗透调节和水分保存的补偿机制［5］；对双峰驼结肠组织的转录组学分析，表明双峰驼在水资源限制的情况下通过降低RNA 和蛋白质的合成，从而减少结肠中的代谢速率［6］；利用RNA-Seq 和蛋白质组相结合的分析表明，在慢性脱水和随后的急性补水条件下，胆固醇水平的调节在单峰驼（Camelus dromedarius）缺水期间肾脏产生高浓度尿液的能力中起着关键作用［7］。

皮肤作为哺乳动物直接同外界环境接触的器官，与动物对环境和温度的适应能力息息相关。研究者在双峰驼、努比亚山羊（Capra nubiana）、乌干达本地山羊（C.hircus）和埃及肥尾羊（Ovis aries）等干旱—沙漠环境生活的有蹄类动物中已筛选出与皮肤保护策略相关的多种候选基因［5，8-11］。此外，齐昱等［12］利用RNA-Seq 技术对蒙古牛（Bos taurus）皮肤组织进行转录组学分析筛选出与该物种抗寒相关的信号通路及候选基因。

塔里木马鹿（Cervus elaphus yarkandensis）是新疆林区及荒漠环境中典型的大型有蹄类动物，是我国珍贵的鹿科（Cervidae）基因资源，2021 年发布的《国家重点保护野生动物名录》，将塔里木马鹿（仅限野外种群）由国家二级升级为国家一级重点保护野生动物。塔里木马鹿对塔里木盆地的极端干旱、炎热、强烈的紫外线辐射和各种应激环境已表现出良好的适应性，这意味着塔里木马鹿需要良好的表皮屏障系统来调节体温，以最大限度地减少皮肤水分流失以及抵御紫外线损伤。通过全基因组重测序的群体遗传学分析研究发现，在塔里木马鹿中受选择的TGM1、EDAR和PRDX6等基因可能与该物种干旱环境适应性进化的皮肤保护策略有关［13］。然而，这些基因在圈养与野生塔里木马鹿皮肤中的表达水平如何，二者皮肤组织中的差异表达基因是否与该物种干旱环境适应性相关等问题，值得更深入地探讨。因此，本研究对野生和圈养塔里木马鹿皮肤组织mRNA 进行测序和比较分析，筛选出野生塔里木马鹿皮肤组织中的差异表达候选基因，并对这些基因进行功能注释和富集分析，尝试寻找二者差异表达基因中是否有与塔里木马鹿干旱环境适应性相关的基因，期望通过转录组分析，揭示野生和圈养塔里木马鹿皮肤组织在基因表达方面的差异，了解其相关的生物学功能和代谢通路，从而为进一步分析和探讨功能基因及其相关适应性表型的试验奠定基础。

1 材料与方法

1.1 样本采集

2019 年10 月，在塔里木胡杨国家级自然保护区和新疆昌吉市盛华商贸有限责任公司养殖场，分别采集3 头成年雄性塔里木马鹿的小块皮肤组织，分别装到已标记的2 mL 冻存管内，并置于液氮管中速冻后带回实验室，超低温（-80 ℃）保存备用。样本地点信息见表1。

表1 塔里木马鹿样本采集地点主要气候数据Tab.1 Main climate data of Tarim red deer sample collection sites

1.2 RNA提取

利用TRIzol（Invitrogen，美国）法提取皮肤组织中的总RNA，用1%琼脂糖凝胶初步检测完整性，使用NanoDrop ND-2000 分光光度计（NanoDrop Lite，Thermo，美国）进行总RNA 的纯度和浓度检测，Qubit 2.0 对总RNA 浓度进行精确定量，Agilent 2100 生物分析仪（Agilent，美国）精确检测RNA 完整度（RNA integrity number，RIN）。一次建库要求：RNA的总量和浓度不得少于15 µL 和50 ng/µL，OD260/280应为1.8～2.2，OD260/230≥2.0，RIN值大于7。

1.3 建库、上机测序、组装和注释

利用Oligo（dT）磁珠富集法将总RNA 中的核糖体RNA 去除，富集mRNA。利用NEB Fragmentation Buffer随机打断mRNA，根据NEB普通建库方法构建cDNA文库。利用Illumina NovaSeq 6000平台边合成边测序（PE 文库，读长2×150 bp）。测序原始数据（raw data/reads）经过滤后得到高质量的测序数据（clean data/reads），并进行统计和质量评估。使用HISAT2软件（http：//ccb.jhu.edu/software/hisat2/index.shtml）将clean data/reads 与马鹿参考基因组CerEla1.0 red deer（Cervus elaphus hippelaphus，Gen-Bank accession GCA_002197005.1）比对，进而获得比对数据（mapped data/reads）。利用已有的参考基因组，使用软件StringTie（http：//ccb.jhu.edu/software/stringtie/）将mapped reads 组装拼接，并与已知的转录本比较，获取没有注释信息的转录本，对其中潜在的新转录本进行功能注释，发掘该物种的新转录本和新基因，从而补充和完善原有的基因组注释信息。

1.4 序列功能注释

获得基因/转录本后，将所有基因/转录本与NR、Swiss-Prot、Pfam、STRING、GO 和KEGG 数据库进行比对，全面获得基因和转录本的功能信息，统计注释情况。

1.5 表达量及表达量差异分析

利用表达定量软件RSEM v1.3.1［14］分别对基因和转录本的表达水平定量分析，FPKM 为定量指标，根据FPKM 分布情况对样本的表达量分布进行统计。用皮尔逊相关系数的平方（r2）表示样本间基因表达水平的相关性。将不同样品间获得的基因/转录本的读段数目（read counts）进行标准化处理，以圈养塔里木马鹿皮肤组织作为对照组，利用R 软件包DEGseq v1.38.0［15］进行组间差异表达分析，参数Padjust＜0.05，且|log2FoldChange|≥1 作为筛选差异基因的阈值［16］，对筛选出的差异表达基因可视化作图，推断差异基因的整体分布情况。

1.6 差异表达基因的功能分类及富集分析

采用软件Goatools v0.6.5［17］对差异表达基因进行GO 功能注释和富集分析。使用Fisher 精确检验，当P-adjust＜0.05 时，认为此GO功能存在显著富集情况。按照GO功能富集分析的计算原理，利用R脚本对基因集中的基因/转录本进行KEGG通路富集分析，当P-adjust＜0.05时，认为显著富集此通路的功能。

1.7 测序结果可靠性验证

为验证测序结果的可靠性，选择6 个差异表达基 因（PRDX2、PRDX6、HSPA4、TRAP1、HSPH1和K2C6A），其中PRDX6和TRAP1基因是在前期研究中通过全基因组重测序筛选出的可能与塔里木马鹿干旱环境适应性相关的受选择候选基因［13］。利用qRT-PCR 方法，以圈养塔里木马鹿皮肤组织作为对照组，以GAPDH作为内参基因进行基因表达量的检测，每个目的基因做3 次重复。具体基因信息见表2。以经检测合格并定量的总RNA 逆转录成cDNA（PrimeScript 1st Strand cDNA Synthesis Kit，大连Ta-KaRa），并将其作为模板，扩增目的基因。qRT-PCR反应体 系：2×SYBR real-time PCR premixture（BioTeke，中国）10.0 µL，10 µmol/L 的PCR 特异引物上游和下游各0.4 µL，cDNA 1.0 µL，Yellow 0.2 µL，加入RNase free dH2O 至20.0 µL。扩增程序：95 ℃预变性5 min；95 ℃变性15 s，各引物退火温度30 s，72 ℃延伸30 s，40 个循环；72 ℃延伸10 min，4 ℃保存。每个样品做平行3 管，以RNase free dH2O 为模板作为阴性对照。采用2-ΔΔCt分析方法计算目的基因的表达差异。

2 结果

2.1 总RNA质检结果

从野生和圈养塔里木马鹿皮肤组织中提取的总RNA 质量检测结果表明，所提取的总RNA无色素、蛋白和糖类等杂质污染，纯度（OD260/280为2.01～2.06，OD260/230为1.80～2.17）和浓度（201.50～393.00 ng/µL）较高。Agilent 2100 检测的RIN 值（8.1～9.3）符合文库构建和测序要求。

2.2 转录组测序与质量评估

2.2.1 测序数据统计与参考基因组比对

6个样本的转录组分析质控后共获得43.368 Gb测序数据，各样本的平均数据量为7.228 Gb，Q30在86.77%以上，GC 比例为50.98%～53.85%。表明低质量读段较少，测序数据质量较高（表3）。

表3 样本测序数据质量统计Tab.3 Sequencing quality data statistics of samples

利用HISAT2 软件将质控后的clean data/reads与参考基因组（CerEla1.0 red deer）进行比对，6 个样本的平均比对率为74.89%，总比对率为73.99%～76.28%，唯一比对率为71.58%～74.46%（表4）。

表4 塔里木马鹿皮肤组织与参考基因组的比对结果Tab.4 The mapping results of skin tissues of Tarim red deer to reference genome

2.2.2 读段在不同区域的分布

比对后的读段在参考基因组不同定位区域的分布情况表明，比对到的reads 主要分布在参考基因组的编码区和基因间区，表明研究用的参考基因组注释较为完整，分布在基因间区的reads 可能转录自新基因或新的非编码RNA（表5）。

表5 读段在基因组上的分布情况Tab.5 Distribution of reads on the genome

2.2.3 转录本覆盖均一性分布评估

转录本覆盖均一性分布评估用于检测转录本内测序reads 是否均一，并且是否存在5'区到3'区的偏差。分析过程中将所有已知转录本归一化为100 nt长度范围的区域，并计算覆盖在每个子区域位置上的reads 数目。所有转录本覆盖情况的叠加结果可以体现在reads 覆盖度分布图上（图1，以TRD_SK1为例）。由图1 可以看出，转录本没有明显的5'或3'偏向性，结果均一。

图1 TRD_SK1的reads覆盖度分布Fig.1 Coverage distribution of reads in TRD_SK1

2.2.4 转录本组装

根据拼接转录本与已知转录本叠加关系进行分类（即class code）统计，获得拼接转录本与已知转录本的位置关系，由图2 可知，新的转录本中，潜在的新转录本或转录片段最多，拼装出来的转录本包含于已知转录本中的转录本最少。

图2 转录本类型分布Fig.2 Pie chart of the distribution of transcript types

2.3 功能注释

鉴别出新基因/转录本后，将所有基因和已知基因、所有转录本和已知转录本分别与6大数据库（NR、Swiss-Prot、Pfam、STRING、GO 和KEGG）进行比对并注释，获得全部基因/转录本的功能信息，统计各数据库注释情况。结果显示，注释到NR数据库的基因/转录本数量最多，而注释到Pfam 数据库的基因数量和注释到STRING数据库的转录本数量最少（图3）。

图3 基因/转录本功能注释统计Fig.3 The summary statistics of functional annotation of genes/transcripts

2.4 样本表达量分布和相关性

从样本表达量的密度分布（图4A）可以看出，样本基因表达量主要分布在1.0～5.0，在一定程度上表明测序得到的基因表达情况较好，可以进行后续的差异表达基因分析。样本表达量间的r2值为0.65～0.93，相关性可用于后续分析（图4B），r2值越接近1，说明样品间表达模式的相似度就越高。

图4 样本中基因表达量密度分布（A）和相关性（B）Fig.4 Expression density distribution（A）and correlation（B）of expressed genes in samples

2.5 差异表达基因的筛选及其聚类分析

对筛选的差异表达基因（DEGs）进行火山图可视化作图，由图5 可以看出，相对于圈养个体，野生塔里木马鹿皮肤组织样本中表达差异的已知基因有3 303 个，其中上调表达基因1 535 个，下调表达基因1 768个。

图5 野生塔里木马鹿组织样本中差异表达基因Fig.5 Differentially expressed genes in tissue samples of wild Tarim red deer

将显著差异表达的mRNA 基因进行表达模式聚类分析（图6），发现圈养和野生塔里木马鹿的表达模式不同，可能是由于生存环境以及受到的压力不同导致。

图6 差异基因的表达模式聚类热图Fig.6 Clustering heat map of expression patterns of differential genes

2.6 差异表达基因的GO富集分析

将野生塔里木马鹿差异表达的基因比对GO 数据库，差异基因注释结果均被分为生物学过程（BP）、分子功能（MF）和细胞组分（CC）三大功能类别，其中前30 条生物学过程功能类别的富集统计如图7 所示。野生塔里木马鹿皮肤组织中的上调表达基因显著富集在基因表达、蛋白质折叠、蛋白质生物合成和代谢相关的生物学过程（BP，n=7 GO terms），翻译（BP，n=7 GO terms），核糖体的生成和大小亚基相关的生物过程及细胞组成（n=20 GO terms），角质化细胞分化及其细胞组成（n=7 GO terms），皮肤、表皮发育及其调控（BP，n=5 GO terms），皮肤屏障的形成，线粒体膜和呼吸链复合体组装相关的细胞组成及生物过程（n=20 GO terms），抗氧化以及氧化还原酶活性（n=5 GO terms）和抗菌肽生产的调控（BP，n=3 GO terms）等221 个GO 功能类别（P-adjust＜0.05，图7A）。下调表达基因显著富集于骨骼系统发育（BP，n=6 GO terms），骨化及其调控（BP，n=10 GO terms），软骨和成骨细胞分化及其调节（BP，n=6 GO terms），血管生成、发育和调控有关的生物过程（BP，n=30 GO terms），眼睛形态发生及视网膜血管形成和发育（BP，n=11 GO terms）等1 104 个GO 功能类别（P-adjust＜0.05，图7B）。

图7 上调表达（A）和下调表达（B）基因显著富集的前30条生物学过程气泡图Fig.7 Bubble chart of the first 30 biological processes with significant enrichment of up-regulated（A）and down-regulated（B）expression genes

2.7 差异表达基因的KEGG富集分析

为查找差异表达基因参与的代谢通路或信号转导过程，对差异表达基因进行KEGG 富集分析，发现野生塔里木马鹿中上调表达的基因显著富集在核糖体（ko03010）、真核生物的核糖体生物发生（ko03008）、剪接体（ko03040）、蛋白酶体（ko03050）、氧化磷酸化（ko00190）、RNA 转运（ko03013）、次级代谢产物的生物合成（ko01110）、酮体的合成与降解（ko00072）和帕金森病（ko05012）等12 个代谢通路上（P-adjust＜0.05，图8A）。下调表达基因主要涉及ECM 受体相 互作用（ko04512）、血小板活化（ko04611）、造血细胞谱系（ko04640）、补体与凝血级联反应（ko04610）、PI3K-Akt 信号通路（ko04151）、原发性免疫缺陷（ko05340）、脂肪细胞脂解的调节（ko04923）、蛋白质消化与吸收（ko04974）和cGMPPKG 信号通路（ko04022）等38 个代谢和信号通路（P-adjust＜0.05，图8B）。

2.8 qRT-PCR验证

为验证RNA-Seq 结果的可靠性，从野生塔里木马鹿皮肤组织上调表达基因中选择了与皮肤发育相关的基因（K2C6A），热应激相关基因（HSPA4、TRAP1和HSPH1）和氧化应激相关基因（PRDX2和PRDX6）进行qRT-PCR 检测。将qRT-PCR 检测到的目的基因的差异倍数与转录组测序结果进行比较，并用柱状图进行可视化（图9），可以看出，目的基因的差异倍数和转录组测序结果有一定的差异，但表达趋势一致，表明RNA-Seq测序结果较可靠，这些差异表达基因在塔里木马鹿干旱环境适应性中可能发挥着潜在的作用。

图9 qRT-PCR和RNA-Seq检测候选基因的相对表达量Fig.9 Comparison of candidate gene expression differences detected by qRT-PCR and RNA-Seq

3 讨论

极端干旱地区的主要特征是太阳辐射较强、极端温度、水和食物资源贫乏［18］。在这些地区，动物长时间受强烈的紫外线辐射，可能会增加皮肤的失水速度或引起皮肤损伤［19］。在双峰驼、埃及肥尾羊、努比亚山羊和乌干达本地山羊等干旱—沙漠环境中生活的有蹄类动物中已筛选出与皮肤保护策略相关的多种候选基因［5，8-11］，塔里木马鹿在塔里木盆地极端干旱环境中，可能也形成了皮肤保护策略来促进其适应性。除此之外，受长期家养驯化影响，圈养和野生塔里木马鹿在饮食和生存环境方面存在较大差异，这种差异使圈养和野生塔里木马鹿的皮肤在体温调节、水分调节和对身体的保护功能方面产生了一定的差异，故而推测其在基因表达方面应该也存在一定的差别，这些差异表达基因是否与该物种干旱环境适应性相关，笔者以此为切入点开展本研究。

3.1 对照组选择

本试验将栖息在塔里木胡杨自然保护区的野生塔里木马鹿作为研究对象，将来自昌吉市马鹿养殖场的圈养塔里木马鹿作为对照组进行RNA-Seq 和qRT-PCR 比较分析，检测目的基因表达量。采集的塔里木马鹿均来自干旱环境，但相对于圈养马鹿的栖息环境，野生塔里木马鹿栖息的塔里木胡杨自然保护区的年均气温较高（10.3 ℃）、年均降水量较低、年均蒸发量特别大，CGIAR-CSI 全球干旱度数据库显示，该地区的干旱指数（aridity index，AI）为0.015～0.023，属于超干旱地区（AI＜0.05），而昌吉市的年均气温较低（6.8 ℃），干旱指数为0.102～0.112，属于干旱地区（0.05＜AI＜0.20）［20］。野生塔里木马鹿世世代代栖息在塔里木盆地恶劣的自然环境中，不仅遭受自然环境中有害因素（如强烈UV 辐射和机械力）的损害，还受夏季干旱少雨、酷热，冬季严寒多风和水资源匮乏的影响［21］。此外，年均气温和年均降水量是影响该物种遗传多样性的主要环境因子［22］。皮肤作为哺乳动物直接同外界环境接触的器官，执行排汗、调节体温和保护身体的功能，与其感受到的外界环境刺激以及水合状态相关［23］。最重要的是，本研究中圈养的塔里木马鹿是从1983 年开始通过人工繁殖方式自繁自养［24］的马鹿后代，其被长期饲养在圈舍内凉棚下，有充足的水源供应，能够有效避免高温和UV辐射对身体的直接照射，从而能够减少皮肤蒸发的水分和强烈UV辐射对皮肤的损伤。以上栖息环境及饲养条件的差别对圈养和野生塔里木马鹿皮肤中相关基因的表达产生一定影响。而且，在生命活动的不同阶段、不同的细胞及不同的环境中，基因的表达均具有较大的差异，因此，转录组包含的信息也存在较大的差异。

3.2 塔里木马鹿干旱环境适应性相关基因筛选

本研究发现，在野生塔里木马鹿皮肤组织中的角蛋白（KRT222）、角蛋白相关蛋白（KRTAP11-1、KRTAP13-2和KRTAP14）、角质形成细胞分化相关蛋白（KRTDAP和KRTCAP3）、角蛋白Ⅰ型表皮（KRT3）、角蛋白Ⅱ型表皮（KRT83和KRT84）、角蛋白Ⅰ型细胞骨架（KRT10、KRT14、KRT17、KRT19、KRT26、KRT28和KRT39）、角蛋白Ⅱ型细胞骨架（KRT1、KRT5、K2C6A、KRT71、KRT72、KRT73、KRT74、KRT75和KRT75）、角蛋白Ⅰ型微纤 维（LOC122695079和LOC122695080）和角蛋白Ⅱ型微纤维（LOC102506783和LOC122681432）等31条基因显著上调表达（P＜0.05）。这些基因在调节毛囊和表皮发育、皮肤角质形成细胞分化和维持分层上皮中发挥重要作用。在机体内，角蛋白相互作用形成角蛋白细丝，组成稠密的网络，为皮肤及其他组织提供支撑，并保护这些组织免受外部摩擦和机械力损伤，减少皮肤损伤，加强保湿滋润，在维持皮肤屏障功能、防止水分流失和体温调节等方面具有重要作用［25］。有些角蛋白（如K2C6A）与其他角蛋白参与皮肤屏障、表皮伤口愈合以及启动和调节皮肤相关免疫反应的过程［26］。除了上述角蛋白及其相关基因外，多种跨 膜蛋白基因（TMEM40、TMEM54、TMEM62、TMEM79、TMEM80、TMEM184A、TMEM254和TMEM256）在野生塔里木马鹿皮肤组织中上调表达，其中TMEM79基因参与多种过程，包括上皮细胞成熟、皮肤屏障的建立和表皮发育的正调控，并在表皮发育、毛囊形态发生以及在角质化的上游或内部起作用［27］。这些角蛋白和跨膜蛋白显著富集于角质化（GO：0031424）、角质形 成细胞分化（GO：0030216）、角化（GO：0070268）、表皮细胞分化（GO：0009913）、皮肤发育（GO：0043588）、表皮发育（GO：0008544）、角蛋白丝（GO：0045095）、表皮lamellar body（GO：0097209）、角质化包膜（GO：0001533）、表皮细胞分化的调控（GO：0045604）、表皮发育的调控（GO：0045682）和皮肤屏障的建立（GO：0061436）等（P-adjust＜0.05）。富集不显著的功能类别中，上调表达基因涉及皮肤表皮发育、伤口愈合、表皮细胞扩散、表皮生长因子受体信号通路及其调控等，与皮肤发育和表皮细胞发育相关的多种生物学过程和分子功能（n=20 GO terms）相关。更重要的是，没有一个角蛋白及其相关蛋白质基因在野生塔里木马鹿皮肤组织中呈现下调表达。以上基因及其参与的生物学功能说明皮肤保护屏障在野生塔里木马鹿极端干旱环境适应过程中可能发挥着重要作用。

本研究中的上调表达基因NUDT5和COL17A1中，COL17A1编码XVI 型胶原的α链，该链是连接基底上皮细胞和基质半桥粒的关键结构成分［28］，该基因突变或其蛋白质的水解能引起皮肤萎缩、脱发和毛囊干细胞的老化［29］。而NUDT5可防止蛋白质的非酶性核糖基化，从核苷酸中去除被氧化的核苷，降低细胞在受氧化应激刺激下出现错误的DNA 复制率，保护细胞免受氧化应激损伤［30-31］。本研究结果在转录水平上初步支持了NUDT5和COL17A1基因在塔里木马鹿干旱环境适应性的潜在作用。

塔里木马鹿在塔里木盆地的极端干旱栖息环境中，遭受UV 辐射、温度和盐度等多种胁迫因素的影响。这些因素会导致塔里木马鹿机体内的活性氧（reactive oxygen species，ROS）大量积累，当ROS的形成超过靶细胞的抗氧化防御能力时，氧化应激会导致氧化损伤。皮肤作为机体第一保护屏障，要保护机体免受这些有害因素引起的氧化损伤。在本研究中，野生塔里木马鹿皮肤组织中抗氧化活性基因（PRDX2、PRDX6和GPX2）、转谷氨酰胺酶编码基因（TGM1和TGM3）、ATP 结合盒亚家族A成员12（ABCA12）和DNA 损伤修复基因（NSMF、HOXC12、GADD45GIP1、GADD45B和GADD45）呈现出明显的上调表达（P＜0.05），其中有些基因在双峰驼（TGM1）［5］、努比亚山羊（ABCA12）［8］、埃及肥尾羊（TGM3）［9］和乌干达本地山羊（HOXC12）［10］等物种中受正选择，并参与这些物种在干旱—炎热及沙漠环境中适应性皮肤保护策略的形成，以保护细胞免受沙漠中强烈的UV辐射损害，并参与炎热环境的体温调节，如HOXC12基因通过调节角蛋白分化特异基因，在毛囊分化、生长发育中发挥作用［11］。TGM3广泛表达于皮肤细胞，特别是角质形成细胞和角朊细胞。在角质形成细胞分化过程中，TGM3交联结构蛋白和脂质形成角质化细胞膜，为表皮提供抵御有害环境刺激（如UV 辐射）的屏障功能［32］。TGM1编码的蛋白质在角质化细胞膜的组装中具有完整的功能，形成稳定的不溶性的保护屏障以防止太阳辐射、水分流失和病原体的侵入，并参与构成表皮屏障功能的修复机制，从而对皮肤起保护作用［33］。ABCA12的脂质转运是角质形成细胞分化和表皮形态发生所必需的，因此在表皮脂质屏障形成和角质形成细胞分化中起重要作用［34］。NSMF基因在DNA 损伤增加的应激条件下，促进复制应激诱导的DNA 损伤反应以维持基因组的完整性［35］。而抗氧化活性基因（PRDX2、PRDX6和GPX2）可以保护角质化细胞免受ROS诱导的凋亡，促进伤口愈合，减轻氧化应激对皮肤造成的损伤［36］。此外，基于全基因组重测序的群体遗传学分析，发现TGM1、PRDX6和NSMF等基因在塔里木马鹿中受选择［13］。GADD45基因能发挥修复DNA 和调控细胞周期等作用［37］。这些DNA 损伤修复相关基因在保持基因组的完整性［9］和维持细胞正常功能方面具有至关重要的作用。因此，本研究结果在转录水平上初步支持了在野生塔里木马鹿皮肤组织中的上调表达基因可能通过参与氧化应激反应，形成皮肤保护屏障和维持皮肤稳态，发挥保护机体的作用。

在漫长的进化过程中，塔里木马鹿不仅对塔里木盆地夏季炎热表现出独特的适应性，还可以适应较低温的环境［38］。细胞对热/冷应激的耐受是由热休克蛋白（HSPs）家族介导的。热休克蛋白是一类高度保守的蛋白质，在逆环境中组成性表达和/或诱导表达，它具有分子伴侣功能，可调节信号转导通路，参与免疫应答和细胞周期，并具有抗氧化应激作用。另一方面，高盐摄入增加了热休克蛋白的表达，从而保护细胞免受高尿素的损伤［39］。本研究发现，热休克蛋白基因（HSPA4、HSPA14、TRAP1、HSPH1、HSPD1、HSP90和HSP110）与热休克蛋白家族具有序列和功能同源性的J 结构域蛋白基因（DNAJC2、DNAJC7、DNAJC19、DNAJC15和DNAJC21）在野生塔里木马鹿皮肤组织中上调表达，显著富集于蛋白 质折叠、蛋白质生物合成和代谢等生物学过程（P-adjust＜0.05）。这些热休克蛋白基因能够帮助蛋白质正确折叠，促进将错误折叠或变性的蛋白质清除，减轻应激对生物体的损害。其中，J 结构域蛋白是HSP40 家族的成员之一［40］，在细胞遭受干旱、高温、低温和强光等条件下能诱导表达，并应答胁迫应激，从而调节热休克蛋白的功能，发挥分子伴侣的作用，参与蛋白质折叠［41］。检测炎热干旱环境中的乌干达本地山羊品种，发现该物种的热休克蛋白与其体温调节适应性相关［10-11］。因此，笔者推测在野生塔里木马鹿皮肤组织中，上调表达的热休克蛋白基因以及J 结构域蛋白基因，在塔里木马鹿干旱环境适应性中的应激耐受性和体温调节方面可能具有重要的作用。