山麥冬果實花青素生物合成中內參基因的篩選與驗證

干思宸，師 悅，梁立軍

（浙江農林大學 風景園林與建筑學院，浙江 杭州 311300）

山麥冬Liriope spicata為百合科Liliaceae多年生草本植物，在園林綠化中多栽培于林下或林緣半陰處，掩飾裸露土壤，起到補充綠地改善不良景觀的作用。山麥冬屬Liriope植物只有8種，中國栽培6種，其中包含3個特有種，但山麥冬屬植物分布廣泛，除極寒地區及高海拔地區外，中國各省均有分布，其地理分布受人為栽培引種因素影響很大，沒有特定的地理分布規律[1]。山麥冬成熟時果實表皮由綠轉黑，9月結果后觀果時期可長達整個冬季，且其花葶較長多矗立于葉子的上方，易于觀察，具有很高的園林應用價值。目前，針對山麥冬成熟過程中呈色物質及調控基因尚未報道，但花青素合成途徑在植物中是保守的，合成途徑中上游合成基因是決定植物組織能否積累花青素的關鍵[2]，而下游修飾基因的表達常與花青素的積累一致，是加深果色花色的關鍵基因[3-5]。此外，花青素的積累還受轉錄因子的調控，其中以MYB轉錄因子與bHLH轉錄因子最為常見[6]。

用于基因表達定量分析的方法比較多，其中實時熒光定量PCR(RT-qPCR)由于定量準確、成本低且高通量，被廣泛應用于基因表達水平研究。但其結果常受RNA質量、反轉錄效率、引物特異性、初始樣品量及擴增效率等因素的影響[7-8]，需要引入1個或多個表達穩定的內參基因(reference genes, RGs)來評估目的基因的相對表達[9]。在植物學研究中，曾以肌動蛋白(actin，ACT)[10-12]、組蛋白(histone)[11]、蛋白磷酸酶(protein phosphatase，PP2A)[13]、甘油醛-3-磷酸-脫氫酶(glyceraldehyde-3-phosphate dehydrogenase，GAPDH)[12]、泛素結合酶(ubiquitin conjugating enzyme, UBC)[14-15]以及18S核糖體RNA(18S ribosomal RNA，18S)[16]等基因作為內參基因。但是常見的內參基因也并非適用于任何研究，且目前還未見山麥冬內參基因的報道。鑒于此，本研究基于山麥冬轉錄組數據，對山麥冬果實發育中穩定表達的內參基因進行研究，為提高果色轉變關鍵基因RT-qPCR分析的準確性提供科學依據。

1 材料與方法

1.1 材料

在浙江農林大學資源圃，選取生長環境相同，且植株生長狀況良好、長勢整齊的山麥冬，隨機均勻采集15～20株山麥冬植株的各一簇花葶的上、中、下部分果實，基于山麥冬果實生長特性，采集山麥冬幼果期(2020年9月)及成熟期(2020年11月) 2個時期樣品，果實從花葶中取下后立即存于-80 ℃冰箱備用。設置3次生物學重復。

1.2 總 RNA 提取及 cDNA 合成

使用天根離心柱型RNA試劑盒(天根生物科技有限公司)從每個時期樣本中提取總RNA。采用質量分數為1%的瓊脂糖凝膠電泳檢測RNA的完整性。總RNA的純度和質量濃度采用NanoDrop ONE微量核酸蛋白濃度測定儀(Therm，美國)測定。總RNA樣本質量濃度均高于4×10-5ng·L-1以上，總RNA純度 [D(260)/D(280)]為 1.9～2.1 。cDNA 的合成使用 PrimerScript? RT Master Mix cDNA (Perfect Real Time)反轉錄試劑盒，所有樣本總RNA加入量按照3×10-5ng·L-1稀釋至同一質量濃度，cDNA置于-20 ℃冰箱保存。

1.3 候選內參基因的篩選及 RT-qRCR

基于已獲得的山麥冬轉錄組數據及京都基因與基因組百科全書(KEGG)注釋，篩選了多條通路的基因作為內參基因參考庫，包括參與山麥冬果實運輸和分解代謝的基因(SLC36等)，參與代謝過程的基因(PP2C、MGL、PDP、G6PD等)，參與信號傳導與轉運的基因(AUX、GPR107、CNNM等)，參與細胞過程的基因(CFL等)，參與植物免疫的基因(Trx等)，參與遺傳信息處理的基因(UGT、PP2A、EF1-α等)共1 648個，參考前人對內參基因的篩選閾值稍作修改后[11-13]，以每千個堿基轉錄每百萬映射讀取的片段(FPKM)高于5的基因(低表達的難以檢測)、變異系數＜0.1、變化倍數＜0.2為篩選條件，得到前15個候選內參基因(表1)。

