丁宇　王馬寅　唐敏強　李子陽　謝尚潛

摘要：棉花（Gossypium hirsutum L.）是一種重要的經(jīng)濟作物，是世界上第二大的天然紡織纖維來源和重要的食用油來源。棉花對生物和非生物脅迫高度敏感，尤其是高溫脅迫極易影響花粉的活力和花藥的開裂。為了解析棉花在高溫脅迫下基因轉(zhuǎn)錄水平的相應(yīng)變化機制，開展了熱敏和耐熱2個棉花品種的全長轉(zhuǎn)錄組響應(yīng)高溫脅迫處理變化的研究。通過轉(zhuǎn)錄本的可變剪切分析發(fā)現(xiàn)，2個棉花品種在高溫脅迫下可變剪切的總數(shù)顯著增加。熱敏品種Che61-72中發(fā)現(xiàn)了2 900個差異表達基因，并且差異基因在加熱前樣本（R0）和加熱12 h后的樣本（R12）中分別識別到了 13 251 個和25 296個可變剪切事件，其中內(nèi)含子保留事件增加得最多，有3 837個。耐熱品種新陸早36號中發(fā)現(xiàn)了 2 437 個差異表達基因，在加熱前樣本（T0）和加熱12 h后的樣本（T12）中鑒定到了11 248個和13 769個可變剪切事件，外顯子跳躍事件變化得最大，增加了4 144個。富集分析發(fā)現(xiàn)，2個品種的差異基因都顯著富集到了光系統(tǒng)Ⅰ的光捕獲、葉綠體類囊體膜和光合作用-天線蛋白通路中，并篩選出5個關(guān)鍵基因（CPB3、A0A1U8IZF2、A0A1U8KCA2、A0A1U8NDW4和A0A1U8NI70），均被注釋為葉綠素a/b結(jié)合蛋白，它們參與了調(diào)控棉花光合作用動態(tài)平衡。本研究為棉花在高溫脅迫的調(diào)節(jié)機制的深入研究提供了理論依據(jù)，為后續(xù)耐高溫的種質(zhì)改良及新品種培育提供了數(shù)據(jù)支持。

關(guān)鍵詞：棉花；高溫脅迫；全長轉(zhuǎn)錄組；納米孔測序；可變剪切

中圖分類號：S562.01文獻標志碼： A

文章編號：1002-1302（2023）05-0001-10

棉花（Gossypium hirsutum L.）是一種重要的經(jīng)濟作物，是世界上第二大的天然紡織纖維來源和重要的食用油來源［1］。近年來，由于溫室效應(yīng)，地表溫度升高對植物生長產(chǎn)生不同程度的影響［2］。棉花在高溫脅迫下花粉活力會受嚴重影響并導(dǎo)致花藥的開裂［3-4］。因此，從不同水平了解棉花抗高溫脅迫分子機理，選育耐高溫棉花材料是十分必要的。

隨著測序技術(shù)的發(fā)展以及基因組研究的興起，棉花抗高溫脅迫已經(jīng)開展了相關(guān)的研究。已有研究表明，棉花的酶活性、萌發(fā)、生長、根系發(fā)育、營養(yǎng)發(fā)育和開花最適溫度在25～31 ℃［5］，高溫環(huán)境通過減弱棉花光合作用和新陳代謝，引起授粉能力減弱［6］，影響棉花花朵的定型、生產(chǎn)率和所產(chǎn)生的纖維的質(zhì)量［7］。Liang等通過納米孔測序轉(zhuǎn)錄組分析了棉花的抗熱脅迫機制，揭示了富集在核糖體、激素信號轉(zhuǎn)導(dǎo)以及內(nèi)質(zhì)網(wǎng)蛋白質(zhì)加工途徑的差異表達基因、葉綠體a/b結(jié)合蛋白基因和lncRNAs對棉花耐熱機制的調(diào)節(jié)作用［8］。Reddy等研究發(fā)現(xiàn)，在20～30 ℃的最適溫度下，棉花生物產(chǎn)量增加，超過最適溫度則導(dǎo)致產(chǎn)量嚴重下降［9］。高于最適溫度會抑制棉花光合作用［7-10］，同時伴隨著碳水化合物的進一步消耗和呼吸作用的增強［11］。Barua等發(fā)現(xiàn)，熱休克蛋白與棉花等作物的光合作用能力密切相關(guān)［12］。Min等發(fā)現(xiàn)，在高溫敏感品系種，棉花花藥在高溫誘導(dǎo)下GhCKI高度表達，該基因可能在高溫脅迫中有重要作用［13］。Wang等研究發(fā)現(xiàn)，miRNAs通過調(diào)節(jié)與生長穩(wěn)態(tài)、活性氧、葉綠體功能、植物-病原菌相互作用和植物激素信號轉(zhuǎn)導(dǎo)途徑相關(guān)的基因來發(fā)揮作用［14］。Ma等揭示了lncRNA對高溫下棉花雄性生殖的作用，并發(fā)現(xiàn)了GhHRK1基因在高溫脅迫反應(yīng)中的負調(diào)控作用［15］。這些研究結(jié)果為理解棉花抗高溫脅迫的基因水平的調(diào)控機制提供了支持，但高溫脅迫條件下關(guān)鍵調(diào)控基因轉(zhuǎn)錄本水平的調(diào)控解析仍有待進一步研究。

