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

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

(1.貴州大學農學院, 貴陽 550025; 2.云南省德宏州農業技術推廣中心, 云南 芒市 678400)

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

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

1 材料與方法

1.1 試驗材料

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

1.2 試驗方法

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

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

1.3 玉米轉錄組測序的SNP 和InDel分析

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

2 結果與分析

2.1 樣品轉錄組測序質量分析

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

表1 樣品轉錄組測序質量統計Table 1 Sample transcriptome sequencing quality statistics

表2 轉錄本和Unigene的組裝結果Table 2 Assembly results of transcripts and Unigene

圖1 Unigene長度統計Fig.1 Unigene size statistics

2.2 SNP特征分析

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

圖2 基于轉錄組測序的玉米不同SNP基因型統計Fig.2 Statistics of different SNP genotypes of maize based on transcriptome sequencing

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

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

2.3 InDel分析

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

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

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

圖4 InDel 類型統計Fig.4 InDel type statistics

2.4 SNP/InDel的區域分布

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

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

2.5 差異表達基因GO和KEGG 功能分類分析

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

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

3 結論與討論

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

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

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

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

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

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

近年來,SNP/InDel 標記越來越多地用于QTL定位研究,因為其在基因組中的豐度避免了高度重復的區域,并且與其他標記系統相比可以提供較高的圖譜分辨率[42-43]。 因此,目前轉錄組數據中已識別的SNP 和InDel標記及其在GO和KEGG富集的差異表達基因將有助于功能多樣性分析,因為它們直接參與代謝和細胞過程的生長和調節。此外,這些標記將豐富玉米中已有的基因組資源,并可用于飽和現有的連鎖和關聯圖譜。