任春艷，魏志遠(yuǎn)，饒軍華，饒義劍，陳建歡

1)江南大學(xué)無(wú)錫醫(yī)學(xué)院，江蘇無(wú)錫 214122；2)江南大學(xué)生物工程學(xué)院，江蘇無(wú)錫214122；3)江南大學(xué)-廣東科學(xué)院動(dòng)物研究所“慢性病靈長(zhǎng)類研究聯(lián)合實(shí)驗(yàn)室”，江蘇無(wú)錫214122；4)廣東科學(xué)院動(dòng)物研究所，廣東廣州510260

核糖核酸(ribonucleic acid, RNA)編輯是在轉(zhuǎn)錄水平上發(fā)生的RNA核苷酸序列的改變[1]．經(jīng)典RNA編輯主要包括由ADAR(adenosine deaminases acting on RNA)蛋白家族介導(dǎo)的腺苷到肌苷(A-I)和APOBEC(apolipoprotein B mRNA editing enzyme, catalytic polypeptide-like)蛋白家族介導(dǎo)的胞苷到尿苷(C-U)的編輯．RNA編輯目前已在動(dòng)物、植物和人等多種生物中被發(fā)現(xiàn)，并參與多種生理及病理學(xué)過(guò)程．近期研究提示，RNA編輯在應(yīng)激、神經(jīng)系統(tǒng)疾病和精神病理性疾病條件下存在明顯變化[2-3]，提示其可能參與這些過(guò)程或疾病的分子機(jī)制．

社交挫敗模型主要用于研究創(chuàng)傷后應(yīng)激障礙(post-traumatic stress disorder，PTSD)、抑郁癥和其他應(yīng)激相關(guān)的疾病．文獻(xiàn)[4]研究表明，社交挫敗應(yīng)激可誘導(dǎo)RNA編輯的變化．DICK等[5]的研究提示，小鼠在慢性社交挫敗應(yīng)激(chronic social defeat stress, CSDS)后，其大腦內(nèi)側(cè)前額葉皮層和基底外側(cè)杏仁核內(nèi)存在A-I編輯的動(dòng)態(tài)調(diào)節(jié)．但是C-U的RNA編輯在情緒應(yīng)激模型中的變化如何，目前尚不明確．有意思的是，單純通過(guò)目擊創(chuàng)傷性事件所引起的情緒應(yīng)激(emotional stress, PS)與通過(guò)身體接觸所引起的生理應(yīng)激(physical stress, ES)誘導(dǎo)的小鼠大腦轉(zhuǎn)錄組表達(dá)譜存在差異，因此，不同應(yīng)激條件下CSDS小鼠模型中C-U的RNA編輯尚待研究．

本研究通過(guò)RNA測(cè)序，在轉(zhuǎn)錄組范圍內(nèi)識(shí)別大腦的腹側(cè)被蓋區(qū)(ventral tegmental area, VTA)中的C-U RNA編輯，研究其分布情況．通過(guò)比較不同條件下的CSDS過(guò)程中C-U RNA編輯情況，分析其對(duì)基因功能和通路的潛在影響，以及參與調(diào)控的潛在分子機(jī)制．

1 材料與方法

1.1 RNA測(cè)序數(shù)據(jù)的獲取

