利用mRNA-seq和Ribo-seq技術研究低氧條件下SW620細胞的基因表達特征

梅 勇，梁綺華，鄧載安，汪勁松

(1.湖北師范大學生命科學學院 食用野生植物保育與利用湖北省重點實驗室，湖北 黃石 435002；2.深圳技術大學健康與環境工程學院，深圳 坪山 518118)

0 引言

核糖體印跡測序(ribosome profiling, Ribo-seq)是由Ingolia等發明的一種對基因翻譯組學研究的手段[1]，可以對正在翻譯的RNA進行選擇性測序，可用于揭示基因翻譯特征。其原理是利用翻譯抑制劑(放線菌酮)阻止核糖體翻譯，然后消化核酸，捕獲和富集核糖體內被保護的免受核酸酶切割的短基因片段，這些小片段又被稱為核糖體足跡(ribosome footprints，RFs)[2]，通過對RFs的高通量測序獲得正在翻譯的基因信息和定量信息。此外，根據數據分析可以獲得核糖體在每個轉錄本上的分布和密度，推斷出起始密碼子位置(包括非ATG起始)、密碼子使用偏愛、上游ORFs (uORFs：開放閱讀框)和翻譯暫停景觀等信息。Ribo-seq還能觀察到單個轉錄物的全局3-nt周期性，可發現小ORFs[3].核酸序列生物信息學分析能選擇與核糖體足跡一致的讀數，將讀數的起點或終點重新校準到核糖體的肽基位點(P位點)，并對轉錄物、密碼子、多密碼子基序的讀數進行定量[4]。

越來越多的實驗和研究表明，低氧環境是實體瘤中心常常出現的現象，特別是結直腸癌[5]。低氧誘導因子1α(HIF1-α)是參與低氧轉錄激活的主要轉錄因子，在低氧條件下，脯氨酰羥化酶和天冬酰胺酰羥化酶被抑制，從而抑制關鍵HIF1-α的羥基化，阻斷泛素-蛋白酶體系統對HIF1-α的降解作用，使得HIF1-α的快速積累以適應低氧環境[6]，HIF1-α通過激活下游靶基因來參與癌細胞的增殖、凋亡、血管生成和轉移以及放化療敏感性等生理病理過程[7]。有臨床數據證明，腫瘤內低氧環境或基因突變引起的HIF1-α高表達與許多類型癌癥患者的死亡率增加有關[8]。然而，低氧條件下的細胞整體基因表達和翻譯特征仍不清楚[9]。

細胞中蛋白質水平取決于基因轉錄水平、mRNA翻譯水平和蛋白質降解水平及其他因素[10]。為了研究低氧條件下結直腸癌細胞基因表達特性和調控機制，我們利用CoCl2化學模擬低氧條件來處理SW620細胞，使用mRNA-seq和Ribo-seq方法檢測CoCl2處理前后轉錄和翻譯水平的變化，并通過比較基因的翻譯效率來尋找低氧條件下調控的基因。本項目的研究為低氧調控結直腸癌細胞的基因表達提供了一個新的角度。

1 材料和方法

1.1 試劑、細胞株和儀器

胎牛血清購自美國Gibco公司(10099141C)；DMEM-高糖培養基(SH30022.01)、PBS(SH30256.01)、青霉素/鏈霉素(SV30010)、胰酶(SH30042.01)均購自美國Hyclone公司；HIF1-α抗體(3716S)購自美國CST公司；GAPDH抗體(60004-1-Ig)、羊抗兔IgG(SA00001-2)、羊抗鼠IgG(SA00001-1)購自武漢三鷹生物技術有限公司；GLUT3抗體(bs-1207R)、AARE抗體(bs-5983R)購買自北京博奧森生物技術有限公司；RIPA細胞裂解液(R0278-500ML)、CoCl2(C8661-100G)購自Sigma公司；RNA提取試劑盒(AG21017)、反轉錄試劑盒(AG11706)均購自湖南艾科瑞生物工程有限公司；引物均購自華大基因；人結腸直腸腺癌SW620細胞系購自南京科佰生物科技有限公司；細胞培養箱購自美國Thermo 公司；凝膠成像系統購自上海天能公司；TAdvanced 96G PCR 擴增儀購自德國Biometra公司；QuantStudio 6 Flex實時熒光定量PCR儀購自美國Thermo公司。