目前，三代測序技術(shù)已經(jīng)被廣泛運用于植物的全長轉(zhuǎn)錄本研究，主要包括鑒定復(fù)雜的可變剪切事件、全長轉(zhuǎn)錄變異體、融合轉(zhuǎn)錄本和選擇性多聚腺苷事件。可變剪切事件的發(fā)生是調(diào)控基因表達的一種重要調(diào)控機制，也是促進轉(zhuǎn)錄組和蛋白質(zhì)組多樣性的動力［16］。三代測序相較于二代測序在序列長度、轉(zhuǎn)錄本識別方面的優(yōu)勢，更有助于識別復(fù)雜的剪接異構(gòu)體，研究植物響應(yīng)各種脅迫的基因調(diào)控機制［17］。前人利用三代測序技術(shù)已經(jīng)在高粱（Sorghum bicolor）［18］、玉米（Zea mays）［19］、異源多倍體棉花（allopolyploid cotton）［20］和花竹（Phyllostachys edulis）［21］中分別鑒定到可觀的可變剪切事件。多倍體棉花中與耐熱機制相關(guān)的可變剪切事件仍有待進一步研究。

本研究基于納米孔測序的全長轉(zhuǎn)錄組數(shù)據(jù)，分析了熱敏（Che61-72）和耐熱（新陸早36號）2個品種的棉花在熱刺激處理前后的可變剪切差異，發(fā)現(xiàn)高溫脅迫增加了可變剪切事件，并進一步分析了在高溫脅迫中差異表達基因的功能，以及可變剪切事件的解析。本研究為棉花高溫脅迫的調(diào)節(jié)機制的深入研究提供了理論依據(jù)，為后續(xù)耐高溫的種質(zhì)改良及新品種培育提供了數(shù)據(jù)支持。

1 材料與方法

1.1 數(shù)據(jù)收集與處理

收集來自NCBI的2個棉花品種（熱敏品種Che61-72和耐熱品種新陸早36號）的納米孔測序數(shù)據(jù)，此數(shù)據(jù)由石河子大學在2021年3月4日上傳發(fā)布，登錄號為PRJNA706603。測序樣品使用光照培養(yǎng)箱（溫度為40 ℃）進行熱應(yīng)激處理，在處理0 h和12 h后測序。使用pipeline-nanopore-ref-isoforms程序（https：//github.com/nanoporetech/ pipeline- nanopore-ref-isoforms）對納米孔測序全長轉(zhuǎn)錄組數(shù)據(jù)進行處理，再用minimap2 v2.24軟件將測序數(shù)據(jù)比對到棉花的參考基因組（https：//www.ncbi.nlm.nih.gov/genome/？term=Gossypium+ hirsutum）［22］，參數(shù)設(shè)置為-ax splice，-secondary=no-C5。

1.2 注釋與可變事件的鑒定

將cDNA_Cupcake（https：//github.com/Magdoll/cDNA_Cupcake）的collapse_isoforms_ by_sam.py腳本用于比對數(shù)據(jù)去除冗余的轉(zhuǎn)錄本和合并同種轉(zhuǎn)錄本，使用SQANTI3 v4.3軟件［23］對處理后的轉(zhuǎn)錄本數(shù)據(jù)結(jié)果與基因組注釋進行比較分析。可變剪切事件通過SUPPA v2.3軟件［24］使用默認參數(shù)進行識別鑒定，SUPPA可以從注釋生成轉(zhuǎn)錄事件和局部選擇性剪切事件，把可變剪切事件分為以下7種類型：外顯子跳躍（SE），外顯子互斥（MX），5′端選擇性剪切（A5），3′端選擇性剪切（A3），內(nèi)含子保留（RI），替代第1個外顯子（AF），替代最后1個外顯子（AL）。

1.3 基因表達水平的量化和差異表達分析

通過Salmon定量軟件［25］將全長reads比對到參考基因組上，并對匹配質(zhì)量大于5的reads進行量化。通過每百萬計數(shù)（CPM）估計基因表達水平（CPM=映射到轉(zhuǎn)錄本的reads/樣本中對齊的總reads×106）。使用DESeq2 R包［26］對2個棉花品種不同時間點的樣品進行差異表達分析。在該分析中，由DESeq2確定的P＜0.01、|foldchange|≥2被認為是差異表達基因（DEG）。

