基于BSA-seq和RNA-seq挖掘水稻株高相關QTL

吳元明 林佳怡 柳雨汐 李丹婷 張宗瓊 鄭曉明 逄洪波

（1. 沈陽師范大學生命科學學院，沈陽 110034；2. 廣西壯族自治區農業科學院水稻研究所 廣西水稻遺傳育種重點實驗室，南寧 530007；3. 中國農業科學院作物科學研究所，北京 100081；4. 海南三亞中國農業科學院國家南繁研究院，三亞 571700；5. 國際水稻研究所，菲律賓馬尼拉 DAPO box 7777）

株高是影響水稻（Oryza sativa L.）產量的重要農藝性狀之一，自“綠色革命”的興起［1-2]以及與株高相關代表基因sd1［3]和Rht［4]被發現并闡明作用機制后，形成了對株高與產量之間關系的正確理解：適度的植株矮化增強了其抗倒伏性，使產量更加穩定。普通野生稻（Oryza rufipogon Griff.）是亞洲栽培稻的祖先，廣泛分布于全球各地，遺傳多樣性豐富，具有較強的適應性和抗逆能力，是栽培稻品種改良和培育的重要種質資源庫。染色體片段代換系（chromosome segment substitution lines, CSSLs）群體作為分子育種研究中常用的試驗材料之一，其主要特點是來源于親本的單個染色體段被替換成輪回親本背景，即每一個CSSL除了攜帶1個或少數幾個來自供體親本的代換片段外，其遺傳背景均和輪回親本一致，能有效簡化遺傳背景，是研究QTL的重要群體之一。

隨著現代生物技術的發展，越來越多與株高性狀相關的基因與數量性狀基因座（quantitative trait locus, QTL）被成功定位。如Liu等［5]在水稻第5染色體定位的抗倒伏基因SBI，該基因編碼GA-2氧化酶，通過抑制赤霉素（gibberellins, GAs）的活性控制水稻株高，遺傳分析表明，SBI是控制水稻抗倒伏能力的半顯性基因，其等位基因表現出一系列不同的產物活性。Tu等［6]通過構建RILs（recombinant inbred lines）在一個主效QTL-qPND1內定位到編碼細胞分裂素（cytokinin, CTK）氧化酶的候選基因OsCKX2/Gn1a，該基因第3外顯子中11 bp的缺失使得帶有不同親本型等位基因個體的抗倒伏能力表現出極顯著的差異，后續研究發現抗倒伏型gn1a與SCM2［7]之間的相互作用可進一步增強水稻的抗倒伏能力。最新發現的OsFBK4［8]和DHT1［9]，前者編碼具有F-box和kelch repeat結構的蛋白，通過抑制GAs活性進而促進矮化，后者的產物參與構成剪接體使獨腳金內酯（strigolactone, SL）受體基因D14正常表達，D14促進SL信號抑制物D53的降解，令SL正常發揮作用，進而抑制分蘗形成，其突變形式dht1使D14的mRNA剪接受阻，SL無法發揮作用，進而促進分蘗形成以促進矮化。這一系列研究發現株高的本質是受一系列基因與環境因素共同制約的數量性狀［10-11]。GA［12]和油菜素內酯（brassinolide,BR）［13]是已知的參與株高調控過程的關鍵成員，后續發現SL［14]和生長素（auxin, IAA）［15]也參與這一過程。

雖然已經在水稻中鑒定到一些與株高相關的QTL和相關基因，但相比于株高這個復雜的數量性狀而言，對于這些基因的功能和調控機制的理解還不充分。需要進行更多的研究來挖掘株高基因，以充分理解它們在株高調控中的作用。本研究使用秈稻油占8號和廣西普通野生稻構建的CSSL群體以及現代分子生物學技術，對廣西普通野生稻進行研究，通過二代測序（next generation sequencing, NGS）和混池分析（bulk segregation analysis, BSA）方法，鑒定廣西普通野生稻基因組中與株高相關的QTL，進而篩選出在株高調控中可能起到關鍵作用的候選基因。同時，運用RNA-seq分析進一步驗證這些候選基因在株高調控過程中的表達量與水平，為進一步探究水稻株高調控機制提供了重要的理論指導和研究思路，為培育理想株型的水稻提供理論依據。