1.2 方法

1.2.1 細胞復蘇及培養 SW620細胞快速從液氮罐取出，放入37℃恒溫水浴鍋中輕輕晃動融化，細胞凍存液轉移到含5 mL(10%胎牛血清和1%雙抗)DMEM-高糖培養基的25 mm2培養瓶中，在37℃含5% CO2細胞培養箱中培養。

1.2.2 CoCl2化學模擬低氧處理 配制200 mM CoCl2母液，SW620細胞生長狀態良好的情況下，分別加入終濃度為0、50、100、200 μM的CoCl2溶液，細胞培養2～3天。

1.2.3 Western blot 加入適量RIPA溶液提取細胞的總蛋白，利用BCA方法測定蛋白濃度；樣品加入上樣緩沖液(Loading Buffer(5×)，在100 ℃金屬浴煮樣5 min；每孔加入15 μg蛋白進行SDS-PAGE電泳；再將蛋白轉移至PVDF膜，在5 %脫脂奶粉溶液中封閉1 h；一抗孵育(GAPDH 1∶20000；HIF1-α 1∶1000；SLC2A3 1∶1000；APEH 1∶1000)，4℃過夜；二抗孵育(羊抗兔IgG 1∶3000；羊抗鼠IgG 1∶5000)，室溫孵育45 min；ECL化學發光顯影并拍照。

1.2.4 mRNA-seq和Ribo-seq 實驗的樣品準備 mRNA-seq和Ribo-seq實驗由廣州表觀生物科技有限公司完成。1)mRNA-seq文庫構建：用Trizol試劑提取總RNA，mRNA-seq根據VAHTSTM Stranded mRNA-seq Library Prep Kit for Illumina v2 進行制備。2)Ribo-seq文庫構建：0 μM和200 μMCoCl2處理的SW620細胞，含放線菌酮(Cycloheximide)PBS溶液洗滌，刮下細胞轉移到離心管中離心5 min，去上清，重復一次。細胞裂解物先經過核酸消化，將ribosome bounding RNA與游離mRNA分離。然后富集核糖體-mRNA復合物，得到RPFs.去除rRNA并用PAGE膠純化后，得到目的RNA片段。在兩端加5’和3’接頭進行末端修飾。反轉錄擴增得到的cDNA進行PAGE膠純化。

1.2.5 測序方法及數據分析 1)將上述方法準備好的樣品，利用Bioptic Qsep100 Analyzer對mRNA-seq和Ribo-seq文庫大小分布進行質檢。2)測序: 在Illumina NovaSeq 6000平臺PE150進行測序。3)數據分析：首先對原始數據進行質量控制處理，獲取Clean data，利用 Fast QC 分析質量分布，堿基含量分布，重復測序片段比例等測序的數據質量；利用 hisat2 軟件將過濾后的mRNA-seq和Ribo-Seq的reads數據與參考基因組比對；使用htseq-count對mRNA-seq上的reads數進行統計，使用feature Counts軟件計算ORF區內Ribo-seq水平的reads數，將其轉換為FPKM (Fragments Per Kilobase Million)值；使用 DEseq2 軟件通過log2FC>0條件篩選差異翻譯效率基因；篩選出的差異翻譯效率基因采用 Fisher 檢驗計算篩選出顯著性GO和KEGG.

1.2.6 RNA提取和RT-qPCR 根據試劑盒說明書，提取總細胞RNA并逆轉錄成cDNA.實時熒光定量PCR儀反應程序為:95 ℃ 10 min，1個循環；95 ℃ 60 s，60 ℃ 30 s，40個循環。GAPDH作為內參，用2-ΔΔCt法計算相對基因表達水平，每個樣品三個重復孔(所有使用的引物如表1所示)。

表1 引物序列信息

1.2.7 統計學分析 通過Origin 9.0統計繪圖軟件進行數據處理，采用獨立樣本t檢驗比較兩組間統計學差異，*p<0.05，**p<0.01，***p<0.001.

2 結果

2.1 CoCl2化學處理SW620細胞模擬低氧條件

