基于RNA-seq技术的单眼形觉剥夺模型大鼠视皮层差异基因筛选、鉴定及功能分析

孙娅玲 严兴科 刘安国

1甘肃中医药大学第一临床医学院，兰州 730000；2甘肃中医药大学针灸推拿学院，兰州 730000

在视觉发育期内，由单眼斜视、屈光参差、高度屈光不正及形觉剥夺引起的最佳矫正视力低于相应年龄或双眼视力相差2行及以上的视觉障碍称为弱视。最新的Meta分析结果显示弱视的综合患病率为1.44%。弱视严重影响患者的学习和生活质量。弱视的发病机制复杂，与视网膜神经节细胞至视觉中枢的视传导通路及相应视皮层的结构和功能损伤有关。当前弱视的发病机制研究主要借助显微成像、电生理和脑功能成像等技术围绕视神经元形态功能以及视觉中枢改变来开展。但是，视功能减退的具体级联性分子机制研究尚不明确，尤其是转录、翻译水平的具体发病机制靶点还有待进一步探索和验证。转录组测序技术(RNA Sequencing，RNA-seq)从RNA水平上研究基因的表达情况，可以从整体水平研究特定时空内生物体某一细胞、组织甚至器官中全部基因的转录本，具有准确性高、通量高和成本低等优点。形觉剥夺是导致弱视发生的常见原因，在出生后视觉发育的关键时期进行单眼形觉剥夺，视皮质的形态结构和功能将发生显著的变化，导致形觉剥夺眼视力显著降低。本研究拟采用单眼眼睑缝合方法制作大鼠形觉剥夺弱视模型，并利用RNA-seq技术筛选形觉剥夺弱视发病相关基因，为弱视的靶向治疗提供新的参考依据。

1 材料与方法

1.1 材料

1.1.1

实验动物 14日龄SD幼鼠18只，体质量20～30 g[由甘肃中医药大学科研动物实验中心提供，使用许可证号：SYXK(甘)2011-0001]。实验动物饲养于甘肃中医药大学SPF动物实验室，环境温度(23±2)℃，相对湿度40%～60%，大鼠饲养期间正常进食和饮水。实验动物的使用遵循ARVO声明，本研究方案符合《甘肃中医药大学动物实验伦理审查标准》，并经甘肃中医药大学伦理委员会审核批准(批文号：2016-58)。

1.1.2

主要试剂及仪器 水合氯醛(天津市光复精细化工研究所)；复方托吡卡胺(邯郸康业制药有限公司)；红霉素眼膏(北京双吉制药有限公司)；TruSeq Stranded mRNA LT Sample Prep Kit试剂盒(美国Illumina公司)；Agencourt AMPure XP试剂盒(美国Beckman Coulter公司)；Qubit dsDNA HS Assay Kit试剂盒、Qubit RNA Assay Kit试剂盒、SuperScript Ⅱ Reverse Transcriptase试剂盒(德国Thermo Fisher公司)；Bioanalyzer 2100 RNA-6000 Nano Kit试剂盒(美国Aglient公司)；mirVanamiRNA ISOlation Kit试剂盒(美国Ambion公司)。RETI-PORT-21 Compact眼电生理诊断系统(德国Roland Consult公司)；Illumina Hiseq 2500型测序仪(美国Illumin公司)；大鼠脑立体定位仪(深圳市瑞沃德生命科技有限公司)；Tannon 2500凝胶成像系统(中国天能公司)；Nano Drop 2000紫外分光光度计(德国Thermo Fisher公司)；Agilent 2100生物分析仪(美国Aglient公司)；台式离心机(美国Beckman Coulter公司)；9700 PCR仪(美国ABI公司)。

1.2 方法

1.2.1