1.4 差異表達基因的功能富集分析

使用eggNOG v2.0.0軟件［27］進行差異基因注釋，從差異基因的GO數(shù)據(jù)庫中獲得注釋背景信息。利用R的clusterProfiler R包［28］對差異基因進行GO富集分析，P值小于0.05的GO條目被認為是差異表達基因顯著富集。進一步使用KOBAS v3.0［29］分析KEGG通路中DEGs的富集情況（α=0.05）。

2 結(jié)果與分析

2.1 數(shù)據(jù)統(tǒng)計分析

本研究從NCBI公共數(shù)據(jù)庫中選取了棉花的耐熱品種（新陸早36號）和熱敏品種（Che61-72）的納米孔測序的轉(zhuǎn)錄組樣本進行分析，其中熱敏品種Che61-72對照組（加熱前樣本，R0）和試驗組（加熱12 h后樣本，R12）總的reads分別為19 301 677和21 456 617，總堿基數(shù)為15.6 G和19.4 G，GC含量為43.2%和42.5%。 序列的平均長度為849 bp和951 bp，N50長度為918 bp和1 064 bp。在耐熱品種新陸早36號的對照組（加熱前樣本，T0）和試驗組（加熱12 h后樣本，T12）中，分別獲得了 19 219 360 和16 564 945條reads，總的堿基數(shù)為14.7 G和13.6 G，GC含量分別為43.3%和42.6%，序列N50長度為868 bp和934 bp，平均長度為801 bp和964 bp（表1）。

2.2 數(shù)據(jù)的處理與注釋

比對2個不同品種棉花獲得轉(zhuǎn)錄本，使用cDNA_Cupcake對全長轉(zhuǎn)錄組序列進一步去除冗余轉(zhuǎn)錄本。4個樣品的轉(zhuǎn)錄組序列比對到基因組上的

比對率分別為99.93%、99.93%、99.91%和99.93%（表2）。

使用SQANTI3v.4.3對轉(zhuǎn)錄組進行分析和注釋，轉(zhuǎn)錄本總數(shù)最多的是熱敏品種Che61-72的試驗組（R12），得到了377 518個轉(zhuǎn)錄本，基因總數(shù)也最多，其中已知基因數(shù)為65 347個，新基因數(shù)量為 12 291 個。而轉(zhuǎn)錄本最少的是耐熱品種新陸早36號的試驗組（T12），有195 797個。基因總數(shù)為 50 437 個，其中已知基因數(shù)為44 711，新基因數(shù)有 5 726 個。對新的基因進一步進行注釋，并用于下游分析。

在熱敏品種Che61-72的R0中每個基因只有一個轉(zhuǎn)錄本的現(xiàn)象占的比重最大，有32.25%，其次是具有2～3個轉(zhuǎn)錄本，占了28.49%。而在R12中具有6個以上轉(zhuǎn)錄本的基因的比重最高，有36.3%。在耐熱品種新陸早36號的T0和T12中，每個基因具有1個轉(zhuǎn)錄本的現(xiàn)象比重最大（圖1）。經(jīng)過熱應(yīng)激處理后，2個品種中6個以上轉(zhuǎn)錄本的數(shù)量增加，可能是生物體內(nèi)基因需要產(chǎn)生更多的可變剪切體來應(yīng)對熱應(yīng)激反應(yīng)。這些結(jié)果表明，使用納米孔測序平臺進行大規(guī)模轉(zhuǎn)錄本測序是重建轉(zhuǎn)錄本的有效方法，有利于復(fù)雜轉(zhuǎn)錄事件的研究。

2.3 熱應(yīng)激處理后可變剪切事件的廣泛動態(tài)變化