表1 山麥冬 15 個候選的內參基因Table 1 15 candidate reference genes of L. spicata

根據轉錄組獲得的核酸序列信息，利用primer 5軟件設計引物，并交由杭州有康生物技術有限公司合成 (表2)。利用 TB Green 染料 (Takara)預反應，體積 20 μL，并使用 LightCycler? 480 Ⅱ型熒光定量PCR 儀 (羅氏，瑞士)進行 RT-qPCR。反應程序：95 ℃ 預變性 5 min；95 ℃ 變性 10 s；60 ℃ 退火延伸 30 s，40個循環。實驗設置3次生物學重復。擴增效率(cDNA稀釋濃度梯度為5-1、5-2、5-3、5-4、5-5)計算公式為E=[10(-1/K)-1]×100%，其中：E為擴增效率，K為斜率。15個候選內參基因的擴增效率為91.7%～108.0%(表2)。

表2 15個候選內參基因的引物序列和擴增子特征Table 2 Primer sequences and amplicon characteristics of 15 candidate reference genes

1.4 內參基因的穩定性分析及驗證

通過4種方法分析內參基因的穩定性：ΔCt值法[17]、geNorm[18]、NormFinder[19]和BestKeeper[20]。利用Excel 2010計算4種方法對候選內參基因幾何平均數的排名，綜合篩選最適的內參基因。同時根據前期轉錄組數據篩選了10種目的基因，涵蓋花青素合成通路上下游基因以及調控基因。這10種基因在轉錄組數據加權共表達分析中屬于中樞基因，表達量高、與花青素相關性強，且在果實成熟過程中顯著上調。目的基因包括C4H、CHS、MT、UFGT、MYB、bHLH，上述基因引物序列及擴增子特征見表3，最后利用 SPSS 19.0 與 Graphpad Prism 8.0 分析及作圖。

表3 10 個目的基因的引物序列和擴增子特征Table 3 Primer sequences and amplicon characteristics of 10 target genes

2 結果與分析

2.1 山麥冬候選內參基因的表達量分析

15個候選內參基因的溶解曲線均為單一峰(圖1)，瓊脂糖凝膠電泳檢測后出現與預期大小一致的單一條帶(圖2)。該結果表明引物具有良好的特異性。

圖1 15 個候選內參基因的溶解曲線Figure 1 Melting curves of fifteen reference genes

圖2 15個候選內參基因PCR擴增產物的瓊脂糖凝膠電泳Figure 2 Agarose gel electrophoresis of the PCR products of the fifteen reference genes

根據原始循環閾值(Ct)分布發現：所有候選內參基因的Ct為15.53～28.81，Ct越高，基因的表達量越低，反之表達量越高。本研究中，EF1-α基因表達量最高，PP2C基因表達量最低，其余基因表達量介于兩者之間。此外，由箱線圖(圖3)跨度可初步判定內參基因的穩定性。PP2C、Trx-1、AUX、PP2A、PDP基因的Ct跨度廣，不穩定，而GPR107、CNNM、EF1-α、G6PD-2、Trx-2基因最為穩定，其中GPR107、CNNM、G6PD-2基因的Ct中位數與平均數接近，即上述基因相對表達量離散程度低，表達更穩定。然而對原始Ct分析內參基因穩定性的不足，還需引入其他的方法。

圖3 15 個候選內參基因的 CtFigure 3 Ct values of the 15 candidate reference genes

2.2 內參基因的穩定性分析

利用ΔCt法、geNorm、NormFider和BestKeeper對15個候選內參基因的穩定性進行分析(表4)。

表4 4種方法評價 15個候選內參基因表達的穩定性Table 4 Expression stability of 15 candidate reference genes evaluated by 4 methods

ΔCt法是在原始Ct值的基礎上，計算每個基因所有樣本與其他基因的Ct值之差，并計算其標準差。一般平均標準差越低，基因穩定性越高。該方法中，EF1-α、PP2C、CFL、CNNM是山麥冬果實發育階段最穩定的內參基因；PDP、G6PD-1、PP2A是最不穩定的內參基因。

geNorm軟件通過平均表達值來描述候選內參基因的穩定性，同時還能計算歸一化因子之間的兩兩變異(Vn/n+1，其中n為可使RT-qPCR結果準確的最少基因數目)。該方法中，所有基因的平均表達值都在1.5以下(穩定內參基因的臨界值)，即該方法判定下的所有基因都可作為內參基因，其中GPR107(0.817)與CNNM(0.847)基因的平均表達值最低，說明最穩定。同時PDP、G6PD-1基因的平均表達值最高，分別為1.390、1.204，最不穩定，這與ΔCt法判定結果一致。此外，利用geNorm計算2個歸一化基因的Vn/n+1，確定適合量化果實生長過程的最優內參基因數目。geNorm首先計算2個最穩定的候選內參基因的歸一化因子值，然后將剩余候選內參基因按其表達穩定性下降的順序依次相加。如果基因之間的Vn/n+1大于或等于0.15，則進行RT-qPCR分析時應該再添加1個基因才能達到可靠的結果，一旦Vn/n+1低于0.15，就不需要添加額外的基因[21]。由圖4可見：從V4/5開始Vn/n+1小于0.15，即需要使用4個內參基因才能得到可靠的RT-qPCR結果。

