日本落葉松SNP和InDel位點(diǎn)及其所在抗生物脅迫基因挖掘

摘要：[目的]挖掘日本落葉松SNP和InDel位點(diǎn)及其所在抗生物脅迫基因，為日本落葉松分子育種提供分子標(biāo)記和候選基因。[方法]使用158份來(lái)自中國(guó)、日本和英國(guó)地區(qū)的日本落葉松轉(zhuǎn)錄組數(shù)據(jù)，首先進(jìn)行分子標(biāo)記位點(diǎn)的鑒定和分類，隨后比較了活動(dòng)期和休眠期基因的表達(dá)水平，最后將帶有可靠分子標(biāo)記的差異表達(dá)基因進(jìn)行功能注釋。[結(jié)果]本研究共鑒定到515 935個(gè)SNP位點(diǎn)和1 056個(gè)InDel位點(diǎn)，它們分布在35 827個(gè)基因上。通過(guò)比較不同地區(qū)的位點(diǎn)數(shù)量，推測(cè)日本地區(qū)的日本落葉松遺傳多樣性較為豐富。至少在50份轉(zhuǎn)錄組中出現(xiàn)的非同義突變SNP位點(diǎn)有6 444個(gè)，InDel位點(diǎn)為10個(gè)，它們分布在3 742個(gè)基因上，可以作為可靠的分子標(biāo)記進(jìn)行利用?；顒?dòng)期和休眠期的轉(zhuǎn)錄組比較后，發(fā)現(xiàn)帶有可靠分子標(biāo)記的2 569個(gè)基因差異表達(dá)；GO注釋后發(fā)現(xiàn)其中101個(gè)基因與植物對(duì)真菌、細(xì)菌、卵菌、病毒、昆蟲和線蟲的抗性反應(yīng)有關(guān)。[結(jié)論]這些結(jié)果不僅為日本落葉松全基因組關(guān)聯(lián)分析和全基因組選擇育種提供了分子標(biāo)記，也為利用轉(zhuǎn)基因和基因編輯手段進(jìn)行遺傳育種提供了候選基因。

關(guān)鍵詞：日本落葉松；分子標(biāo)記；抗性基因；轉(zhuǎn)錄組

中圖分類號(hào)：S722;S791.223 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1001-1498（2024）05-0054-11

單核苷酸多態(tài)性標(biāo)記（single nucleotidepolymorphism，SNP）是由單堿基的轉(zhuǎn)換及顛倒等突變而產(chǎn)生的不同個(gè)體之間堿基變異，是物種可遺傳變異中最常見的一種，目前已被應(yīng)用于植物抗性育種研究中，是全基因組關(guān)聯(lián)分析的基礎(chǔ)。Calic等人在山毛櫸（Fagus grandifolia Ehrh.）中檢測(cè)出4個(gè)與抗病性狀極顯著相關(guān)的SNP。插入缺失性標(biāo)記（insertion and deletion，InDel）是指同一物種不同個(gè)體之間基因組同一位點(diǎn)的序列發(fā)生不同大小核苷酸片段的插入或缺失，目前在植物抗性育種研究中報(bào)道較少。

日本落葉松（Larix kaempferi （Lamb.） Carr.）原產(chǎn)日本，被引種到世界各地作為重要造林樹種。目前，關(guān)于日本落葉松林染病的案例越來(lái)越多，例如因松材線蟲（Bursaphelenchus xylophilus）引起的松材線蟲病，因疫霉（Phytophthoraramorum）引發(fā)的疫霉病。為了適應(yīng)生存環(huán)境，在一個(gè)生長(zhǎng)周期中日本落葉松需要經(jīng)歷兩個(gè)時(shí)期：活動(dòng)期和休眠期。在活動(dòng)期，落葉松生長(zhǎng)旺盛，能夠?yàn)椴≡锾峁┴S富的資源和營(yíng)養(yǎng)，且生長(zhǎng)環(huán)境也更適宜病原物生長(zhǎng)和繁殖；而在休眠期，落4p/YAJ1nS+ij3rSQpcwiEw==葉松停止生長(zhǎng)，為病原物提供的資源和營(yíng)養(yǎng)有限，再加上環(huán)境的影響，病原物的活動(dòng)和繁殖受到限制。在這兩個(gè)時(shí)期中有很多基因差異表達(dá)，構(gòu)成了日本落葉松適應(yīng)生存環(huán)境的分子基礎(chǔ)，因此“在活動(dòng)期和休眠期差異表達(dá)的部分基因是否在日本落葉松防御病原物攻擊的過(guò)程中發(fā)揮重要作用”以及“這些基因中存在的序列變異是否可能造成不同基因型擁有不同抗性”值得探討。

