玉米灰斑病病原菌脅迫下玉米轉(zhuǎn)錄組測序的SNP和InDel分析

蔣 滔, 張曉梅, 楊 梅, 張 志, 邱紅波

(1.貴州大學(xué)農(nóng)學(xué)院, 貴陽 550025; 2.云南省德宏州農(nóng)業(yè)技術(shù)推廣中心, 云南 芒市 678400)

玉米灰斑病(Gray leaf spot,GLS) 是玉米十大常見病害之一,它是由病菌Cercosporazeae-maydis和Cercosporazeina引起的,是全球玉米種植區(qū)最嚴(yán)重的玉米葉病害之一[1-2]。2009年玉米全基因組測序的完成促進(jìn)了玉米分子標(biāo)記輔助育種研究,也推動了玉米新品種的培育。同時,玉米灰斑病抗性屬于數(shù)量性狀,Pedro等[1]、Nsibo等[2]通過不同類型的分子標(biāo)記對其展開了研究,其中圖位克隆的方式克隆未知抗病基因成為當(dāng)下研究熱點(diǎn),而SNP和InDel分子標(biāo)記作為圖位克隆的重要工具,通過對轉(zhuǎn)錄組測序結(jié)果進(jìn)行SNP和InDel標(biāo)記分析,可為開發(fā)篩選玉米抗病種質(zhì)資源標(biāo)記奠定基礎(chǔ)。

SNP/InDel作為新基因或等位基因變異驅(qū)動因素的寶貴來源,當(dāng)結(jié)果表型表現(xiàn)出有利特征時,可以通過自然或人工方式選擇這些優(yōu)良變異[3]。利用下一代測序技術(shù)進(jìn)行基于轉(zhuǎn)錄組測序的標(biāo)記搜索,可以快速發(fā)現(xiàn)專注于編碼區(qū)域的多態(tài)性,避免了高度重復(fù)的基因組區(qū)域,使其成為檢測變異/突變的強(qiáng)大工具[4-5]。盡管從轉(zhuǎn)錄組測序數(shù)據(jù)中獲得了大量的研究成果和豐富的信息,但其在識別SNP和InDel等分子標(biāo)記方面的應(yīng)用仍然很少[6]。SNP/InDel是植物基因組中最豐富的DNA變異,因其高效性而被廣泛應(yīng)用于植物改良和作物育種[7-8]。與其他分子標(biāo)記相比,SNP和InDel便于基因組序列或轉(zhuǎn)錄組序列的比較[9-10]。單核苷酸多態(tài)性 (SNP) 標(biāo)記是單個堿基發(fā)生變異的多態(tài)性序列,具有突變率低、遺傳穩(wěn)定性高及可以自動化檢測等特點(diǎn),使標(biāo)記在廣泛應(yīng)用時更為有效[11-12]。插入/缺失多態(tài)性(InDel)標(biāo)記是根據(jù)核苷酸片段的插入或缺失而開發(fā)的,有分布廣、可重復(fù)性高、密度高、成本較低、變異率低、多態(tài)性強(qiáng)且易于檢測等優(yōu)點(diǎn)[13-14]。目前,SNP和InDel已應(yīng)用于研究許多優(yōu)良物種的遺傳多樣性分析、遺傳圖譜構(gòu)建、關(guān)聯(lián)作圖分析和標(biāo)記輔助選擇(MAS)育種[15-16]。隨著使用新一代測序技術(shù)收集越來越多的基因組或轉(zhuǎn)錄組數(shù)據(jù),開發(fā)SNP和InDel變得越來越容易。到目前為止,已經(jīng)使用新一代測序技術(shù)獲得了很多與玉米抗病過程相關(guān)的轉(zhuǎn)錄組數(shù)據(jù),并在NCBI數(shù)據(jù)集中發(fā)布[17-18]。本研究使用可用的玉米轉(zhuǎn)錄組數(shù)據(jù)集識別和評估玉米灰斑病病原菌脅迫下SNP和InDel標(biāo)記,以玉米灰斑病抗、感自交系T 32和J 51為材料,通過高通量測序技術(shù),獲得了大量轉(zhuǎn)錄組信息,對SNP和InDel標(biāo)記進(jìn)行分析,旨在開發(fā)出更豐富的分子標(biāo)記,為抗灰斑病的SNP 和InDel標(biāo)記的后續(xù)開發(fā)利用,為玉米抗灰斑病分子標(biāo)記輔助育種以及候選基因功能分析提供理論參考。

