基于RNA-seq技术筛选低钙培养奶牛乳腺上皮细胞差异表达基因

王 娟，姜忠玲，丛 霞，曹荣峰，田文儒，丰艳妮，李华涛

(青岛农业大学动物医学院，山东 青岛 266109)

钙是奶牛体内不可或缺的常量无机元素，参与体内多种生理过程[1]，如凝血、骨生成[2]、肌肉收缩、神经传导以及能量代谢[3]等，并发挥着十分重要的作用。当奶牛受到各种因素的影响，致使其体内的钙浓度变化超出自身的调节范围时，就会引起钙代谢紊乱。尤其是围产期奶牛因妊娠和泌乳导致钙大量流失[4]，极易引发低钙血症。目前，临床上对奶牛围产期低钙血症的检测多依赖于血清生化指标，但由于个体差异的原因，不同奶牛对血钙浓度变化范围耐受程度不同，血清生化指标并不能准确的预测低钙血症的发生与发展，也无法准确判定亚临床型低钙血症。因此，找出更为灵敏且精准判定奶牛围产期低钙血症的生物标志物将有重要的应用价值。随着测序技术的发展，转录组测序(RNA sequencing,RNA-seq)技术作为一种新兴技术被广泛应用于各个领域，其测序结果也在后续研究中相继被证实。虽然目前RNA-seq技术仍存在一些不足，但其具有筛选范围广、精度高、成本低、无损检测等特点，仍是临床发掘疾病生物标志物不可取代的手段。

本研究拟筛选低钙培养奶牛乳腺上皮细胞(Cow mammary epithelial cells, CMEC)差异表达基因，为临床寻找更为准确、灵敏的奶牛低钙血症生物标志物提供方法与思路。

1 材料与方法

1.1 细胞来源及处置 本实验室冻存的自行培养并且已经鉴定的CMEC细胞系，将细胞分为低钙培养组(DE)及空白对照组(DP)，提取总RNA后用于后续试验。转录组测序每组设3个重复，定量反转录聚合酶连锁反应(qRT-PCR)验证每组，设9个重复。

1.2 试剂 TRIzol、反转录试剂盒、DL2 000 DNA marker、TB GreenTMpremix ExTMTaqⅡ荧光染料，均购自TaKaRa公司；焦碳酸二乙酯(Diethyl pyrocarbonate，DEPC)、乙二醇双(2-氨基乙基醚)四乙酸[Ethylenebis(oxyethylenenitrilo)tetraacetic acid, EGTA]、青链霉素混合液(双抗，100×)、50×TAE缓冲液,均购自Solarbio公司；D-MEM/F-12培养基，购自Gibco公司(使用时添加10%胎牛血清和1%双抗)；胎牛血清，购自ExCell Bio公司；钙试剂盒(微板法)及CCK-8细胞活力检测试剂盒，均购自南京建成生物研究所；PBS缓冲液用DEPC处理过的水配制；EASY spin组织/细胞RNA快速提取试剂盒，购自艾德莱生物有限公司；Kadaq 2×PCR Master Mix，购自Abm公司；引物由青岛派诺森基因科技有限公司合成；其余未特别指明的生化试剂，均购自Sigma公司。

1.3 方法

1.3.1 依他酸(EGTA)浓度筛选及CMEC低钙培养 将经复苏后培养的CMEC随机分为空白对照(DP)、0.5 mmol/L EGTA(DE1)、1.0 mmol/L EGTA(DE2)和2.0 mmol/L EGTA(DE3)组，分别培养24、48 h和72 h后镜检观察细胞生长状态，并用Toup View软件拍照记录。按照钙试剂盒说明书对4个分组的培养基进行Ca2+浓度测定。按照CCK-8检测试剂盒说明书分别于培养后24～144 h检测细胞活力(每隔24 h检测1次)，以造成细胞生长明显抑制且不造成细胞大量凋亡为标准选择低钙培养浓度。