日本落葉松基因組序列[9-10]以及大量轉(zhuǎn)錄組的測(cè)定，為鑒定序列變異提供了豐富的序列資源，也為推動(dòng)分子標(biāo)記輔助育種技術(shù)在日本落葉松抗性育種的應(yīng)用奠定基礎(chǔ)。本研究利用日本落葉松158份轉(zhuǎn)錄組數(shù)據(jù)對(duì)SNP和InDel位點(diǎn)進(jìn)行鑒定和分析，并比較了這些位點(diǎn)所在基因在活動(dòng)期和休眠期的表達(dá)，以期為日本落葉松抗性育種提供分子標(biāo)記和候選基因。

1 材料和方法

本研究的全流程如圖1所示。

1.1 轉(zhuǎn)錄組數(shù)據(jù)收集

78份轉(zhuǎn)錄組測(cè)序材料采自中國(guó)林業(yè)科學(xué)研究院，液氮冷凍后-80℃超低溫冰箱保存，干冰運(yùn)至北京諾禾致源科技股份有限公司進(jìn)行轉(zhuǎn)錄組測(cè)序，轉(zhuǎn)錄組數(shù)據(jù)已提交到國(guó)家生物信息中心（https：//www.cncb.ac.cn）數(shù)據(jù)庫(kù)（登錄號(hào)：PRJCA022796）。并且在公共數(shù)據(jù)庫(kù)NCBI（www.ncbi.nlm.nih.gov）上，篩選和下載測(cè)序平臺(tái)盡量一致、Strategy為RNA-Seq、Layout為PAIRED的日本落葉松轉(zhuǎn)錄組數(shù)據(jù)（表1），共同用于鑒定分子標(biāo)記。參考轉(zhuǎn)錄本來(lái)自日本落葉松的50 690條編碼序列（Codingsequence， CDS）（http：//btg.kazusa.or.jp/blast.html）。

1.2 轉(zhuǎn)錄組數(shù)據(jù)處理

利用fastp軟件（v0.32）進(jìn)行數(shù)據(jù)過(guò)濾，去掉測(cè)序接頭，保留長(zhǎng)度在20 bp以上的序列；將過(guò)濾后的數(shù)據(jù)按照長(zhǎng)度進(jìn)行篩選，去掉長(zhǎng)度小于30bp，或者只有一端的read。將質(zhì)控后得到的cleanreads，使用BWA軟件比對(duì)到參考轉(zhuǎn)錄本上，獲取reads在參考轉(zhuǎn)錄本上的定位信息。使用Salmon軟件統(tǒng)計(jì)比對(duì)到每個(gè)轉(zhuǎn)錄本上的read數(shù)量，并計(jì)算每個(gè)轉(zhuǎn)錄本的表達(dá)量，使用Transcripts perkilobase million（TPM）表示。

1.3 基因組變異檢測(cè)和分類

使用變異檢測(cè)軟件bcftools，從1 58份轉(zhuǎn)錄組數(shù)據(jù)中分別檢測(cè)SNP和InDel。隨后對(duì)原始結(jié)果進(jìn)行過(guò)濾。其中mpileup設(shè)定參數(shù)：-C50-d 80 000 -L 10 000-q 20；call設(shè)定參數(shù)-vMm；filter設(shè)定參數(shù)-i＇%QUAL＞10＆（DP4[2]+DP4[3]）＞2‘。最后使用SnpE仟工具對(duì)檢測(cè)到的變異進(jìn)行功能分類。使用TBtools軟件，對(duì)來(lái)自不同地區(qū)（中國(guó)黑龍江，中國(guó)遼寧，日本和英國(guó)）轉(zhuǎn)錄組中的SNP和InDel進(jìn)行韋恩分析。

1.4 差異基因分析

158份轉(zhuǎn)錄組數(shù)據(jù)包含本課題組之前測(cè)定的12份來(lái)自活動(dòng)期和休眠期的日本落葉松轉(zhuǎn)錄組數(shù)據(jù)。本研究使用TBtools軟件對(duì)這12份轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行了差異基因分析（p＜0.05），隨后利用eggNOG和SwissProt數(shù)據(jù)庫(kù)對(duì)得到的差異基因進(jìn)行注釋。

2 結(jié)果

2.1 轉(zhuǎn)錄組數(shù)據(jù)統(tǒng)計(jì)