1 材料與方法

1.1 試驗(yàn)材料

選用玉米(ZeamaysL.)灰斑病抗、感自交系T 32和J 51為試驗(yàn)材料,均由貴州大學(xué)玉米研究所提供,于2020年在云南省芒市(98.58°E,24.43°N)玉米灰斑病重發(fā)區(qū)種植玉米材料作為玉米灰斑病病菌自然接菌處理,在芒市無灰斑病發(fā)生地種植作為不接菌處理。在抽雄期進(jìn)行穗位葉取樣并送往深圳華大生物公司進(jìn)行轉(zhuǎn)錄組測序。

1.2 試驗(yàn)方法

使用擎科生物科技有限公司的植物總RNA提取試劑盒提取RNA,詳細(xì)步驟參照說明書。對用于轉(zhuǎn)錄組測序的6組樣本材料(T 32的3個生物學(xué)重復(fù),命名為T 32_Repeat_1、T 32_Repeat_2和T 32_Repeat_3;J 51的3個生物學(xué)重復(fù),命名為J 51_Repeat_1、J 51_Repeat_2、J 51_Repeat_3)的總RNA進(jìn)行mRNA純化、反轉(zhuǎn)錄、轉(zhuǎn)錄組測序(在深圳華大基因有限公司完成)等。

測序后的數(shù)據(jù)通過測序平臺的軟件將其轉(zhuǎn)換為FASTQ格式的原始數(shù)據(jù),再對每個樣品的原始數(shù)據(jù)分別進(jìn)行整理統(tǒng)計,包括樣品名稱、Q 20、Q 30和GC等。使用過濾軟件SOAPnuke對原始數(shù)據(jù) (Raw reads) 進(jìn)行過濾處理,把獲得的高質(zhì)量數(shù)據(jù) (Clean reads) 比對到參考基因組 (B 73_Ref Gen_V 4) 上,使用GATK檢測樣品中的SNP和InDel信息,利用HISAT 2軟件進(jìn)行數(shù)據(jù)預(yù)處理;利用Haplotype Caller對SNP和InDel變異位點(diǎn)進(jìn)行檢測[19]。最后統(tǒng)計SNP和InDel的變異類型及分布情況,采用Excel 2010軟件進(jìn)行數(shù)據(jù)分析與作圖。

1.3 玉米轉(zhuǎn)錄組測序的SNP 和InDel分析

經(jīng)過處理后的轉(zhuǎn)錄組數(shù)據(jù)通過Varscan軟件獲取SNP和InDel位點(diǎn),再對得到的SNP和InDel進(jìn)行統(tǒng)計分析,包括: 1) 樣品轉(zhuǎn)錄組測序質(zhì)量分析; 2) SNP特征分析; 3) InDel分析; 4) SNP/InDel在玉米基因組中Up 2 k,Exon,Intron,Down 2 k,Intergenic 等5種功能元件上的區(qū)域分布; 5) 差異表達(dá)基因GO 和 KEGG 功能分類分析。

2 結(jié)果與分析

2.1 樣品轉(zhuǎn)錄組測序質(zhì)量分析