1.3.2 CMEC的转录组测序及数据分析 选取生长状态良好，细胞量>5×106个的DP组和DE2组CMEC，用预冷的PBS洗涤2次，加入1 mL TRIzol试剂反复吹打，至溶液澄清度均一后，移入冻存管中，分别标记为DP和DE，密封后于-80 ℃保存，将制备好的样品保存于干冰中送公司进行转录组测序。将测序获得的原始数据经碱基识别分析转化为原始测序序列(Raw reads)，过滤后得到Clean reads，后续的数据分析都基于Clean reads。将Clean reads与参考基因组的序列进行对比，得到基因的表达量，经分析后获得差异表达基因。最后对基因功能进行注释及GO、京都基因与基因组百科全书(Kyoto encyclopedia of genes and genomes,KEGG)分析。

1.3.3 奶牛乳腺上皮细胞的转录组测序的qRT-PCR验证

1.3.3.1 测序样品制备 按照RNA提取试剂盒说明书从低钙组及对照组细胞中提取总RNA。用Eppendrof 蛋白质核酸定量仪检测提取的总RNA浓度，A260/A280和A260/A230值均在2.0以上时可用于后续试验。将提取的总RNA于-80 ℃保存备用。

1.3.3.2 引物合成 使用牛β-actin为参考基因，从测序结果中筛选8个差异表达显著性基因，参照美国国家生物技术信息中心(National Center for Biotechnology Information,NCBI)公布的基因序列，采用NCBI官网的Primer-BLAST进行引物设计，序列见表1。

表1 qRT-PCR验证基因的引物序列Table 1 Primers used for qRT-PCR in the study

1.3.3.3 反转录及核酸电泳验证 按照TaKaRa反转录试剂盒说明书配制10 μL gDNA去除体系，根据RNA提取浓度，保证体系中含有1 μg RNA，震荡混匀后离心5～10 s，于42 ℃孵育2 min。配制反转录扩增体系，加入gDNA去除后的RNA，漩涡混匀、离心后，于PCR仪中37 ℃孵育15 min、85 ℃孵育5 s 后得到cDNA，于-20 ℃保存备用。使用β-actin引物，配制10 μL PCR扩增体系，对cDNA进行PCR扩增。PCR扩增程序：94 ℃预变性3 min；94 ℃ 变性30 s，60 ℃退火30 s，72 ℃延伸30 s，共29个循环；72 ℃再延伸7 min，最后4 ℃孵育。核酸电泳对扩增的cDNA进行验证，每组设9个重复。电泳结束后，在凝胶成像系统下进行观察并保存图像。

1.3.3.4 qRT-PCR验证 配制20 μL qRT-PCR扩增体系，按照Light Cycler®96 SW荧光定量PCR仪操作说明书，对所选8个差异基因进行qRT-PCR验证。以β-actin为参考基因，采用2-ΔΔCt方法对所测基因mRNA的相对表达情况进行分析。

1.4 统计学分析 使用GraphPad Prism 5软件对数据进行统计学分析，组间采用单因素方差分析，计算Log2FC值后与RNA-seq测序结果进行差异倍数比较。

2 结果

2.1 EGTA浓度筛选及CMEC低钙培养 采用不同Ca2+浓度培养基培养奶牛乳腺上皮细胞24 h后，DE1、DE2组与DP组间无明显差异，而DE3组明显可见大量细胞皱缩成圆形亮点，细胞质空泡化，核质凝缩。培养48 h及72 h后，DE1组较DP组无明显肉眼可见变化，DE2组细胞较DP组细胞生长缓慢，并有少量细胞凋亡，DE3组细胞开始大量凋亡脱落(中插彩版图1)。

图1 不同Ca2+浓度培养基培养后24、48 h和72 h的奶牛乳腺上皮细胞 (40×)Fig.1 Bovine mammary epithelial cells cultured within different concentrations of Ca2+ for 24, 48 h and 72 h

