基于转录组测序的铝胁迫下甘蓝型油菜新内参基因的发掘与引物开发

丁 戈，黄 杨，陈伦林，李书宇，宋来强，熊 洁

(江西省农业科学院 作物研究所，国家油料改良中心南昌分中心，江西省油料作物生物学重点实验室，农业农村部长江中下游作物生理生态与耕作重点实验室，江西 南昌 330200)

酸性土壤(pH 值<5.5)的溶解铝(Al3+)通过抑制根系生长，影响水分和养分吸收，限制了作物的生长与产量[1-3]。植物根系是铝毒害的主要作用部位，根系通过诱导或增强相关铝抗性基因表达，如转录、RNA剪接、蛋白质翻译或翻译后修饰等调控，合成功能性基因产物来应对铝毒害[4-9]。实时荧光相对定量PCR(Real-time fluorescence relative quantitative PCR)是研究基因转录水平表达量的常用技术，具有高通量、特异性强、灵敏度高、重复性高等特点。其对目标基因的相对表达量分析需要选择高度保守且稳定表达的基因作为内参基因。常用内参基因通常是在所有细胞中表达的持家基因(Housekeeping gene)，但持家基因的表达稳定性是相对的，已经不能满足多样化的试验体系中相对定量分析对内参基因的要求[10]。内参基因需要根据物种、组织器官、生长发育阶段和试验条件等因素进行筛选、评估和验证[11]。根据具体的样品类型选择合适的内参基因作为参照物校正样品间的试验误差，通过对结果数据的标准化处理，确保实时荧光相对定量PCR分析结果的准确性[12]。随着转录组测序技术(RNA-seq)[13-15]的发展，基于高通量大数据的新内参基因发掘可以不再局限于常见的持家基因。本研究采用有参考基因组的转录组测序技术筛选候选内参基因，随后采用实时荧光相对定量PCR分析、基因表达稳定性分析及相对表达量分析等技术，评估、验证候选内参基因的表达稳定性，进行铝胁迫下甘蓝型油菜苗期根组织新内参基因的发掘以及实时荧光相对定量PCR引物的开发。

1 材料和方法

1.1 试验材料

植物材料培养于人工气候温室(昼夜时间14 h/10 h，昼夜温度25 ℃/20 ℃)。供试材料为前期研究筛选获得1对苗期耐铝性差异显著的油菜品种：耐铝品种赣油杂7号(Ganyouza 7)和铝敏感品种蓉油18(Rongyou 18)。将2个品种分别播种在装有纯净水的培养箱中，7 d后每品种选择长势一致的40株幼苗，移栽于装有改良型Hoagland(pH值 5.8～6.0)营养液[16]的培养箱中，依次1/4、1/2、完全营养液各培养7 d。在苗龄28 d时，分为6个试验组，设置3次重复，对照组为未施加AlCl3时的0 h取样，处理组为施加100 μmol/L AlCl3(pH值 4.5)后分别在3，24 h取样。每个试验组各取样3个生物重复样品，共18个样品。每个样品根部取样大于4 g，液氮速冻后保存于-70 ℃冰柜。后续的文库制备及有参考基因组的转录组测序，总RNA提取，cDNA逆转录反应及实时荧光相对定量PCR由广州基迪奥生物科技有限公司完成。

1.2 试验方法

1.2.1 文库制备及转录组测序分析 使用TRIzolTMReagent(Invitrogen，美国)提取样品的总RNA，质检合格后，真核生物转录组文库构建使用NEB#7530试剂盒(New England Biolabs，美国)，文库质检(DNA 1000 assay Kit，Agilent，美国)合格后，每个样品构建一个300～500 bp片段长度的双端(Paired-end，PE)测序文库，在Illumina HiSeq2500测序仪(Illumina，美国)上进行PE150测序。经测序原始数据质控，基于基因组注释的定量质控、基因分析、表达量统计等步骤，获得基因表达的FPKM值与组间差异倍数|log2FC|值。本研究中使用的甘蓝型油菜基因序列均参考甘蓝型油菜Darmor-bzh参考基因组[17-18](Darmor-bzhv4.1，下载自http：//www.genoscope.cns.fr/brassicanapus)。

1.2.2 总RNA提取与质量检测 用于cDNA逆转录反应及实时荧光相对定量PCR分析的总RNA提取参照赛默飞世尔科技(中国)有限公司InvitrogenTMTRIzolTMReagent产品说明书。提取完成后存于-70 ℃备用。微量分光光度计(型号NanoDrop 2000，Thermo Scientific，美国)检测RNA的质量与浓度。RNA完整性检测采用1%(0.5 μg/mL EB)非变性琼脂糖凝胶(1×TAE缓冲液，pH值 8.5，40 mmol/L Tris-乙酸，2 mmol/L EDTA)，RNA原液上样1 μL。电泳仪参数设置为恒压5 V/cm，电泳时间20～30 min或溴酚蓝染料迁移至胶下游3/4处时停止电泳。EB染色的凝胶置于凝胶成像系统(型号Gel DocTMXR+，BIO-RAD，美国)观察并照相。

