基于RNA-Seq技術(shù)的野生和圈養(yǎng)塔里木馬鹿皮膚組織轉(zhuǎn)錄組比較分析

布威海麗且姆·阿巴拜科日，單文娟，馬合木提·哈力克，鐘林強(qiáng)

（1.新疆維吾爾醫(yī)學(xué)專科學(xué)校，和田，848000；2.新疆大學(xué)生命科學(xué)與技術(shù)學(xué)院，烏魯木齊，830017）

物種適應(yīng)性進(jìn)化的遺傳基礎(chǔ)可以通過基因序列或有助于適應(yīng)性表型的基因表達(dá)的變化來表現(xiàn)［1］。物種基因組序列的變異可能通過順式或反式作用引起基因表達(dá)譜的變化，或改變蛋白質(zhì)編碼序列從而改變相應(yīng)的適應(yīng)性表型［2］。相對(duì)于基因組的變異，基因表達(dá)的變化更短暫，能夠反映相應(yīng)適應(yīng)性表型的變化，可以增強(qiáng)表型的靈活性，因此基因的表達(dá)調(diào)控是一種特別有用的環(huán)境適應(yīng)機(jī)制［3］。在物種環(huán)境適應(yīng)性研究中，轉(zhuǎn)錄組學(xué)（RNA-Seq）可作為比較基因組學(xué)和群體遺傳學(xué)的補(bǔ)充，能夠有效地識(shí)別那些作為強(qiáng)烈自然選擇目標(biāo)的潛在特征的適應(yīng)性基因［2］。利用轉(zhuǎn)錄組學(xué)對(duì)同一物種的不同組織或不同處理?xiàng)l件下的同一個(gè)物種或組織進(jìn)行特異性差異表達(dá)分析，可以揭示物種環(huán)境適應(yīng)性機(jī)制［3］。此外，轉(zhuǎn)錄組學(xué)只需對(duì)少數(shù)個(gè)體的血液或組織取樣就能夠識(shí)別與干旱—沙漠適應(yīng)表型相關(guān)的基因［3］。因此，在物種環(huán)境適應(yīng)性機(jī)制的研究中轉(zhuǎn)錄組學(xué)被廣泛應(yīng)用，如基于RNA-Seq 技術(shù)分析干旱環(huán)境適應(yīng)性旗尾更格蘆鼠（Dipodomys spectabilis）的腎臟和脾臟組織表達(dá)譜，揭示該物種腎臟組織中高表達(dá)的基因與腎臟濾液中的水重吸收作用有關(guān)［4］；對(duì)雙峰駝（Camelus bactrianus）沙漠環(huán)境適應(yīng)性研究發(fā)現(xiàn)，在禁水條件下飼養(yǎng)的雙峰駝腎皮質(zhì)和髓質(zhì)中上調(diào)表達(dá)的基因參與調(diào)節(jié)鈉和水分的重吸收、葡萄糖代謝和能量產(chǎn)生、滲透保護(hù)和糖酵解等重要的過程，而下調(diào)表達(dá)的基因主要參與腎髓質(zhì)的滲透調(diào)節(jié)，進(jìn)一步揭示雙峰駝在高血糖水平下獨(dú)特的滲透調(diào)節(jié)和水分保存的補(bǔ)償機(jī)制［5］；對(duì)雙峰駝結(jié)腸組織的轉(zhuǎn)錄組學(xué)分析，表明雙峰駝在水資源限制的情況下通過降低RNA 和蛋白質(zhì)的合成，從而減少結(jié)腸中的代謝速率［6］；利用RNA-Seq 和蛋白質(zhì)組相結(jié)合的分析表明，在慢性脫水和隨后的急性補(bǔ)水條件下，膽固醇水平的調(diào)節(jié)在單峰駝（Camelus dromedarius）缺水期間腎臟產(chǎn)生高濃度尿液的能力中起著關(guān)鍵作用［7］。

皮膚作為哺乳動(dòng)物直接同外界環(huán)境接觸的器官，與動(dòng)物對(duì)環(huán)境和溫度的適應(yīng)能力息息相關(guān)。研究者在雙峰駝、努比亞山羊（Capra nubiana）、烏干達(dá)本地山羊（C.hircus）和埃及肥尾羊（Ovis aries）等干旱—沙漠環(huán)境生活的有蹄類動(dòng)物中已篩選出與皮膚保護(hù)策略相關(guān)的多種候選基因［5，8-11］。此外，齊昱等［12］利用RNA-Seq 技術(shù)對(duì)蒙古牛（Bos taurus）皮膚組織進(jìn)行轉(zhuǎn)錄組學(xué)分析篩選出與該物種抗寒相關(guān)的信號(hào)通路及候選基因。