CoCl2化學處理是一種常用的化學模擬低氧處理試劑，CoCl2是羥化酶活性的競爭性抑制劑，可以抑制低氧因子HIF1蛋白的羥基化，減少HIF1與pVHL結合來達到穩定HIF1蛋白目的[11,12]。在本實驗中，我們分別利用0 μM, 50 μM, 100 μM, 200 μM的CoCl2處理SW620細胞，處理后的細胞利用Western blot技術檢測不同條件下HIF1-α的蛋白含量，實驗結果顯示SW620細胞隨CoCl2濃度升高，HIF1-α蛋白的表達逐漸增加，在200 μM表達量最高(圖1)。因此，我們選用了200 μMCoCl2處理的SW620細胞進行mRNA-seq和Ribo-seq實驗。

圖1 CoCl2處理SW620細胞HIF1-α蛋白含量

2.2 mRNA-seq和Ribo-seq整體數據特征分析

為了比較低氧條件對SW620細胞的轉錄水平和翻譯水平的影響，我們將200 μM的CoCl2處理的SW620細胞分別分成兩組，一組提取總RNA，用于mRNA-seq采集數據，一組用放線菌酮處理，按照Ribo-seq實驗流程獲得采集數據。數據顯示，在mRNA-seq中，0 μM和200 μMCoCl2處理的樣本的Total Read分別為86226028和93874676，Q30_rate分別為0.978和0.978；Ribo-seq實驗中，Total Read分別為95491377和92705635，Q30_rate分別為0.954和0.967(表2)。這些數據顯示mRNA-seq和Ribo-seq采集的數據總量和質量可滿足下游分析。

表2 mRNA-seq和Ribo-seq的樣本序列質控統計結果

在Ribo-seq數據中，我們發現0 μM和200 μMCoCl2處理的樣本RFs主要分布在25～34 nt(圖2a, 2c)。利用riboWaltzR對25～34 nt的評估RFs在RNA上的3堿基的周期性分布，如圖2b, 2d看出，橫坐標為距離start-codon or stop-codon 的距離，縱坐標為比對上的reads.0 μM和200 μMCoCl2處理的SW620兩組的25～34 nt的RFs符合3 nt特征，這也說明本研究中Ribo-seq數據質量可靠。

圖2 SW620 在0 μM和200 μM CoCl2處理下核糖體的RFs長度分布(a，c)以及25～34 nt長度的RFs在翻譯起始位點后50 nt和終止密碼子前50 nt的3 nt特征(b，d)

2.3 CoCl2化學模擬低氧處理對基因轉錄水平及翻譯水平的影響

為了分析每個樣本中的基因轉錄水平與翻譯水平之間的關系，我們對每個樣本的mRNAs和核糖體結合mRNAs的豐度采用TPM (Transcripts Per Kilobase Million)方法來進行均一化處理，然后計算每個基因的log10(TPM+1)值，通過相關性分析發現，兩個樣本中基因的轉錄水平和翻譯水平的皮爾森相關系數(Pearson correlation)分別為0.77和0.76(圖3a,3b)，說明在我們這兩個樣本基因的表達水平和翻譯水平具有強的正相關性。

圖3 0 μM和200 μM CoCl2處理下轉錄和翻譯水平上log10(FTM+1)關系(a,b)以及200 μM CoCl2處理后轉錄或翻譯水平的差異表達基因數量(c)和轉錄和翻譯水平上的基因關系(d)

根據mRNA-seq數據顯示在0 μM CoCl2處理檢測到24 613個基因，200 μM CoCl2處理檢測到24 321個基因。與0 μM CoCl2處理相比，200 μM CoCl2處理后的細胞有118個基因上調，119個基因下調。根據Ribo-seq數據我們發現0 μM CoCl2處理下檢測到6 829個基因，200 μM CoCl2處理后檢測到4 941個基因。與0 μM CoCl2處理相比，200 μM CoCl2處理后的細胞中有254個基因翻譯水平上調，646個基因翻譯水平下調(圖3c)。其中，4個基因在轉錄水平和翻譯水平都上調，16個基因在轉錄水平和翻譯水平都下調(圖3d)。上調的基因分別為NDRG1、CXCL8、MIR22HG、SLC2A3，下調的基因分別為PTK7FIBP、PRSS2、CHCHD10、APEH、BEX3、RAB13、DNPH1、ATF6B、FAM210B、TSTD1、STRA6、SCRN2、IGFBP6、VWA1、TMEM176A.KEGG分析結果顯示，低氧條件下，翻譯效率差異基因主要富集在10個信號通路：癌癥中蛋白聚糖，中性粒細胞胞外陷阱形成，酗酒，病毒致癌作用，系統性紅斑狼瘡，細胞衰老，補體系統，FOXO信號通路，ECM受體相互作用，膀胱癌(圖4d)。

