不同果肉厚度苦瓜轉錄組比較分析

裘波音，鄭 旋，李大忠，林 琿，張前榮，劉建汀，朱海生，溫慶放

(1.福建省蔬菜遺傳育種重點實驗室,福州 350013；2.福建省農業科學院作物研究所,福州 350013；3.福建省蔬菜工程技術研究中心,福州 350013；4.福建省種子總站(福建種子質量檢驗站)，福州 350003)

【研究意義】苦瓜(MomordicacharantiaL.)屬于葫蘆科苦瓜屬一年生蔓性草本植物，起源于熱帶，喜光喜濕，對肥料要求較高，花、果期5—10月[1]。苦瓜因含有羅漢果苷和苦瓜素等成分而呈苦味，是一種重要的食藥兼用園藝作物，具有重要經濟價值。果肉厚度是瓜果類作物重要的品質性狀，也是產量的決定性因素，即果肉越厚，果實的可食用部分就越大[2]。果肉發育是果實發育的重要階段，也是關系果實整體品質的核心因子[3]。果實發育過程受一系列基因表達引起的生理生化變化調控，最終形成不同的顏色、質地、香氣、營養成分等，導致產量及品質上的差異[4-5]。研究果肉發育一方面有利于深化植物發育相關機制的探討，另一方面也對推進理論成果轉化應用起到積極作用。【前人研究進展】高通量轉錄組測序作為后基因時代率先發展起來的技術，已經在生物學領域得到了廣泛的應用。在蔬菜相關研究中，轉錄組測序技術主要用于分子標記開發，如劉峰等[6]基于5份辣椒的轉錄組數據篩選獲得18 159個SNP，且選取的15對驗證引物中有8對表現出多態性；林琿等[7]在24份花椰菜中挖掘出10 715個轉錄組SSR位點，分析獲得具有潛在高多態性的引物6119對；劉建汀等[8]在印度南瓜11 871個Unigenes中發現其中5391個含有SSR位點，且含有至少2個SSR位點的有319個；遺傳育種和進化分析，如Koenig等[9]利用轉錄組分析了番茄栽培品種和5個近緣野生品種的基因序列和表達差異，發現序列變異方面有50多個基因出現了正向選擇，而在表達水平方面，有數千個基因出現了變化，其中很多包括果色、耐熱、耐鹽等在內的相關基因其表達差異是由選擇壓力造成的，表明人工馴化在栽培植物的形態、生理和生活史方面產生了重要影響；差異表達基因分析和功能注釋，如王明霞等[10]對安徽烏菜與普通白菜共3份材料進行測序，其中200個KEGG注釋基因中參與植物病原菌互作通路的占11%，位居第一，參與植物激素信號轉導通路和類苯基丙烷生物合成通路的占9.5%，排第二；吳萍等[11]對來源于6個產地的不同種白芷進行分析共獲得15 939個差異基因，其功能主要與細胞組分、催化活性、代謝過程等相關，其中在川白芷中上調最高的10個基因多與蛋白質加工、激素信號傳導等相關，而下調最多的10個基因則與RNA運輸、果糖代謝等相關；功能基因挖掘，如基于轉錄水平獲得的BvM14-SAMS2基因在甜菜M14品系根中表達最高，并受鹽脅迫誘導表達[12]；而WRKY家族各成員在甜菜不同組織間表達差異大，其中BvWRKY33在熱脅迫和鹽脅迫下均呈上調表達[13]，表明上述基因對甜菜逆境生理調控具有重要作用。【本研究切入點】近年來，苦瓜營養價值被不斷開發和利用，促使其生產在全世界逐漸流行，這對苦瓜種質創新提出了挑戰。了解不同果肉厚度苦瓜果實成熟過程中的基因表達情況，有利于探究果實發育調控機理，推動良種培育，滿足人類對其持續消費的需求。【擬解決的關鍵問題】本研究立足果肉厚度，對不同苦瓜種質進行統計分析，篩選果肉厚度差異顯著材料，再利用轉錄組測序技術對其進行比較分析，挖掘差異表達基因，預測轉錄因子，并對差異表達基因進行功能注釋、代謝通路富集等分析，為進一步揭示苦瓜果肉發育分子調控機制提供理論依據。

1 材料與方法

1.1 試驗材料