塔里木馬鹿（Cervus elaphus yarkandensis）是新疆林區(qū)及荒漠環(huán)境中典型的大型有蹄類動(dòng)物，是我國珍貴的鹿科（Cervidae）基因資源，2021 年發(fā)布的《國家重點(diǎn)保護(hù)野生動(dòng)物名錄》，將塔里木馬鹿（僅限野外種群）由國家二級(jí)升級(jí)為國家一級(jí)重點(diǎn)保護(hù)野生動(dòng)物。塔里木馬鹿對(duì)塔里木盆地的極端干旱、炎熱、強(qiáng)烈的紫外線輻射和各種應(yīng)激環(huán)境已表現(xiàn)出良好的適應(yīng)性，這意味著塔里木馬鹿需要良好的表皮屏障系統(tǒng)來調(diào)節(jié)體溫，以最大限度地減少皮膚水分流失以及抵御紫外線損傷。通過全基因組重測序的群體遺傳學(xué)分析研究發(fā)現(xiàn)，在塔里木馬鹿中受選擇的TGM1、EDAR和PRDX6等基因可能與該物種干旱環(huán)境適應(yīng)性進(jìn)化的皮膚保護(hù)策略有關(guān)［13］。然而，這些基因在圈養(yǎng)與野生塔里木馬鹿皮膚中的表達(dá)水平如何，二者皮膚組織中的差異表達(dá)基因是否與該物種干旱環(huán)境適應(yīng)性相關(guān)等問題，值得更深入地探討。因此，本研究對(duì)野生和圈養(yǎng)塔里木馬鹿皮膚組織mRNA 進(jìn)行測序和比較分析，篩選出野生塔里木馬鹿皮膚組織中的差異表達(dá)候選基因，并對(duì)這些基因進(jìn)行功能注釋和富集分析，嘗試尋找二者差異表達(dá)基因中是否有與塔里木馬鹿干旱環(huán)境適應(yīng)性相關(guān)的基因，期望通過轉(zhuǎn)錄組分析，揭示野生和圈養(yǎng)塔里木馬鹿皮膚組織在基因表達(dá)方面的差異，了解其相關(guān)的生物學(xué)功能和代謝通路，從而為進(jìn)一步分析和探討功能基因及其相關(guān)適應(yīng)性表型的試驗(yàn)奠定基礎(chǔ)。

1 材料與方法

1.1 樣本采集

2019 年10 月，在塔里木胡楊國家級(jí)自然保護(hù)區(qū)和新疆昌吉市盛華商貿(mào)有限責(zé)任公司養(yǎng)殖場，分別采集3 頭成年雄性塔里木馬鹿的小塊皮膚組織，分別裝到已標(biāo)記的2 mL 凍存管內(nèi)，并置于液氮管中速凍后帶回實(shí)驗(yàn)室，超低溫（-80 ℃）保存?zhèn)溆谩颖镜攸c(diǎn)信息見表1。

表1 塔里木馬鹿樣本采集地點(diǎn)主要?dú)夂驍?shù)據(jù)Tab.1 Main climate data of Tarim red deer sample collection sites

1.2 RNA提取

利用TRIzol（Invitrogen，美國）法提取皮膚組織中的總RNA，用1%瓊脂糖凝膠初步檢測完整性，使用NanoDrop ND-2000 分光光度計(jì)（NanoDrop Lite，Thermo，美國）進(jìn)行總RNA 的純度和濃度檢測，Qubit 2.0 對(duì)總RNA 濃度進(jìn)行精確定量，Agilent 2100 生物分析儀（Agilent，美國）精確檢測RNA 完整度（RNA integrity number，RIN）。一次建庫要求：RNA的總量和濃度不得少于15 μL 和50 ng/μL，OD260/280應(yīng)為1.8～2.2，OD260/230≥2.0，RIN值大于7。

1.3 建庫、上機(jī)測序、組裝和注釋

利用Oligo（dT）磁珠富集法將總RNA 中的核糖體RNA 去除，富集mRNA。利用NEB Fragmentation Buffer隨機(jī)打斷mRNA，根據(jù)NEB普通建庫方法構(gòu)建cDNA文庫。利用Illumina NovaSeq 6000平臺(tái)邊合成邊測序（PE 文庫，讀長2×150 bp）。測序原始數(shù)據(jù)（raw data/reads）經(jīng)過濾后得到高質(zhì)量的測序數(shù)據(jù)（clean data/reads），并進(jìn)行統(tǒng)計(jì)和質(zhì)量評(píng)估。使用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）比對(duì)，進(jìn)而獲得比對(duì)數(shù)據(jù)（mapped data/reads）。利用已有的參考基因組，使用軟件StringTie（http：//ccb.jhu.edu/software/stringtie/）將mapped reads 組裝拼接，并與已知的轉(zhuǎn)錄本比較，獲取沒有注釋信息的轉(zhuǎn)錄本，對(duì)其中潛在的新轉(zhuǎn)錄本進(jìn)行功能注釋，發(fā)掘該物種的新轉(zhuǎn)錄本和新基因，從而補(bǔ)充和完善原有的基因組注釋信息。

1.4 序列功能注釋

獲得基因/轉(zhuǎn)錄本后，將所有基因/轉(zhuǎn)錄本與NR、Swiss-Prot、Pfam、STRING、GO 和KEGG 數(shù)據(jù)庫進(jìn)行比對(duì)，全面獲得基因和轉(zhuǎn)錄本的功能信息，統(tǒng)計(jì)注釋情況。

1.5 表達(dá)量及表達(dá)量差異分析