圖4 精確歸一化的最佳內參基因數量Figure 4 Optimal number of control genes for accurate normalization

NormFinder軟件可分析候選內參基因的兩兩變異性，其中穩定值越小，候選內參基因越穩定。CNNM與GPR107基因的穩定值最小，分別為0.157、0.167，即CNNM與GPR107基因最穩定，這與geNorm分析結果一致；此外，對最差的內參基因評價也與上述2種方法一致：PDP、G6PD-1、AUX是量化果實發育階段最不適合的內參基因。

Bestkeeper與geNorm、NormFinder軟件不同，需導入原始Ct值平均數，計算候選內參基因在所有樣品中的標準差、變異系數、相關系數。一般地，穩定的內參基因擁有低的標準差、變異系數及高的相關系數。在Bestkeeper評價中，與geNorm、NormFinder分析結果一致，CNNM與PDP基因分別還是最穩定與最不穩定的內參基因。除此之外，還發現G6PD-2為該方法中最穩定的內參基因，其標準差與變異系數最低，分別為0.290、1.323，相關系數為0.750。

最后通過幾何平均數對這4種方法的分析結果進行綜合性排序(表5)。根據表5的排名與geNorm推薦的內參基因數目，篩選CNNM、GPR107、EF1-α、G6PD-2作為標準化山麥冬果實RT-qPCR的最優內參組合，PDP為最差內參基因，通過4種算法得出的結果也與最初候選內參基因原始Ct值分布箱線圖分析結果一致。

表5 15 個候選內參基因的綜合排名Table 5 Comprehensive ranking of reference genes for normalization

2.3 內參基因穩定性的驗證

為驗證內參基因的有效性，選擇10種花青素合成結構基因與調控基因作為目的基因。用單一內參基因：最優內參(CNNM)、最差內參(PDP)，及2種內參組合：排名前2位的內參基因(CNNM、GPR107)和排名前4位的內參基因(CNNM、GPR107、EF1-α、G6PD-2)進行歸一化。從圖5可見：在山麥冬果實花青素合成過程中，使用4種內參方式歸一化時，所有的目的基因都上調表達，但變化倍數稍有不同。在山麥冬果實成熟期，使用PDP基因作為內參時，所有目的基因相對表達量均顯著高于其他3類，特別是對轉錄因子bHLH基因的量化時產生嚴重偏差，使用PDP基因與CNNM+GPR107+EF1-α+G6PD-2基因組合作為內參，bHLH基因的相對表達量分別為6.28與15.70，兩者差異高達2.5倍。然而，當使用最優內參基因CNNM進行標準化時，除UFGT基因外，CNNM、GPR107、EF1-α、G6PD-2內參組合無顯著差異，使用CNNM基因標準化時，UFGT相較幼果期上調表達50.71倍，使用4種內參組合時，UFGT上調72.49倍。此外，本研究還分析了候選內參排名前2位的基因(CNNM、GPR107)作為目的基因的表達量，發現選用2種內參基因與geNorm軟件推薦使用4種內參基因，在10個目的基因中均無顯著差異。

圖5 不同內參基因歸一化后10個目的基因的相對表達量Figure 5 Relative expression levels of ten target genes after normalized by different reference genes

從圖6可見：利用最差內參PDP得到的目的基因表達量與4種內參基因組合得到的目的基因表達量相關系數為0.868 6 (P＜0.01)，當使用最優基因CNNM作為內參時，與4種內參組合相關系數可達0.991 6 (P＜0.01)。對2種內參組合與geNorm推薦的4個數目內參組合比較發現：通過這2種方法標準化得到的目的基因相關性可達0.999 9 (P＜0.01)，即僅使用CNNM、GPR107基因作為雙內參也可達到geNorm軟件推薦的4個內參數目組合的效果。

圖6 不同內參基因標準化后10種目的基因表達量的相關性分析Figure 6 Correlation analysis for relative expression levels of ten target genes after normalized by different reference genes

3 討論