1 材料與方法

1.1 材料

以廣西普通野生稻（Oryza rufipogon Griff.）（Guangxi common wild Rice, GCWR）和秈稻油占8號（Oryza sativa L.）（Youzhan8, YZ8）分別作為供體和受體，按照圖1所示流程，采用單粒傳法構建得到285個CSSL群體，由廣西壯族自治區農業科學院水稻研究所提供。

圖1 本研究所用的回交群體的構建過程Fig. 1 Construction process of the backcross populations used in this study

1.2 方法

1.2.1 株高表型鑒定與極端性狀混池的構建 每個CSSL群體隨機選取20粒種子種植于海南省三亞市國家南繁研究院的試驗田中，親本與子代各群體均隨機選取5株用于測量抽穗期株高。分別選取最高株和最矮株的20個系構建高株池（H-pool）和矮株池（L-pool），收集其葉鞘和節間組織，液氮快速冷凍后，提取其基因組DNA和總RNA。

1.2.2 總DNA的提取及NGS測序 采用CTAB法從每個極端群體提取總量為1.5 μg的基因組DNA作為構建文庫的原料，按照生產商提供的標準使用Truseq Nano DNA HT Sample preparation Kit（Illumina USA）構建測序文庫。大致流程為：總DNA由超聲波處理形成長度為350 bp的片段，對其進行末端修復、3'末端添加A尾、連接接頭等處理用于后續的Illumina測序和PCR擴增。PCR產物在AMPure XP system中進行純化，文庫經Agilent2100 Bioanalyzer和real-time PCR完成定量分析。經上述方法構建完成的文庫在Illumina HiSeq 4000 platform中完成測序，最終生成全長為350 bp的雙端測序片段。

1.2.3 測序數據質控及比對 原始測序數據需要經過一系列質量控制程序的處理，去除人為誤差后可用于參考基因組比對和后續分析。質量控制標準包含以下3個方面：（1）去除未知堿基含量≥10%的片段；（2）去除堿基質量值＜5的片段；（3）去除與接頭互補核苷酸＞10 nt的片段，該過程允許的最高錯誤匹配幾率為10%。

將上述質控處理后的樣品測序數據用Burrows-Wheeler Aligner（BWA）［16]軟件對比至Nipponbare（NIP）基因組中（settings: mem -t 4 -k 32 -M -R），參考基因組MSU 7.0下載自http://rice.plantbiology.msu.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/。SAMtools software［17]用于將比對結果轉化為bam文件（settings:-bS -t）。

1.2.4 基于NGS測序的BSA-seq分析與候選基因挖掘 GATK［18]軟件中的Unified Genotyper function用于完成對各樣品SNP和InDel的檢測，并分別使用GATK中不同的參數對SNP（settings: --cluster-WindowSize 4, --filter Expression “QD＜4.0||FS＞60.0||MQ＜40.0”,-G_filter“GQ＜20”）和InDel（settings:--filter Expression “QD＜4.0||FS＞200.0||Read PosRankSum＜-20.0||Inbreeding Coeff＜-0.8”）進行過濾以防止低質量和無效的分子標記對后續分析造成影響。用ANNOVA［19]將過濾后的SNP和InDel在參考基因組文件中進行注釋。

以日本晴基因組MSU 7.0作為參考，比對子代池中與參考基因組在特定堿基位點相同或不同的reads數量，計算不相同reads數量占總數的比例，即為SNP-index。InDel-index的計算方法與SNPindex基本相同，只是InDel涉及長度不等的DNA片段而非單個堿基。若兩個子代池中出現SNP/InDelindex＜0.3的位點，則將其過濾掉，以減少測序和比對誤差造成的影響［20]。選用1 Mb為窗口大小、1 kb為步長作為滑動窗口方法的標準參數，并用該方法對SNP/InDel-index在染色體上的分布進行作圖。計算2個子代池中的SNP/InDel-index的差值，結果即為Δ（SNP/InDel-index）。SNP/InDel-index關聯分析所需的統計學置信區間按照標準方法設置為95%和99%。根據計算的SNP/InDel-index和Δ SNP/InDelindex，挑選出SNP/InDel-index在2個池中差異較大，Δ SNP/InDel-index在一定置信水平下的窗口即為候選區間。在候選區間中挑選出SNP/InDel index在一個池中為1，另一個池中為0的SNP/InDel作為原始的候選SNP，根據這些SNP的注釋信息，挑選在外顯子上，且是非同義突變、stop gain、stop loss等重要突變的SNP，這些SNP作為候選SNP，這些SNP所在的基因作為候選基因。