利用表達(dá)定量軟件RSEM v1.3.1［14］分別對(duì)基因和轉(zhuǎn)錄本的表達(dá)水平定量分析，F(xiàn)PKM 為定量指標(biāo)，根據(jù)FPKM 分布情況對(duì)樣本的表達(dá)量分布進(jìn)行統(tǒng)計(jì)。用皮爾遜相關(guān)系數(shù)的平方（r2）表示樣本間基因表達(dá)水平的相關(guān)性。將不同樣品間獲得的基因/轉(zhuǎn)錄本的讀段數(shù)目（read counts）進(jìn)行標(biāo)準(zhǔn)化處理，以圈養(yǎng)塔里木馬鹿皮膚組織作為對(duì)照組，利用R 軟件包DEGseq v1.38.0［15］進(jìn)行組間差異表達(dá)分析，參數(shù)Padjust＜0.05，且|log2FoldChange|≥1 作為篩選差異基因的閾值［16］，對(duì)篩選出的差異表達(dá)基因可視化作圖，推斷差異基因的整體分布情況。

1.6 差異表達(dá)基因的功能分類及富集分析

采用軟件Goatools v0.6.5［17］對(duì)差異表達(dá)基因進(jìn)行GO 功能注釋和富集分析。使用Fisher 精確檢驗(yàn)，當(dāng)P-adjust＜0.05 時(shí)，認(rèn)為此GO功能存在顯著富集情況。按照GO功能富集分析的計(jì)算原理，利用R腳本對(duì)基因集中的基因/轉(zhuǎn)錄本進(jìn)行KEGG通路富集分析，當(dāng)P-adjust＜0.05時(shí)，認(rèn)為顯著富集此通路的功能。

1.7 測序結(jié)果可靠性驗(yàn)證

為驗(yàn)證測序結(jié)果的可靠性，選擇6 個(gè)差異表達(dá)基 因（PRDX2、PRDX6、HSPA4、TRAP1、HSPH1和K2C6A），其中PRDX6和TRAP1基因是在前期研究中通過全基因組重測序篩選出的可能與塔里木馬鹿干旱環(huán)境適應(yīng)性相關(guān)的受選擇候選基因［13］。利用qRT-PCR 方法，以圈養(yǎng)塔里木馬鹿皮膚組織作為對(duì)照組，以GAPDH作為內(nèi)參基因進(jìn)行基因表達(dá)量的檢測，每個(gè)目的基因做3 次重復(fù)。具體基因信息見表2。以經(jīng)檢測合格并定量的總RNA 逆轉(zhuǎn)錄成cDNA（PrimeScript 1st Strand cDNA Synthesis Kit，大連Ta-KaRa），并將其作為模板，擴(kuò)增目的基因。qRT-PCR反應(yīng)體 系：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。擴(kuò)增程序：95 ℃預(yù)變性5 min；95 ℃變性15 s，各引物退火溫度30 s，72 ℃延伸30 s，40 個(gè)循環(huán)；72 ℃延伸10 min，4 ℃保存。每個(gè)樣品做平行3 管，以RNase free dH2O 為模板作為陰性對(duì)照。采用2-ΔΔCt分析方法計(jì)算目的基因的表達(dá)差異。

2 結(jié)果

2.1 總RNA質(zhì)檢結(jié)果

從野生和圈養(yǎng)塔里木馬鹿皮膚組織中提取的總RNA 質(zhì)量檢測結(jié)果表明，所提取的總RNA無色素、蛋白和糖類等雜質(zhì)污染，純度（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）符合文庫構(gòu)建和測序要求。

2.2 轉(zhuǎn)錄組測序與質(zhì)量評(píng)估

2.2.1 測序數(shù)據(jù)統(tǒng)計(jì)與參考基因組比對(duì)

6個(gè)樣本的轉(zhuǎn)錄組分析質(zhì)控后共獲得43.368 Gb測序數(shù)據(jù)，各樣本的平均數(shù)據(jù)量為7.228 Gb，Q30在86.77%以上，GC 比例為50.98%～53.85%。表明低質(zhì)量讀段較少，測序數(shù)據(jù)質(zhì)量較高（表3）。

表3 樣本測序數(shù)據(jù)質(zhì)量統(tǒng)計(jì)Tab.3 Sequencing quality data statistics of samples

利用HISAT2 軟件將質(zhì)控后的clean data/reads與參考基因組（CerEla1.0 red deer）進(jìn)行比對(duì)，6 個(gè)樣本的平均比對(duì)率為74.89%，總比對(duì)率為73.99%～76.28%，唯一比對(duì)率為71.58%～74.46%（表4）。

表4 塔里木馬鹿皮膚組織與參考基因組的比對(duì)結(jié)果Tab.4 The mapping results of skin tissues of Tarim red deer to reference genome

2.2.2 讀段在不同區(qū)域的分布

比對(duì)后的讀段在參考基因組不同定位區(qū)域的分布情況表明，比對(duì)到的reads 主要分布在參考基因組的編碼區(qū)和基因間區(qū)，表明研究用的參考基因組注釋較為完整，分布在基因間區(qū)的reads 可能轉(zhuǎn)錄自新基因或新的非編碼RNA（表5）。

表5 讀段在基因組上的分布情況Tab.5 Distribution of reads on the genome