山麥冬作為一種優良的地被園林植物及藥用植物，研究多集中于提高栽培技術及塊莖產量，而針對園林觀賞應用的研究較少。在本研究之前沒有山麥冬內參的研究報道，作為沿階草族植物，其近源種也僅有麥冬Ophiopogon japonicus抗逆性研究中曾以微管蛋白基因(tubulin)[22]及Actin[23]作為參考基因。但這2類基因在前期轉錄組篩選中由于變異系數及變化倍數在候選內參中就已經被排除。本研究根據幾何平均數的綜合排名，推薦使用內參基因CNNM、GPR107、EF1-α、G6PD-2作為研究山麥冬花青素合成的最優內參組合。EF1-α、G6PD-2屬于常見的內參基因，在植物生長發育、抗逆反應、代謝合成中已被廣泛應用[24-25]。基于前期轉錄組數據，新型內參基因CNNM、GPR107也可作為RT-qPCR分析的內參基因，CNNM編碼過渡金屬轉運蛋白，可參與多種金屬吸收、排除及區分化[26]，GPR107編碼G蛋白偶聯受體107，廣泛存在于細胞表面的膜蛋白，可參與植物體多種細胞信號轉導及調控機制保守[27]。上述2種基因在山麥冬果實中表達穩定，其相對表達量平均值與中位數相近，離散程度低，且表達量適中，符合內參基因的標準。在觀賞植物中，由于新型內參基因穩定性強于傳統內參基因，常被選用標準化目的基因的表達。例如，在異型花柱連翹Forsythia suspensa中，轉錄組中變化微小的未知基因是研究花開放最適合的內參基因[28]；太行花Taihangia rupestris花器官有復雜的性別決定機制，鑒定兩性花與雄性花的內參基因是編碼鐵硫簇組裝蛋白、3-巰基丙酮酸硫轉移酶與跨膜蛋白50的新型內參基因[11]。SmDnaJ基因在旱柳Salix matsudana各種非生物脅迫下表達最為穩定[29]。bHLH在觀賞百合Lilium oriental×Trumpet hybrid體胚誘導、體胚發育中表達最穩定[30]，但bHLH是植物顏色育種中的重要靶基因，并不適合作為本研究的內參基因候選，這也證實了不同目標性狀需采用不同的內參基因，沒有一種內參基因是普適的。

花青素合成路徑在植物中是保守的，其中MYB轉錄因子與bHLH轉錄因子可形成二元復合體，激活花青素合成酶基因[31-32]。大量研究表明：MYB、bHLH轉錄因子基因與花青素合成酶基因在紫色系植物組織發育過程中協同上調[3, 33-34]。為驗證內參基因的結果，挑選了10個在山麥冬花青素合成調控網絡的中樞基因(相關性強且表達量高)作為驗證，其中包括轉錄因子與結構基因(C4H、CHS、MT、UFGT、MYB、bHLH)，這10種基因在4種歸一化方法下表達模式均顯著上調，但趨勢稍有不同，選用較差內參PDP標準化結果偏差最大，在山麥冬成熟黑果中所有基因都顯著高于其他基因。盡管最優內參基因CNNM對目的基因的歸一化可以達到與4種內參組合很高的相關系數，但對UFGT基因的量化存在顯著差異，而UFGT基因作為花青素合成通路的下游修飾，對花青素積累至關重要，特別是在山麥冬這類組織顏色深即富含花青素的類型[2, 35]，例如在葡萄Vitis vinifera果皮[36]、玫瑰Rosa rugosa[37]、紫皮石刁柏Asparagus officinalis[33]中UFGT都被驗證為關鍵基因，因此僅選用單一基因作為研究山麥冬果皮花青素積累的內參是不合適的，繼而在CNNM基因基礎上又引入GPR107來規避單內參基因的誤差，該內參組合與geNorm推薦的內參組合相關系數最高，在10種目的基因的驗證結果中與4種內參組合均無顯著差異，且選用雙內參組合比4種內參組合可操作性強，因此判定使用CNNM、GPR107作為雙內參即可得到可靠的RT-qPCR結果。雙內參組合聯合使用可以減少實驗因素對基因表達的影響，且結果更為準確。暴露于UV-B輻射下的番茄Lycopersicon esculentum幼苗不同組織都應選用特定的內參組合，例如葉中選用肌動蛋白基因與微管蛋白基因，而根中選用微管蛋白與UV-B抗性位點基因更加適合[38]；UBQ和EF1-α基因由于表達穩定，可作為內參基因用于鵝掌草Anemone flaccida各器官的不同發育階段[39]。

4 結論

本研究基于轉錄組數據篩選了15個候選內參基因，分析其在山麥冬果實不同時期的表達穩定性。經過10種目的基因驗證后，表明以CNNM、GPR107基因作為組合是山麥冬果實花青素生物合成研究的最佳內參基因，而常用的內參基因卻并不適用于本研究，這為篩選新型內參基因提供了新思路。