在高等動植物基因組中，大部分基因在轉(zhuǎn)錄過程中存在可變剪切，從而引起同個基因的多個轉(zhuǎn)錄本。植物在脅迫條件下會產(chǎn)生更多的剪切變異來微調(diào)其基因的表達模式［30］。本研究使用了SUPPA軟件去識別可變剪切事件，可變剪切事件分為7種類型：外顯子跳躍（SE），外顯子互斥（MX），5′端選擇性剪切（A5），3′端選擇性剪切（A3），內(nèi)含子保留（RI），替代第1個外顯子（AF），替代最后1個外顯子（AL）（圖2-A）。在熱敏品種Che61-72的R0樣品中，從35 162個基因中檢測到了66 638個可變剪切事件，在R12樣品中，從61 097個基因中檢測到了137 248可變剪切事件（圖2-B）。研究結(jié)果表明，經(jīng)過40 ℃熱應(yīng)激的外部條件刺激12 h后，R12與R0相比表達的基因變多，可變剪切事件增加，特別是內(nèi)含子保留（RI）增加了27 591，占了總的可變剪切數(shù)量的30.06%，其次是3′端選擇性剪切（A3），占了28.46%（圖2-C）。在耐熱品種新陸早36號的T0樣品中，從31 687個基因識別到了 59 639 個可變剪切事件，其中內(nèi)含子保留（RI）事件最多。在T12樣品中，從34 977個基因中識別到了69 483個可變剪切事件，而3′端選擇性剪切事件最多，占了總事件的30.73%（圖2-B和圖2-C）。在40 ℃熱應(yīng)激刺激后， T12與T0相比表達的基因和可變剪切事件都有所增加，但耐熱品種對熱應(yīng)激的反應(yīng)沒有熱敏品種這般強烈，最多的3′端選擇性剪切事件增加了4 023個。3′端選擇性剪切（A3）是棉花中最普遍的形式［31］。已有研究表明，內(nèi)含子保留（RI）可能是植物應(yīng)對脅迫的重要分子機制［32］。棉花在熱應(yīng)激處理后，內(nèi)含子保留（RI）事件發(fā)生了顯著變化，這表明內(nèi)含子保留（RI）可能在響應(yīng)高溫脅迫的轉(zhuǎn)錄后轉(zhuǎn)錄組調(diào)整中發(fā)揮作用。

2.4 響應(yīng)熱應(yīng)激的差異表達基因與可變剪切的變化

為了探究熱應(yīng)激條件下基因表達的差異，分析了熱敏品種Che61-72和耐熱品種新陸早36號在溫度為40 ℃條件刺激12 h后的基因表達情況。在Che61-72中發(fā)現(xiàn)了2 900個差異基因（DEGs），其中1 592個基因表達上調(diào)，1 308個基因表達下調(diào)。差異基因在R0樣本中共識別到了13 251個可變剪切事件，在R12樣本差異基因中，可變剪切事件大量增加，一共識別到了25 296個可變剪切事件，其中RI事件增加得最多，有3 837個（圖3-A）。在新陸早36號中發(fā)現(xiàn)了2 437個差異表達基因，其中1 369個上調(diào)基因，1 068個下調(diào)基因。在T0樣本中鑒定到了11 248個可變剪切事件，在T12樣本中鑒定到了13 769個可變剪切事件，在熱應(yīng)激后，可變剪切事件增加了2 521個，其中外顯子跳躍事件變化得最大，增加了4 144個。

從2個不同品種的棉花樣本的可變剪切變化中研究發(fā)現(xiàn)，熱敏品種Che61-72對熱應(yīng)激反應(yīng)非常強烈，不同類別的可變剪切事件都呈上升趨勢，而耐熱品種新陸早36號的不同類別的可變剪切事件除了外顯子跳躍顯著增加和外顯子互斥少量增加之外，其他5個類型都有減少的現(xiàn)象，這可能是與耐熱品種中耐熱機制對熱應(yīng)激的自身調(diào)控有關(guān)。

為了防止因為材料本身的差異導(dǎo)致的結(jié)果誤差，本研究繼續(xù)分析了熱敏品種（R0）和耐熱品種（T0）之間的基因表達的差異情況。一共發(fā)現(xiàn)了318個差異基因，其中206個基因表達上調(diào)，112個基因表達下調(diào)。其與熱敏品種的熱應(yīng)激的表達差異基因中有37個共有基因（圖3-C），與耐熱品種的熱應(yīng)激的表達差異基因中有38個共有基因（圖3-D）。在后續(xù)的功能注釋分析時，剔除共有差異基因。

2.5 差異表達基因的GO與KEGG富集分析

進一步分析差異表達基因的功能，本研究對Che61-72和新陸早36號的差異基因進行了GO富集分析。在Che61-72的差異基因中，一共富集到了973條GO通路。在生物過程通路中，差異基因主要富集到了光合作用、光系統(tǒng)Ⅰ的光捕獲、對光刺激的反應(yīng)、蛋白質(zhì)折疊等；在細胞成分通路中，主要富集在胞質(zhì)、葉綠體類囊體膜和質(zhì)體球等；在分子功能類別中，主要富集在DNA結(jié)合轉(zhuǎn)錄因子活性、折疊蛋白的結(jié)合和水通道活動等。在新陸早36號的差異基因被856條GO通路注釋，在生物過程中，差異基因主要涉及到了核糖體的結(jié)構(gòu)成分、核糖核酸綁定和葉綠體類囊體膜等條目；在細胞成分中，主要涉及了質(zhì)體球、蛋白質(zhì)泛素化和核斑點中；在分子功能類別中，主要富集在血紅素結(jié)合、細胞質(zhì)大核糖體亞基和細胞質(zhì)小核糖體亞基等條目之中（圖4）。