2.2.3 轉(zhuǎn)錄本覆蓋均一性分布評(píng)估

轉(zhuǎn)錄本覆蓋均一性分布評(píng)估用于檢測轉(zhuǎn)錄本內(nèi)測序reads 是否均一，并且是否存在5'區(qū)到3'區(qū)的偏差。分析過程中將所有已知轉(zhuǎn)錄本歸一化為100 nt長度范圍的區(qū)域，并計(jì)算覆蓋在每個(gè)子區(qū)域位置上的reads 數(shù)目。所有轉(zhuǎn)錄本覆蓋情況的疊加結(jié)果可以體現(xiàn)在reads 覆蓋度分布圖上（圖1，以TRD_SK1為例）。由圖1 可以看出，轉(zhuǎn)錄本沒有明顯的5'或3'偏向性，結(jié)果均一。

圖1 TRD_SK1的reads覆蓋度分布Fig.1 Coverage distribution of reads in TRD_SK1

2.2.4 轉(zhuǎn)錄本組裝

根據(jù)拼接轉(zhuǎn)錄本與已知轉(zhuǎn)錄本疊加關(guān)系進(jìn)行分類（即class code）統(tǒng)計(jì)，獲得拼接轉(zhuǎn)錄本與已知轉(zhuǎn)錄本的位置關(guān)系，由圖2 可知，新的轉(zhuǎn)錄本中，潛在的新轉(zhuǎn)錄本或轉(zhuǎn)錄片段最多，拼裝出來的轉(zhuǎn)錄本包含于已知轉(zhuǎn)錄本中的轉(zhuǎn)錄本最少。

圖2 轉(zhuǎn)錄本類型分布Fig.2 Pie chart of the distribution of transcript types

2.3 功能注釋

鑒別出新基因/轉(zhuǎn)錄本后，將所有基因和已知基因、所有轉(zhuǎn)錄本和已知轉(zhuǎn)錄本分別與6大數(shù)據(jù)庫（NR、Swiss-Prot、Pfam、STRING、GO 和KEGG）進(jìn)行比對(duì)并注釋，獲得全部基因/轉(zhuǎn)錄本的功能信息，統(tǒng)計(jì)各數(shù)據(jù)庫注釋情況。結(jié)果顯示，注釋到NR數(shù)據(jù)庫的基因/轉(zhuǎn)錄本數(shù)量最多，而注釋到Pfam 數(shù)據(jù)庫的基因數(shù)量和注釋到STRING數(shù)據(jù)庫的轉(zhuǎn)錄本數(shù)量最少（圖3）。

圖3 基因/轉(zhuǎn)錄本功能注釋統(tǒng)計(jì)Fig.3 The summary statistics of functional annotation of genes/transcripts

2.4 樣本表達(dá)量分布和相關(guān)性

從樣本表達(dá)量的密度分布（圖4A）可以看出，樣本基因表達(dá)量主要分布在1.0～5.0，在一定程度上表明測序得到的基因表達(dá)情況較好，可以進(jìn)行后續(xù)的差異表達(dá)基因分析。樣本表達(dá)量間的r2值為0.65～0.93，相關(guān)性可用于后續(xù)分析（圖4B），r2值越接近1，說明樣品間表達(dá)模式的相似度就越高。

圖4 樣本中基因表達(dá)量密度分布（A）和相關(guān)性（B）Fig.4 Expression density distribution（A）and correlation（B）of expressed genes in samples

2.5 差異表達(dá)基因的篩選及其聚類分析

對(duì)篩選的差異表達(dá)基因（DEGs）進(jìn)行火山圖可視化作圖，由圖5 可以看出，相對(duì)于圈養(yǎng)個(gè)體，野生塔里木馬鹿皮膚組織樣本中表達(dá)差異的已知基因有3 303 個(gè)，其中上調(diào)表達(dá)基因1 535 個(gè)，下調(diào)表達(dá)基因1 768個(gè)。

圖5 野生塔里木馬鹿組織樣本中差異表達(dá)基因Fig.5 Differentially expressed genes in tissue samples of wild Tarim red deer

將顯著差異表達(dá)的mRNA 基因進(jìn)行表達(dá)模式聚類分析（圖6），發(fā)現(xiàn)圈養(yǎng)和野生塔里木馬鹿的表達(dá)模式不同，可能是由于生存環(huán)境以及受到的壓力不同導(dǎo)致。

圖6 差異基因的表達(dá)模式聚類熱圖Fig.6 Clustering heat map of expression patterns of differential genes

2.6 差異表達(dá)基因的GO富集分析