采用钙试剂盒分别对DP、DE1、DE2组和DE3组培养基进行Ca2+浓度检测发现，4个组钙离子浓度分别为1.02、0.87、0.43 mmol/L和0.03 mmol/L。CCK-8结果显示，与DP组相比，DE1、DE2、DE3组对CMEC的生长均有抑制作用，DE1组与DE2组细胞在48 h存活率明显上升，72 h存活率下降，之后逐渐趋于平稳(图2)。DE3组细胞的生长受到显著抑制，且随着培养时间的延长，细胞逐渐凋亡。据此，将选取DE2组细胞培养72 h后进行后续试验。

图2 不同Ca2+浓度培养基培养后奶牛乳腺上皮细胞存活率Fig.2 Effects of different concentrations of Ca2+ on survival rate of bovine mammary epithelial cells

2.2 CMEC的转录组测序及数据分析

2.2.1 测序数据质量分析 根据DP组和DE组的碱基质量分数和碱基质量分布，不论是正向测序还是反向测序，其碱基质量分数90%以上集中于30%～40%，说明其错误率控制在0.1%～0.01%。DP组和DE组的GC含量大多集中在40%～50%，说明在建库过程中序列是随机打断的且测序文库的质量良好。通过Hisat 2软件将所测序列与参考基因组Bos taurus进行比对，其比对率均达96%以上(表2)，说明所选参考基因组合适，且样品可满足后续分析的需求。

表2 参考基因组比对统计分析Table 2 Statistical analysis of comparison with reference genome

2.2.2 差异表达基因数据分析 以差异倍数∣log2(Fold change,FC)∣≥1 及显著水平P<0.05为差异显著性标准筛选DE组与DP组之间的差异表达基因，共计得到显著差异表达基因116个，其中52个表达上调基因，64个表达下调基因(见中插彩版图3A)。对2组的差异表达基因进行GO富集分析，发现其主要集中于三大类30个小类中，主要参与了染色质构型、细胞内组成、染色质结合等(中插彩版图3B)。差异表达基因大多富集于硫代谢、硫中继系统、病毒致癌作用、半胱氨酸和蛋氨酸代谢等信号通路(中插彩版图3C)，根据富集结果，以与钙平衡调节相关的甲状腺素信号通路和甲状腺素合成通路为例，查找到差异表达上调基因KAT2A和差异表达下调基因DUOXA2。甲状腺素相关通路对机体钙水平调节具有重要意义，该差异基因可能与细胞受低钙环境培养有关。

图3 两组间差异表达基因分析Fig.3 Analysis of differentially expressed genes between the two groupsA:火山图，图中一个点代表一个基因，红点为上调基因，蓝点为下调基因；B:GO富集分析图; C:KEGG通路富集分析图，图中点越大代表所富集的基因数量越多，颜色越深代表P值越小A: Volcano plot, a dot in the diagram represents a gene, red dot represents up-regulated gene, blue dot represents down-regulated gene;B: GO enrichment analysis diagram; C: KEGG pathway enrichment analysis diagram. The larger the dot in the picture, the more genes are enriched. The darker color, the lower P value.

2.3 CMEC转录组测序的qRT-PCR验证

2.3.1 样品的qRT-PCR检测 经蛋白质核酸定量检测仪检测，提取的DP组和DE组总RNA浓度均在440 ng/μL以上，且A260/A280、A260/A230值均>2.0，提取质量优良。以参考基因β-actin为引物进行cDNA的PCR扩增反应，核酸电泳对扩增产物进行验证。结果可得，产物均在250 bp左右，可初步判断产物为β-actin片段。2组及其重复均可见较亮的条带(图4)，可用于后续试验。

图4 cDNA扩增结果Fig.4 Amplification results of cDNAM：DL-2 000 DNA相对分子质量标准； DP-W：空白对照； DE-W：低钙培养组M: DL-2 000 DNA marker; DP-W: Blank control; DE-W: Low calcium medium group