利用KOBAS v3.0分析了2個樣本差異基因的KEGG富集通路。在Che61-72的差異基因中，一共被122條通路注釋到，分別是光合作用-天線蛋白、晝夜節(jié)律-植物、代謝途徑和次生代謝物的生物合成等通路；新陸早36號的差異基因一共富集到了113條通路之中，其中有光合作用-天線蛋白、核糖體、代謝途徑和次生代謝物的生物合成等通路（圖5）。通過富集分析發(fā)現(xiàn)，在GO條目中2個品種的差異基因都顯著富集在光系統(tǒng)Ⅰ的光捕獲、葉綠體類囊體膜中，在KEGG通路中最顯著富集于光合作用-天線蛋白中，可見熱應(yīng)激對棉花的光合作用過程產(chǎn)生了影響。

2.6 響應(yīng)熱應(yīng)激的信號通路和基因分析

在本研究中，2個品種的差異基因都顯著富集到了光系統(tǒng)Ⅰ的光捕獲、葉綠體類囊體膜和光合作用-天線蛋白中。在光系統(tǒng)I的光捕獲之中，Che61-72 富集上了39個差異基因，其中2個基因上調(diào)，37個基因下調(diào)。使用string蛋白質(zhì)數(shù)據(jù)庫（https：//cn.string-db.org/）對基因進行注釋，發(fā)現(xiàn)有11個基因被注釋成葉綠素a/b結(jié)合蛋白。新陸早36號富集到了35個差異基因，表達下調(diào)，有9個基因被注釋成葉綠素a/b結(jié)合蛋白。在葉綠體類囊體膜之中，Che61-72富集上了52個差異基因，有6個基因表達上調(diào)，46個基因表達下調(diào)。在新陸早36號富集到了48個差異基因，其中11個基因表達上調(diào)，37個基因表達下調(diào)。富集到的基因分別被注釋成葉綠素a/b結(jié)合蛋白、核糖核蛋白A、伴侶蛋白、光系統(tǒng)Ⅰ反應(yīng)中心亞基和蛋白質(zhì)曲率類囊體1A。在KEGG最顯著富集的通路光合作用-天線蛋白中，Che61-72富集到了41個差異基因，2個基因表達上調(diào)，39個基因表達下調(diào)，新陸早36號在通路中富集到了37個差異基因，表達下調(diào)。此通路的基因都屬于捕光葉綠素a/b結(jié)合蛋白（LHC）家族。

在熱敏品種所富集上的基因中，在熱應(yīng)激前（R0）的轉(zhuǎn)錄本總數(shù)為952個，可變剪切的總數(shù)為735個，其中內(nèi)含子保留的類型最多有231個，其次是3′端可變剪切有221個。在熱應(yīng)激12 h后（R12），轉(zhuǎn)錄本數(shù)量為565個，比R0減少了387個，可變剪切總數(shù)有571個，最多的類型是3′端可變剪切有167個，其次為內(nèi)含子保留152個。與R0相比，內(nèi)含子保留的數(shù)量減少的最多，減少了79個。在耐熱品種所富集上的基因里，熱應(yīng)激前（T0）的轉(zhuǎn)錄本總數(shù)為776個，可變剪切總數(shù)為609個，3′端可變剪切和5′端可變剪切的數(shù)量最多，分別為171個和178個。熱應(yīng)激12 h后（T12）轉(zhuǎn)錄本數(shù)量和T0相比，減少了383個，可變剪切的數(shù)量減少了281個，其中3′端可變剪切減少的數(shù)量最多，減少了77個（表3）。

在自然條件下，強烈的陽光會削弱光合作用，導(dǎo)致光合活性受損。本研究中，在熱應(yīng)激12 h的條件下處于這3條通路上的基因幾乎都呈現(xiàn)下調(diào)表達（圖6），導(dǎo)致轉(zhuǎn)錄本數(shù)目和可變剪切數(shù)量減少，其中CPB3、A0A1U8IZF2、A0A1U8KCA2、A0A1U8NDW4和A0A1U8NI70同時參與了這3條通路的調(diào)控。

3 討論與結(jié)論

納米孔測序技術(shù)作為一種長讀長測序技術(shù)，并且可以表征全長轉(zhuǎn)錄本以及量化轉(zhuǎn)錄表達，它能夠在轉(zhuǎn)錄和轉(zhuǎn)錄后水平上剖析植物轉(zhuǎn)錄本的復(fù)雜性［33］。本研究分析了納米孔測序平臺的2個棉花品種的全長轉(zhuǎn)錄本數(shù)據(jù)。分析測序結(jié)果得到了超過63 Gbp的有效數(shù)據(jù)，并且讀取的平均長度超過800個堿基。超過 99% 的reads成功比對到棉花基因組，這表明納米孔測序reads的高質(zhì)量。大規(guī)模、高質(zhì)量的納米孔數(shù)據(jù)保證了轉(zhuǎn)錄本鑒定和定量的準確性，這表明納米孔測序是一種具有成本效益的植物全長轉(zhuǎn)錄本分析方法。

