引用本文： 楊琪翔, 史平玲, 盧聰, 宋昊, 宋宗明. 光感受器661細胞系早期低氧損傷轉錄組測序的生物信息學分析. 中華眼底病雜志, 2021, 37(3): 214-223. doi: 10.3760/cma.j.cn511434-20201009-00483 復制
光感受器細胞低氧損傷是高海拔視網膜病變、視網膜中央或分支動脈阻塞、黃斑視網膜脫離、糖尿病視網膜病變和青光眼等缺血缺氧性眼底疾病的共同致病過程[1-2]。661W細胞系因具有視錐細胞標記，廣泛應用于光感受器細胞的相關研究[3]。目前研究已證實，661W細胞在低氧條件下會發生時間依賴性死亡[4]。另有研究表明，凋亡參與了低氧因素介導的細胞損傷[5-8]，而自噬的增加起到了保護作用[4,9]。但目前低氧對661W細胞相關分子機制的研究尚少，且缺乏高通量轉錄組測序數據。本研究利用轉錄組測序技術和生物信息學分析初步探討了低氧對661W細胞基因表達水平和相關通路的早期影響，以期為低氧相關性眼底疾病的研究提供數據支持，也為早期干預缺血缺氧相關眼底疾病提供潛在分子基礎及藥物靶點。現將結果報道如下。
1 材料和方法
1.1 細胞培養
小鼠661W細胞系由溫州醫科大學附屬眼視光醫院谷峰教授饋贈。細胞培養在含質量濃度為5.5 g/L葡萄糖、100 U/ml青霉素、100 μg/ml鏈霉素和含10%胎牛血清（杭州四季青生物工程材料有限公司）的改良Eagle高糖培養基（北京索寶生物科技有限公司）中。從液氮中取出第6代661W細胞系后經過消化傳代4次恢復細胞活性，用于后續實驗。將細胞分為低氧處理組和常氧對照組。低氧處理組細胞置于37 ℃、體積分數為1%O2和5% CO2培養箱中培養；常氧對照組細胞置于37 ℃、體積分數為5% CO2培養箱中培養。
1.2 轉錄組測序
低氧處理組、常氧對照組每組設3個生物學重復，為反映樣本間基因表達的相關性，計算每兩個樣本之間所有基因表達量的Pearson相關系數，以熱圖表示，相關系數越高，基因表達水平越相似。樣本收集后由華大基因公司武漢分部進行RNA提取、質檢和文庫構建。采用Trizol法提取細胞總RNA。使用全自動核酸片段分析系統（Fragment Analyzer）和NanoDropTM分光光度計檢測樣品質量。滿足建庫要求后，利用BGISEQ-500平臺進行轉錄組測序。測序儀單次測序所得到的堿基序列記為原始數據。通過去除測序原始數據中低質量（質量值＜10的堿基占該原始數據總堿基數的比例＞20%）、接頭污染（含接頭的原始數據）以及未知堿基N含量過高（未知堿基N含量＞5%）的原始數據得到過濾數據（測序儀單次測序所得到的高質量堿基序列）。使用HISAT軟件將過濾數據比對到參考基因組序列GCF_000001635.26_GRCm38.p6（http://www.ncbi.nlm.nih.gov/assembly/GCF_000001635.26），使用Bowtie2將過濾數據比對到參考基因序列GCF_000001635.26，得到比對結果及各組樣品中每個基因的表達量。
1.3 差異表達基因（DEG）篩選
使用Dr.Tom在線分析軟件篩選低氧處理組、常氧對照組顯著DEG；篩選標準為差異表達倍數（FC）≥ 1.5且校正后q值≤0.05[10]。基因表達水平以每百萬原始數據中來自于某基因每千堿基長度的原始數據（FPKM）表示。FPKM≥1認為該基因有表達；篩選常氧對照組、低氧處理組表達的基因目錄，繪制韋恩圖和火山圖。
1.4 生物信息學分析方法
使用R軟件對數據進行生物信息學分析，包括樣本間的相關性分析、DEG篩選、基因注釋（GO）分析、信號通路分析等。使用Bioconducter軟件包pheatmap繪制DEG的聚類熱圖，紅色越重表示DEG表達水平越高，綠色越重表示DEG表達水平越低。利用Dr.Tom平臺對DEG進行GO功能富集分析和基于京都基因與基因組百科全書（KEGG）的信號通路富集分析。對測序結果進行基因集富集分析（GSEA）時，保留在GO條目中基因在15～500的閉區間，并以校正后的歸一化富集指數的絕對值≥1.00、名義P值≤0.05、矯正多重假設檢驗≤0.25的標準篩選有意義的基因集；選擇富集程度最高的20個條目繪制富集柱狀圖。DEG行蛋白-蛋白互作網絡（PPI）分析時，使用STRING 11.0數據庫查找506個DEG關聯蛋白，根據現有數據庫和實驗驗證的證據篩選出中度可信度≥0.40的192個關聯節點，應用Cytoscape 3.5.1軟件網絡分析工具繪制PPI圖（圖1）；并采用CytoHubba插件最大集團中心度算法篩選得到關鍵基因，根據計算結果對節點進行排序。
圖1
差異表達基因蛋白-蛋白互作網絡圖。節點越大表示關聯節點數量越多
1.5 逆轉錄聚合酶鏈反應（RT-PCR）檢測mRNA表達水平
篩選有意義的GO條目中與低氧相關的基因并按Log2（FC）值排序。FPKM值較大代表該DEG mRNA表達水平較高。根據排序結果選擇排名靠前且FPKM值較大的DEG進行RT-PCR驗證。低氧處理4、8、12 h后采用RT-PCR檢測661W細胞mRNA表達水平。根據試劑操作說明，使用RNAiso Plus（日本TaKaRa公司）提取兩組細胞總RNA，并用PrimeScriptTM RT reagent kit with gDNA 試劑盒（日本TaKaRa公司）逆轉錄得到cDNA。使用Primer premier 5.0軟件設計相關引物，內參照為β-肌動蛋白（β-actin）（表1）。使用TB Green Premix EX Taq Ⅱ（日本TaKaRa公司）試劑盒進行RT-PCR檢測。每個擴增體系體積為20 μl，其中包括10 μl 2×SYBR Premix Ex Taq Ⅱ、2 μl 正向引物、2 μl反向引物、5.6 μl cDNA模板、0.4 μl 50×ROX Reference Dye。擴增反應條件為95 ℃ 30 s；95 ℃ 5 s，60 ℃ 30 s，40個循環；95 ℃ 15 s，60 ℃ 1 min，95 ℃ 15 s。采用2-ΔΔCT法計算目的基因mRNA相對表達情況。
表1
基因擴增和測序引物
1.6 流式細胞術檢測細胞凋亡
將661W細胞置于低氧和常氧條件下分別培養24、48、72 h。收集細胞培養基中的未貼壁細胞，磷酸鹽緩沖液（PBS）洗滌細胞2次，每次都收集洗滌過的PBS。用不含乙二胺四乙酸的胰酶（中國索萊寶公司）消化細胞并進行細胞計數，每個生物學重復孔分別收集約1~6×105個細胞，加入500 μl 緩沖液懸浮細胞后，加入5 μl Annexin Ⅴ FITC-A和5 μl DNA 7-AAD-A混勻（凋亡試劑盒購于武漢伊萊瑞特生物科技股份有限公司）。室溫避光孵育30 min后用流式細胞儀檢測不同時間點細胞的凋亡情況，將早期、晚期凋亡細胞合并統計凋亡率。
1.7 倒置活細胞拍照及活/死細胞熒光染色檢測細胞狀態
將661W細胞均勻鋪在直徑35 mm的培養皿中，每個皿約有1~6×105個細胞。放入細胞培養箱靜置12 h讓細胞牢固貼壁。隨機分組并將細胞按分組放入低氧培養箱和常氧培養箱中，培養條件如前所述。觀察24、48、72 h時細胞的生長狀態，用倒置顯微鏡拍攝圖片。在每個時間點用Hoechst 33342、鈣黃綠素AM和碘化丙啶（美國invitrogen生命技術公司）分別熒光標記細胞核（藍色）、活細胞（綠色）和死細胞（紅色），并用熒光顯微鏡拍攝圖片。
1.8 統計學方法
采用SPSS 26.0軟件行統計學分析處理。低氧處理組和常氧對照組測序的組間基因表達差異水平行t檢驗；樣本間重復性行Pearson相關性系數分析；低氧處理組和常氧對照組細胞的DEG mRNA表達水平比較和流式細胞術檢測細胞凋亡情況行曼-惠特尼檢驗。檢驗標準α=0.05。P＜0.05為差異有統計學意義。
2 結果
與常氧對照組比較，低氧處理2 h時，低氧處理組細胞中N-myc下游調控基因-1（Ndrg1）mRNA相對表達水平升高，差異有統計學意義（Z =－2.309，P=0.029）；而血管內皮生長因子（VEGF）A、Mt1 mRNA相對表達水平差異無統計學意義（Z =－1.852、－1.414，P=0.133、0.229）。低氧處理4 h時，低氧處理組細胞中Ndrg1、VEGFA、Mt1 mRNA相對表達水平較常氧對照組升高2倍，差異有統計學意義（Z =－2.121、－2.309、－2.121，P=0.034、0.021、0.034）。低氧處理8 h時，兩組細胞中Ndrg1、VEGFA、Mt1 mRNA相對表達水平的差異進一步擴大（Z =－1.964、－1.964、－2.121，P=0.050、0.050、0.034）（圖2）。本次測序選擇4 h作為低氧處理時間，收集樣本后送往華大基因公司武漢分部進行轉錄組測序。
圖2
低氧處理2、4、8 h后低氧處理組、常氧對照組細胞中血管內皮生長因子A（VEGFA）、Mt1、N-myc下游調控基因-1（Ndrg1）mRNA相對表達量比較。2A~2C分別示VEGFA、Mt1、Ndrg1。*P＜0.05
661W細胞RNA建庫檢測結果顯示，樣本片段大小和濃度合格（圖3，表2），滿足建庫要求。
圖3
低氧處理組、常氧對照組樣本間基因表達量Pearson相關系數
表2
661W細胞RNA建庫測序結果
DEG分析，共篩選出506個DEG，其中上調基因459個，下調基因47個（圖4A，4B）。常氧對照組、低氧處理組同時表達的基因為10 947個，僅在常氧對照組組、低氧處理組表達的基因分別為126、118個（圖4C）。
圖4
差異表達基因（DEG）分析結果。4A示聚類熱圖，紅色越深表示DEG表達水平越高，綠色越深表示DEG表達水平越低。4B示火山圖，紅色和橘紅色表示上調基因，綠色表示下調基因，深灰色表示非DEG；紅色和深綠色表示差異表達倍數≥2，淺綠色和橘紅色表示差異表達倍數≥1.5。4C示韋恩圖
GO功能富集分析結果顯示，細胞成分方面，細胞質、胞質溶液、磷酸丙酮酸水合酶復合物、6-磷酸果糖激酶復合物、前膠原-脯氨酸4-雙加氧酶復合物在胞質區的表達、神經元投射和有被小窩被顯著富集（q＜0.05）；分子功能方面，顯著富集的是葡萄糖氧化還原反應相關激酶催化活性、DNA結合轉錄阻遏物活性、RNA聚合酶Ⅱ特異性轉錄調節活性、L-抗壞血酸及激酶結合（q＜0.05）；生物過程方面，細胞對低氧反應、糖酵解、細胞增生負調控、節律及凋亡等被顯著富集（q＜0.05）（圖5）。GSEA分析共富集17 978個基因在GO條目中，富集程度較高的集中于細胞對缺氧的反應、糖酵解、糖原代謝和合成、肽交聯和鈣離子內流的正調控（表3）。
圖5
基因注釋（GO）功能顯著富集分析。條柱代表富集GO條目的q值，柱線越高表示q值越小，富集越顯著；折線代表富集GO條目下的靶基因個數。1. 胞質溶液；2. 細胞質；3. 磷酸丙酮酸水合酶復合物；4. 6磷酸果糖激酶復合物；5. 前膠原脯氨酸4雙加氧酶復合物；6. 神經元投射；7. 有被小窩；8. 谷氨酸能突觸；9. 細胞質核周區域；10. 突觸膜；11. 細胞皮層區；12. 核；13. 細胞表面；14. 軸突；15. 髓鞘；16. 海馬苔蘚纖維對CA3突觸；17. 神經元投射終點；18. 糖原顆粒；19. 晚期內體；20. VCB復合體；21. 雙加氧酶活性；22. 蛋白質結合；23. RNA聚合酶特異性胞皮層區；24. L抗壞血酸結合；25. 催化活性；26. 葡萄糖結合；27. 蛋白質均二聚活性；28. 磷酸丙酮酸水合酶活性；29. 氧化還原酶活性；30. CXC趨化因子受體活性； 31. 激酶結合；32. 轉錄共受體結合；33. 轉錄因子結合；34. 6磷酸果糖激酶活性；35. 肽脯氨酸雙加氧酶活性；36. 組蛋白脫甲基酶活性；37. 磷酸酶活性；38. RNA聚合酶調控特異性DNA結合；39. β-catenin結合；40. 酶結合；41. 細胞對缺氧的反應；42. 糖酵解過程；43. 對缺氧的反應；44. 細胞增生的負調控；45. 常規糖酵解；46. 肽基脯氨酸羥基化；47. 冷致產熱；48. 調節細胞增生；49. 糖原生物合成過程；50. 1,6-二磷酸果糖代謝過程；51. RNA聚合酶Ⅱ對轉錄的負調控；52. 基因表達的晝夜調節；53. 葡萄糖穩態；54. 凋亡過程；55. 碳水化合物代謝；56. 甲基乙二醛的生物合成過程；57. 6-磷酸果糖的糖酵解過程；58. 紅細胞穩態；59. 光周期晝夜節律；60. 節律
表3
基因集富集分析結果
KEGG信號通路富集分析結果顯示，糖酵解與糖異生、果糖代謝、磷酸戊糖途徑以及淀粉和蔗糖代謝等糖代謝途徑被顯著富集。主要富集的通路為缺氧誘導因子（HIF）-1α信號通路、晝夜節律、FoxO信號通路、胰島素信號通路和腺苷酸活化蛋白激酶（AMPK）信號通路（q＜0.05）（圖6）。
圖6
京都基因與基因組百科全書（KEGG）的信號通路富集分析結果。條柱為富集的KEGG 通路q值，柱線越高表示q值越小，富集越顯著；折線為代表富集KEGG通路下靶基因個數。1. 缺氧誘導因子-1信號通路；2. 糖酵解和糖異生；3. 碳代謝；4. 甲烷代謝；5. 果糖和甘露糖代謝；6. 晝夜節律；7. 次生代謝產物的合成；8. 抗生素的生物合成；9. 不同環境微生物代謝；10. 氨基酸的生物合成；11. 磷酸戊糖途徑；12. 半乳糖代謝；13.淀粉和蔗糖代謝；14. FoxO信號通路；15. 代謝途徑；16. 胰島素信號通路；17. 腺苷酸活化蛋白激酶信號通路；18. 鏈霉素的生物合成；19. 氨基糖和核苷酸糖代謝；20. 新霉素卡那霉素和慶大霉素的生物合成
PPI分析結果顯示，網絡中富集程度最高的前10位關鍵基因分別是Aldoa、Aldoc、Gpi1、Hk2、Hk1、Pfkl、Pfkp、Vhl、Fbxo10、Fbxo27（圖7）。
圖7
蛋白互作網絡圖。以紅色至黃色的顏色變化代表排名位次列出10個關鍵基因節點，均一的粉紅色表示與之關聯的節點基因
RT-PCR檢測結果顯示，低氧處理組細胞中15個DEG mRNA相對表達水平較常氧對照組升高，差異有統計學意義（圖8，表4）。
圖8
低氧處理組、常氧對照組細胞中15個差異基因 mRNA表達水平比較。*P＜0.05
表4
低氧處理組、常氧對照組檢測差異表達基因mRNA表達水平
流式細胞術檢測結果顯示，低氧處理24、48、72 h時，低氧處理組細胞的凋亡率分別為19.7%、42.9%和64.2%，明顯高于常氧對照組的4.5%、7.9%和25.3%（P＜0.05），并具有一定的時間相關性（圖9）。
圖9
流式細胞術檢測低氧處理組和常氧對照組細胞24、48、72 h凋亡率。9A示流式細胞術檢測圖。Q1（Annexin V-FITC）-/7-AAD+區域為壞死細胞，也包括少數晚期凋亡細胞和機械損傷細胞；Q2（Annexin V-FITC）+/7-AAD+區域為晚期凋亡細胞；Q3（Annexin V-FITC）-/7-AAD-區域為活細胞；Q4（Annexin V-FITC）+/7-AAD-區域為早期凋亡細胞。9B示低氧處理組和常氧對照組細胞24、48、72 h凋亡率比較
活細胞倒置圖像及活/死細胞熒光染色結果顯示，與常氧對照組比較，低氧處理組細胞在低氧處理24 h時細胞突觸減少，細胞變圓，開始出現死細胞；低氧處理48 h時細胞形態仍偏圓，死亡細胞增多；低氧處理72 h時死亡細胞數量進一步增多，細胞增生速度較慢（圖10）。
圖10
低氧處理組和常氧對照組細胞低氧處理24、48、72 h時活細胞倒置圖像及活/死細胞熒光染色像。10A~10C示常氧對照組活細胞倒置圖像；10D~10F示常氧對照組活/死細胞熒光染色像；10G~10I示低氧處理組活細胞倒置圖像；10J~10L示低氧處理組活/死細胞熒光染色像。與常氧對照組比較，低氧處理組細胞在低氧處理24 h時細胞突觸減少，細胞變圓，開始出現死細胞；低氧處理48 h時細胞形態仍偏圓，死亡細胞增多；低氧處理72 h時死亡細胞數量進一步增多，細胞增生速度較慢 標尺：200 μm
3 討論
本研究首次對低氧處理的661W細胞進行轉錄組測序分析，共篩選得到506個661W細胞低氧相關DEG，其中表達上調基因459個，表達下調基因47個。KEGG通路富集分析和GO功能富集分析結果顯示，糖代謝相關通路、HIF-1α通路、晝夜節律、FoxO通路和AMPK通路參與661W細胞對低氧的早期反應；缺氧反應、細胞增生負調控、鈣離子內流正調控和細胞凋亡也參與了這一過程。AMPK通路、FoxO通路和鈣離子內流的正調控可能參與了低氧下細胞的損傷過程[11-12]。這為研究低氧對光感受器細胞糖代謝改變、細胞增生和損傷提供了理論依據。
目前研究表明，N-myc下游調控基因-1（Ndrg1）是一種有效的轉移抑制因子，受缺氧、鐵等金屬離子、NO自由基等的調控[13]。Takita等[14]研究發現，Ndrg1b是鯉魚的錐特異性基因，并和其同源基因共同參與斑馬魚光感受器外節的形成，且參與了斑馬魚視錐和視桿細胞外節的正常發育。Hu等[15]的研究進一步證實了Ndrg1家族基因參與保護了甲基亞硝基脲誘導的神經節細胞損傷。本研究結果顯示，低氧使該基因表達水平升高，可能發揮了對神經節細胞的保護作用，但仍需進一步的研究探索。Ⅰ類膜型基質金屬蛋白酶（MT1-MMP）是基質金屬蛋白酶家族成員的關鍵成員，可以通過降解腫瘤周圍間質屏障，增強腫瘤的侵襲力，低氧條件下可大量表達[16-17]。大量研究已經表明，VEGFA表達升高參與低氧下視網膜新生血管的形成[18]，與視網膜靜脈阻塞、糖尿病視網膜病變和視網膜脫離等疾病密切相關。而MT1-MMP的表達也與VEGF的自分泌相關[19]。低氧時Mt1和VEGFA表達水平升高，說明小鼠光感受器細胞啟動了對周圍間質的酶解和新生血管的生成。本研究驗證了低氧處理4 h時 VEGFA、Mt1和Ndrg1 mRNA表達水平已明顯增高。倒置活細胞圖像及活/死細胞熒光染色也顯示早期低氧處理48 h即出現了神經節細胞的大量死亡以及細胞增生能力的明顯下降，由此得到了低氧早期損傷出現的時間節點。
視網膜光感受器細胞是視網膜數量最多且代謝旺盛的細胞，營養和氧主要來源于脈絡膜毛細血管的血液供應[20]。目前研究認為，視網膜缺血、缺氧等疾病情況下，光感受器細胞會出現退化，造成視力損害[21]。本研究分析發現，低氧條件下，糖代謝能量代謝的改變是最主要的改變，其次是鈣離子內流正調控、細胞增生的負調控、細胞凋亡等。HIF-1α通路、FoxO通路和AMPK通路主要介導了這些改變。借助PPI蛋白互作分析，得到光感受器細胞早期缺氧的關鍵基因Aldoa、Aldoc、Gpi1、Hk2、Hk1、Pfkl、Pfkp、Vhl、Fbxo10、Fbxo27。其中，前8個關鍵基因功能集中于糖代謝、中央碳代謝等代謝相關生物過程，可能參與了糖尿病視網膜病變等能量代謝性眼底疾病的發生發展；后2個關鍵基因Fbxo10、Fbxo27與細胞凋亡及泛素化生物學過程相關，并與Rnf19a、Rnf217、Rnf126、Slah2、Egln1、Neurl2形成蛋白交互網絡。Fbxo10是UvrD家族DNA解旋酶中成員，也是進化保守的泛素連接酶基因，蛋白結構包含一個F-box域[22]，通過參與細胞凋亡的關鍵蛋白Bcl-2的降解，誘導細胞死亡過程[23-26]。另有文獻報道，Fbxo27可泛素化損傷溶酶體表面糖蛋白，通過誘導自噬來減少凋亡，產生一定的保護作用[27]。但目前關于兩者的功能研究都集中于腫瘤方向，視網膜組織或神經系統方面的研究則鮮見報道。Fbxo10和Fbxo27作為661W細胞低氧條件下的關鍵基因，可能與缺血缺氧發生時神經節細胞的損傷和保護相關，有待進一步研究探索。
本研究結果顯示，選取的15個低氧相關DEG mRNA表達水平差異均有統計學意義，表達趨勢與測序結果一致，驗證了測序結果的準確性。應當注意的是，磷酸甘油醛脫氫酶基因受低氧影響較大，低氧相關實驗時不應選擇其作為內參基因。
本研究存在的不足之處：（1）缺血缺氧性眼底疾病情況復雜，影響疾病發生發展的因素較多，單一實驗干預因素無法完全模擬體內情況，還需要進一步動物實驗結果加以驗證結論；（2）低氧干預時間4 h時已檢測到轉錄組水平的改變，生物信息學分析了早期低氧相關基因及信號通路的改變，但未觀察這些基因及信號通路在更長時間低氧干預下的情況，仍需要進一步實驗探索。
本研究結果表明，糖代謝改變、鈣離子內流正調控、細胞增生的負調控、細胞凋亡、HIF-1α通路、FoxO通路和AMPK通路參與了661W細胞的早期低氧反應。更進一步分析出Aldoa、Aldoc、Gpi1、Hk2、Hk1、Pfkl、Pfkp、Vhl、Fbxo10、Fbxo27可能在661W細胞應對低氧反應中起到關鍵作用。表達上調基因Fbxo10和Fbxo27可能參與了光感受器細胞的凋亡和自噬過程。Ndrg1、Mt1及VEGFA可作為研究缺血缺氧相關眼底疾病的潛在功能性靶基因。
光感受器細胞低氧損傷是高海拔視網膜病變、視網膜中央或分支動脈阻塞、黃斑視網膜脫離、糖尿病視網膜病變和青光眼等缺血缺氧性眼底疾病的共同致病過程[1-2]。661W細胞系因具有視錐細胞標記，廣泛應用于光感受器細胞的相關研究[3]。目前研究已證實，661W細胞在低氧條件下會發生時間依賴性死亡[4]。另有研究表明，凋亡參與了低氧因素介導的細胞損傷[5-8]，而自噬的增加起到了保護作用[4,9]。但目前低氧對661W細胞相關分子機制的研究尚少，且缺乏高通量轉錄組測序數據。本研究利用轉錄組測序技術和生物信息學分析初步探討了低氧對661W細胞基因表達水平和相關通路的早期影響，以期為低氧相關性眼底疾病的研究提供數據支持，也為早期干預缺血缺氧相關眼底疾病提供潛在分子基礎及藥物靶點。現將結果報道如下。
1 材料和方法
1.1 細胞培養
小鼠661W細胞系由溫州醫科大學附屬眼視光醫院谷峰教授饋贈。細胞培養在含質量濃度為5.5 g/L葡萄糖、100 U/ml青霉素、100 μg/ml鏈霉素和含10%胎牛血清（杭州四季青生物工程材料有限公司）的改良Eagle高糖培養基（北京索寶生物科技有限公司）中。從液氮中取出第6代661W細胞系后經過消化傳代4次恢復細胞活性，用于后續實驗。將細胞分為低氧處理組和常氧對照組。低氧處理組細胞置于37 ℃、體積分數為1%O2和5% CO2培養箱中培養；常氧對照組細胞置于37 ℃、體積分數為5% CO2培養箱中培養。
1.2 轉錄組測序
低氧處理組、常氧對照組每組設3個生物學重復，為反映樣本間基因表達的相關性，計算每兩個樣本之間所有基因表達量的Pearson相關系數，以熱圖表示，相關系數越高，基因表達水平越相似。樣本收集后由華大基因公司武漢分部進行RNA提取、質檢和文庫構建。采用Trizol法提取細胞總RNA。使用全自動核酸片段分析系統（Fragment Analyzer）和NanoDropTM分光光度計檢測樣品質量。滿足建庫要求后，利用BGISEQ-500平臺進行轉錄組測序。測序儀單次測序所得到的堿基序列記為原始數據。通過去除測序原始數據中低質量（質量值＜10的堿基占該原始數據總堿基數的比例＞20%）、接頭污染（含接頭的原始數據）以及未知堿基N含量過高（未知堿基N含量＞5%）的原始數據得到過濾數據（測序儀單次測序所得到的高質量堿基序列）。使用HISAT軟件將過濾數據比對到參考基因組序列GCF_000001635.26_GRCm38.p6（http://www.ncbi.nlm.nih.gov/assembly/GCF_000001635.26），使用Bowtie2將過濾數據比對到參考基因序列GCF_000001635.26，得到比對結果及各組樣品中每個基因的表達量。
1.3 差異表達基因（DEG）篩選
使用Dr.Tom在線分析軟件篩選低氧處理組、常氧對照組顯著DEG；篩選標準為差異表達倍數（FC）≥ 1.5且校正后q值≤0.05[10]。基因表達水平以每百萬原始數據中來自于某基因每千堿基長度的原始數據（FPKM）表示。FPKM≥1認為該基因有表達；篩選常氧對照組、低氧處理組表達的基因目錄，繪制韋恩圖和火山圖。
1.4 生物信息學分析方法
使用R軟件對數據進行生物信息學分析，包括樣本間的相關性分析、DEG篩選、基因注釋（GO）分析、信號通路分析等。使用Bioconducter軟件包pheatmap繪制DEG的聚類熱圖，紅色越重表示DEG表達水平越高，綠色越重表示DEG表達水平越低。利用Dr.Tom平臺對DEG進行GO功能富集分析和基于京都基因與基因組百科全書（KEGG）的信號通路富集分析。對測序結果進行基因集富集分析（GSEA）時，保留在GO條目中基因在15～500的閉區間，并以校正后的歸一化富集指數的絕對值≥1.00、名義P值≤0.05、矯正多重假設檢驗≤0.25的標準篩選有意義的基因集；選擇富集程度最高的20個條目繪制富集柱狀圖。DEG行蛋白-蛋白互作網絡（PPI）分析時，使用STRING 11.0數據庫查找506個DEG關聯蛋白，根據現有數據庫和實驗驗證的證據篩選出中度可信度≥0.40的192個關聯節點，應用Cytoscape 3.5.1軟件網絡分析工具繪制PPI圖（圖1）；并采用CytoHubba插件最大集團中心度算法篩選得到關鍵基因，根據計算結果對節點進行排序。
圖1
差異表達基因蛋白-蛋白互作網絡圖。節點越大表示關聯節點數量越多
1.5 逆轉錄聚合酶鏈反應（RT-PCR）檢測mRNA表達水平
篩選有意義的GO條目中與低氧相關的基因并按Log2（FC）值排序。FPKM值較大代表該DEG mRNA表達水平較高。根據排序結果選擇排名靠前且FPKM值較大的DEG進行RT-PCR驗證。低氧處理4、8、12 h后采用RT-PCR檢測661W細胞mRNA表達水平。根據試劑操作說明，使用RNAiso Plus（日本TaKaRa公司）提取兩組細胞總RNA，并用PrimeScriptTM RT reagent kit with gDNA 試劑盒（日本TaKaRa公司）逆轉錄得到cDNA。使用Primer premier 5.0軟件設計相關引物，內參照為β-肌動蛋白（β-actin）（表1）。使用TB Green Premix EX Taq Ⅱ（日本TaKaRa公司）試劑盒進行RT-PCR檢測。每個擴增體系體積為20 μl，其中包括10 μl 2×SYBR Premix Ex Taq Ⅱ、2 μl 正向引物、2 μl反向引物、5.6 μl cDNA模板、0.4 μl 50×ROX Reference Dye。擴增反應條件為95 ℃ 30 s；95 ℃ 5 s，60 ℃ 30 s，40個循環；95 ℃ 15 s，60 ℃ 1 min，95 ℃ 15 s。采用2-ΔΔCT法計算目的基因mRNA相對表達情況。
表1
基因擴增和測序引物
1.6 流式細胞術檢測細胞凋亡
將661W細胞置于低氧和常氧條件下分別培養24、48、72 h。收集細胞培養基中的未貼壁細胞，磷酸鹽緩沖液（PBS）洗滌細胞2次，每次都收集洗滌過的PBS。用不含乙二胺四乙酸的胰酶（中國索萊寶公司）消化細胞并進行細胞計數，每個生物學重復孔分別收集約1~6×105個細胞，加入500 μl 緩沖液懸浮細胞后，加入5 μl Annexin Ⅴ FITC-A和5 μl DNA 7-AAD-A混勻（凋亡試劑盒購于武漢伊萊瑞特生物科技股份有限公司）。室溫避光孵育30 min后用流式細胞儀檢測不同時間點細胞的凋亡情況，將早期、晚期凋亡細胞合并統計凋亡率。
1.7 倒置活細胞拍照及活/死細胞熒光染色檢測細胞狀態
將661W細胞均勻鋪在直徑35 mm的培養皿中，每個皿約有1~6×105個細胞。放入細胞培養箱靜置12 h讓細胞牢固貼壁。隨機分組并將細胞按分組放入低氧培養箱和常氧培養箱中，培養條件如前所述。觀察24、48、72 h時細胞的生長狀態，用倒置顯微鏡拍攝圖片。在每個時間點用Hoechst 33342、鈣黃綠素AM和碘化丙啶（美國invitrogen生命技術公司）分別熒光標記細胞核（藍色）、活細胞（綠色）和死細胞（紅色），并用熒光顯微鏡拍攝圖片。
1.8 統計學方法
采用SPSS 26.0軟件行統計學分析處理。低氧處理組和常氧對照組測序的組間基因表達差異水平行t檢驗；樣本間重復性行Pearson相關性系數分析；低氧處理組和常氧對照組細胞的DEG mRNA表達水平比較和流式細胞術檢測細胞凋亡情況行曼-惠特尼檢驗。檢驗標準α=0.05。P＜0.05為差異有統計學意義。
2 結果
與常氧對照組比較，低氧處理2 h時，低氧處理組細胞中N-myc下游調控基因-1（Ndrg1）mRNA相對表達水平升高，差異有統計學意義（Z =－2.309，P=0.029）；而血管內皮生長因子（VEGF）A、Mt1 mRNA相對表達水平差異無統計學意義（Z =－1.852、－1.414，P=0.133、0.229）。低氧處理4 h時，低氧處理組細胞中Ndrg1、VEGFA、Mt1 mRNA相對表達水平較常氧對照組升高2倍，差異有統計學意義（Z =－2.121、－2.309、－2.121，P=0.034、0.021、0.034）。低氧處理8 h時，兩組細胞中Ndrg1、VEGFA、Mt1 mRNA相對表達水平的差異進一步擴大（Z =－1.964、－1.964、－2.121，P=0.050、0.050、0.034）（圖2）。本次測序選擇4 h作為低氧處理時間，收集樣本后送往華大基因公司武漢分部進行轉錄組測序。
圖2
低氧處理2、4、8 h后低氧處理組、常氧對照組細胞中血管內皮生長因子A（VEGFA）、Mt1、N-myc下游調控基因-1（Ndrg1）mRNA相對表達量比較。2A~2C分別示VEGFA、Mt1、Ndrg1。*P＜0.05
661W細胞RNA建庫檢測結果顯示，樣本片段大小和濃度合格（圖3，表2），滿足建庫要求。
圖3
低氧處理組、常氧對照組樣本間基因表達量Pearson相關系數
表2
661W細胞RNA建庫測序結果
DEG分析，共篩選出506個DEG，其中上調基因459個，下調基因47個（圖4A，4B）。常氧對照組、低氧處理組同時表達的基因為10 947個，僅在常氧對照組組、低氧處理組表達的基因分別為126、118個（圖4C）。
圖4
差異表達基因（DEG）分析結果。4A示聚類熱圖，紅色越深表示DEG表達水平越高，綠色越深表示DEG表達水平越低。4B示火山圖，紅色和橘紅色表示上調基因，綠色表示下調基因，深灰色表示非DEG；紅色和深綠色表示差異表達倍數≥2，淺綠色和橘紅色表示差異表達倍數≥1.5。4C示韋恩圖
GO功能富集分析結果顯示，細胞成分方面，細胞質、胞質溶液、磷酸丙酮酸水合酶復合物、6-磷酸果糖激酶復合物、前膠原-脯氨酸4-雙加氧酶復合物在胞質區的表達、神經元投射和有被小窩被顯著富集（q＜0.05）；分子功能方面，顯著富集的是葡萄糖氧化還原反應相關激酶催化活性、DNA結合轉錄阻遏物活性、RNA聚合酶Ⅱ特異性轉錄調節活性、L-抗壞血酸及激酶結合（q＜0.05）；生物過程方面，細胞對低氧反應、糖酵解、細胞增生負調控、節律及凋亡等被顯著富集（q＜0.05）（圖5）。GSEA分析共富集17 978個基因在GO條目中，富集程度較高的集中于細胞對缺氧的反應、糖酵解、糖原代謝和合成、肽交聯和鈣離子內流的正調控（表3）。
圖5
基因注釋（GO）功能顯著富集分析。條柱代表富集GO條目的q值，柱線越高表示q值越小，富集越顯著；折線代表富集GO條目下的靶基因個數。1. 胞質溶液；2. 細胞質；3. 磷酸丙酮酸水合酶復合物；4. 6磷酸果糖激酶復合物；5. 前膠原脯氨酸4雙加氧酶復合物；6. 神經元投射；7. 有被小窩；8. 谷氨酸能突觸；9. 細胞質核周區域；10. 突觸膜；11. 細胞皮層區；12. 核；13. 細胞表面；14. 軸突；15. 髓鞘；16. 海馬苔蘚纖維對CA3突觸；17. 神經元投射終點；18. 糖原顆粒；19. 晚期內體；20. VCB復合體；21. 雙加氧酶活性；22. 蛋白質結合；23. RNA聚合酶特異性胞皮層區；24. L抗壞血酸結合；25. 催化活性；26. 葡萄糖結合；27. 蛋白質均二聚活性；28. 磷酸丙酮酸水合酶活性；29. 氧化還原酶活性；30. CXC趨化因子受體活性； 31. 激酶結合；32. 轉錄共受體結合；33. 轉錄因子結合；34. 6磷酸果糖激酶活性；35. 肽脯氨酸雙加氧酶活性；36. 組蛋白脫甲基酶活性；37. 磷酸酶活性；38. RNA聚合酶調控特異性DNA結合；39. β-catenin結合；40. 酶結合；41. 細胞對缺氧的反應；42. 糖酵解過程；43. 對缺氧的反應；44. 細胞增生的負調控；45. 常規糖酵解；46. 肽基脯氨酸羥基化；47. 冷致產熱；48. 調節細胞增生；49. 糖原生物合成過程；50. 1,6-二磷酸果糖代謝過程；51. RNA聚合酶Ⅱ對轉錄的負調控；52. 基因表達的晝夜調節；53. 葡萄糖穩態；54. 凋亡過程；55. 碳水化合物代謝；56. 甲基乙二醛的生物合成過程；57. 6-磷酸果糖的糖酵解過程；58. 紅細胞穩態；59. 光周期晝夜節律；60. 節律
表3
基因集富集分析結果
KEGG信號通路富集分析結果顯示，糖酵解與糖異生、果糖代謝、磷酸戊糖途徑以及淀粉和蔗糖代謝等糖代謝途徑被顯著富集。主要富集的通路為缺氧誘導因子（HIF）-1α信號通路、晝夜節律、FoxO信號通路、胰島素信號通路和腺苷酸活化蛋白激酶（AMPK）信號通路（q＜0.05）（圖6）。
圖6
京都基因與基因組百科全書（KEGG）的信號通路富集分析結果。條柱為富集的KEGG 通路q值，柱線越高表示q值越小，富集越顯著；折線為代表富集KEGG通路下靶基因個數。1. 缺氧誘導因子-1信號通路；2. 糖酵解和糖異生；3. 碳代謝；4. 甲烷代謝；5. 果糖和甘露糖代謝；6. 晝夜節律；7. 次生代謝產物的合成；8. 抗生素的生物合成；9. 不同環境微生物代謝；10. 氨基酸的生物合成；11. 磷酸戊糖途徑；12. 半乳糖代謝；13.淀粉和蔗糖代謝；14. FoxO信號通路；15. 代謝途徑；16. 胰島素信號通路；17. 腺苷酸活化蛋白激酶信號通路；18. 鏈霉素的生物合成；19. 氨基糖和核苷酸糖代謝；20. 新霉素卡那霉素和慶大霉素的生物合成
PPI分析結果顯示，網絡中富集程度最高的前10位關鍵基因分別是Aldoa、Aldoc、Gpi1、Hk2、Hk1、Pfkl、Pfkp、Vhl、Fbxo10、Fbxo27（圖7）。
圖7
蛋白互作網絡圖。以紅色至黃色的顏色變化代表排名位次列出10個關鍵基因節點，均一的粉紅色表示與之關聯的節點基因
RT-PCR檢測結果顯示，低氧處理組細胞中15個DEG mRNA相對表達水平較常氧對照組升高，差異有統計學意義（圖8，表4）。
圖8
低氧處理組、常氧對照組細胞中15個差異基因 mRNA表達水平比較。*P＜0.05
表4
低氧處理組、常氧對照組檢測差異表達基因mRNA表達水平
流式細胞術檢測結果顯示，低氧處理24、48、72 h時，低氧處理組細胞的凋亡率分別為19.7%、42.9%和64.2%，明顯高于常氧對照組的4.5%、7.9%和25.3%（P＜0.05），并具有一定的時間相關性（圖9）。
圖9
流式細胞術檢測低氧處理組和常氧對照組細胞24、48、72 h凋亡率。9A示流式細胞術檢測圖。Q1（Annexin V-FITC）-/7-AAD+區域為壞死細胞，也包括少數晚期凋亡細胞和機械損傷細胞；Q2（Annexin V-FITC）+/7-AAD+區域為晚期凋亡細胞；Q3（Annexin V-FITC）-/7-AAD-區域為活細胞；Q4（Annexin V-FITC）+/7-AAD-區域為早期凋亡細胞。9B示低氧處理組和常氧對照組細胞24、48、72 h凋亡率比較
活細胞倒置圖像及活/死細胞熒光染色結果顯示，與常氧對照組比較，低氧處理組細胞在低氧處理24 h時細胞突觸減少，細胞變圓，開始出現死細胞；低氧處理48 h時細胞形態仍偏圓，死亡細胞增多；低氧處理72 h時死亡細胞數量進一步增多，細胞增生速度較慢（圖10）。
圖10
低氧處理組和常氧對照組細胞低氧處理24、48、72 h時活細胞倒置圖像及活/死細胞熒光染色像。10A~10C示常氧對照組活細胞倒置圖像；10D~10F示常氧對照組活/死細胞熒光染色像；10G~10I示低氧處理組活細胞倒置圖像；10J~10L示低氧處理組活/死細胞熒光染色像。與常氧對照組比較，低氧處理組細胞在低氧處理24 h時細胞突觸減少，細胞變圓，開始出現死細胞；低氧處理48 h時細胞形態仍偏圓，死亡細胞增多；低氧處理72 h時死亡細胞數量進一步增多，細胞增生速度較慢 標尺：200 μm
3 討論
本研究首次對低氧處理的661W細胞進行轉錄組測序分析，共篩選得到506個661W細胞低氧相關DEG，其中表達上調基因459個，表達下調基因47個。KEGG通路富集分析和GO功能富集分析結果顯示，糖代謝相關通路、HIF-1α通路、晝夜節律、FoxO通路和AMPK通路參與661W細胞對低氧的早期反應；缺氧反應、細胞增生負調控、鈣離子內流正調控和細胞凋亡也參與了這一過程。AMPK通路、FoxO通路和鈣離子內流的正調控可能參與了低氧下細胞的損傷過程[11-12]。這為研究低氧對光感受器細胞糖代謝改變、細胞增生和損傷提供了理論依據。
目前研究表明，N-myc下游調控基因-1（Ndrg1）是一種有效的轉移抑制因子，受缺氧、鐵等金屬離子、NO自由基等的調控[13]。Takita等[14]研究發現，Ndrg1b是鯉魚的錐特異性基因，并和其同源基因共同參與斑馬魚光感受器外節的形成，且參與了斑馬魚視錐和視桿細胞外節的正常發育。Hu等[15]的研究進一步證實了Ndrg1家族基因參與保護了甲基亞硝基脲誘導的神經節細胞損傷。本研究結果顯示，低氧使該基因表達水平升高，可能發揮了對神經節細胞的保護作用，但仍需進一步的研究探索。Ⅰ類膜型基質金屬蛋白酶（MT1-MMP）是基質金屬蛋白酶家族成員的關鍵成員，可以通過降解腫瘤周圍間質屏障，增強腫瘤的侵襲力，低氧條件下可大量表達[16-17]。大量研究已經表明，VEGFA表達升高參與低氧下視網膜新生血管的形成[18]，與視網膜靜脈阻塞、糖尿病視網膜病變和視網膜脫離等疾病密切相關。而MT1-MMP的表達也與VEGF的自分泌相關[19]。低氧時Mt1和VEGFA表達水平升高，說明小鼠光感受器細胞啟動了對周圍間質的酶解和新生血管的生成。本研究驗證了低氧處理4 h時 VEGFA、Mt1和Ndrg1 mRNA表達水平已明顯增高。倒置活細胞圖像及活/死細胞熒光染色也顯示早期低氧處理48 h即出現了神經節細胞的大量死亡以及細胞增生能力的明顯下降，由此得到了低氧早期損傷出現的時間節點。
視網膜光感受器細胞是視網膜數量最多且代謝旺盛的細胞，營養和氧主要來源于脈絡膜毛細血管的血液供應[20]。目前研究認為，視網膜缺血、缺氧等疾病情況下，光感受器細胞會出現退化，造成視力損害[21]。本研究分析發現，低氧條件下，糖代謝能量代謝的改變是最主要的改變，其次是鈣離子內流正調控、細胞增生的負調控、細胞凋亡等。HIF-1α通路、FoxO通路和AMPK通路主要介導了這些改變。借助PPI蛋白互作分析，得到光感受器細胞早期缺氧的關鍵基因Aldoa、Aldoc、Gpi1、Hk2、Hk1、Pfkl、Pfkp、Vhl、Fbxo10、Fbxo27。其中，前8個關鍵基因功能集中于糖代謝、中央碳代謝等代謝相關生物過程，可能參與了糖尿病視網膜病變等能量代謝性眼底疾病的發生發展；后2個關鍵基因Fbxo10、Fbxo27與細胞凋亡及泛素化生物學過程相關，并與Rnf19a、Rnf217、Rnf126、Slah2、Egln1、Neurl2形成蛋白交互網絡。Fbxo10是UvrD家族DNA解旋酶中成員，也是進化保守的泛素連接酶基因，蛋白結構包含一個F-box域[22]，通過參與細胞凋亡的關鍵蛋白Bcl-2的降解，誘導細胞死亡過程[23-26]。另有文獻報道，Fbxo27可泛素化損傷溶酶體表面糖蛋白，通過誘導自噬來減少凋亡，產生一定的保護作用[27]。但目前關于兩者的功能研究都集中于腫瘤方向，視網膜組織或神經系統方面的研究則鮮見報道。Fbxo10和Fbxo27作為661W細胞低氧條件下的關鍵基因，可能與缺血缺氧發生時神經節細胞的損傷和保護相關，有待進一步研究探索。
本研究結果顯示，選取的15個低氧相關DEG mRNA表達水平差異均有統計學意義，表達趨勢與測序結果一致，驗證了測序結果的準確性。應當注意的是，磷酸甘油醛脫氫酶基因受低氧影響較大，低氧相關實驗時不應選擇其作為內參基因。
本研究存在的不足之處：（1）缺血缺氧性眼底疾病情況復雜，影響疾病發生發展的因素較多，單一實驗干預因素無法完全模擬體內情況，還需要進一步動物實驗結果加以驗證結論；（2）低氧干預時間4 h時已檢測到轉錄組水平的改變，生物信息學分析了早期低氧相關基因及信號通路的改變，但未觀察這些基因及信號通路在更長時間低氧干預下的情況，仍需要進一步實驗探索。
本研究結果表明，糖代謝改變、鈣離子內流正調控、細胞增生的負調控、細胞凋亡、HIF-1α通路、FoxO通路和AMPK通路參與了661W細胞的早期低氧反應。更進一步分析出Aldoa、Aldoc、Gpi1、Hk2、Hk1、Pfkl、Pfkp、Vhl、Fbxo10、Fbxo27可能在661W細胞應對低氧反應中起到關鍵作用。表達上調基因Fbxo10和Fbxo27可能參與了光感受器細胞的凋亡和自噬過程。Ndrg1、Mt1及VEGFA可作為研究缺血缺氧相關眼底疾病的潛在功能性靶基因。