2.3.2 筛选基因的qRT-PCR验证 将qRT-PCR结果经计算后与RNA-seq测序结果的差异倍数Log2FC进行比较，7个基因与RNA-seq测序结果相符，有1个 基因与测序结果不符(图5)。说明RNA-seq测序结果较为可靠，准确性较高。

图5 RNA-seq与qRT-PCR差异倍数比较Fig.5 Comparison of differential multiples between RNA-seq and qRT-PCR

3 讨论

围产期低钙血症是奶牛常见的钙代谢紊乱性疾病，处于低钙血症状态的奶牛常因免疫力低下而诱发子宫复旧不全、子宫内膜炎、胎衣不下、乳房炎[5-8]和真胃变位[8]，同时还会影响奶牛的生产性能[9-12]，对于分娩后的奶牛来说，大量的钙随乳汁流逝更易引发低钙血症。泌乳期CMEC在调节钙的转运以及维持胞内钙平衡等方面起重要作用[13]，但其具体机制尚不明确。目前，已有研究人员利用RNA-seq技术对动物[14]、植物[15-16]等在遭受外界刺激后的转录组进行测序，从而发现了相关差异表达基因，为动物、植物遭受外界应激后相关分子机制的研究奠定了基础。人医中也有利用该技术用于疾病诊断[17]以及寻找疾病生物标志物的相关报道[18-20]。对于低钙血症的检测临床多依赖于血清生化指标，但取血时操作不当易引发应激、细菌感染等问题，十分不便。故本研究选用CMEC为研究对象，利用RNA-seq技术研究低钙培养对CMEC基因表达水平的影响，筛选差异表达基因，以研究该细胞应对低钙环境的分子机制，以期寻找到可经乳汁中乳腺上皮细胞检测围产期奶牛低钙血症的生物标志物。

经差异表达基因GO富集分析及KEGG通路分析显示，低钙培养可影响细胞的基因表达，且2组间差异表达基因富集于多个通路之中，以与钙调节相关的甲状腺通路为例查找到富集基因KAT2A和DUOXA2。甲状腺激素主要作用于骨骼，能够调节骨骼的新陈代谢，影响骨骼的重建过程，促进破骨细胞的活性发生改变，使骨吸收过程大于骨形成过程[21]，可通过增加尿液中的钙含量和粪钙的排出量来调节机体的钙平衡。KAT2A和DUOXA2基因分别指导赖氨酸乙酰基转移酶和双氧化酶成熟因子的合成。赖氨酸乙酰基转移酶是组蛋白乙酰转移酶的一种，其可作为转录激活因子将乙酰基从乙酰辅酶A中转移到特定组蛋白氨基末端的赖氨酸上，使与该段结合的DNA双螺旋结构变疏松，利于转录因子的结合。有研究表明，KAT2A基因可通过经典Wnt通路调节牙周膜干细胞的成骨分化作用，且与之呈正相关，KAT2A表达水平降低可使细胞的成骨分化能力减弱，使牙齿骨量减少[22-23]。甲状腺素是由TPO在H2O2的参与下经催化碘的有机化、酪氨酸碘化以及碘化酪氨酸耦联等过程生成的，双氧化酶2与H2O2的生成有关，双氧化酶成熟因子2则主要起协助双氧化酶2蛋白转运、成熟以及定位的作用[24]。本研究中KAT2A和DUOXA2基因表达分别上、下调的结果说明，低钙环境下培养可能会影响细胞甲状腺素的合成，致使甲状腺素含量异常，影响骨的代谢过程[25]，进而影响钙水平的调节。因此，KAT2A和DUOXA2基因差异表达可能间接与该细胞应对低钙培养环境有关，这为寻找围产期奶牛低钙血症生物标志物提供了研究方向，但其具体参与机制和临床应用还需进一步的研究。