可變剪切會產(chǎn)生多個轉(zhuǎn)錄本，并且可能會從同一基因產(chǎn)生不止一種蛋白質(zhì)，它會因?qū)χ参锷L和發(fā)育產(chǎn)生負面影響的環(huán)境壓力而改變。植物脅迫相關(guān)基因特別容易發(fā)生選擇性剪切事件，這通常會調(diào)節(jié)活性和非活性異構(gòu)體之間的比例以響應(yīng)非生物脅迫，從而微調(diào)關(guān)鍵脅迫調(diào)節(jié)因子的表達［32］。本研究在熱應(yīng)激的脅迫條件下，從全轉(zhuǎn)錄水平展示了棉花的熱敏品種Che61-72和耐熱品種新陸早36號的可變剪切的變化。2個品種的棉花在高溫脅迫下，可變剪切的總數(shù)大量增加，特別是熱敏品種Che61-72，在被熱刺激的12 h之后，可變剪切的數(shù)量增加了2倍多，并且內(nèi)含子保留的事件最為顯著。然后分析了2個棉花品種在熱應(yīng)激12 h前后基因表達模式的變化。在熱敏品種Che61-72中發(fā)現(xiàn)了2 900個差異基因，并且差異基因在R0樣本中識別到了13 251個可變剪切事件，在R12樣本中識別到了25 296個可變剪切事件，其中內(nèi)含子保留事件增加得最多。在新陸早36號中發(fā)現(xiàn)了2 437個差異表達基因，在T0樣本中鑒定到了11 248個可變剪切事件，在T12樣本中鑒定到了13 769個可變剪切事件，與熱敏品種不同的是，耐熱品種的外顯子跳躍事件變化得最大。

高溫會導(dǎo)致光合作用裝置的直接損害，而光系統(tǒng)Ⅱ通常被認為是熱誘導(dǎo)的光合作用失活的主要目標，光系統(tǒng)Ⅰ被認為可以確定生命系統(tǒng)中的全球焓量，因此溫度升高的影響關(guān)于光系統(tǒng)Ⅰ可能對于調(diào)節(jié)變化環(huán)境中所有光合自養(yǎng)生物的光合反應(yīng)至關(guān)重要［34］。葉綠體的宏觀結(jié)構(gòu)會因熱刺激下被暴露而改變，類囊體膜的完整性受到損害，在35～45 ℃ 的溫度范圍內(nèi)可以觀察到膜堆積減少和類囊體膜的一般重組［35］。差異基因GO和KEGG富集分析發(fā)現(xiàn)，2個品種的差異基因都顯著富集到了光系統(tǒng)Ⅰ的光捕獲、葉綠體類囊體膜和光合作用-天線蛋白通路中，并且篩選出5個關(guān)鍵基因（CPB3、A0A1U8IZF2、A0A1U8KCA2、A0A1U8NDW4和A0A1U8NI70），均被注釋為葉綠素a/b結(jié)合蛋白（LHC），屬于光捕獲葉綠素a/b結(jié)合（LHC）蛋白家族，是光捕獲復(fù)合體的組成部分，作為一個光受體，它捕獲并將激發(fā)能量傳遞給與其密切相關(guān)的光系統(tǒng)。Liang等在棉花納米孔測序轉(zhuǎn)錄組分析中也發(fā)現(xiàn)葉綠體a/b結(jié)合蛋白基因（LHC）有利于耐熱品種的光合穩(wěn)態(tài)，參與棉花抗熱調(diào)節(jié)機制［8］。在高溫處理下的大麥（Hordeum vulgare）和葡萄（Vitis vinifera）中， 葉綠體a/b結(jié)合蛋白（LHC）的蛋白豐度和轉(zhuǎn)錄本均下降［36］。紫花苜蓿（Medicago sativa）葉片中的3種葉綠a/b結(jié)合蛋白在24 h的高溫處理下，表達豐度提高，但在延長的高溫脅迫過程中（72 h）蛋白豐度下降［37］，表明葉綠體素結(jié)合蛋白的豐度變化可能與耐熱性有關(guān)，其可能取決于高溫脅迫的強度、持續(xù)時間和植物種類。此外，在低溫處理的芒屬植物中，葉綠素a/b結(jié)合蛋白基因上調(diào)表達，表明可能有助于保護光合作用［38］。本研究中，高溫脅迫是導(dǎo)致棉花葉綠素a/b結(jié)合蛋白基因下調(diào)的根本原因，5個LHC基因與調(diào)控棉花光合作用動態(tài)平衡密切相關(guān)。