圖4 差異翻譯效率基因GO細胞組分(a)、分子功能(b)和生物過程(c)分析結果以及差異翻譯效率基因KEGG排名前10的富集結果(d)

2.4 RT-qPCR和WB驗證基因和蛋白表達

為了探究預測的mRNA-seq和Ribo-seq數據是否與蛋白質水平相匹配。我們從轉錄水平和翻譯水平變化一致的基因中，挑選了部分基因進行了RT-qPCR(NDRG1、SLC2A3、APEH、TMEM176A、FAM210B)和Western blot(SLC2A3、APEH)實驗。RT-qPCR結果顯示，在低氧條件下，NDRG1和SLC2A3在轉錄水平上調(圖5a, b)，APEH、TMEM176A、FAM210B在低氧條件下轉錄水平下調，但SLC2A3不顯著。Western blot結果顯示SLC2A3蛋白質表達水平沒有顯著變化，APEH蛋白質水平降低(圖5c)。這些結果顯示暗示有的基因的蛋白水平可能與基因的轉錄水平和翻譯水平不一致。

圖5 0 μM和200 μM CoCl2處理下HIF1-a、NDRG1、SLC2A3基因表達(a),APEH、TMEM176A、FAM210B基因表達(b)和APEH、SLC2A3蛋白表達(c)情況

3 討論

低氧是腫瘤組織的重要特征之一，但是低氧如何調控基因表達的研究仍需要深入研究。本項目中，我們利用mRNA-seq和Ribo-seq技術分析了CoCl2化學模擬低氧處理SW620細胞后基因在轉錄水平和翻譯水平變化，并對基因的翻譯效率進行了分析。據我們所知，這項研究是率先利用mRNA-seq和Ribo-seq技術研究低氧條件下SW620基因表達特征的報告，這些數據揭示了低氧條件下SW620的基因表達模式，將有助于進一步認識結直腸癌細胞的特征。

多種方法可以實現細胞的體外低氧培養，包括低氧培養箱，低氧培養裝置，化學處理等。CoCl2還可以模擬建立體外Caco-2結直腸癌細胞和MCF-7乳腺癌細胞的缺氧模型[13]。在本研究中，對HIF1蛋白的Western blot分析發現CoCl2處理SW620細胞也可使HIF1-a蛋白含量升高，也證明了CoCl2處理SW620細胞模擬低氧環境是可行的。

本項目中，我們分別采集了CoCl2處理的SW620細胞的轉錄組學(mRNA-seq)和翻譯組學(Ribo-seq)數據。每個處理組的轉錄組學totalread超過100M，QC30超過了95%，可以判斷兩組轉錄組學數據量和質量可滿足下游分析的。在翻譯組學數據中，RFs長度主要在25～34 nt之間，符合預期的核糖體保護的RNA片段(RFs)的長度在28 nt左右的特征。翻譯過程中，核糖體在可翻譯的RNA以密碼子長度(3 nt)為單位移動，完成一次肽段延伸[14]。因為核糖體在密碼子第一堿基位置停留時間最長，因此RFs 比對位置對應密碼子第一個堿基的比例通常最高，而第二第三個堿基都偏低，RFs片段應該在翻譯的RNA的開放閱讀框(ORF)中呈現數據呈現“高-低-低”三堿基周期模式，即3 nt特征。我們的Ribo-seq實驗中的25～34 nt 的RFs也符合這些特征。由此，我們認為本項目中的轉錄組學(mRNA-seq)和翻譯組學(Ribo-seq)數據質量好。