70份苦瓜品系由福建省農科院作物研究所蔬菜研究室提供，于2020—2021年春季種植于福州市試驗基地(119.28°E, 26.08°N)。

1.2 果肉厚度差異苦瓜篩選

于成熟期采樣，每個品系采長勢相似的4個瓜，在離瓜蒂2/3處橫切，用游標卡尺測量橫徑和瓤腔寬度，計算果肉厚度，公式如下：果肉厚度=(橫徑-瓤腔寬度)×0.5。根據果肉厚度的顯著性差異篩選厚肉和薄肉材料各1份進行后續試驗。

1.3 RNA提取與檢測

于樣品成熟期對厚肉型苦瓜和薄肉型苦瓜各1份進行取樣，每份材料取9個瓜，設3個生物學重復，即每個重復為3個瓜混樣。采摘后立即置于液氮中保存，用于總RNA的提取。利用OmniPlant RNA(DNase I)試劑盒(康為世紀，北京)進行RNA提取；1%瓊脂糖凝膠電泳分析樣品RNA完整性及是否存在DNA污染；利用NanoPhotometer spectrophotometer檢測RNA純度(OD260/280及OD260/230比值)；利用Agilent 2100 bioanalyzer精確檢測RNA完整性。樣品送福建華諾通達公司進行轉錄組測序。

1.4 文庫構建與質檢

使用Illumina的NEBNext?UltraTM RNA Library Prep Kit對總RNA進行建庫。文庫構建完成后，先使用Qubit2.0 Fluorometer進行初步定量，隨后使用Agilent 2100 bioanalyzer對文庫的insert size進行檢測，qRT-PCR對文庫有效濃度進行準確定量(文庫有效濃度高于2 nmol/L)，以保證文庫質量。

1.5 上機測序及功能注釋

采用Illumina HiSeq2500進行測序，對測序所獲得的原始數據，按照華諾通達公司的數據處理方法進行過濾，去除帶接頭(adapter)、含N(N表示無法確定堿基信息)和低質量reads(Qphred≤20的堿基數占整個read長度的50%以上的reads)，再對clean data進行Q20，Q30和GC含量計算。后續所有分析均是基于clean data進行的高質量分析。

本研究采用有參轉錄組測序(Momordica charantia ASM199503v1，https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/995/035/GCF_001995035.1_ASM199503v1/)，使用HISAT2 v2.0.5構建參考基因組的索引，并將配對末端clean reads與參照基因組比對。通過BlastX軟件，與NCBI的non-redundant(NR)、UniProt/Swiss-Prot、Gene Ontology(GO)、Kyoto Encyclopedia of Genes and Gnomes(KEGG)等數據庫進行比對，以獲得其功能注釋信息。

1.6 差異基因鑒定及富集分析

通過FPKM法來計算基因的表達量，即每百萬fragments中來自某一基因每千堿基長度的fragments數目。差異基因篩選基于“有生物學重復”類型，利用DESeq2軟件和負二項分布的P-value計算模型，以“|log2(Fold Change)| > 0，padj <0.05”為標準進行。用BLAST2GO軟件進行GO注釋，用KOBAST軟件進行KEGG注釋。并采用clusterProfiler軟件對差異基因集進行GO功能富集分析和KEGG通路富集分析。

1.7 數據分析

數據采用SPSS 21.0進行統計分析，采用WPS Office進行作圖。

2 結果與分析

2.1 不同果肉厚度苦瓜篩選

如表1所示，測定的70個品系材料果肉厚度分布范圍較廣，最大值為12.85 mm，最小值為8.65 mm，平均厚度為9.84 mm，整體呈正態分布，其中大部分厚度在9～11 mm，小于9 mm的有7個，大于11 mm的有3個。基于一般線性模型的單因素分析進行多重比較，方法選用Tamhane’s T2(P≤0.01)，發現LX1-3(12.85 mm)和ZK54(8.83 mm)果肉厚度具有極顯著差異，前者是后者的1.46倍，因此，將LX1-3作為厚肉型材料，ZK54作為薄肉型材料進行后續分析。

表1 苦瓜不同果肉厚度統計Table 1 Summary of flesh thickness of different bitter gourds

2.2 轉錄組測序統計