本研究分析了2個棉花品種在高溫脅迫下轉(zhuǎn)錄組的可變剪切事件，發(fā)現(xiàn)在熱處理后可變剪切事件的總數(shù)量增多，其中熱敏品種增加了70 610個，耐熱品種增加了9 844個。這表明轉(zhuǎn)錄后剪接在棉花應(yīng)對高溫脅迫的生物調(diào)節(jié)中起重要作用。富集分析發(fā)現(xiàn)，差異基因顯著富集到了光系統(tǒng)Ⅰ的光捕獲、葉綠體類囊體膜和光合作用-天線蛋白通路中，并篩選出5個關(guān)鍵基因，它們均為葉綠素a/b結(jié)合蛋白基因，參與高溫脅迫下調(diào)控棉花光合作用動態(tài)平衡。本研究探索了棉花在高溫脅迫條件下的調(diào)節(jié)機制，為后續(xù)耐高溫的種質(zhì)選育奠定基礎(chǔ)、提供重要了的數(shù)據(jù)支持。

參考文獻：

［1］Han Z，Wang C，Song X，et al. Characteristics，development and mapping of Gossypium hirsutum derived EST-SSRs in allotetraploid cotton ［J］. Theoretical and Applied Genetics，2006，112（3）：430-439.

［2］Grover A，Mittal D，Negi M，et al. Generating high temperature tolerant transgenic plants：achievements and challenges ［J］. Plant Science，2013，205：38-47.

［3］Zahid K R，Ali F，Shah F，et al. Response and tolerance mechanism of cotton Gossypium hirsutum L. to elevated temperature stress：a review ［J］. Frontiers in Plant Science，2016，7：937.

［4］Pettigrew W T. The effect of higher temperatures on cotton lint yield production and fiber quality ［J］. Crop Science，2008，48（1）：278-285.

［5］Burke J J，Wanjura D F. Plant responses to temperature extremes［M］//Physiology of cotton. Berlin：Springer，2010：123-128.

［6］Oosterhuis D M，Snider J L. High temperature stress on floral development and yield of cotton ［M］//Oosterhuis D M. Stress physiology in cotton.? Cordova，Tennessee：The Cotton Foundation，2011：1-24.

［7］Reddy K R，Kakanl V G，Zhao D，et al. Interactive effects of ultraviolet-B radiation and temperature on cotton physiology，growth，development and hyperspectral reflectance ［J］. Photochemistry and Photobiology，2004，79（5）：416-427.

［8］Liang Y，Gong Z，Wang J，et al. Nanopore-based comparative transcriptome analysis reveals the potential mechanism of high-temperature tolerance in cotton （Gossypium hirsutum L.） ［J］. Plants，2021，10（11）：2517.

［9］Reddy K，Hodges H，Reddy V. Temperature effects on cotton fruit retention ［J］. Agronomy Journal，1992，84（1）：26-30.

［10］Reddy V，Reddy K，and Baker D. Temperature effect on growth and development of cotton during the fruiting period ［J］. Agronomy Journal，1991，83（1）：211-217.

［11］Soliz L M A，Oosterhuis D M，Coker D L，et al. Physiological response of cotton to high night temperature ［J］. Am J Plant Sci Biotechnol，2008，2：63-68.

［12］Barua D，Downs C A，Heckathorn S A. Variation in chloroplast small heat-shock protein function is a major determinant of variation in thermotolerance of photosynthetic electron transport among ecotypes of Chenopodium album［J］. Functional Plant Biology，2003，30（10）：1071-1079.

［13］Min L，Zhu L，Tu L，et al. Cotton Gh CKI disrupts normal male reproduction by delaying tapetum programmed cell death via inactivating starch synthase ［J］. The Plant Journal，2013，75（5）：823-835.

［14］Wang Q，Liu N，Yang X，et al. Small RNA-mediated responses to low-and high-temperature stresses in cotton ［J］. Scientific Reports，2016，6（1）：1-14.

［15］Ma Y，Min L，Wang J，et al. A combination of genome-wide and transcriptome-wide association studies reveals genetic elements leading to male sterility during high temperature stress in cotton ［J］. New Phytologist，2021，231（1）：165-181.

［16］Barash Y，Calarco J A，Gao W，et al. Deciphering the splicing code ［J］. Nature，2010，465（7294）：53-59.

［17］Zhao L，Zhang H，Kohnen M V，et al. Analysis of transcriptome and epitranscriptome in plants using PacBio Iso-Seq and nanopore-based direct RNA sequencing ［J］. Frontiers in Genetics，2019，10：253.