使用BLAST方法將候選區間內的編碼基因在GO、KEGG、Uniprot、NOG數據庫中進行詳細注釋，以注釋結果為基礎對候選基因進行篩選，初步確定參與水稻株高調控的候選基因。

1.2.5 轉錄組數據分析 根據水稻GH998已有的水稻抽穗期株高轉錄組數據，對高株和矮株的差異表達基因進行分析，所用的轉錄組數據庫從NCBI數據庫中下載（GH998 MU, SRR16686933-SRR16686935和GH998 WT, SRR16686936, SRR16686938,SRR16686939），不同版本的基因ID編號借助The Rice Annotation Project Database［21]中的ID Converter功能完成，利用TBtools［22]軟件從RNA-seq原始數據中提取株高差異表達基因并繪制熱圖。

1.2.6 總RNA的提取和候選基因表達模式的驗證 TRIzol法提取葉鞘和節間組織總RNA，按照如下步驟完成質量檢測：（1）1%瓊脂糖凝膠電泳檢測總RNA是否存在降解和污染；（2）NanoPhotometer? spectrophotometer（IMPLEN, CA,USA）檢測純度；（3）Qubit? RNA Assay Kit in Qubit? 2.0 Flurometer（Life Technologies, CA, USA）測定RNA濃度；（4）RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system（Agilent Technologies,CA, USA）檢測RNA樣品的完整性。

隨機六堿基引物和M-MuLV反轉錄酶（New England Biolabs, Ipswich, MA）用于合成cDNA的第一鏈，DNA聚合酶I和RNase H合成cDNA的第二鏈。Invitrogen試劑盒完成雙鏈cDNA的末端修復、加尾和接頭添加。選用肌動蛋白作為內參基因。所有用于RT-qPCR的引物均采用Primer 3（https://primer3.ut.ee）軟件設計（表1）。各基因設置3次生物學重復以減少誤差。相對表達模式經由2-ΔΔCt方法［23]計算獲得。

表1 用于RT-qPCR分析的引物列表Table 1 List of primers for RT-qPCR analysis

1.2.7 數據分析 將測量得到的CSSL群體株高求平均值分析分布頻率，用Student's t-test檢驗RTqPCR結果的差異顯著性，用GraphPad Prism v8.4.2軟件繪制圖形。

2 結果

2.1 株高性狀評估

通過對285個CSSL群體平均株高表型進行鑒定（圖2）。所有CSSL群體的平均抽穗期株高分布范圍為101.9-128.2 cm，GCWR抽穗期株高為127.4 cm，YZ8抽穗期株高為103 cm。頻率分布表明，株高在CSSL群體中的分布基本符合正態分布的特征，即滿足數量性狀的定義，因此，這些CSSL群體可用于BSA-seq分析。

圖2 CSSL群體的株高表型分析Fig. 2 Plant height phenotype analysis of CSSL populations

2.2 NGS數據質控與比對

親本與極端性狀子代池的基因組DNA經Illumina平臺測序后（表2），供體親本GCWR和受體親本YZ8分別產生11.92 Gb（Q20≥97.75%，Q30≥93.72%，GC含量為44.23%）和11.13 Gb（Q20≥97.57%，Q30≥93.31%，GC含量為44.6%）的原始數據。2個子代池H-pool和L-pool分別產生14.38 Gb（Q20≥97.34%，Q30≥92.92%，GC含量為46.54%）和14.27 Gb（Q20≥97.1%，Q30≥92.45%，GC含量為44.38%）的原始數據。每組樣品的測序數據經整理和過濾后至少有97.89%成功比對至日本晴（Oryza sativa var Nipponbare）基因組MSU 7.0中。各樣本的測序數據量足夠，質量合格，所有樣本的比對結果正常，可用于后續分析。