測序后6組樣品共獲得39.68 G的可用數(shù)據(jù),經(jīng)過濾得到T 32和J 51最終可用的reads 數(shù)量分別是130 369 594個和134 221 272個,占原始數(shù)據(jù)的95.8%和96.7%。堿基質(zhì)量值Q 20大于94.5%,Q 30大于88.15%(表1)。兩組樣本通過triniy組裝得到61 947條轉(zhuǎn)錄本,平均長度為1 639個核苷酸,用Tgicl 軟件去除冗余得到37 268條Unigene,平均長度為1 728個核苷酸。將所有的Unigene 按照從長到短排序并依次相加,當(dāng)相加長度等于 Unigene總長度的一半時,所對應(yīng)的Unigene的長度即為N 50值。本實(shí)驗(yàn)所得轉(zhuǎn)錄本N 50等于2 307 nt,Unigene的N 50等于2 040 nt(表2)。Unigene長度在500～2 100 nt范圍的數(shù)量最多,占69%,隨著Unigene長度的增加,數(shù)量呈梯度下降,小于300 nt的數(shù)量最少(圖1)。

表1 樣品轉(zhuǎn)錄組測序質(zhì)量統(tǒng)計Table 1 Sample transcriptome sequencing quality statistics

表2 轉(zhuǎn)錄本和Unigene的組裝結(jié)果Table 2 Assembly results of transcripts and Unigene

圖1 Unigene長度統(tǒng)計Fig.1 Unigene size statistics

2.2 SNP特征分析

使用GATK軟件對兩組樣品中的SNP進(jìn)行分析,共獲得109 977個SNP。為保證SNP位點(diǎn)的準(zhǔn)確性,篩選SNP應(yīng)保證兩個轉(zhuǎn)錄本的覆蓋率大于20個重疊群之和,候選SNP位點(diǎn)兩側(cè)至少有5 bp的保守序列。根據(jù)上述條件,篩選出A/C,A/G,A/T,C/T,C/G,G/T六種堿基突變類型,其中轉(zhuǎn)換類型71 352個,顛換類型38 625個,轉(zhuǎn)換類型是顛換類型的1.85倍。在轉(zhuǎn)換位點(diǎn)中,C/T 和A/G轉(zhuǎn)換分別有37 894個和33 458個位點(diǎn);在顛換位點(diǎn)中,A/T、C/G、G/T和A/C顛換分別有8 125個、10 652個、10 896個和8 952個位點(diǎn)(圖 2)。

圖2 基于轉(zhuǎn)錄組測序的玉米不同SNP基因型統(tǒng)計Fig.2 Statistics of different SNP genotypes of maize based on transcriptome sequencing

SNP跨染色體和基因的分布對于評估其基因組覆蓋度和標(biāo)記密度特別重要。 因此,分析了所有玉米染色體的SNP分布(圖3)。結(jié)果表明,在轉(zhuǎn)錄組數(shù)據(jù)中鑒定出的109 977個SNP位點(diǎn)分布在10條染色體中,每間隔585 bp出現(xiàn)1個SNP。位于第1染色體上的SNP數(shù)量最多,占16.11%,其次是第5染色體上,占12.78%,第10染色體分布最少,占7%,平均每條染色體上有10 997.7個SNP。

圖3 不同染色體轉(zhuǎn)錄組 SNP 位點(diǎn)密度分布頻率Fig.3 Distribution frequency of SNP site density in different chromosome transcriptomes

2.3 InDel分析

T 32共有12 147個插入位點(diǎn),13 141個缺失位點(diǎn),J 51共有9 964個插入位點(diǎn),11 130個缺失位點(diǎn)。在T 32的10條染色體上,InDel位點(diǎn)在第1染色體最多,為4 214個,第9染色體最少,為1 620個,平均每條染色體有2 528個。在J 51的10條染色體上,InDel

注:A為J 51的SNP位點(diǎn)在玉米基因組上的區(qū)域分布情況;B為T 32的SNP位點(diǎn)在玉米基因組上的區(qū)域分布情況;C為J 51的InDel位點(diǎn)在玉米基因組上的區(qū)域分布情況;D為T 32的InDel位點(diǎn)在玉米基因組上的區(qū)域分布情況(Up 2 k表示基因上游2 kbp;Down 2 k表示基因下游2 kbp;Intron表示內(nèi)含子區(qū);Exon表示外顯子區(qū);Intergenic表示基因間隔區(qū))。圖5 SNP/InDel在玉米基因組中的區(qū)域分布Fig.5 Regional distribution of SNP/InDel in the maize genome