將野生塔里木馬鹿差異表達(dá)的基因比對(duì)GO 數(shù)據(jù)庫，差異基因注釋結(jié)果均被分為生物學(xué)過程（BP）、分子功能（MF）和細(xì)胞組分（CC）三大功能類別，其中前30 條生物學(xué)過程功能類別的富集統(tǒng)計(jì)如圖7 所示。野生塔里木馬鹿皮膚組織中的上調(diào)表達(dá)基因顯著富集在基因表達(dá)、蛋白質(zhì)折疊、蛋白質(zhì)生物合成和代謝相關(guān)的生物學(xué)過程（BP，n=7 GO terms），翻譯（BP，n=7 GO terms），核糖體的生成和大小亞基相關(guān)的生物過程及細(xì)胞組成（n=20 GO terms），角質(zhì)化細(xì)胞分化及其細(xì)胞組成（n=7 GO terms），皮膚、表皮發(fā)育及其調(diào)控（BP，n=5 GO terms），皮膚屏障的形成，線粒體膜和呼吸鏈復(fù)合體組裝相關(guān)的細(xì)胞組成及生物過程（n=20 GO terms），抗氧化以及氧化還原酶活性（n=5 GO terms）和抗菌肽生產(chǎn)的調(diào)控（BP，n=3 GO terms）等221 個(gè)GO 功能類別（P-adjust＜0.05，圖7A）。下調(diào)表達(dá)基因顯著富集于骨骼系統(tǒng)發(fā)育（BP，n=6 GO terms），骨化及其調(diào)控（BP，n=10 GO terms），軟骨和成骨細(xì)胞分化及其調(diào)節(jié)（BP，n=6 GO terms），血管生成、發(fā)育和調(diào)控有關(guān)的生物過程（BP，n=30 GO terms），眼睛形態(tài)發(fā)生及視網(wǎng)膜血管形成和發(fā)育（BP，n=11 GO terms）等1 104 個(gè)GO 功能類別（P-adjust＜0.05，圖7B）。

圖7 上調(diào)表達(dá)（A）和下調(diào)表達(dá)（B）基因顯著富集的前30條生物學(xué)過程氣泡圖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 差異表達(dá)基因的KEGG富集分析

為查找差異表達(dá)基因參與的代謝通路或信號(hào)轉(zhuǎn)導(dǎo)過程，對(duì)差異表達(dá)基因進(jìn)行KEGG 富集分析，發(fā)現(xiàn)野生塔里木馬鹿中上調(diào)表達(dá)的基因顯著富集在核糖體（ko03010）、真核生物的核糖體生物發(fā)生（ko03008）、剪接體（ko03040）、蛋白酶體（ko03050）、氧化磷酸化（ko00190）、RNA 轉(zhuǎn)運(yùn)（ko03013）、次級(jí)代謝產(chǎn)物的生物合成（ko01110）、酮體的合成與降解（ko00072）和帕金森病（ko05012）等12 個(gè)代謝通路上（P-adjust＜0.05，圖8A）。下調(diào)表達(dá)基因主要涉及ECM 受體相 互作用（ko04512）、血小板活化（ko04611）、造血細(xì)胞譜系（ko04640）、補(bǔ)體與凝血級(jí)聯(lián)反應(yīng)（ko04610）、PI3K-Akt 信號(hào)通路（ko04151）、原發(fā)性免疫缺陷（ko05340）、脂肪細(xì)胞脂解的調(diào)節(jié)（ko04923）、蛋白質(zhì)消化與吸收（ko04974）和cGMPPKG 信號(hào)通路（ko04022）等38 個(gè)代謝和信號(hào)通路（P-adjust＜0.05，圖8B）。

2.8 qRT-PCR驗(yàn)證

為驗(yàn)證RNA-Seq 結(jié)果的可靠性，從野生塔里木馬鹿皮膚組織上調(diào)表達(dá)基因中選擇了與皮膚發(fā)育相關(guān)的基因（K2C6A），熱應(yīng)激相關(guān)基因（HSPA4、TRAP1和HSPH1）和氧化應(yīng)激相關(guān)基因（PRDX2和PRDX6）進(jìn)行qRT-PCR 檢測。將qRT-PCR 檢測到的目的基因的差異倍數(shù)與轉(zhuǎn)錄組測序結(jié)果進(jìn)行比較，并用柱狀圖進(jìn)行可視化（圖9），可以看出，目的基因的差異倍數(shù)和轉(zhuǎn)錄組測序結(jié)果有一定的差異，但表達(dá)趨勢一致，表明RNA-Seq測序結(jié)果較可靠，這些差異表達(dá)基因在塔里木馬鹿干旱環(huán)境適應(yīng)性中可能發(fā)揮著潛在的作用。

圖9 qRT-PCR和RNA-Seq檢測候選基因的相對(duì)表達(dá)量Fig.9 Comparison of candidate gene expression differences detected by qRT-PCR and RNA-Seq

3 討論

極端干旱地區(qū)的主要特征是太陽輻射較強(qiáng)、極端溫度、水和食物資源貧乏［18］。在這些地區(qū)，動(dòng)物長時(shí)間受強(qiáng)烈的紫外線輻射，可能會(huì)增加皮膚的失水速度或引起皮膚損傷［19］。在雙峰駝、埃及肥尾羊、努比亞山羊和烏干達(dá)本地山羊等干旱—沙漠環(huán)境中生活的有蹄類動(dòng)物中已篩選出與皮膚保護(hù)策略相關(guān)的多種候選基因［5，8-11］，塔里木馬鹿在塔里木盆地極端干旱環(huán)境中，可能也形成了皮膚保護(hù)策略來促進(jìn)其適應(yīng)性。除此之外，受長期家養(yǎng)馴化影響，圈養(yǎng)和野生塔里木馬鹿在飲食和生存環(huán)境方面存在較大差異，這種差異使圈養(yǎng)和野生塔里木馬鹿的皮膚在體溫調(diào)節(jié)、水分調(diào)節(jié)和對(duì)身體的保護(hù)功能方面產(chǎn)生了一定的差異，故而推測其在基因表達(dá)方面應(yīng)該也存在一定的差別，這些差異表達(dá)基因是否與該物種干旱環(huán)境適應(yīng)性相關(guān)，筆者以此為切入點(diǎn)開展本研究。