实验分组及单眼形觉剥夺大鼠模型建立 选取14日龄大鼠，按照体质量大小依次排序编号，按照随机数字表法随机等分为空白对照组与单眼形觉剥夺组，每组9只。参照文献[10]方法，选取单眼形觉剥夺组幼鼠，采用体积分数10%水合氯醛(0.3 ml/100 g)腹腔注射麻醉；将大鼠右眼周围皮肤严格消毒，剪去右眼眼睑上、下缘各0.5～1.0 mm；对齐创缘，用6/0缝线小心缝合上下眼睑3～5针，缝合手术后，在创面涂少量红霉素眼膏预防感染，并将模型大鼠置于自然光线、室温条件下饲养，连续饲养14 d。每日早、晚检查动物的眼睑缝合是否裂开、是否有感染等；有裂缝者给予及时修补，否则不纳入实验。

1.2.2

大鼠视皮质图形视觉诱发电位检测 2个组大鼠在同等生活环境下连续饲养14 d，将单眼形觉剥夺组大鼠右眼缝合的部位剪开，暴露其角膜，各大鼠均暗适应30 min后称量动物体质量，并按体质量给予体积分数10%水合氯醛(0.3 ml/100 g)腹腔注射，使大鼠处于浅麻醉状态，然后左右眼依次给予复方托吡卡胺滴眼液各1滴扩瞳；采用RETI-PORT-21 Compact眼电生理诊断系统记录大鼠右侧视皮质图形视觉诱发电位(pattern visual evoked potential，P-VEP)P波的潜伏期和振幅，电极内置电极膏，连接13 mm长毫针，参考电极于大鼠内眼角连线中点刺入，记录电极于大鼠枕骨粗隆上0.5 cm刺入。参数设置：屏幕刺激方式选用黑白棋盘方格，刺激频率为2.00 Hz，上限频率为50 Hz，下限频率为1 Hz，叠加次数为128次，放大倍数为200 000，采样时程为300 ms。检查时用不透明黑色眼罩遮盖其左侧眼,每只大鼠右眼P-VEP数据连续测3次，取其平均值。

1.2.3

RNA提取 收集各组大鼠双侧视皮质组织样本，每3个样本合并成1个样本进行分析，参照mirVanamiRNA ISOlation Kit试剂盒提取视皮质总RNA。加入600 μl Lysis/Binding Buffer，匀浆后加入30 μl miRNA Homogenate Additive，冰浴10 min；加入630 μl现配等体积酚(酸酚)氯仿溶液，高速离心(13 000 r/min，30 s)取上清，加乙醇，高速离心弃上清，加入350 μl miRNA Wash Solution，高速离心弃上清，将离心柱重新放置到收集管中；取沉淀混合10 μl DNaseⅠ和70 μl Buffer RDD QIAGEN上样至离心柱中的膜上；加350 μl miRNA Wash Solution 1，高速离心弃上清，将离心柱重新放置到收集管中；500 μl Wash Solution 2/3洗脱2次，高速离心弃上清，将离心柱重新放置到收集管中；空柱离心1 min，将离心柱放置到新的收集管中，柱中心加入100 μl 95 ℃洗脱液，放置2 min后离心，收集洗脱液，即为总RNA。使用NanoDrop微量分光光度计和Agilent 2100生物分析仪检测总RNA的完整性及质量，测定浓度与纯度后-80 ℃保存。

1.2.4

链特异性转录组建库及测序 使用TruSeq Stranded mRNA LT Sample Prep Kit试剂盒进行转录组建库，主要步骤如下：样品总RNA使用DNase消化DNA后，用带有Oligo(dT)的磁珠富集mRNA；加入打断试剂将mRNA打断成短片段，以打断后的mRNA为模板，用六碱基随机引物合成一链cDNA，然后配制二链合成反应体系合成二链cDNA，并使用试剂盒纯化双链cDNA；纯化的双链cDNA再进行末端修复、加A尾并连接测序接头，然后进行片段大小选择，最后进行PCR扩增；构建好的文库用Agilent 2100生物分析仪质检合格后，采用Illumina Hiseq 2500型测序仪进行RNA-seq测序。

1.2.5