位點(diǎn)在第1染色體最多,為3 467個,第10染色體最少,為1 389個,平均每條染色體有2 111個(表3)。InDel插入/缺失片段以1～3 bp長度為主,片段從小到大數(shù)量依次減少,均僅在6 bp長度時有小幅度上升,在大于11 bp長度之后逐漸下降,長度大于20 bp的數(shù)目占比不超過1%(圖4)。

圖4 InDel 類型統(tǒng)計Fig.4 InDel type statistics

2.4 SNP/InDel的區(qū)域分布

SNP/InDel 位點(diǎn)在基因組的5個位置區(qū)域(Up 2 k,Exon,Intron,Down 2 k,Intergenic)的分布(圖5),在外顯子(Exon)區(qū)域分布最多,占比50.7%～56%,超過了SNP/InDel位點(diǎn)總數(shù)的一半,其次是內(nèi)含子(Intron)區(qū)域,占比分別為32.6%、32.3%、38.8%、39.2%;在基因的間隔區(qū)(Intergenic),兩材料SNP位點(diǎn)的占比分別為8.7%、8.6%%,InDel位點(diǎn)的占比均為5.9%;其余基因上下游(Up 2 k,Down 2 k)區(qū)域內(nèi)SNP/InDel位點(diǎn)數(shù)量最少,不超過總量的5%。

表3 T 32、J 51兩自交系InDel在染色體上的分布Table 3 Chromosome distribution of InDel,T 32 and J 51 inbred lines 單位:個

2.5 差異表達(dá)基因GO和KEGG 功能分類分析

通過使用BLAST針對NCBI參考序列蛋白質(zhì)數(shù)據(jù)庫的序列相似性搜索,對得到的最終組裝轉(zhuǎn)錄本進(jìn)行驗(yàn)證和注釋?；谛蛄型葱垣@得的61 947條轉(zhuǎn)錄本的功能注釋,在分析中轉(zhuǎn)錄本被分為3個標(biāo)準(zhǔn)類別:生物過程通路、細(xì)胞成分通路和分子功能通路。在生物過程中,差異表達(dá)基因富集在代謝過程和細(xì)胞過程,而免疫系統(tǒng)過程、細(xì)胞增殖、碳利用、細(xì)胞殺傷、氮利用和硫的利用等沒有差異表達(dá)基因富集。細(xì)胞成分通路分為14個小類,在細(xì)胞成分中,差異表達(dá)基因也富集在細(xì)胞和細(xì)胞部分,而超分子復(fù)合物、細(xì)胞外區(qū)域部分和類核素等只有幾個或沒有差異表達(dá)基因富集。在分子功能中,差異表達(dá)基因與催化活性和黏合相關(guān)性顯著,而與營養(yǎng)庫活性、分子載體活性、蛋白標(biāo)記和翻譯調(diào)節(jié)活性相關(guān)性不顯著(圖6)。玉米轉(zhuǎn)錄組數(shù)據(jù)庫中含SNP/InDel標(biāo)記基因經(jīng)KEGG數(shù)據(jù)庫注釋后,差異表達(dá)基因被定位到KEGG中的規(guī)范參考通路途徑。所有獨(dú)特的差異表達(dá)基因被分配到5個第一層級通路途徑和第二層級34個KEGG通路途徑。其中全局和概覽圖在代謝相關(guān)通路途徑中占主導(dǎo)地位;翻譯和折疊、分類和降解在遺傳信息過程中占主導(dǎo)地位。此外,大多數(shù)差異表達(dá)基因由環(huán)境信息處理途徑中的信號轉(zhuǎn)導(dǎo)表示,而運(yùn)輸和分解代謝、細(xì)胞生長和死亡在細(xì)胞過程途徑中占主導(dǎo)地位,有機(jī)系統(tǒng)途徑中環(huán)境適應(yīng)占主導(dǎo)地位(圖7)。