組學技術發展有助于我們對基因表達在不同水平的進行研究，為我們了解基因表達調控提供思路。本項目中，我們利用兩組轉錄組學和翻譯組學中的RNA豐度進行了相關性分析，研究發現兩組中的基因在轉錄水平和翻譯水平都呈強相關性。基因差異表達分析發現翻譯水平上的變化基因數目大于轉錄水平上變化的基因。同時，我們發現翻譯水平被抑制的基因多于增強的基因。這暗示著低氧條件可抑制基因的表達，這與Lei等人在玉米中的研究發現低氧抑制基因翻譯的現象相符合[15]。

我們對翻譯水平和轉錄水平上變化一致的基因進行了統計，這暗示著低氧條件可以對這些基因的轉錄水平和翻譯水平都存在著調控。NDRG1 (N-myc下游調節基因1)是一種與腫瘤侵襲和遷移密切相關的多功能基因，其通過G1/S期阻滯細胞周期來抑制Caco2細胞的增殖，且NDRG1過度表達時，侵襲和遷移的強度降低[16]。多個研究組表明該蛋白受HIF1-α的上調[17,18]。SLC2A3是編碼葡萄糖轉運蛋白的溶質載體家族，也稱為GLUT3，介導葡萄糖從細胞外轉運到細胞內，在癌細胞的代謝重編程中起關鍵作用[19]。SLC2A3基因的上調與結直腸癌患者的OS和DFS降低有關[20]。Taniguchi-Ponciano等在新冠肺炎重癥患者中發現HIF1a及其下游與碳水化合物代謝相關的SLC2A3基因表達都增加[21]。但是在本研究中，我們發現SLC2A3的轉錄水平受到低氧誘導表達但不顯著，蛋白水平也沒有顯著變化。這可能是由于不同的細胞組織低氧對SLC2A3誘導表達的程度不同導致的，也有可能是由于不同細胞對SLC2A3蛋白穩定性調控不同。APEH(N-酰肽水解酶)是一種細胞溶質酶，又稱氨基酸釋放酶(AARE)，屬于絲氨酸肽酶的脯氨酰寡肽酶家族，在泛素系統以及蛋白質分解過程中起重要作用[22]。Palmieri等發現，與健康對照相比，在阿茨海默癥(AD)樣品中APEH活性顯著降低[23]。TMEM176A(人跨膜蛋白176A)的表達誘導H1299和H23細胞G2/M期阻滯導致細胞凋亡，并抑制集落形成、細胞增殖、遷移和侵襲，其通過抑制ERK1/2信號來抑制肺癌生長[24]。FAM210B一種線粒體外膜蛋白的癌癥進展抑制基因，FAM210B的低表達導致線粒體呼吸能力增加和糖酵解減少，從而激活EMT并增強遷移和侵襲特性，FAM210B缺失與體內和體外存活率降低和轉移增強顯著相關[25]。

翻譯效率可以反映基因在翻譯水平上的精細調節。為了研究低氧對基因表達調控機制，我們分析了低氧條件下的基因翻譯效率的差異基因。KEGG分析發現低氧條件下最主要影響就是癌癥中蛋白聚糖，性粒細胞胞外陷阱形成，酗酒，病毒致癌作用，系統性紅斑狼瘡，細胞衰老，補體系統，FOXO信號通路，ECM-受體相互作用，膀胱癌信號通路中的基因翻譯效率。其中蛋白聚糖信號通路是富集程度最高的一個通路，該通路中的26個基因的翻譯效率受低氧調控，之前的研究表明相關的基因有助于腫瘤生長基質的形成，影響細胞—細胞和細胞—基質相互作用以及腫瘤細胞的信號轉導，還調節腫瘤細胞的表型和腫瘤基質血管生成[26，27]。Zhang等人在卵巢癌細胞中發現低氧可調節蛋白聚糖信號通路[28]。本研究揭示了低氧調控相關的基因的表達可能是涉及翻譯水平的調控。

本文聯合了轉錄組學和翻譯組學對低氧條件下結直腸癌細胞SW620的轉錄水平進行了分析，發現了一系列在轉錄水平和翻譯水平變化一致的基因，更重要的是通過比較不同條件下基因翻譯效率發現，低氧可調控多個信號通路的基因翻譯效率，為闡述低氧調控基因表達機制提供新的角度。