158份轉(zhuǎn)錄組數(shù)據(jù)統(tǒng)計(jì)顯示，每份數(shù)據(jù)大小在17 292 864～159 873 888 bp之間；質(zhì)控后得到的數(shù)據(jù)大小在16 053 498～155 706 686 bp之間，Q30在80.76%～97.3g%之間，高質(zhì)量數(shù)據(jù)比例在80.89%～99.74%之間。將質(zhì)控后得到clean reads Lk對(duì)到參考轉(zhuǎn)錄本上，發(fā)現(xiàn)共有48 094個(gè)轉(zhuǎn)錄本得到比對(duì)，占總轉(zhuǎn)錄本數(shù)量的94.88%，每份轉(zhuǎn)錄組中有55.72%～88.00%的clean readsLC對(duì)到參考轉(zhuǎn)錄本上。

2.2 SNP、InDel鑒定及分類

每份轉(zhuǎn)錄組中SNP的數(shù)量在50 520～154 343個(gè)之間，InDel的數(shù)量在3 045～9 457個(gè)之間。158份轉(zhuǎn)錄組共有SNP位點(diǎn)515 935個(gè)（表2），共有InDel位點(diǎn)1 056個(gè)（表3），這些位點(diǎn)分布在35 827個(gè)轉(zhuǎn)錄本上，占總轉(zhuǎn)錄本數(shù)量的70.68%。在SNP中，注釋為錯(cuò)義突變的有302 087個(gè)，占比最大，為58.55%，注釋為同義突變的有203 529個(gè)，占比39.45%；注釋為終止密碼子獲得的有9 091個(gè)，注釋為起始密碼子丟失的有318個(gè)；同時(shí)注釋為兩種類型的有861個(gè)。在InDel中，注釋為移碼突變的最多，有742個(gè)，占比70.27%；同時(shí)注釋為兩種類型的有34個(gè)，且其中一種類型都為終止密碼子獲得。

分別從日本、英國(guó)、中國(guó)遼寧和中國(guó)黑龍江來(lái)源的轉(zhuǎn)錄組中隨機(jī)選擇6份進(jìn)行SNP和InDel位點(diǎn)分析。其中，日本獨(dú)有的位點(diǎn)有52 798個(gè)，英國(guó)獨(dú)有的有25 785個(gè)，中國(guó)黑龍江獨(dú)有的有46 1 70個(gè)，中國(guó)遼寧獨(dú)有的有24 056個(gè)，這4個(gè)地區(qū)共有的有12 879個(gè)（圖2）。

2.3 位點(diǎn)篩選

CDS內(nèi)的非同義突變可能會(huì)導(dǎo)致基因功能的改變。本研究篩選出注釋為非同義突變的SNP位點(diǎn)有312 399個(gè)，InDel位點(diǎn)全部為非同義突變，有1 056個(gè)，這些位點(diǎn)所對(duì)應(yīng)的基因共有33 101個(gè)。在158份轉(zhuǎn)錄組中都出現(xiàn)的位點(diǎn)極少，SNP位點(diǎn)僅有755個(gè)（0.24%），InDel位點(diǎn)僅有1個(gè)（0.09%）。在158份轉(zhuǎn)錄組中，至少在25份里出現(xiàn)的SNP和InDel分別為9 642個(gè)（3.09%）和16個(gè)（1.52%）；至少在50份里出現(xiàn)的SNP和InDel分別為6 444個(gè)（2.06%）和10個(gè)（0.95%）；至少在75份里出現(xiàn)的SNP和InDel分別為5 281個(gè)（1.69%）和6個(gè)（0.57%）（圖3）。

基于此，本研究選擇至少在50份里出現(xiàn)的SNP和InDel可作為可靠分子標(biāo)記，它們所對(duì)應(yīng)的基因共有3 742個(gè)。

2.4 抗性相關(guān)位點(diǎn)與基因分析

本研究利用日本落葉松活動(dòng)期和休眠期的轉(zhuǎn)錄組進(jìn)行差異基因分析，鑒定到1 7 828個(gè)差異表達(dá)基因（p＜0.05）（圖1）。將這些差異表達(dá)基因與帶有可靠分子標(biāo)記的3 742個(gè)基因進(jìn)行韋恩分析（圖1），發(fā)現(xiàn)有2 569個(gè)基因既在活動(dòng)期和休眠期差異表達(dá)，又?jǐn)y帶可靠分子標(biāo)記（圖1、4）。