表2 用于BSA-seq分析的NGS測序數據統計概況Table 2 Summary of NGS data for BSA-seq analysis

2.3 基于NGS的BSA-seq數據分析

SNP-index和InDel-index兩種算法被用于鑒定與株高相關的基因組區域。來自H-pool與L-pool的測序數據分別經BWA、GATK和ANNOVAR比對、分子標記檢測和注釋后，共得到3955 417個SNP和571046個InDels，過濾后得到13362個SNP（其中包含1088個非同義SNP）和1746個InDels。

利用Δ（SNP-index）算法在Chr.7和Chr.10上分別鑒定到3.205（1-3205 000 bp）和1.311 Mb（3782 001-5093 000 bp）的候選基因組區域（圖3）。根據Δ（InDel-index）算法，在Chr.7和Chr.10上分別鑒定到大小為2.848（1-2848 000 bp）和1.292 Mb（3796 001-5088 000 bp）的候選基因組區域（圖4）。值得注意的是，利用Δ（InDel-index）算法得到的候選區間完全包含在Δ（SNP-index）算法鑒定到的候選區間內。上述所有候選基因組區域的置信閾值均大于99%（P＜0.01）。為了縮小定位株高相關的基因組區域，將這兩種算法所得到的候選區間的共有部分作為最終的關聯區段，相關信息如表3所示，分布于第7和第10染色體上的最終區段長度分別為2.848和1.292 Mb，共包含276個已注釋的基因。值得注意的是，在Chr.7的候選區段中包含有已報道的與水稻株高調控相關的基因LOC_Os07g05900（PROG1）和LOC_Os07g05720（OsTCP21）。

表3 Δ（SNP-index）和Δ（InDel-index）定位與株高相關的候選基因組區段Table 3 Candidate genome regions related to plant height mapping by Δ（SNP-index）and Δ（InDel-index）

圖3 基于BSA-seq分析的矮株（A）和高株（B）的SNP-index圖以及Δ（SNP-index）圖（C）Fig. 3 SNP-index graphs of bulk-L（A）and bulk-H（B）, as well as Δ（SNP-index）graph（C）based on BSA-seq analysis

圖4 基于BSA-seq分析的矮株（A）和高株（B）的InDel-index圖以及Δ（InDel-index）圖（C）Fig. 4 InDel-index graphs of bulk-L（A）and bulk-H（B）, as well as Δ（InDel -index）graph（C）based on BSA-seq analysis

2.4 候選區間內基因的功能注釋與富集分析

根據候選基因組區間內SNP、InDel調用和注釋的結果，共有413個高質量的非同義突變SNP（4個nsSNP被定位在日本晴參考基因組的Chr.10內，其余nsSNP位于Chr.7內）包含在92個基因內（3個基因位于參考基因組的Chr.10內，其余基因位于Chr.7內），139個高質量InDels（125個屬于上游類別，其余被歸類于移碼突變）被定位到85個基因內（除了位于Chr.10內的1個基因外，剩余基因全部位于Chr.7內）。將利用兩種分子標記所關聯到的可能的候選基因列表取交集，得到78個共有基因，這些基因即為與株高潛在相關的候選基因。

將78個基因在GO、KEGG、Uniprot、eggNOG數據庫中進行功能注釋以縮小候選基因的選區范圍（表4）。在所有注釋到的高質量nsSNPs和upstream類SNP中分別有15個nsSNP和20個upstream類SNP都位于相同的5個基因內，這些位點的Δ（SNPindex）均等于1。此外，這5個基因還包含有6個Δ（InDel-index）=1的upstream類InDels，其具體的功能注釋信息如表5所示。LOC_Os07g02770在參考基因組中被注釋為含有F-box結構域的蛋白，這是一類廣泛分布于真核生物中的蛋白，具有豐富的多樣性，以參與泛素化蛋白酶體途徑為主要功能。LOC_Os07g02850被注釋為含DUF538結構域的蛋白，該類蛋白為植物所獨有，在雙子葉植物與單子葉植物中均有分布。目前，已鑒定到相當數量的DUF538類蛋白，但其具體的生物學功能仍然未知。LOC_Os07g04220被注釋為與水稻中唯一的紅光和遠紅光感受器——光敏色素信號通路相關的受體類激酶。LOC_Os07g05050被注釋為含有DEAD-box結構域的蛋白，其廣泛分布于動物與植物中，在植物中的已知主要功能是參與對多種脅迫因素的響應過程。LOC_Os07g05720的產物是功能和作用機制已知的OsTCP21轉錄因子蛋白。