經轉錄組測序分析統計，LX1-3共計得到241 306 304條初始序列，過濾后得到234 825 400條干凈序列，17.61 G干凈堿基對，平均Q30值為92.74%；比對到參考基因組上的序列平均為34 455 982條，占總序列數的88.07%，其中，單一比對序列為34 027 357條，占比86.98%，多重比對序列為428 625條，占比1.10%，GC含量為46.54%；ZK54共計得到250 983 052條初始序列，過濾后得到241 318 020條干凈序列，18.1 G干凈堿基對，平均Q30值為92.26%；比對到參考基因組上的序列平均為35 079 431條，占總序列數的87.23%，其中，單一比對序列為34 607 181條，占比86.06%，多重比對序列為472 250條，占比1.17%，GC含量為46.75%(表2)。

表2 樣品轉錄組數據統計Table 2 Summary of transcriptome sequencing results of different samples

對序列長度統計分析(圖1)顯示，長度小于200 bp的有54條，占0.28%；長度為200～500 bp的有684條，占3.49%；長度為501～1000 bp的有3672條，占18.73%；長度為1001～2000 bp的合計共有8578條，占43.76%；長度為2001～3000 bp的有4169條，占21.27%；長度為3001～4000的有1555條，占7.93%；長度超過4000 bp的有890條，占4.54%。總體呈“兩頭低中間高”的趨勢，數量最多的序列長度集中在1001～2000 bp。

圖1 苦瓜轉錄組序列的長度分布Fig.1 The length distribution of unigenes of transcriptome in bitter gourd

2.3 差異基因鑒定及轉錄因子預測

經篩選共得到1492個差異基因，其中796個基因上調表達，696個基因下調表達，分別占差異基因數量的53.35%和46.65%，從圖2中可以看岀,2組材料之間的基因表達水平呈顯著差異。本次測序共預測到598個轉錄因子，分屬于216類，│log2(Fold Change)│值在0.62～11.24。其中，Pkinase數目最多(45個，7.53%)，其次是p450(29個，4.85%)，接著為Pkinase_Tyr(21個，3.51%)，2OG-FeII_Oxy(17個，2.84%)，UDPGT(16個，2.68% )，Myb_DNA-binding(14個，2.34%)，AP2(12個，2.01%)和HLH(10個，1.67%)，其他家族共有基因434個，占72.58%(圖3)。

圖2 不同果肉厚度苦瓜差異基因Fig.2 Differential gene of bitter gourd with different sarcocarp thickness

圖3 苦瓜序列的轉錄因子預測 Fig.3 Transcription factor prediction of Unigenes in bitter gourd

2.4 GO功能富集分析

共有1452個差異基因被注釋到了GO數據庫中，得到616條功能注釋，共分為三大功能，分別是生物過程、細胞組分和分子功能，依次得到308、93和215條功能注釋，占比分別為50.00%、15.10%、34.90%。每項功能中重要注釋(P<0.05)如表3所示。生物過程含DEGs 556個，共獲得308條注釋。最主要的11條功能注釋共包含151個DEGs，占27.16%。從差異基因數量上可看出，“碳水化合物代謝過程”條目居于第一位，含47個DEGs，其次為“細胞碳水化合物代謝過程”和“多糖代謝過程”，分別含18和13個DEGs，條目“細胞葡聚糖代謝過程”，“葡聚糖代謝過程”，“細胞多糖代謝過程”和“抗氧化反應”并居第四，分別含12個DEGs。從差異基因的上下調變化上得出，11個條目中上調基因數目多于下調基因數的有7條，且其中6條位于前六位。細胞組分含DEGs 180個，功能注釋93條。最主要的前9條功能注釋(P<0.05)共含75個DEGs，占41.67%。其中，“細胞周邊”條目含DEGs數目最多，為16個，其次為“細胞壁”和“外部封裝結構”，均含有11個。9個條目共計含有上調基因55個，下調基因20個，上調基因數目是下調基因的2.275倍。分子功能含DEGs 716個，功能注釋215條。前27條主要注釋條目含有587個DEGs，占81.98%。“血紅素結合”和“四吡咯結合”最多，均含有43個DEGs，其次為“轉移酶的活性/轉移糖基” “DNA結合轉錄因子活性” “轉錄因子活性”，各含有41個DEGs。27個條目共計含有上調基因328個，下調基因259個，上調基因數目多于下調基因。從上述結果看出，生物過程，細胞成分和分子功能三方面的上調基因數目均大于下調基因數，說明果肉厚度的發育受上調基因的表達影響較大。