［18］Abdel-Ghany S E，Hamilton M，Jacobi J L，et al. A survey of the sorghum transcriptome using single-molecule long reads ［J］. Nature Communications，2016，7（1）：1-11.

［19］Wang B，Tseng E，Regulski M，et al. Unveiling the complexity of the maize transcriptome by single-molecule long-read sequencing ［J］. Nature communications，2016，7（1）：1-13.

［20］Wang Q，Rio D C. JUM is a computational method for comprehensive annotation-free analysis of alternative pre-mRNA splicing patterns ［J］. Proceedings of the National Academy of Sciences，2018，115（35）：EB181-EB190.

［21］Wang X，Xu Y，Zhang S，et al. Genomic analyses of primitive，wild and cultivated citrus provide insights into asexual reproduction ［J］. Nature Genetics，2017，49（5）：765-772.

［22］Li H. Minimap2：pairwise alignment for nucleotide sequences ［J］. Bioinformatics，2018，34（18）：3094-3100.

［23］Tardaguila M，de la Fuente L，Marti C，et al. SQANTI：extensive characterization of long-read transcript sequences for quality control in full-length transcriptome identification and quantification ［J］. Genome Research，2018，28（3）：396-411.

［24］Trincado J L，Entizne J C，Hysenaj G，et al. SUPPA2：fast，accurate，and uncertainty-aware differential splicing analysis across multiple conditions ［J］. Genome Biology，2018，19（1）：40.

［25］Patro R，Duggal G，Love M I，et al. Salmon provides fast and bias-aware quantification of transcript expression ［J］. Nature Methods，2017，14（4）：417-419.

［26］Love M I，Huber W，Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 ［J］. Genome Biology，2014，15（12）：550.

［27］Cantalapiedra C P，Hernández-Plaza A，Letunic L，et al. eggNOG-mapper v2：functional annotation，orthology assignments，and domain prediction at the metagenomic scale ［J］. Molecular Biology and Evolution，2021，38（12）：5825-5829.

［28］Yu G，Wang L G，Han Y，et al. clusterProfiler：an R package for comparing biological themes among gene clusters ［J］. Omics，2012，16（5）：284-287.

［29］Xie C，Mao X，Huang J，et al. KOBAS 2.0：a web server for annotation and identification of enriched pathways and diseases ［J］. Nucleic Acids Res，2011，39（Web Server issue）：W316-W322.

［30］Chaudhary S，Jabre I，Reddy A S N，et al. Perspective on alternative splicing and proteome complexity in plants ［J］. Trends in Plant Science，2019，24（6）：496-506.

［31］Kim S，Kim T H. Alternative splicing for improving abiotic stress tolerance and agronomic traits in crop plants ［J］. Journal of Plant Biology，2020，63（6）：409-420.

［32］Laloum T，Martín G，Duque P. Alternative splicing control of abiotic stress responses ［J］. Trends in Plant Science，2018，23（2）：140-150.

［33］Cui J W，Shen N，Lu Z G，et al. Analysis and comprehensive comparison of PacBio and nanopore-based RNA sequencing of the Arabidopsis transcriptome ［J］. Plant Methods，2020，16（1）：85.

［34］Ivanov A G，Velitchkova M Y，Allakhverdiev S I，et al. Heat stress-induced effects of photosystem Ⅰ：an overview of structural and functional responses ［J］. Photosynthesis Research，2017，133（1/2/3）：17-30.

［35］Allakhverdiev S I，Kreslavski V D，Klimov V V，et al. Heat stress：an overview of molecular responses in photosynthesis ［J］. Photosynthesis Research，2008，98（1/2/3）：541-550.

［36］Shakeel S N，Aman S，Haq N U，et al. Proteomic and transcriptomic analyses of Agave americana in response to heat stress ［J］. Plant Molecular Biology Reporter，2013，31（4）：840-851.

［37］Li W M，Wei Z W，Qiao Z H，et al. Proteomics analysis of alfalfa response to heat stress ［J］. PLoS One，2013，8（12）：e82725.

［38］Spence A K，Boddu J，Wang D F，et al. Transcriptional responses indicate maintenance of photosynthetic proteins as key to the exceptional chilling tolerance of C4 photosynthesis in Miscanthus×giganteus［J］. Journal of Experimental Botany，2014，65（13）：3737-3747.

收稿日期：2022-05-17

基金項目：國家自然科學基金（編號：32060149、31760316）；海南省自然科學基金（編號：320RC500、321RC469）。

作者簡介：丁 宇（1996—），女，海南?？谌?，碩士研究生，主要從事轉(zhuǎn)錄組、基因組等生物信息學研究。E-mail：woshidingyuya@163.com。

通信作者：謝尚潛，博士，教授，主要從事生物信息學和基于三代測序的基因組學和轉(zhuǎn)錄組學研究。E-mail：sqianxie@foxmail.com。