3.1 對(duì)照組選擇

本試驗(yàn)將棲息在塔里木胡楊自然保護(hù)區(qū)的野生塔里木馬鹿作為研究對(duì)象，將來自昌吉市馬鹿養(yǎng)殖場的圈養(yǎng)塔里木馬鹿作為對(duì)照組進(jìn)行RNA-Seq 和qRT-PCR 比較分析，檢測目的基因表達(dá)量。采集的塔里木馬鹿均來自干旱環(huán)境，但相對(duì)于圈養(yǎng)馬鹿的棲息環(huán)境，野生塔里木馬鹿棲息的塔里木胡楊自然保護(hù)區(qū)的年均氣溫較高（10.3 ℃）、年均降水量較低、年均蒸發(fā)量特別大，CGIAR-CSI 全球干旱度數(shù)據(jù)庫顯示，該地區(qū)的干旱指數(shù)（aridity index，AI）為0.015～0.023，屬于超干旱地區(qū)（AI＜0.05），而昌吉市的年均氣溫較低（6.8 ℃），干旱指數(shù)為0.102～0.112，屬于干旱地區(qū)（0.05＜AI＜0.20）［20］。野生塔里木馬鹿世世代代棲息在塔里木盆地惡劣的自然環(huán)境中，不僅遭受自然環(huán)境中有害因素（如強(qiáng)烈UV 輻射和機(jī)械力）的損害，還受夏季干旱少雨、酷熱，冬季嚴(yán)寒多風(fēng)和水資源匱乏的影響［21］。此外，年均氣溫和年均降水量是影響該物種遺傳多樣性的主要環(huán)境因子［22］。皮膚作為哺乳動(dòng)物直接同外界環(huán)境接觸的器官，執(zhí)行排汗、調(diào)節(jié)體溫和保護(hù)身體的功能，與其感受到的外界環(huán)境刺激以及水合狀態(tài)相關(guān)［23］。最重要的是，本研究中圈養(yǎng)的塔里木馬鹿是從1983 年開始通過人工繁殖方式自繁自養(yǎng)［24］的馬鹿后代，其被長期飼養(yǎng)在圈舍內(nèi)涼棚下，有充足的水源供應(yīng)，能夠有效避免高溫和UV輻射對(duì)身體的直接照射，從而能夠減少皮膚蒸發(fā)的水分和強(qiáng)烈UV輻射對(duì)皮膚的損傷。以上棲息環(huán)境及飼養(yǎng)條件的差別對(duì)圈養(yǎng)和野生塔里木馬鹿皮膚中相關(guān)基因的表達(dá)產(chǎn)生一定影響。而且，在生命活動(dòng)的不同階段、不同的細(xì)胞及不同的環(huán)境中，基因的表達(dá)均具有較大的差異，因此，轉(zhuǎn)錄組包含的信息也存在較大的差異。

3.2 塔里木馬鹿干旱環(huán)境適應(yīng)性相關(guān)基因篩選