圖7 玉米轉(zhuǎn)錄組序列中SNP/InDel基因的KEGG功能分類Fig.7 KEGG functional classes of SNP/InDel genes in the maize transcriptome sequence

3 結(jié)論與討論

在下一代測序數(shù)據(jù)中發(fā)現(xiàn)的SNP和InDel已廣泛用于人類和許多其他動物物種以及植物的基因組和轉(zhuǎn)錄組分析, 如在棉花[20]、水稻[21]和大麥[22]等植物中大量應(yīng)用。然而,在玉米病原菌脅迫中的應(yīng)用很少,以往的應(yīng)用研究主要局限于SSR標(biāo)記的開發(fā)[23],SNP 和 InDel 標(biāo)記很少研究,特別是在玉米灰斑病病原菌脅迫下,玉米中用于轉(zhuǎn)錄組測序的SNP和InDel分子標(biāo)記的開發(fā)未有報道。

本研究探索了6組樣品在玉米灰斑病病原菌脅迫下獲得的39.68 Gb高質(zhì)量轉(zhuǎn)錄組序列數(shù)據(jù),并對轉(zhuǎn)錄組測序數(shù)據(jù)進(jìn)行SNP和InDel分析。收集了玉米灰斑病抗病自交系T 32和感病自交系J 51的轉(zhuǎn)錄組數(shù)據(jù)集,并得到61 947條轉(zhuǎn)錄本,平均長度為1 639個核苷酸。從轉(zhuǎn)錄組數(shù)據(jù)集中生成了37 268條非冗余轉(zhuǎn)錄本,結(jié)果表明,單個轉(zhuǎn)錄組不包含玉米中所有基因的表達(dá)圖譜,這在小麥[24]、油菜[25]和水稻[26]中也觀察到。此外,在探索基因表達(dá)圖譜的過程中,可能不會檢測低表達(dá)水平的轉(zhuǎn)錄本,一些不足以產(chǎn)生全長組裝的轉(zhuǎn)錄本也會影響分析[26]。由于轉(zhuǎn)錄組代表基因的時間和空間表達(dá)圖譜,因此參考轉(zhuǎn)錄組由不同的基因型對于使用轉(zhuǎn)錄組數(shù)據(jù)開發(fā)SNP和InDel 至關(guān)重要[27]。同時,在開發(fā)標(biāo)記之前評估轉(zhuǎn)錄數(shù)據(jù)的質(zhì)量是必要的。

SNP是植物基因組中最豐富的DNA標(biāo)記,已廣泛用于遺傳研究和育種計劃。在本實(shí)驗(yàn)中,兩種玉米材料樣品統(tǒng)計分析共檢測到17 499個Unigenes有SNP位點(diǎn),共檢索到109 977個SNP位點(diǎn),平均每條Unigene大約有2.9個SNP,每間隔585 bp 出現(xiàn)1個SNP;其中轉(zhuǎn)換類型71 352種,占64.88%,顛換類型38 625種,占35.12%,轉(zhuǎn)換類型約為顛換類型的1.85倍,轉(zhuǎn)換類型明顯大于顛換類型,因?yàn)閴A基替換是產(chǎn)生遺傳變異和推動進(jìn)化的機(jī)制之一,表明轉(zhuǎn)錄組序列堿基的轉(zhuǎn)換變異存在特定的作用機(jī)制[28]。此外,本研究共有6種SNP類型,并且2種轉(zhuǎn)換類型的SNP位點(diǎn)數(shù)量遠(yuǎn)大于4種顛換類型,這是由于基因組中的不同堿基突變頻率存在差異造成的[29]。高頻率的轉(zhuǎn)換可能反映了甲基化后高水平的C/T突變[30];同時,轉(zhuǎn)換/顛換的高比率表明基因組比較中的遺傳差異水平較低[31-32]。

