侯穎輝, 王少銘, 羅莉斯, 李晉華, 冷家歸, 汪志燚, 李德文
(1.貴州省農(nóng)業(yè)科學院香料研究所, 貴陽 550006; 2.貴州省農(nóng)業(yè)科學院油料研究所, 貴陽 550006)
小黃姜(Dioscoreazingiberensis)是生姜的一種,因其根狀莖較小,顏色金黃而得名,與普通生姜相比,其姜辣素、姜黃素等活性成分含量更高,更具藥食及研究價值[1]。貞豐小黃姜、水城小黃姜已成為地理標志性產(chǎn)品。由于小黃姜為無性繁殖,播種需要消耗大量姜塊,加上生姜連作障礙嚴重,易發(fā)姜瘟病,小黃姜的產(chǎn)量相對較低,為了提高產(chǎn)量和減少成本,小黃姜莖尖離體培養(yǎng)技術(shù)被廣泛研究[2]。據(jù)報道,莖尖離體培養(yǎng)技術(shù)能夠有效降低植株體內(nèi)病毒數(shù)量,降低病毒病的發(fā)生率[3-4]。通過小黃姜莖尖離體培養(yǎng)獲得的小黃姜組培苗栽培后產(chǎn)生了一定頻率的表型突變,突變率為10%左右,突變株植株矮小、分枝數(shù)多、產(chǎn)量較低。本研究主要從品質(zhì)表現(xiàn)及轉(zhuǎn)錄組學兩個方面分析突變株與正常株的主要差異,揭示突變產(chǎn)生的機制,以期在育苗期篩選出突變株,從而減少后續(xù)實驗投入和栽培成本。同時,對小黃姜種質(zhì)資源創(chuàng)新具有重要理論意義。
以貴州地理標志性產(chǎn)品鎮(zhèn)寧小黃姜為母本,通過莖尖離體培養(yǎng)技術(shù)得到小黃姜突變株和正常株原種,按照每小區(qū)30株定植,株距15 cm,行距30 cm。使用采收后的鮮姜進行品質(zhì)檢測,在根狀莖膨大期后期選擇健康且長勢一致的脫毒正常株和突變株根狀莖,在北京百邁客生物科技有限公司進行轉(zhuǎn)錄組分析。
1.2.1品質(zhì)檢測
采用恒溫干燥法[5]檢測樣品中干物質(zhì)含量,水蒸氣及GC-MS技術(shù)[6]分析樣品中精油含量,試劑盒法檢測纖維素和木質(zhì)素含量,有機溶劑萃取法和HPLC技術(shù)[7]分析姜辣素(以6-姜酚計)含量。
1.2.2轉(zhuǎn)錄組學測序
Trizol法提取樣品RNA,利用邊合成邊測序技術(shù)與Illumina Hiseq高通量測序平臺[8]和分析cDNA文庫,獲得的高質(zhì)量Reads再利用生物信息學分析,通過與NR、KEGG、Swiss-Prot、COG、GO和Pfam數(shù)據(jù)庫[9]的比對Unigene進行功能注釋。
通過品質(zhì)檢測(圖1)發(fā)現(xiàn),小黃姜突變株干物質(zhì)、纖維素含量顯著低于正常株,分別為13%和0.6%,而木質(zhì)素、精油及姜辣素含量則顯著高于正常株,分別為0.63%、0.043 mL/g和 8.5 mg/g,分別為正常株的1.5倍、1.54倍和1.26倍。
圖1 小黃姜突變株與正常株品質(zhì)比較Fig.1 Comparison of quality between mutant and normal ginger of Dioscorea zingiberensis
各樣品Q30堿基百分比均在93.79%及以上,Clean Data均達到6.01 Gb,6個樣品共獲得46.33 Gb Clean Data。組裝后共獲得Unigene 100 280條,長度在1 kb以上的有30 687條。
2.2.1測序數(shù)據(jù)產(chǎn)出統(tǒng)計
如表1 所示,不同樣品的Clean Data與組裝的Unigene庫進行比對得知突變株映射比率為59.54%~64.52%;正常株的映射比率為56.74%~64.55%。
表1 Clean Data與組裝的Unigene庫的比對結(jié)果Table 1 Comparison results of Clean Dat and assembly Unigene
2.2.2Unigene功能注釋及表達量分析
使用BLAST[10]軟件(E-value≤1e-5)將獲得的Unigene序列分別與GO[13]、COG[14]、KOG[15]、eggNOG4.5[16]等多個數(shù)據(jù)庫進行比對分析, 并使用KOBAS2.0[17]分析Unigene序列在KEGG中的功能注釋,氨基酸序列預(yù)測之后再通過HMMER[18]軟件(E-value≤1e-10)在Pfam[19]數(shù)據(jù)庫分析比對[9],從而獲得不同序列的注釋信息。最終獲得有注釋信息的Unigene 43 989個(表2)。
表2 Unigene注釋統(tǒng)計Table 2 Unigene annotation statistics
將測序得到的Reads與Unigene庫通過Bowtie[20]軟件進行比對,并結(jié)合RSEM[21]估計其表達量水平。Unigene的表達豐度[22]利用FPKM值[23](Fragments Per Kilobase of transcript per Million mapped reads)表示。由圖2可知,蛋白質(zhì)編碼基因表達水平FPKM值橫跨10-2到104六個數(shù)量級。另外,從圖3中可以得知單個樣品基因表達水平分布的離散程度和整體基因表達豐度。
注:a為FPKM density distribution;b為FPKM box diagram。圖2 各樣品的密度分布對比以及FPKM箱線圖Fig.2 Comparison of density distribution among different samples and FPKM box diagram and FPKM box diagram
圖3 SSR密度分布Fig.3 SSR density profile
2.2.3SSR及SNP分析
利用MISA(MIcroSAtellite identification tool)軟件對1 kb 以上的 Unigene序列進行分析,鑒定出6種類型的SSR[24](表3),合計13 342個,并對不同類型的SSR進行密度分布統(tǒng)計(圖3)。
表3 SSR分析結(jié)果統(tǒng)計Table 3 Statistical table of SSR analysis results
利用軟件STAR[25]和GATK[26]針對RNA-Seq識別單核苷酸多態(tài)性(Single Nucleotide Polymorphism,SNP)位點。根據(jù)SNP位點的等位(Allele)數(shù)目將SNP位點分為純合型SNP位點(HomoSNP,只有一個等位)和雜合型SNP位點(HeteSNP,兩個或多個等位)。各樣品SNP位點數(shù)目統(tǒng)計見表4,并對SNP在Unigene上的分布密度進行統(tǒng)計(圖4)。
圖4 SNP密度分布Fig.4 SNP density profile
表4 SNP數(shù)量統(tǒng)計Table 4 SNP number statistics
2.2.4DEG分析
采用DESeq2[27]軟件結(jié)合FDR(False Discovery Rate) 與FC(Fold Change)對樣品進行組間差異表達分析,共得到DEG(Differentially Expressed Genes)338條,其中164條基因表達量上調(diào),174條基因表達量下調(diào)(圖5)。一個點代表一個基因,橫坐標為基因在兩樣品中表達量差異倍數(shù)的對數(shù)值,兩樣品間的表達量倍數(shù)差異越大,對應(yīng)的絕對值越大。并對篩選出的DEG做層次聚類分析(圖6)。通過與不同功能數(shù)據(jù)庫比對,獲得注釋DEG 239條,分別在COG、GO、KEGG等數(shù)據(jù)庫[9]注釋到88條、132條、81條、130條、186條、171條、224條、238條(表5)。
注:綠色為下調(diào)DEG,紅色為上調(diào)DEG,黑色為無顯著性表達差異基因。圖5 差異表達基因MA圖Fig.5 MA map of DEGs
圖6 差異表達基因表達模式聚類Fig.6 Clustering map of differential expression gene expression patterns
圖7 DEG的KOG富集分析Fig.7 KOG enrichment analysis of DEG
圖8 DEG的eggNOG富集分析Fig.8 eggNOG enrichment analysis of DEG
圖9 DEG 的GO二級節(jié)點注釋統(tǒng)計Fig.9 Annotated statistical map of GO secondary nodes of DEGs
圖10 DEG的KEGG分類Fig.10 KEGG classification of DEGs
表5 注釋的差異表達基因數(shù)量統(tǒng)計Table 5 Annotated statistical on the number of differentially expressed genes
通過不同功能數(shù)據(jù)庫的富集分析發(fā)現(xiàn),DEG主要集中于次生代謝合成、信號傳導機制、類黃酮代謝、光合作用碳的固定等過程。主要DEG的功能分析見表6,其中植物激素信號傳導相關(guān)基因4個(AXU/IAA、B-ARR、ERFL1/2、MYC2),光合作用相關(guān)基因1個(GOT1),植物生理節(jié)律相關(guān)基因2個(CRY、CHS),類黃酮生物合成相關(guān)基因1個(CHS)。除B-ARR、CRY為下調(diào)外,其余基因均上調(diào)。
圖11 DEG 的KEGG通路富集散點圖Fig.11 Enrichment scatter plot of KEGG pathway of DEGs
轉(zhuǎn)錄學分析是在mRNA水平上研究表型差異的有效方法,轉(zhuǎn)錄組測序技術(shù)能夠準確獲得不同組織,不同時期相關(guān)基因的表達水平??椎抡娴萚28]利用轉(zhuǎn)錄組學分析鑒定出鹽脅迫、堿脅迫相關(guān)的DEG 900多個,陳浩等[29]基于轉(zhuǎn)錄組測序獲得SSR引物10 320對,鄭世茂等[30]基于轉(zhuǎn)錄組測序獲得SSR引物37 173對,成功率為70.39%。楊博涵等[31]利用轉(zhuǎn)錄組測序技術(shù)分析了陽春砂花絲、花柱,得到62 174條Unigene注釋,涉及25個功能分類。
本研究中,小黃姜突變株與正常株來源于同一植株的莖尖組織,通過組織培養(yǎng)后產(chǎn)生不同的表型,經(jīng)檢測發(fā)現(xiàn),其品質(zhì)與正常株和母體均具有顯著差異。通過轉(zhuǎn)錄組分析發(fā)現(xiàn),注釋到的差異表達基因239條,分別在COG、GO、KEGG、KOG、Pfam、Swiss-Prot、eggNOG功能數(shù)據(jù)庫注釋到88條、132條、 81條、130條、186條、171 條和224條,主要涉及光合作用碳的固定、次生代謝物合成、植物信號傳導等20余個功能分類。其中,光合作用碳的固定與干物質(zhì)含量、纖維素和木質(zhì)素含量關(guān)系密切[32],根據(jù)光合作用KEGG途徑推測可能GOT1基因的高表達會促進草酰乙酸轉(zhuǎn)化為天冬氨酸,減少了蘋果酸的形成,從而降低碳水化合物的積累。生姜精油、姜辣素等次生代謝物含量差異主要由于次生代謝物合成相關(guān)基因的差異表達導致[33]。植物激素信號傳導調(diào)控植物的整個生命過程,包括細胞分裂和分化,是植株分枝/分蘗的主要內(nèi)在因素[34],該通路相關(guān)基因的差異表達導致了突變株與正常株之間分枝數(shù)與產(chǎn)量的差異,尤其AXU/IAA,作為生長素應(yīng)答因子ARF的抑制因子[35],其高表達抑制了植株根狀莖的生長,從而導致產(chǎn)量下降。獲得DEG以及SNP序列有助于在栽培前期篩選出突變株,以減少栽培的時間和物質(zhì)成本。
另外,由于小黃姜為無性繁殖作物,種質(zhì)創(chuàng)新較難,該突變株的產(chǎn)生對于小黃姜綜合利用和種質(zhì)資源創(chuàng)新具有重要理論意義。