差异表达基因的筛选 采用Hisat2软件将Clean reads与大鼠的参考基因组进行比对，用htseq-count软件获取每个样本中比对到基因上的Reads数。各样本中能够定位到基因组上的测序序列Reads数与参考基因组匹配率均维持在97%左右；各样本测序获得的Reads数中，能够唯一匹配到标准基因库中的Reads数目，维持在90%以上。采用Cufflinks软件来计算基因的表达量FPKM值，计算公式如下：FPKM (

A

)=。利用RNA-seq数据比较分析2个组样品中同一个基因是否存在差异表达，以样品中同一个基因表达水平的变化倍数(fold change，FC)>2且

P

<0.05为差异表达基因。用错误发现率(false discovery rate，FDR)<0.05进行差异表达基因过滤。

1.2.6

差异表达基因的功能分析 利用基因本体论(gene ontology，GO)富集分析对与弱视发病相关基因进行生物信息学分析，统计每个GO条目中所包括的差异基因个数，并用超几何分布检验方法计算每个GO条目中差异基因富集的显著性，筛选富集数大于2的GO条目。利用京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes,KEGG)数据库对相关差异基因进行代谢通路分析。

表1 2个组间大鼠右侧视皮质术前及术后P100潜伏期和波幅比较(x±s)Table 1 Comparison of preoperative P100 latency and amplitude of N45-P100 for rats between two groups (x±s)组别样本量潜伏期(ms)振幅(μV)术前术后术前术后空白对照组9106.64±5.52105.97±6.289.42±0.4310.48±0.96单眼形觉剥夺组9107.22±4.35121.34±3.64a10.21±1.555.17±1.44a 注:P100潜伏期:F分组=5.452,P=0.02;F时间=3.382,P=0.01.振幅:F分组=4.342,P=0.01;F时间=4.811,P=0.01.与各自术后空白对照组比较,aP<0.05(重复测量两因素方差分析,LSD-t检验) Note:P100 latency:Fgroup=5.452,P=0.02;Ftime=3.382,P=0.01.Amplitude of N45-P100:Fgroup=4.342,P=0.01;Ftime=4.811,P=0.01.Compared with respective blank control group,aP<0.05 (Two-way repeated measures ANOVA,LSD-t test)

表2 左侧视皮质差异表达基因的GO富集分析Table 2 Go enrichment analysis of DEGs in the left visual cortex基因功能分类功能描述P值LOC100910771、Mast4、Sgk2、Aurkb生物过程蛋白磷酸化0.007Foxr1、Sp7、Mafa生物过程RNA聚合酶Ⅱ启动子的转录调控0.019Grm2生物过程谷氨酸分泌0.031Sgk2、Aurkb分子功能蛋白丝氨酸/苏氨酸激酶活性0.005Tshz2、Sp7、Pdcd5、Mafa分子功能DNA结合0.034LOC100910771、Sgk2、Aurkb、Yme1l1、Katnal1分子功能ATP结合0.064 注:P值表示富集显著性 Note:P value indicated the significance of enrichment

1.3 统计学方法

采用SPSS 19.0统计学软件进行统计分析，计量资料经Shapiro-Wilk检验证实符合正态分布，以表示。本研究中各组手术前后P潜伏期和波幅总体比较采用重复测量两因素方差分析，组间多重比较采用LSD-

t

检验。

P

<0.05为差异有统计学意义。

2 结果

2.1 各组P100波形比较

与空白对照组相比，单眼形觉剥夺组大鼠P波潜伏期显著延长，振幅显著降低，差异有统计学意义(均

P

<0.05)，P波形发生了显著变化，表明造模成功(表1)。

2.2 差异表达基因筛选结果

单眼形觉剥夺组与空白对照组比较，左侧视皮质差异表达基因数为40个，其中表达上调基因19个，下调基因21个；右侧视皮质差异表达基因数为63个，其中上调基因39个，下调基因24个(图1)。