表3 不同果肉厚度苦瓜差異表達基因的GO富集分析Table 3 GO enrichment analysis of differentially expressed genes of bitter gourd with different sarcocarp thickness

續表3 Continued table 3

2.5 KEGG通路富集分析

共有308個差異基因富集到了KEGG通路上，包括188個上調基因和120個下調基因，結果如圖4所示。上調基因有178個富集到20條主要的代謝通路。富集因子最高的為“卟啉和葉綠素代謝”通路，其次為“光合作用生物的碳固定作用”和“光合作用”通路。而富集數目最多的通路為“核糖體”，其次為“碳代謝”，兩者所占DEGs比列分別為12.77%和9.04%。下調基因有111個富集到20條主要的代謝通路。其中，富集因子最高的通路為“苯丙素生物合成”，其次為“二萜生物合成”通路。富集下調差異基因最多的通路為“苯丙素生物合成”，占差異基因比例為10.00%，其次為“淀粉和蔗糖代謝”“氨基酸的生物合成”和“碳代謝”，三者數量相當，所占比列為7.5%。從富集DEGs的數量角度也可以看出，厚肉苦瓜材料的代謝途徑中上調基因數量高于下調基因數量，是下調基因數量的1.56倍，說明果肉厚度發育可能主要受上調基因表達調控。

圓點的位置代表代謝途徑，圓點的大小代表基因的數量The position of the dot represents the enriched pathway，the size of the dot indicates the number of differential genes圖4 差異表達基因的KEGG分析Fig.4 The KEGG analysis of the differentially expressed genes

3 討 論

苦瓜是重要的蔬菜作物，其果實生長發育受細胞分裂、激素水平、信號傳導、糖碳基礎代謝、環境條件等多種因素的影響。果肉厚度作為一個重要的性狀，對其產量和質量起決定性作用。立足果肉厚度的早期研究報道不多，何曉明[14]、徐強[15]等對黃瓜產量和果實性狀進行調查，得出果肉厚度對黃瓜產量影響重大，基于相關分析也得出果實的鮮重與果肉厚度的遺傳呈正相關，楊濤等[16]也發現辣椒果肉厚度與果實橫徑、單果重也極顯著正相關。本研究通過對70份苦瓜的果肉厚度測定分析，發現不同材料的果肉厚度差異較大，最大值與最小值差值占比平均果肉厚度超四成，說明果肉厚度受基因型調控表現出種系間差異，進而可能對產量及品質造成影響。LX1-3和ZK54兩者果肉厚度具有極顯著差異，前者是后者的1.46倍，但其它生長性狀基本相似，因此將兩者作為后續轉錄組測序的研究材料對保證結果真實性和準確性具有有利的一面，對深入研究其次生代謝、產量形成、品質進化等具有重要意義。

轉錄組測序是目前生物學研究中的主流技術，通過分析不同材料在不同時期或環境下的轉錄組差異，有助于了解它們的差異表達在不同生物學過程中的作用[17]。早期研究發現，果實品質在成熟過程中會發生一系列生理生化變化，這些變化與各種基因的表達調控有關，如黃瓜基因Csa2GO58670與兩親本(厚果肉和薄果肉)果肉厚度均顯著相關，果實發育不同時間該基因在兩個品種中相對表達量存在差異[18]。本研究基于不同果肉厚度苦瓜進行轉錄組分析，發現多個差異表達基因，其中上調基因數目明顯多于下調基因數，表明其果肉發育可能受上調基因影響較大。作為轉錄過程的重要元件，轉錄因子能與5’端上游特定序列專一性結合，從而保證目的基因以特定的強度在特定的時間與空間表達，在真核生物中，轉錄因子能與RNA聚合酶Ⅱ形成轉錄起始復合體，共同協助其完成復雜的轉錄過程。 本研究預測得到598個轉錄因子，但分屬于多個家族，說明參與苦瓜果肉發育調控的轉錄因子種類繁多，再次證實植物發育是一個極其復雜的過程，最終果實的質量和品質是多種因素綜合調控的結果。