RNA測(cè)序原始數(shù)據(jù)從美國(guó)國(guó)家生物信息中心(National Center for Biotechnology Information, NCBI)的Gene Expression Omnibus(GEO)數(shù)據(jù)庫(kù)中(https://www.ncbi.nlm.nih.gov/geo/)下載獲得．ES組、PS組以及對(duì)照組(每組設(shè)平行實(shí)驗(yàn)N=3)的成年雄性小鼠大腦腹側(cè)被蓋區(qū)(ventral tegmental area, VTA)RNA測(cè)序數(shù)據(jù)(GSE36005)[4]．RNA編輯位點(diǎn)的驗(yàn)證數(shù)據(jù)包含44個(gè)成年小鼠大腦VTA RNA測(cè)序(GSE89692)[6-7]．

1.2 RNA測(cè)序數(shù)據(jù)的比對(duì)

利用FASTQC軟件對(duì)以上獲得的原始測(cè)序數(shù)據(jù)進(jìn)行質(zhì)量控制分析．質(zhì)控合格的測(cè)序讀段(reads)，利用RNA STAR (2.7.0e版)[8]軟件與加利福尼亞大學(xué)圣克魯茲分校(University of California, Santa Cruz, UCSC)小鼠的基因組序列(mm10版)進(jìn)行比對(duì)，分析其RNA剪接接頭，將其定位(map)到基因組上，并生成二進(jìn)制比對(duì)映射(binary alignment map, BAM)格式的比對(duì)結(jié)果文件．利用Samtools(1.9版)軟件將BAM文件中的讀段去重復(fù)(deduplication)[9]，只保留有單一定位的讀段，并使用GATK(4.1.3版)軟件對(duì)BAM文件進(jìn)行堿基質(zhì)量評(píng)分重校準(zhǔn)[10]．

1.3 RNA編輯位點(diǎn)的鑒定及注釋

利用VarScan(2.4.3版)軟件對(duì)校準(zhǔn)處理后的BAM文件進(jìn)行RNA單核苷酸變異(single nucleotide variant, SNV)識(shí)別[11]．識(shí)別標(biāo)準(zhǔn)為堿基質(zhì)量≥25，測(cè)序深度≥10，U堿基(RNA測(cè)序中被轉(zhuǎn)換為T堿基)深度≥2且頻率≥1%，并用VarScan的 fpfilter命令和默認(rèn)參數(shù)過(guò)濾假陽(yáng)性變異．對(duì)達(dá)到以下任一條件的SNV均予以濾除[12]：① 位于簡(jiǎn)單重復(fù)序列或5個(gè)連續(xù)相同堿基內(nèi)；② 位于線粒體內(nèi)；③ 距離RNA剪接接頭6個(gè)核苷酸以內(nèi)；④ 距離插入或刪除變異1個(gè)核苷酸以內(nèi)；⑤ 在dbSNP數(shù)據(jù)庫(kù)中的已知變異相吻合．為進(jìn)一步提高RNA編輯位點(diǎn)識(shí)別的可信度，只保留在樣本或驗(yàn)證數(shù)據(jù)樣本中都被檢測(cè)到的，即可重復(fù)驗(yàn)證的編輯水平≥1%的C-T SNV，將其定義為高置信度變異．SNV的注釋使用Ensembl Variant Effect Predictor (VEP， https://www.ensembl.org/vep)軟件進(jìn)行[13]．

1.4 基因表達(dá)的定量及差異分析

采用FeatureCounts軟件對(duì)RNA STAR生成的比對(duì)文件進(jìn)行分析，獲得RNA表達(dá)的計(jì)數(shù)，并用EdgeR軟件(3.7版)計(jì)算標(biāo)準(zhǔn)化的基因表達(dá)水平(transcripts per kilobase million，TPM)[14]．

1.5 主成分分析

利用1.4節(jié)獲得的基因表達(dá)的定量數(shù)據(jù)，使用R軟件(3.6.3版)中的函數(shù)Prcomp進(jìn)行主成分分析(principal component analysis, PCA)，并使用ggplot2(2.2.1版)軟件包進(jìn)行數(shù)據(jù)可視化．

1.6 基因功能和通路的富集分析

為了理解RNA編輯的潛在生物學(xué)效應(yīng)，利用Enrichr軟件(https://maayanlab.cloud/Enrichr/)[15]分析發(fā)生RNA編輯的基因所富集的gene ontology(GO)和京都基因與基因組百科全書(shū)(Kyoto encyclopedia of genes and genomes，KEGG)通路等，以錯(cuò)誤發(fā)現(xiàn)率(false discovery rate, FDR) < 0.05作為篩選條件．

1.7 統(tǒng)計(jì)學(xué)分析

使用Kruskal-Wallis(KW)非參數(shù)檢驗(yàn)，比較樣本組間RNA編輯水平或基因表達(dá)水平．組間的兩兩比較用Dunn校驗(yàn)進(jìn)行，并用Benjamini-Hochberg法對(duì)多重比較的P值進(jìn)行校正．以P<0.05作為顯著差異篩選標(biāo)準(zhǔn)．采用Spearman相關(guān)分析計(jì)算編輯位點(diǎn)間的相關(guān)系數(shù)R． 在有多個(gè)編輯位點(diǎn)的基因中，差異位點(diǎn)的富集分析使用Fisher精確檢驗(yàn)進(jìn)行．

2 結(jié)果與分析

2.1 小鼠VTA的C-U RNA 編輯

在3組成年小鼠VTA的RNA-Seq中，共在5 866個(gè)基因中發(fā)現(xiàn)了16 198個(gè)高置信度C-U RNA編輯位點(diǎn)，見(jiàn)圖1(C-U RNA編輯位點(diǎn)的具體信息請(qǐng)掃描論文末頁(yè)右下角二維碼查看表S1)．這些位點(diǎn)廣泛存在于各染色體中，其樣本間最高的RNA編輯頻率為1%～100%．從每個(gè)基因的編輯位點(diǎn)數(shù)量來(lái)看，部分基因位點(diǎn)密度較亮，這些位點(diǎn)涵蓋了多種功能類型，按數(shù)量多少依次為3′-非編輯區(qū)(untranslated region, UTR)變異(5 616個(gè))、錯(cuò)義變異(4 953個(gè))、同義變異(4 172個(gè))、無(wú)義變異(716個(gè))、5′-UTR變異(594個(gè))和非編碼轉(zhuǎn)錄本的外顯子變異(147個(gè))， 如圖2．

圖1 成年小鼠VTA轉(zhuǎn)錄組中基因表達(dá)(外圈)、C-U RNA編輯位點(diǎn)(中圈)和位點(diǎn)間相關(guān)性的Circos圖Fig.1 Circos plot showing the gene expression (outer circle), distribution of C-U RNA editing events (middle circle)and their correlation in adult mouse VTA

圖2 成年小鼠VTA轉(zhuǎn)錄組中C-U RNA編輯SNV的功能類型Fig.2 Functional categories of C-U RNA editing variants in the mouse VTA

由圖2可見(jiàn)，部分基因檢測(cè)到多個(gè)編輯位點(diǎn)．編輯位點(diǎn)數(shù)量排前10位的基因請(qǐng)掃描論文末頁(yè)右下角二維碼查看表S2．

圖3 48個(gè)位點(diǎn)在組間存在顯著差異RNA編輯(所有位點(diǎn)的KW檢驗(yàn)均有P<0.05)Fig.3 Forty-eight sites with significant differential RNA editing among the three groups RNA editing(all KW test P<0.05)

2.2 不同應(yīng)激條件下的差異C-U RNA編輯

為分析不同應(yīng)激條件下RNA編輯的改變，通過(guò)KW檢驗(yàn)比較位點(diǎn)在不同分組間的RNA編輯水平，發(fā)現(xiàn)共有48個(gè)基因中的48個(gè)位點(diǎn)存在顯著差異(圖3)．這些基因/位點(diǎn)有22個(gè)在ES中， 9個(gè)在PS組中，15個(gè)在對(duì)照組中有顯著較高的編輯水平，另有2個(gè)在ES中有較低的編輯水平．PCA結(jié)果表明，這些差異編輯可以很好地將3組小鼠分開(kāi)，其中，PC1和PC2的貢獻(xiàn)率分別為51.78%和30.74%(圖4)．Fisher檢驗(yàn)分析發(fā)現(xiàn)，這48個(gè)呈差異編輯的基因中44個(gè)(91.7%)具有1個(gè)以上的編輯位點(diǎn)，顯示差異編輯顯著富集于具多個(gè)編輯位點(diǎn)的基因中(優(yōu)勢(shì)比=8.5，95%置信區(qū)間為3.1～32.5，F(xiàn)isher檢驗(yàn)P=1.74×10-7)．其中，差異C-U RNA編輯位點(diǎn)在有多個(gè)編輯位點(diǎn)基因中的富集情況請(qǐng)掃描論文末頁(yè)右下角二維碼查看表S3．

圖4 48個(gè)差異編輯位點(diǎn)的主成分分析Fig.4 PCA analysis using forty-eight sites with differential RNA editing

2.3 RNA編輯的功能富集

圖5 各組中在2個(gè)或以上樣本中檢測(cè)到的C-U RNA編輯位點(diǎn)比較的維恩圖Fig.5 Venn plot comparing the C-U editing sites detected in two or more samples among groups

圖6 GO分析顯示各組獨(dú)有的編輯位點(diǎn)的基因呈現(xiàn)出差異化的基因功能富集Fig.6 GO analysis showing differential enrichment of gene function in genes with editing sites unique to groups

進(jìn)一步比較不同組中RNA編輯的存在情況，發(fā)現(xiàn)ES、PS和對(duì)照組中分別有307、223和301個(gè)編輯位點(diǎn)特異地存在于組內(nèi)至少2個(gè)樣本中，而在另外2組未檢測(cè)到(圖5)．Enrichr基因列表富集分析結(jié)果表明，這些各組獨(dú)有的RNA編輯呈現(xiàn)出不同的基因功能特點(diǎn)(圖6)．ES組特有的RNA編輯主要富集在軸突發(fā)生、神經(jīng)系統(tǒng)發(fā)育、細(xì)胞內(nèi)信號(hào)轉(zhuǎn)導(dǎo)的負(fù)調(diào)控和軸突發(fā)育等生物學(xué)過(guò)程．這些基因主要參與Ras(rat sarcoma)等多種GTP(guanosine-5′-triphosphate)酶的結(jié)合活性、GTP酶調(diào)控因子的活性和鈣黏蛋白的結(jié)合活性等分子功能，與樹(shù)突等細(xì)胞組分有關(guān)．PS組獨(dú)有的RNA編輯主要富集在淀粉樣蛋白及其原纖維形成的調(diào)控和液泡介導(dǎo)運(yùn)輸?shù)壬飳W(xué)過(guò)程，與細(xì)胞皮層區(qū)的細(xì)胞組分有關(guān)．而對(duì)照組中獨(dú)有的RNA編輯主要富集在神經(jīng)系統(tǒng)發(fā)育、神經(jīng)元發(fā)育和具細(xì)胞質(zhì)膜的細(xì)胞投射組裝的正調(diào)控等生物學(xué)過(guò)程．這些基因主要參與小GTP酶結(jié)合活性等分子功能，與裂解性小泡膜、晚期內(nèi)體、溶酶體、溶酶體膜、基于肌動(dòng)蛋白的細(xì)胞投射、SWI/SNF(SWItch/sucrose non-fermentable)復(fù)合物、高爾基體和反面高爾基網(wǎng)絡(luò)等細(xì)胞組分有關(guān)．

2.4 ES組中血清和糖皮質(zhì)激素依賴性激酶1基因的改變

通過(guò)比較差異表達(dá)的基因和特異編輯位點(diǎn)發(fā)現(xiàn)，與對(duì)照相比，血清和糖皮質(zhì)激素依賴性激酶1(serum/glucocorticoid regulated kinase 1,Sgk1)基因是ES組表達(dá)差異最為顯著上調(diào)的基因之一，而在PS組中改變不明顯(圖7(a))，同時(shí)，在10號(hào)染色體的第21 999 752位點(diǎn)的3′-UTR位置，Sgk1基因也僅在ES組中有特異性編輯(圖7(b))．

圖7 小鼠VTA中Sgk1 mRNA的表達(dá)存在顯著的組間差異及其在ES組中的特異性C-U RNA編輯Fig.7 Significantly differential mRNA expression and ES-specific C-U RNA editing of Sgk1 mRNA in the mouse VTA

3 討 論

CSDS的嚙齒類模型被廣泛用于抑郁癥等相關(guān)疾病的研究，為理解這部分疾病提供了重要的科學(xué)參考和依據(jù)．本研究首次在轉(zhuǎn)錄組范圍內(nèi)對(duì)CSDS小鼠模型中大腦VTA內(nèi)存在的C-U RNA編輯進(jìn)行了系統(tǒng)性研究．

VTA是大腦獎(jiǎng)賞系統(tǒng)的關(guān)鍵區(qū)域，與應(yīng)激相關(guān)行為關(guān)系密切[16]．因此，在VTA中廣泛存在的C-U RNA編輯可能與VTA的重要生物學(xué)功能有密切關(guān)系．這些編輯位點(diǎn)多集中于RNA的3′-UTR ，這與現(xiàn)有的C-U RNA編輯的研究結(jié)果一致．這些RNA編輯可能通過(guò)影響微小RNA與3′-UTR 的結(jié)合，從而調(diào)控轉(zhuǎn)錄后RNA的穩(wěn)定性．此外，有研究還提示3′-UTR RNA編輯可能影響RNA翻譯，從而作用于蛋白水平[17]．C-U RNA編輯在不同基因間的分布并不均勻，部分基因存在數(shù)量較多的編輯位點(diǎn)，可能是編輯熱點(diǎn)．編輯位點(diǎn)數(shù)量最多的10個(gè)基因多與中樞神經(jīng)系統(tǒng)的功能及精神疾病有關(guān)，可能參與相關(guān)的生物學(xué)過(guò)程．例如，在ATP酶 Na+/K+轉(zhuǎn)運(yùn)亞其Atp1a3基因發(fā)現(xiàn)了47個(gè)C-U RNA編輯位點(diǎn)．ATP1A3基因編碼神經(jīng)元特異性ATP依賴的跨膜鈉鉀泵的α亞基，有研究報(bào)導(dǎo)ATP1A3的錯(cuò)義變異與兒童期精神分裂癥相關(guān)[18]．此外，CSDS差異編輯基因91%以上集中于這些有1個(gè)以上編輯位點(diǎn)的基因中，提示這些編輯熱點(diǎn)所在的基因與CSDS之間存在潛在的重要聯(lián)系．與之相關(guān)的是，小膠質(zhì)細(xì)胞中C-U編輯功能的缺陷會(huì)加劇與年齡相關(guān)的中樞神經(jīng)系統(tǒng)病理生理學(xué)改變[19]．值得注意的是，已發(fā)現(xiàn)的RNA編輯變異中，編碼變異占了相當(dāng)大的比例(56.3%)，而編碼變異中的錯(cuò)義和無(wú)義變異加起來(lái)占編碼變異的54.3%．這可能意味著C-U RNA編輯可能對(duì)VTA中表達(dá)的蛋白序列有潛在影響．這種影響的結(jié)果及其產(chǎn)生的原因是什么，尚待進(jìn)一步研究．

與已報(bào)導(dǎo)的CSDS小鼠VTA中的A-I位點(diǎn)相比，C-U RNA編輯位點(diǎn)編輯水平普遍較低[5]．這可能意味著與A>I RNA編輯相比，從RNA測(cè)序數(shù)據(jù)中識(shí)別C-U RNA編輯可能面臨的難度更大．本研究中采用獨(dú)立樣本進(jìn)行編輯位點(diǎn)的驗(yàn)證，并針對(duì)多等位基因位點(diǎn)進(jìn)行了處理，從而提高了分析流程的可靠性．

在CSDS小鼠模型中，發(fā)現(xiàn)48個(gè)基因中的48個(gè)位點(diǎn)的編輯水平有顯著的組間差異．同時(shí)PCA結(jié)果顯示，這些位點(diǎn)PC1和PC2加起來(lái)可以解釋超過(guò)81%的樣本間變異．這些結(jié)果提示在ES、PS和對(duì)照組的中存在較明顯的差異編輯．這些存在差異編輯水平的基因部分為已知的精神類疾病易感基因，如谷氨酸離子受體GRIA4被報(bào)導(dǎo)是尼古丁依賴和重度抑郁共同的風(fēng)險(xiǎn)基因[20]；也有部分為編碼RNA結(jié)合蛋白，如ELAV樣RNA結(jié)合蛋白1 I(Elavl1)基因編碼的Hu抗原R，能結(jié)合mRNA的3′-UTR并增加其穩(wěn)定性，并與RNA編輯的功能調(diào)節(jié)相關(guān)[21]．

谷氨酸受體已知是A>I RNA編輯重要的靶基因[22]．本研究分析結(jié)果顯示，離子型谷氨酸受體Gria4(glutamate ionotropic receptor AMPA type subunit 4)中存在顯著的差異C-U RNA，提示谷氨酸受體相關(guān)基因可能也是重要靶點(diǎn)，并與應(yīng)激反應(yīng)相關(guān)．

特別重要的是，在這些特異性位點(diǎn)中， ES組中存在Sgk1基因同時(shí)存在特異性RNA編輯和表達(dá)水平改變．Sgk1編碼的血清和糖皮質(zhì)激素調(diào)節(jié)激酶1屬于絲氨酸/蘇氨酸蛋白激酶的一種，在應(yīng)激和抑郁癥等疾病中起重要作用[23]．在慢性、輕度、不可預(yù)見(jiàn)應(yīng)激或生命早期應(yīng)激誘導(dǎo)的抑郁癥大鼠模型中，海馬的Sgk1 mRNA水平顯著增加[24]．在未用藥抑郁癥患者的外周血中同樣檢測(cè)到Sgk1基因表達(dá)量增加[24]．首次發(fā)現(xiàn)的與ES相關(guān)的Sgk1 mRNA可能與其表達(dá)水平的增加有關(guān)，具體生物學(xué)功能有待進(jìn)一步研究．

結(jié) 語(yǔ)

通過(guò)首次識(shí)別成年雄性小鼠大腦VTA轉(zhuǎn)錄組中存在的C-U RNA編輯，將不同CSDS模型與對(duì)照進(jìn)行對(duì)比．結(jié)果揭示了CSDS模型特異性的C-U RNA編輯改變，為進(jìn)一步理解RNA編輯參與CSDS的分子機(jī)制提供了依據(jù)．