1.2.3 cDNA逆转录反应体系 用于实时荧光相对定量PCR分析的cDNA逆转录反应体系的配制参照南京诺唯赞生物科技股份有限公司HiScript®Ⅱ Q RT SuperMix for qPCR(+gDNA wiper)产品说明书。使用1 μg总RNA按两步骤依次建立16 μL基因组DNA去除反应体系及20 μL cDNA逆转录反应体系。反应完成后存于-70 ℃备用。

1.2.4 实时荧光相对定量PCR引物设计与合成、反应体系及程序设置 遵循实时荧光定量PCR的引物设计原则，采用Primer Premier 5设计，上海捷瑞生物工程有限公司合成及PAGE方式纯化引物。Actin7上游引物序列为5′-GTTACGCTCTTCCT CACGCT-3′，下游引物序列为5′-CAGGGGCTTCCAT TCCG-3′，扩增子碱基对数为315 bp；ARFA1E上游引物序列为5′-TCGTCACAACCATTCCCACC-3′，下游引物序列为5′-AGATAAGCCCCTGCGTGTTC-3′，扩增子碱基对数为144 bp；SAP5上游引物序列为5′-CGCAGAGAACGGAGAAGGAA-3′，下游引物序列为5′-ACGACTTTGGCTGGTGATGG-3′，扩增子碱基对数为238 bp；UXS3上游引物序列为5′-GACCC TCTCATCCACCCTCA-3′，下游引物序列为5′-CACC ATCTTTATCTCTATGCTCGG-3′，扩增子碱基对数为432 bp；UBC9上游引物序列为5′-GCCCTTACTCTG GTGGTGTT-3′，下游引物序列为5′-CGAACAGAT AGATAGCAGCACC-3′，扩增子碱基对数为197 bp。反应体系的配制参照南京诺唯赞生物科技股份有限公司ChamQTMSYBRTMqPCR Master Mix(High ROX Premixed)产品说明书，使用4 μL cDNA稀释液(原液稀释10×)为模板配制20 μL反应体系。SYBRTMGreen Ⅰ嵌合荧光实时定量PCR反应在StepOne PlusTMReal-Time PCR System(Applied Biosystems，美国)上运行，步骤依次为，步骤1预变性：95 ℃预变性90 s；步骤2循环反应(三步法)：进行40个循环，95 ℃变性5 s，60 ℃退火15 s，72 ℃延伸20 s；步骤3熔解曲线采集：升温(速率100%)至95 ℃变性15 s，降温(速率100%)至65 ℃退火60 s，缓慢升温(0.3 ℃/s)至95 ℃变性15 s。分别在循环反应延伸阶段、熔解曲线缓慢升温及后续变性阶段采集实时荧光信号。