本研究發(fā)現(xiàn)，在野生塔里木馬鹿皮膚組織中的角蛋白（KRT222）、角蛋白相關(guān)蛋白（KRTAP11-1、KRTAP13-2和KRTAP14）、角質(zhì)形成細(xì)胞分化相關(guān)蛋白（KRTDAP和KRTCAP3）、角蛋白Ⅰ型表皮（KRT3）、角蛋白Ⅱ型表皮（KRT83和KRT84）、角蛋白Ⅰ型細(xì)胞骨架（KRT10、KRT14、KRT17、KRT19、KRT26、KRT28和KRT39）、角蛋白Ⅱ型細(xì)胞骨架（KRT1、KRT5、K2C6A、KRT71、KRT72、KRT73、KRT74、KRT75和KRT75）、角蛋白Ⅰ型微纖 維（LOC122695079和LOC122695080）和角蛋白Ⅱ型微纖維（LOC102506783和LOC122681432）等31條基因顯著上調(diào)表達(dá)（P＜0.05）。這些基因在調(diào)節(jié)毛囊和表皮發(fā)育、皮膚角質(zhì)形成細(xì)胞分化和維持分層上皮中發(fā)揮重要作用。在機(jī)體內(nèi)，角蛋白相互作用形成角蛋白細(xì)絲，組成稠密的網(wǎng)絡(luò)，為皮膚及其他組織提供支撐，并保護(hù)這些組織免受外部摩擦和機(jī)械力損傷，減少皮膚損傷，加強(qiáng)保濕滋潤，在維持皮膚屏障功能、防止水分流失和體溫調(diào)節(jié)等方面具有重要作用［25］。有些角蛋白（如K2C6A）與其他角蛋白參與皮膚屏障、表皮傷口愈合以及啟動(dòng)和調(diào)節(jié)皮膚相關(guān)免疫反應(yīng)的過程［26］。除了上述角蛋白及其相關(guān)基因外，多種跨 膜蛋白基因（TMEM40、TMEM54、TMEM62、TMEM79、TMEM80、TMEM184A、TMEM254和TMEM256）在野生塔里木馬鹿皮膚組織中上調(diào)表達(dá)，其中TMEM79基因參與多種過程，包括上皮細(xì)胞成熟、皮膚屏障的建立和表皮發(fā)育的正調(diào)控，并在表皮發(fā)育、毛囊形態(tài)發(fā)生以及在角質(zhì)化的上游或內(nèi)部起作用［27］。這些角蛋白和跨膜蛋白顯著富集于角質(zhì)化（GO：0031424）、角質(zhì)形 成細(xì)胞分化（GO：0030216）、角化（GO：0070268）、表皮細(xì)胞分化（GO：0009913）、皮膚發(fā)育（GO：0043588）、表皮發(fā)育（GO：0008544）、角蛋白絲（GO：0045095）、表皮lamellar body（GO：0097209）、角質(zhì)化包膜（GO：0001533）、表皮細(xì)胞分化的調(diào)控（GO：0045604）、表皮發(fā)育的調(diào)控（GO：0045682）和皮膚屏障的建立（GO：0061436）等（P-adjust＜0.05）。富集不顯著的功能類別中，上調(diào)表達(dá)基因涉及皮膚表皮發(fā)育、傷口愈合、表皮細(xì)胞擴(kuò)散、表皮生長因子受體信號(hào)通路及其調(diào)控等，與皮膚發(fā)育和表皮細(xì)胞發(fā)育相關(guān)的多種生物學(xué)過程和分子功能（n=20 GO terms）相關(guān)。更重要的是，沒有一個(gè)角蛋白及其相關(guān)蛋白質(zhì)基因在野生塔里木馬鹿皮膚組織中呈現(xiàn)下調(diào)表達(dá)。以上基因及其參與的生物學(xué)功能說明皮膚保護(hù)屏障在野生塔里木馬鹿極端干旱環(huán)境適應(yīng)過程中可能發(fā)揮著重要作用。

本研究中的上調(diào)表達(dá)基因NUDT5和COL17A1中，COL17A1編碼XVI 型膠原的α鏈，該鏈?zhǔn)沁B接基底上皮細(xì)胞和基質(zhì)半橋粒的關(guān)鍵結(jié)構(gòu)成分［28］，該基因突變或其蛋白質(zhì)的水解能引起皮膚萎縮、脫發(fā)和毛囊干細(xì)胞的老化［29］。而NUDT5可防止蛋白質(zhì)的非酶性核糖基化，從核苷酸中去除被氧化的核苷，降低細(xì)胞在受氧化應(yīng)激刺激下出現(xiàn)錯(cuò)誤的DNA 復(fù)制率，保護(hù)細(xì)胞免受氧化應(yīng)激損傷［30-31］。本研究結(jié)果在轉(zhuǎn)錄水平上初步支持了NUDT5和COL17A1基因在塔里木馬鹿干旱環(huán)境適應(yīng)性的潛在作用。