表4 候選區間內基因的功能注釋結果統計Table 4 Statistics of functional annotation results for genes in candidate intervals

表5 候選基因功能注釋結果Table 5 Results of functional annotation of candidate genes

2.5 關聯區間內候選基因表達模式分析與驗證

對已有的株高轉錄組數據（GH998MU_S和GH998WT_S）進行分析，結果（圖5）表明，兩種分子標記篩選到的78個共有基因中，有54個在轉錄組數據中可以找到。其中，編碼F-box蛋白的LOC_Os07g02770（Os07g0118900）和編碼OsTCP21轉錄因子的LOC_Os07g05720（Os07g0152000）在轉錄組數據中呈現差異表達，而剩余3個基因在轉錄組數據中表達量較低且無明顯差別。

為了檢驗候選基因的表達模式，利用兩種極端性狀個體的總RNA進行RT-qPCR試驗。結果（圖6）顯示，5個基因的表達模式在兩種個體中均呈差異性表達，表明它們在株高調控過程中均具有一定的功能。3個基因（LOC_Os07g02770、LOC_Os07g04220和LOC_Os07g05720）在矮株中的表達量上調，說明這些基因促進了矮株的形成。2個基因（LOC_Os07g02850和LOC_Os07g05050）則在高株個體中表達量更高，表明其參與株高的正向調控。此外，在5個用于RT-qPCR驗證試驗的基因中，LOC_Os07g05720的產物為功能已知的TCP家族的轉錄因子，與miRNA互作以負向調控株高。該基因的RT-qPCR驗證結果與已有研究相一致，證明本研究中BSA-seq結果的可靠性。

圖6 RT-qPCR驗證候選基因Fig. 6 RT-qPCR validation of selected candidate genes

3 討論

3.1 InDels與SNPs在BSA分析方法中的應用

與質量性狀相比，數量性狀受到更多基因的共同制約，各基因之間相互作用的復雜程度也遠大于前者。隨著測序通量的不斷提升，下一代測序技術的關注點由目標物種的全基因組序列轉變為與數量性狀潛在關聯的多種分子標記。在過去的研究中，SNP標記已經成功地應用于許多研究中［24-30]，因為它們易于獲取和分離，并且通常具有較高的多態性。然而，隨著研究深入，InDels標記作為另一種遺傳變異類型，也逐漸引起了人們的關注和重視，其在BSA-seq中的應用已經逐步得到了證明，如番茄（Solanum lycopersicum L.）［31]和木豆（Cajanus cajan L.）［32]中通過關聯InDels初步定位到與疾病抗性相關的基因組區域。與SNP標記相比，InDels標記可以提供更多的遺傳信息，因為它們不僅包括堿基的改變，還包括DNA序列的插入和缺失。因此，InDels標記可以擴大連鎖不平衡區域的范圍，增加標記密度，并提高遺傳位點的鑒定準確性。另外，由于InDels位于遺傳連鎖的邊界區域，其類型和數量的變化也較為明顯，有助于精細地調控目標基因，并為研究數量性狀的遺傳基礎提供更全面、詳盡的信息［33]。與此同時，InDels標記也存在一些局限性，比如其檢測和分離的難度較大，特別是對于長度不同的InDels，并且對復雜和大基因組物種的數據集的分析需要對用于標記性狀關聯分析的InDels選擇標準進行一些額外的修改等。