GO和KEGG分析是基因注釋的重要手段，GO的功能分類有利于了解基因功能，而KEGG富集分析除了可以推測基因功能，還可以研究基因在不同代謝通路中的位置及作用[19]。本研究GO注釋表明，生物過程在三大功能中占主要地位，深入挖掘發現糖類代謝相關注釋所含差異基因最多且上調基因數均多余下調基因數，同時，KEGG富集分析也發現，上調基因與葉綠素、碳固定以及光合作用相關途徑富集因子最高，而核糖體及碳代謝途徑富集的差異基因數量最多，表明在維持苦瓜果肉正常發育中，基礎糖碳代謝和光合作用以及源庫關系在整個生物過程中至關重要。光合作用是植物進行物質生產的基本生理過程，植物絕大部分干物質是由光合作用提供的，光合作用的重要產物及其運輸的主要形式是蔗糖[20]，光合組織中蔗糖不僅影響著光合同化產物的分配運輸，其含量的多少直接影響籽實中糖分含量，對糖類整體代謝、果實的膨大與生長速率均造成干擾，進而影響作物產量與品質[21-23]。同時，各類酶在保障葉綠體內光合碳代謝、蔗糖運輸、多糖積累平衡等生物過程中起著重要的調節作用，并且與多種調控作物生長和產量的指標密切相關[24-25]。本研究基于轉錄組測序發現參與苦瓜果肉發育生物過程的重要調控因子與早期研究結果類似，再次證明糖碳代謝是保證作物生長發育的關鍵。早期研究表明，源庫器官的數量及其機能的相互協調，是作物生產力的重要決定因素。在本次試驗中，厚肉型苦瓜相較薄肉型苦瓜在糖代謝和光合作用相關基因表達水平更高，表明其具有較高的“源”性能，為多糖形成提供基礎，而“淀粉和蔗糖代謝”基因下調，表明“庫”性能在化合物分配時將較多的基底物質傾向于葡聚糖等多糖的形成，因此推測，改善光合作用性能、調節源庫平衡可能是提高苦瓜產量及品質的重要途徑。

眾所周知，生長素、赤霉素、細胞分裂素等植物內源激素對果實的發育和成熟也發揮著重要作用，不同的激素在不同品種以及同一品種的不同發育階段作用不一。IAA含量的提高能促進發育早期的番茄子房細胞的分裂和膨大以及誘導果實發育[26]；GA種類眾多，在幼果中含量較高，能刺激果實和種子發育[27]；茄子子房(果實)發育前期的IAA、ZR含量明顯高于發育后期，其中單性結實品系的IAA、ABA含量和ABA/(IAA+GA4+ZR)均高于非單性結實品系[28]；黃瓜Q30(CsSUN/CsLNG1)的BR/ABA與GA3/ABA在開花前6 d顯著低于其近等基因系材料[27]。本研究對下調基因進行KEGG分析發現，“苯丙素生物合成”是富集因子最高且基因數最多的通路，這個通路主要涉及一些過氧化物酶，其次為“二萜生物合成”，以赤霉素加雙氧酶為主。植物體內的過氧化物酶對生長素有破環作用，對調控生長素水平，維持生長素庫源的平衡非常關鍵。而赤霉素加雙氧酶參與赤霉素代謝，對維持植物體內活性赤霉素的水平及含量具有重要作用。上述結果表明植物激素參與苦瓜果肉發育，再次證實激素平衡在果實發育中的重要地位。

細胞壁組分和代謝酶對植物生長發育具有重要影響，如基于甜瓜和黃瓜材料在成熟過程的轉錄組測序數據，發現了果實膨大和發育與細胞壁相關基因包括擴張蛋白基因(Csa6G014540)、木葡聚糖內轉葡糖基酶基因(MelO3C014468，Csa1G422480)等關系密切[29-30]；而隨著果實的不斷成熟，砂梨中β-半乳糖苷酶基因PpGal1和PpGal4表達量呈逐漸上升趨勢[31]，鱷梨中PaGal2表達量增高但PaGal4則不表達[32]。本次研究發現細胞組分條目中涉及細胞周邊、細胞壁以及外部封裝結構的基因如果膠酯酶基因、木葡聚糖內糖基轉移酶/水解酶基因等數目占比較大，表明細胞壁相關結構在苦瓜果肉發育中也起到了較大的作用。

4 結 論

通過研究獲得的果肉厚度差異苦瓜種質為培育高產優質苦瓜提供潛力材料；高通量轉錄組測序挖掘的差異基因以及相關功能注釋和富集分析結果，為苦瓜果實發育分子調控機制的探討提供理論基礎。