塔里木馬鹿在塔里木盆地的極端干旱棲息環(huán)境中，遭受UV 輻射、溫度和鹽度等多種脅迫因素的影響。這些因素會(huì)導(dǎo)致塔里木馬鹿機(jī)體內(nèi)的活性氧（reactive oxygen species，ROS）大量積累，當(dāng)ROS的形成超過靶細(xì)胞的抗氧化防御能力時(shí)，氧化應(yīng)激會(huì)導(dǎo)致氧化損傷。皮膚作為機(jī)體第一保護(hù)屏障，要保護(hù)機(jī)體免受這些有害因素引起的氧化損傷。在本研究中，野生塔里木馬鹿皮膚組織中抗氧化活性基因（PRDX2、PRDX6和GPX2）、轉(zhuǎn)谷氨酰胺酶編碼基因（TGM1和TGM3）、ATP 結(jié)合盒亞家族A成員12（ABCA12）和DNA 損傷修復(fù)基因（NSMF、HOXC12、GADD45GIP1、GADD45B和GADD45）呈現(xiàn)出明顯的上調(diào)表達(dá)（P＜0.05），其中有些基因在雙峰駝（TGM1）［5］、努比亞山羊（ABCA12）［8］、埃及肥尾羊（TGM3）［9］和烏干達(dá)本地山羊（HOXC12）［10］等物種中受正選擇，并參與這些物種在干旱—炎熱及沙漠環(huán)境中適應(yīng)性皮膚保護(hù)策略的形成，以保護(hù)細(xì)胞免受沙漠中強(qiáng)烈的UV輻射損害，并參與炎熱環(huán)境的體溫調(diào)節(jié)，如HOXC12基因通過調(diào)節(jié)角蛋白分化特異基因，在毛囊分化、生長發(fā)育中發(fā)揮作用［11］。TGM3廣泛表達(dá)于皮膚細(xì)胞，特別是角質(zhì)形成細(xì)胞和角朊細(xì)胞。在角質(zhì)形成細(xì)胞分化過程中，TGM3交聯(lián)結(jié)構(gòu)蛋白和脂質(zhì)形成角質(zhì)化細(xì)胞膜，為表皮提供抵御有害環(huán)境刺激（如UV 輻射）的屏障功能［32］。TGM1編碼的蛋白質(zhì)在角質(zhì)化細(xì)胞膜的組裝中具有完整的功能，形成穩(wěn)定的不溶性的保護(hù)屏障以防止太陽輻射、水分流失和病原體的侵入，并參與構(gòu)成表皮屏障功能的修復(fù)機(jī)制，從而對(duì)皮膚起保護(hù)作用［33］。ABCA12的脂質(zhì)轉(zhuǎn)運(yùn)是角質(zhì)形成細(xì)胞分化和表皮形態(tài)發(fā)生所必需的，因此在表皮脂質(zhì)屏障形成和角質(zhì)形成細(xì)胞分化中起重要作用［34］。NSMF基因在DNA 損傷增加的應(yīng)激條件下，促進(jìn)復(fù)制應(yīng)激誘導(dǎo)的DNA 損傷反應(yīng)以維持基因組的完整性［35］。而抗氧化活性基因（PRDX2、PRDX6和GPX2）可以保護(hù)角質(zhì)化細(xì)胞免受ROS誘導(dǎo)的凋亡，促進(jìn)傷口愈合，減輕氧化應(yīng)激對(duì)皮膚造成的損傷［36］。此外，基于全基因組重測序的群體遺傳學(xué)分析，發(fā)現(xiàn)TGM1、PRDX6和NSMF等基因在塔里木馬鹿中受選擇［13］。GADD45基因能發(fā)揮修復(fù)DNA 和調(diào)控細(xì)胞周期等作用［37］。這些DNA 損傷修復(fù)相關(guān)基因在保持基因組的完整性［9］和維持細(xì)胞正常功能方面具有至關(guān)重要的作用。因此，本研究結(jié)果在轉(zhuǎn)錄水平上初步支持了在野生塔里木馬鹿皮膚組織中的上調(diào)表達(dá)基因可能通過參與氧化應(yīng)激反應(yīng)，形成皮膚保護(hù)屏障和維持皮膚穩(wěn)態(tài)，發(fā)揮保護(hù)機(jī)體的作用。

在漫長的進(jìn)化過程中，塔里木馬鹿不僅對(duì)塔里木盆地夏季炎熱表現(xiàn)出獨(dú)特的適應(yīng)性，還可以適應(yīng)較低溫的環(huán)境［38］。細(xì)胞對(duì)熱/冷應(yīng)激的耐受是由熱休克蛋白（HSPs）家族介導(dǎo)的。熱休克蛋白是一類高度保守的蛋白質(zhì)，在逆環(huán)境中組成性表達(dá)和/或誘導(dǎo)表達(dá)，它具有分子伴侶功能，可調(diào)節(jié)信號(hào)轉(zhuǎn)導(dǎo)通路，參與免疫應(yīng)答和細(xì)胞周期，并具有抗氧化應(yīng)激作用。另一方面，高鹽攝入增加了熱休克蛋白的表達(dá)，從而保護(hù)細(xì)胞免受高尿素的損傷［39］。本研究發(fā)現(xiàn)，熱休克蛋白基因（HSPA4、HSPA14、TRAP1、HSPH1、HSPD1、HSP90和HSP110）與熱休克蛋白家族具有序列和功能同源性的J 結(jié)構(gòu)域蛋白基因（DNAJC2、DNAJC7、DNAJC19、DNAJC15和DNAJC21）在野生塔里木馬鹿皮膚組織中上調(diào)表達(dá)，顯著富集于蛋白 質(zhì)折疊、蛋白質(zhì)生物合成和代謝等生物學(xué)過程（P-adjust＜0.05）。這些熱休克蛋白基因能夠幫助蛋白質(zhì)正確折疊，促進(jìn)將錯(cuò)誤折疊或變性的蛋白質(zhì)清除，減輕應(yīng)激對(duì)生物體的損害。其中，J 結(jié)構(gòu)域蛋白是HSP40 家族的成員之一［40］，在細(xì)胞遭受干旱、高溫、低溫和強(qiáng)光等條件下能誘導(dǎo)表達(dá)，并應(yīng)答脅迫應(yīng)激，從而調(diào)節(jié)熱休克蛋白的功能，發(fā)揮分子伴侶的作用，參與蛋白質(zhì)折疊［41］。檢測炎熱干旱環(huán)境中的烏干達(dá)本地山羊品種，發(fā)現(xiàn)該物種的熱休克蛋白與其體溫調(diào)節(jié)適應(yīng)性相關(guān)［10-11］。因此，筆者推測在野生塔里木馬鹿皮膚組織中，上調(diào)表達(dá)的熱休克蛋白基因以及J 結(jié)構(gòu)域蛋白基因，在塔里木馬鹿干旱環(huán)境適應(yīng)性中的應(yīng)激耐受性和體溫調(diào)節(jié)方面可能具有重要的作用。