1.2.5 基因表达稳定性分析的算法 geNorm[19](下载地址https：//genorm.cmgg.be)采用成对比较法(Pairwise comparison approach)，计算出基因表达稳定性M值(Average pairwise variationV，平均配对变异V值)。当M值低于1.5时，M值越低，稳定性越高。当配对变异值Vn/Vn+1<0.15时，内参基因组合的最适基因是高稳定性排名的前n个。geNorm算法的数据输入格式为，每个样品(含生物学重复)相对定量EΔCt值，其中样本i中基因j的Ct(i，j)值为其技术重复Ct值的算数平均数(Arithmetic mean)，ΔCt(i，j)=minCt(m，j)-Ct(i，j)。minCt(m，j)为该基因j在所有样本m中最小Ct值,E为扩增效率(效率为100%时，E=2)。NormFinder[20](下载地址https：//www.moma.dk/normfinder-software)采用基于模型的方法(The Model-Based Approach)，计算出基因组内方差(Intra-group variances)与组间方差(Inter-group variances)，再计算出基因表达稳定性ρ值，组间方差越接近0，同时平均组内方差越小，ρ值越低，稳定性越高。NormFinder算法的数据输入格式为，每个样品(含生物学重复)的技术重复Ct值的中位数(Median)，并勾选"log transform data"选项进行自然对数log(x)转换，或数据输入相对定量EΔCt值。BestKeeper[21](下载地址https：//www.gene-quantification.de/bestkeeper.html)采用重复配对相关分析(Repeated pair-wise correlation analysis)、相关分析(Correlation analysis)和回归分析(Regression analysis)，计算Ct值的标准偏差(Standard deviation，s)、变异系数(Coefficient of variation，CV)和配对相关系数(Coefficient of correlation，r)。当s<1时，r值越大，CV值和s值越低，稳定性越高。BestKeeper算法的数据输入格式为，每个样品(含生物学重复)的技术重复Ct值的算数平均数以及基因的扩增效率。

2 结果与分析

2.1 转录组测序的表达量分析

有参考基因组的转录组测序的基因表达量计算使用FPKM法(Fragments per kilobase of transcript per million mapped fragments，每百万片段中来自于某基因每千碱基长度的片段)[22]，FPKM法适用于双端测序文库，其公式能够抵消基因编码区碱基数和样品测序片段数差异对基因表达量运算的影响，计算得到的FPKM值可直接用于比较不同组间的基因表达差异，其比值即为2组间基因表达差异倍数(Fold Change，FC)，对FC取以2为底的对数之后即为log2FC，在|log2FC|>1时，表达差异2倍以上。持家基因具有2个特征，即在试验组织中稳定表达，且表达水平高于背景值，因此，以FPKM值(大于10，表达水平高于背景值)和log2FC绝对值(小于0.1，表达变化较小)作为双层过滤筛选条件，从6个试验组的转录组测序数据中筛选出符合持家基因2个特征的5个候选内参基因(Actin7、ARFA1E、SAP5、UXS3、UBC9)。此外，由于它们属于不同的功能分类，在理论上降低了这些基因会受到共调控的可能性，因此，它们可作为不同类别的候选内参基因(图1-A、表1)。传统持家基因β-肌动蛋白基因(β-actin)的异型体ACT2[23]、磷酸化的甘油醛-3-磷酸脱氢酶基因GAPDH的细胞质同工型GAPC1(Glyceraldehyde-3-phosphate dehydrogenase C subunit 1)[24]、α-微管蛋白基因(α-tubulin)的亚型TUA4(Tubulin alpha-4 chain)[25]在2个品种的铝胁迫时间梯度(3，24 h)处理组中表达量发生变化，不适宜作为铝胁迫分析的内参基因(图1-B、表1)。基于转录组测序数据的候选内参基因的表达稳定性尚不确定，需要后续结合实时荧光相对定量PCR分析、基因表达稳定性分析及相对表达量分析等方法进行评价。

表1 油菜候选内参基因与拟南芥同源基因的蛋白质序列比对分析Tab.1 Alignment analysis of the protein sequences of rape candidate reference genes and Arabidopsis thaliana homologous genes

2.2 总RNA质量检测

提取的总RNA进行琼脂糖凝胶电泳和紫外分光光度检测，6个试验组共18个样品的大片段28S rRNA和18S rRNA条带明亮清晰，泳道无弥散区，加样孔无基因组及蛋白质杂质残留(图2)，样品浓度在183.9～462.7 ng/μL，OD260/280均值在1.92～2.06，OD260/230均值在1.89～1.96，提取的样品RNA完整性好，质量较好，可用于后续分析。

2.3 扩增曲线及熔解曲线分析

SYBR®Green I嵌合荧光实时定量PCR分析表明，各候选内参基因扩增子(Amplicon)的扩增曲线清晰可辨认(图3)，其所有PCR复孔(包括每个样品的3个技术性重复，n=54)扩增曲线重复性高，未出现平台期明显降低的引物二聚体扩增曲线。在三步法循环反应完成后，进行熔解曲线分析。熔解曲线分析是一种简单、直接地检查反应中引物二聚体和扩增子的方法。其原理是核酸的熔解温度(Tm)受碱基对数目、GC碱基含量以及碱基错配等因素影响，因此，可根据熔解特性区分扩增产物。熔解曲线分析能够确保反应的特异性，可以替代扩增产物的凝胶电泳检测步骤。熔解曲线(图3、表2)分析表明，各候选内参基因所有PCR复孔(n=54)的扩增产物熔解曲线均为显著的尖锐单一熔解峰，峰的Tm值大于80 ℃，且起峰到落峰温度跨度不高于7 ℃，表明引物特异性高，扩增产物中未出现非特异性扩增子及引物二聚体，实时荧光相对定量PCR的结果可用于后续分析。

表2 扩增子的熔解温度与扩增效率Tab.2 Tm and amplification efficiency of amplicons

2.4 扩增效率分析

在完成高质量的实时荧光PCR后，需要运用合适的程序[26-28]对原始荧光值(Raw fluorescence)或阈值循环(Threshold cycle，Ct)数据进行正确分析。程序LinRegPCR[29-30](Version 2020.0.0.3，下载地址https：//www.genetargetsolutions.com.au/qpcr/linregpcr-and-factor-qpcr)，将每个PCR反应孔(包括技术性重复的复孔)在第1～40个循环扩增的荧光信号实时检测绝对值(Non-baseline corrected data)经对数log10(x)转换，通过线性回归分析计算线性区段(扩增曲线指数期)的斜率，能够准确获得每个反应孔的具体扩增效率(E)。LinRegPCR分析表明，各候选内参基因的所有PCR复孔(n=54)的扩增线性区段的相关系数r值为0.997 42～0.999 67，表明该区段的循环数与对数转换后的荧光信号值之间的线性相关程度高。各候选内参基因扩增子的扩增效率为87.3%～92.1%(表2)。

2.5 基因表达稳定性分析

采用常用的比较基因表达变化及稳定性的geNorm、NormFinder和BestKeeper等3种Microsoft Excel 2003 宏或加载项对候选内参基因进行表达稳定性分析。因算法差异，各候选内参基因在3个软件中的排序略有差异，对各候选内参基因的3种排序值进行几何平均值(Geometric mean)运算后进行综合排名，优先选择综合排名靠前的候选内参基因，其中表达稳定性综合排名前两位的候选内参基因分别是SAP5和UXS3(表3)。

2.6 基因相对表达量分析

表3 候选内参基因表达稳定性分析Tab.3 Expression stability analysis of candidate reference genes

表4 样品组候选内参基因相对表达量Tab.4 Relative expression of candidate reference genes in sample group

3 结论与讨论

实时荧光PCR的相对定量数据分析方法主要有2-ΔΔCt法[34]、双标准曲线法[35]、线性回归分析法等。相对定量2-ΔΔCt法不需要计算扩增效率，但需要对试验条件进行严格优化，使目标基因和内参基因的扩增效率接近100%且偏差在5%以内。由于实际扩增效率会不同于理想扩增效率(100%)，为了保证结果的准确性，需要检测引物对的扩增效率。线性回归分析法采用的数据全部来自指数期，计算得到的理论扩增效率会低于真实值，可能会引入一定的误差[36]，但是相对定量数据分析的最终数据处理是目的基因和内参基因的相对表达量比率，可以抵消一部分试验误差。本研究采用LinRegPCR线性回归分析法作为实时荧光PCR扩增效率分析的方法。与双标准曲线法相比，LinRegPCR法可以进行自动化分析，不需要用稳定的标准品构建多条标准曲线，分析效率高。

新内参基因的发掘需要综合衡量基因表达稳定性及相对表达量分析的结果。持家基因细胞骨架结构蛋白Actin基因家族成员属于传统的内参基因，在多个物种及重金属胁迫条件下能够稳定表达[37-38]。本研究通过基因表达稳定性及相对表达量分析，同样证实Actin基因家族成员Actin7在甘蓝型油菜受到铝胁迫时能够稳定表达，可以作为内参基因。在基因表达稳定性分析中，UBC9的综合排名高于Actin7，但是基因相对表达量分析表明样品组UBC9的相对表达量上调极显著，不适宜选作内参基因。此外，Actin7和新发掘的内参基因ARFA1E、SAP5、UXS3分别属于4类不同的生物学功能范畴，受到共调控的可能性低，针对它们开发的4对新内参基因引物可以提高铝胁迫试验体系实时荧光PCR试验的准确性，为甘蓝型油菜遗传学研究和分子育种研究提供更多种类的内参基因资源。

基因表达的稳定性是相对的，例如组织器官、发育阶段等因素可能会引起基因的特异性表达，但特定试验系统都存在着相对稳定表达的特定基因[10]，因此，需要根据具体的样品类型，通过系统的试验及分析，确定候选的功能基因作为内参基因的适用性。Actin7、ARFA1E、SAP5、UXS3在耐铝性差异显著的2个甘蓝型油菜品种的对照组和处理组中均稳定表达，暗示这些基因没有参与铝毒害的根部响应过程，这些试验结果为这些基因在该试验体系的内参基因适用性提供了科学依据，表明这些基因可以作为不同甘蓝型油菜品种根组织铝胁迫条件下的内参基因。这些基因在根部和地上部(茎、叶、花、角果)的时空表达调控特性以及作为地上部内参基因的适用性，需要后续试验数据的分析和验证。

综上所述，本研究根据候选内参基因序列开发的4对新内参基因引物特异性高，在铝胁迫试验体系的实时荧光PCR试验中有效实用。新发掘的4个内参基因在2个铝胁迫耐受差异性品种苗期根组织的未处理组(0 h)和铝胁迫(100 μmol/L AlCl3，pH值 4.5)的2个时间梯度(3，24 h)处理组中均稳定表达，按表达稳定性综合排名由高至低排序依次为UXS3、SAP5、ARFA1E、Actin7，其中最优组合为SAP5和UXS3。