InDel 是動物、植物和細(xì)菌中密切相關(guān)的DNA序列之間序列差異的主要來源。本研究InDel分析結(jié)果表明,在11 346條Unigene上發(fā)現(xiàn)29 849個InDel位點(diǎn),InDel在每條Unigene上的出現(xiàn)頻率為48%,其中InDel缺失位點(diǎn)數(shù)量略大于插入位點(diǎn)數(shù)量,大部分InDel插入/缺失片段以1～3 bp長度為主,長度大于20 bp的數(shù)目占比不超過1%。在一定程度上,較長的 InDel 在基因組上的分布密度相對較低,從而導(dǎo)致較低的多態(tài)性產(chǎn)生,這在番茄[33]中得到了證實(shí)。此外,InDel可能會導(dǎo)致移碼突變,使mRNA 在翻譯過程中出現(xiàn)錯誤的終止密碼子造成遺傳變異。

SNP和InDel的分布也因序列區(qū)域類型而異,在基因間區(qū)域的分布頻率相對低于基因區(qū)域[34-35]。劉小紅[35]研究表明,超過60%的總SNP和InDel存在于玉米的基因區(qū)域,其中編碼區(qū)的外顯子區(qū)分布超過50%。本研究中,篩選了基因組的所有區(qū)域也得到相同結(jié)果,發(fā)現(xiàn)兩材料SNP/InDel 位點(diǎn)在外顯子區(qū)域分布最多(超過50%),其次是內(nèi)含子區(qū)域(超過35%),其他區(qū)域內(nèi)分布不超過總量的5%。此外,在包括外顯子和內(nèi)含子序列在內(nèi)的基因區(qū)的SNP和InDel分布數(shù)量,內(nèi)含子區(qū)域約占小麥[36]和大豆[37]等作物中總SNP和InDel的一半。

在本研究中,61 947個轉(zhuǎn)錄本被分配到3個主要的GO類別。Hufford等[38]研究表明,大多數(shù) GO 代表的轉(zhuǎn)錄與催化活性和黏合(在分子功能中)、細(xì)胞和細(xì)胞部分(在細(xì)胞成分中)以及代謝和細(xì)胞過程相關(guān)的轉(zhuǎn)錄物(在生物過程中)在玉米基因型中顯示相似的基因功能(圖 6)。通過KEGG通路的分析有助于進(jìn)一步了解基因功能的生物學(xué)相關(guān)性[39],基于KEGG通路數(shù)據(jù)庫,遺傳信息過程和代謝被表示為優(yōu)勢通路(圖7)。遺傳信息處理途徑由折疊、分類和降解、翻譯、轉(zhuǎn)錄、復(fù)制和修復(fù)等4個亞類組成,同時RNA家族在生物活動中起著最重要的作用。在轉(zhuǎn)錄水平、遺傳信息過程和代謝途徑確保其組成基因精確同步操作,以最大限度地減少錯誤,這個過程對細(xì)胞生長和增殖很重要[40]。此外,代謝途徑包含13個不同的亞類,包括全局和概覽圖、碳水化合物代謝、能量代謝和聚糖生物合成和代謝等,并且參與玉米生長發(fā)育[41]。因此,可以強(qiáng)調(diào)遺傳信息處理和代謝過程在玉米生長發(fā)育中的重要性。

近年來,SNP/InDel 標(biāo)記越來越多地用于QTL定位研究,因?yàn)槠湓诨蚪M中的豐度避免了高度重復(fù)的區(qū)域,并且與其他標(biāo)記系統(tǒng)相比可以提供較高的圖譜分辨率[42-43]。 因此,目前轉(zhuǎn)錄組數(shù)據(jù)中已識別的SNP 和InDel標(biāo)記及其在GO和KEGG富集的差異表達(dá)基因?qū)⒂兄诠δ芏鄻有苑治?因?yàn)樗鼈冎苯訁⑴c代謝和細(xì)胞過程的生長和調(diào)節(jié)。此外,這些標(biāo)記將豐富玉米中已有的基因組資源,并可用于飽和現(xiàn)有的連鎖和關(guān)聯(lián)圖譜。