對(duì)這2 569個(gè)基因進(jìn)行GO注釋：注釋到“對(duì)真菌的反應(yīng)（G0：0009620）”，“對(duì)真菌的防御反應(yīng)（GO：0050832）”，“對(duì)共生真菌的反應(yīng)（GO：0009610）”和“對(duì)真菌防御反應(yīng)的調(diào)控（GO：1900150）”的基因有26個(gè)，其中在活動(dòng)期上調(diào)的有17個(gè)（表4）；注釋到“對(duì)細(xì)菌的反應(yīng)（GO：0009617）”，“對(duì)細(xì)菌的防御反應(yīng)（GO：0042742）”和“對(duì)細(xì)菌防御反應(yīng)的調(diào)控（GO：1900424）”的基因有50個(gè)，其中在活動(dòng)期上調(diào)的有31個(gè)（表5）；注釋到“對(duì)卵菌的反應(yīng)（GO：0002239）”，“對(duì)卵菌的防御反應(yīng)（GO：0002229）”和“對(duì)卵菌防御反應(yīng)的調(diào)控（GO：1902288）”的基因有5個(gè)（表6）；注釋到“對(duì)線蟲的反應(yīng)（GO：0009624）”，“對(duì)線蟲的防御反應(yīng)（GO：0002215）”，“對(duì)昆蟲的反應(yīng)（GO：0009625）”和“對(duì)昆蟲的防御反應(yīng)（GO：0002213）”的基因有8個(gè)（表7）；注釋到“對(duì)病毒的反應(yīng)（GO：0009615）”，“細(xì)胞對(duì)病毒的反應(yīng)（GO：0098586）”，“病毒的傳播（GO：0046794）”，“對(duì)病毒的防御反應(yīng)（GO：0051607）”，“病毒誘導(dǎo)的基因沉默（GO：0009616）”和“病毒在多細(xì)胞宿主中的轉(zhuǎn)運(yùn)（GO：0046739）”的基因有11個(gè)（表8）。

3 討論

本研究在CDS區(qū)共鑒定到515 935個(gè)SNP位點(diǎn)。注釋后發(fā)現(xiàn)，39.45%的SNP位點(diǎn)為同義突變，即這些突變不引起氨基酸改變，這與芒果（Mangifera odorata （Huani））和木瓜（Caricapapaya L．）中同義突變SNP的分布頻率一致。同義突變SNP本身不會(huì)改變蛋白質(zhì)序列，因此編碼蛋白的穩(wěn)定性仍可以維持，這可能構(gòu)成了日本落葉松、芒果和木瓜等植物適應(yīng)復(fù)雜生存環(huán)境的遺傳基礎(chǔ)。

本研究共鑒定到1 056個(gè)InDel位點(diǎn)。與SNP相比，InDel發(fā)生在CDS區(qū)更可能會(huì)導(dǎo)致蛋白質(zhì)結(jié)構(gòu)和功能的改變，從而影響生物體的性狀。所以為了維護(hù)蛋白質(zhì)序列和結(jié)構(gòu)的穩(wěn)定，InDel發(fā)生在CDS區(qū)上的數(shù)量比SNP少，這與在哺乳動(dòng)物和果蠅（Drosophila melanogaster）中的研究結(jié)果一致。

分子標(biāo)記位點(diǎn)可以揭示林木種群間及種群內(nèi)的遺傳多樣性。李培等人利用SRAP標(biāo)記對(duì)來(lái)自中國(guó)的29個(gè)紅椿（Toona ciliata Roem.）種源及1個(gè)澳大利亞種源進(jìn)行遺傳多樣性分析，發(fā)現(xiàn)在總的遺傳變異中，種源間分化占79.26%，種源內(nèi)分化僅占20.74%。使用來(lái)自不同地區(qū)的日本落葉松的轉(zhuǎn)錄組數(shù)據(jù)，本研究發(fā)現(xiàn)日本地區(qū)的分子標(biāo)記位點(diǎn)最多（78 219個(gè)，41 .42%），說(shuō)明日本種群的遺傳多樣性較為豐富。日本是日本落葉松的原生地，隨后被引種到不同地區(qū)，包括中國(guó)、英國(guó)。引種后，和原生地存在時(shí)空上的隔離，使得不同地區(qū)日本落葉松之間的基因交換頻率降低。因此，這四個(gè)地區(qū)共有的位點(diǎn)相對(duì)較少，僅有12 879個(gè)（6.82%）（圖2）。