本研究使用了SNP和InDels兩種分子標記共同進行BSA-seq分析。這兩種分子標記分別涉及基因編碼區堿基類型和數量的變化，由于密碼子的簡并性，SNP造成的影響可能會被抵消，而InDels所在DNA區域內堿基數量發生的變化，極易使基因功能發生改變。其次，在已有的SNP標記基礎上進行InDels標記，有利于增加標記密度、合并候選區間并為開發新的InDels標記提供理論基礎，而本研究結果也表明了InDels對于性狀定位和功能標記開發的重要性。

3.2 候選區域中與株高相關的基因

本研究利用SNP和InDel共同關聯到5個與水稻株高相關的候選基因，且全部位于Chr.7內。其中OsTCP21、LOC_Os07g05050和LOC_Os07g02850在矮株中呈下調表達，LOC_Os07g04220和LOC_Os07g02770在高株中呈下調表達。其中OsTCP21隸屬于TCP轉錄因子，該基因已被證實與促進水稻矮化相關［34]，其轉錄產物可被miR319靶標進而導致功能受到抑制。此外，通過深入分析發現，Chr.7的候選基因組區域中含有已報道的基因——LOC_Os07g05900（PROG1），其產物為C2H2型鋅指蛋白轉錄因子，PROG1功能的喪失是栽培稻株型改良、得以直立生長的關鍵［35-37]。

株高作為決定水稻產量的重要因素，一直是遺傳學和育種學研究的焦點。目前，對水稻株高的遺傳調控機制的研究主要集中在幾個方向。首先是關于植物激素生長素調控水稻株高的研究。近年來，越來越多的研究表明生長素識別和轉導信號通路對水稻株高起著重要作用［38-39]，生長素水平的調節對植物細胞伸長和分化等過程有很大影響。其次，還有一些研究致力于探究一些轉錄因子調控水稻株高的機制，比如TCP類［40]和bZIP類［41]轉錄因子等。這類轉錄因子可以與生長素信號通路進行交叉調控，并對細胞伸長、分化等發育方式進行調節。另外，一些關于miRNA在水稻株高調節中的作用研究表明［42]，miRNA以通過調節特定基因表達活性的方式對水稻的株高進行調控。但是，以上研究大多只停留在基因表達水平的探究，對于這些基因如何影響水稻株高的機制還有待進一步探索。

本研究利用高通量測序及生物信息分析方法，篩選出與水稻株高關聯的5個候選基因，這些基因均位于Chr.7內。LOC_Os07g05050編碼隸屬于RNA解旋酶超家族II［43]的DEAD-box蛋白，這類蛋白與真核生物的生長發育密切相關，如擬南芥（Arabidopsis thaliana）中與低溫相關的LOS4［44]，水稻中提升高溫耐性的TOGR1［45]以及葡萄（Vitis vinifera L.）中與干旱耐性相關的VviDEADRH25a［46]。LOC_Os07g02850產物具有植物特有的DUF538結構域，已報道的參與光合作用的水溶性葉綠素結合蛋白CaWSCP［47]和水解葉綠素的羧酸酯酶［48]均具有該結構域。LOC_Os07g04220編碼與光敏色素信號相關的受體類激酶，光敏色素C已被研究與玉米（Zea mays L.）株高調控相關［49]，近期的一項研究中發現，光敏色素C編碼基因PHYC的功能缺失突變體參與水稻株高的調節［50]。LOC_Os07g02770的產物是具有F-box結構域的蛋白，該類蛋白通過介導泛素化蛋白酶體途徑來調節植物的多種生命活動［51]，如擬南芥中的F-box蛋白——SAGL1突變體表現出明顯的株型矮化［52]。本研究結果為進一步探究水稻株高的遺傳調控機制提供了新的思路。本研究發現這些基因都位于Chr.7上，這是一個重要的發現。因為Chr.7上有很多其他具有基因功能的序列，這些序列可能也與水稻株高相關。在未來的研究中，可以將這些序列與發現的基因進行比較，以確定它們之間的關系。此外，還需進一步研究這些基因與生長素、RNA等方面的交叉互作關系，揭示更深層次的水稻株高遺傳調控機制。

4 結論

位于Chr.7中的5個候選基因的差異表達可能與水稻株高性狀密切相關，其中，OsTCP21可能是株高的重要調節基因，其表達受miR319靶向和下調的調控。