图1 2个组大鼠左侧、右侧视皮质差异表达基因火山图 A：2个组大鼠左侧视皮质差异表达基因火山图 B：2个组大鼠右侧视皮质差异表达基因火山图 FC:变化倍数Figure 1 Volcano plot of DEGs in bilateral visual corte of rats A:Volcano plot of DEGs in the left visual cortex of rats between two groups B:Volcano plot of DEGs in the right visual cortex of rats between two groups FC:fold change

2.3 差异表达基因的生物学功能分析

对左侧视皮质40个差异表达基因进行GO功能富集分析，富集数大于2的GO条目共6个。其中参与生物过程的3个GO条目为蛋白磷酸化、RNA聚合酶Ⅱ启动子的转录调控和蛋白质水解，参与分子功能的3个GO条目为蛋白丝氨酸/苏氨酸激酶活性、DNA结合和ATP结合(表2、图2)。

图2 左侧视皮质差异表达基因的生物学过程及分子功能分析Figure 2 Biological process and molecular function of DEGs in the left visual cortex

对筛选得到的右侧视皮质63个差异表达基因进行GO富集分析，富集数大于2的GO条目共8个。其中参与生物过程的2个GO条目为转录/DNA模板化、转录的负调控/DNA模板化；参与细胞组分的2个GO条目为细胞内、内质网的膜；参与分子功能的4个GO条目为核酸绑定、钙离子结合、磷脂酶A活性和锌离子结合(表3、图3)。

图3 右侧视皮质差异表达基因的生物学过程、细胞组分及分子功能分析Figure 3 Biological process,cellular component and molecular function of DEGs in the right visual cortex

2.4 发病基因的KEGG通路分析结果

KEGG基因组数据库的查询显示

Grm2

和

Pla2g2a

基因与视功能异常改变密切相关(表4)。

表3 右侧视皮质差异表达基因的GO富集分析Table 3 Go enrichment analysis of DEGs in the right visual cortex 基因功能分类功能描述P值LOC100910771、Brd2、LOC102552640、Atf3生物过程转录,DNA模板化0.003Lhx9、LOC102552640、Atf3生物过程转录的负调控,DNA模板化0.010Siglec8、LOC100910771、LOC103690038、Zbbx细胞组分细胞内0.005Fmo2、Oas1f、Fkbp6细胞组分内质网的膜0.007LOC103690038、Nol8分子功能核酸绑定0.030LOC100911363、LOC100911486、LOC100910178分子功能钙离子结合0.030Pla2g2a分子功能磷脂酶A2活性0.047Lhx9、Apobec1、Zbbx分子功能锌离子结合0.033 注:P值表示富集显著性 Note:P value indicated the significance of enrichment

表4 Grm2和Pla2g2a基因相关代谢通路Table 4 Metabolic pathway of Grm2 and Pla2g2agenes 基因物种代谢通路Grm2大鼠钙离子信号通路、FoxO信号通路、磷脂酶D信号通路、脑组织神经活动配体-受体相互作用、长时程增强、谷氨酸能突触、长时程抑制Pla2g2a大鼠甘油磷脂代谢、醚脂代谢、花生四烯酸代谢、亚油酸代谢、亚麻酸代谢、新陈代谢通路、维生素消化吸收

3 讨论

本研究结果显示，视觉发育关键期内右眼形觉剥夺模型大鼠两侧视皮质差异表达基因的GO生物学功能虽然不尽相同，但绝大部分代谢通路都与视功能的异常改变密切相关，最终筛选了

Grm2

和

Pla2g2a

这2个弱视相关发病基因，并对其进行KEGG通路分析，推测其调控的生物学过程改变可能是导致视功能下降的主要病理机制之一。

Grm2