本研究通過(guò)比較活動(dòng)期與休眠期轉(zhuǎn)錄組，發(fā)現(xiàn)在3 742個(gè)可靠位點(diǎn)所在基因中有2 569個(gè)差異表達(dá)，其中有101個(gè)與真菌、細(xì)菌、卵菌、病毒、昆蟲和線蟲抗性反應(yīng)相關(guān)。這些基因不僅構(gòu)成了日本落葉松適應(yīng)生存環(huán)境的分子基礎(chǔ)，可能也在日本落葉松防御病原物攻擊的過(guò)程中發(fā)揮重要作用；另外，這些基因中存在的SNP和InDel可以作為分子標(biāo)記在日本落葉松抗性育種中進(jìn)一步研究和利用。

在真菌感染植物細(xì)胞時(shí)，植物細(xì)胞會(huì)通過(guò)激活特定的防御機(jī)制。本研究共鑒定到17個(gè)與真菌相關(guān)且在活動(dòng)期表達(dá)量高的基因（表4），其中PUX2存在著錯(cuò)義突變SNP。擬南芥（Arabidopsisthaliana（L．）Heynh．）PUX2的一種突變會(huì)降低白粉菌（Golovinomyces orontii）在擬南芥上的繁殖。因此，日本落葉松PUX2的SNP突變也可能影響日本落葉松對(duì)真菌的抗性。

本研究共鑒定到31個(gè)與細(xì)菌相關(guān)且在日本落葉松活動(dòng)期表達(dá)量高的基因（表5）。有研究發(fā)現(xiàn)，EIN3和EIL1突變后擬南芥增強(qiáng)了對(duì)丁香假單胞菌（Pseudomonas syringae）的抵抗力。具有類似情況的基因還有NFXL1，PUB13等。這些結(jié)果表明在活動(dòng)期表達(dá)量高的基因在日本落葉松對(duì)細(xì)菌的抗性反應(yīng)中可能發(fā)揮負(fù)調(diào)控作用。針對(duì)這些負(fù)調(diào)控抗性的基因，可以通過(guò)人工合成miRNA的方式干擾它們的表達(dá)以增強(qiáng)日本落葉松的抗性，這種方法在擬南芥、煙草（Nicotianatabacum L．）及番茄（Solanum lycopersicum L．）等植物上的抗病毒研究中已有應(yīng)用。

本研究鑒定到8個(gè)與昆蟲和線蟲相關(guān)且在活動(dòng)期和休眠期差異表達(dá)的基因（表7），其中包括和線蟲相關(guān)的兩個(gè)基因ZIF/和MTPc2。已有研究發(fā)現(xiàn)，接種根結(jié)線蟲（Meloidogyne spp.）后擬南芥ZIF/的表達(dá)量降低，而MTPc2的表達(dá)量在接種后第一周表達(dá)先升高，隨后降低。這些數(shù)據(jù)表明，ZIF1和MTPc2在日本落葉松對(duì)線蟲的防御中發(fā)揮一定的作用。

4 結(jié)論

通過(guò)分析來(lái)自不同地區(qū)的158份日本落葉松轉(zhuǎn)錄組，共獲取了515 935個(gè)SNP位點(diǎn)和1 056個(gè)InDel位點(diǎn)，它們分布在35 827個(gè)基因上。根據(jù)不同地區(qū)轉(zhuǎn)錄組中的位點(diǎn)數(shù)量，推測(cè)日本地區(qū)的日本落葉松遺傳多樣性較為豐富。至少在50份樣品中出現(xiàn)的非同義突變SNP位點(diǎn)有6 444個(gè)，InDel位點(diǎn)為10個(gè)，它們分布在3 742個(gè)基因上，可以作為可靠的分子標(biāo)記進(jìn)行利用?；顒?dòng)期和休眠期的轉(zhuǎn)錄組比較后，發(fā)現(xiàn)帶有可靠分子標(biāo)記的2 569基因差異表達(dá)；GO注釋后發(fā)現(xiàn)其中101個(gè)基因與植物對(duì)真菌、細(xì)菌、卵菌、病毒、昆蟲和線蟲的抗性反應(yīng)有關(guān)。這些結(jié)果不僅為日本落葉松全基因組關(guān)聯(lián)分析和全基因組選擇育種提供了分子標(biāo)記，也為利用轉(zhuǎn)基因和基因編輯手段進(jìn)行遺傳育種提供了候選基因。

（責(zé)任編輯：張研）

基金項(xiàng)目：科技創(chuàng)新2030 -重大項(xiàng)目（2022ZD0401602; 2022ZD0401705）