基因主要参与调控谷氨酸突触、长时程增强(long-term potentiation，LTP)、长时程抑制(long-term depression，LTD)等视觉信号通路过程。谷氨酸是视网膜到外侧膝状体再到视皮层透射轴突末梢主要的兴奋性神经递质，与视觉通路的兴奋性传导和视功能的正常发挥密切相关。代谢型谷氨酸受体(metabotropic glutamate receptors，mGluRs)是谷氨酸的一类神经递质类受体，可调控神经兴奋性信号传导，改变突触膜及神经可塑性，诱导产生LTP和LTD。有研究证实，给予mGluR的激动剂3,5-二羟基苯甘氨酸可诱导海马区抑制性中间神经元上的α-氨基-3-羟基-5-甲基-4-异恶唑丙酸受体激活，从而介导LTP的产生；而给予mGluR的抑制剂可阻断LTP效应的产生。另外，mGluR被激活后，可刺激细胞产生活化磷脂酶C和三磷酸肌醇通路从而释放更多的Ca，再通过激活Ca依赖的蛋白激酶而增加Ca的内流，从而产生LTD效应。可见，mGluR被激活将引起相关视觉信号传递及其介导的细胞内级联反应，从而产生LTP和LTD现象。而本研究证实形觉剥夺模型大鼠视皮层

Grm2

基因表达显著下调，

Grm2

基因介导的mGluRs信号通路及下游相关生物学活性物质的表达受限，可能是产生视觉剥夺效应的一条重要的分子生物学途径。

Pla2g2a

基因主要参与α-亚麻酸代谢通路和花生四烯酸代谢通路。研究证实，上述2个信号通路与正常视功能的发挥密切相关。α-亚麻酸对大脑神经元发育和视觉系统发育都具有较高的营养性作用，α-亚麻酸能够促进大鼠的脑发育和视功能发育，提高子代大鼠的学习记忆能力和视网膜电位水平，α-亚麻酸缺乏将对中枢神经系统生理功能，尤其是对视觉功能造成不利的影响。膳食摄入α-亚麻酸能够在基因调控下经脱氢与碳链延长，生成一系列代谢产物，转变为长链多不饱和脂肪酸(polyunsaturated fatty acids，PUFAs)。而二十二碳六烯酸(docosahexaenoic acid，DHA)作为PUFAs的一种，是视网膜视杆细胞和视锥细胞外部的主要不饱和脂肪酸，能够影响视网膜发育中的神经节苷脂代谢，提高二唾液酸神经节苷脂3的水平，从而促进视网膜的发育成熟，同时，保护视觉光感受器免受氧化应激损伤。另外，视网膜感光细胞膜和神经细胞膜中含有高浓度DHA，表明了DHA对于细胞膜相关的生理功能具有重要作用，其参与调节了视神经信号传导、神经传递效率和神经元突触的可塑性修复。研究发现，在孕期的最后3个月胎儿视网膜磷脂酰乙醇胺中的DHA含量显著增加，证实了胎儿在后期视网膜感光器的快速发育需要依赖大量的DHA。本研究结果同样证实，形觉剥夺模型大鼠视皮层存在

Pla2g2a

基因表达下调，可能导致α-亚麻酸及其下游的DHA转化及吸收异常，进一步导致视觉光转导异常，这可能是弱视形成视觉剥夺效应的另一条重要的分子生物学途径。

以往在弱视发病机制的研究中发现谷氨酸及其离子型受体N-甲基-D-天冬氨酸在视皮层结构和功能可塑性中发挥重要的作用。本研究又从基因表达水平揭示了谷氨酸及其不同亚型的mGluRs也参与了神经元突触结构和功能可塑性的改变过程，为今后弱视发病机制研究和视觉系统的可塑性研究提供了新的参考依据。在今后的研究中将结合现有基因在cDNA的表达水平，同时对差异基因调控的上游、下游的代谢产物进行分析，进一步验证本研究结论。

利益冲突

所有作者均声明不存在利益冲突

作者贡献声明

孙娅玲：直接参与起草文章、文章撰写、分析/解释数据；刘安国：直接参与酝酿和设计实验、实施研究、采集数据、对文章的知识性内容作批评性审阅；严兴科：统计分析、获取研究经费、指导实验