油茶sRNA分析及其對光合油脂代謝的調控

關鍵詞：油茶；sRNA；miRNA；靶基因；調控；光合油脂代謝

中圖分類號：S601；S794.4 文獻標志碼：A 文章編號：1003—8981（2024）03—0025—11

sRNA 是植物發育和種子形成的重要調節因子。植物中的sRNA 包括short interfering RNA（siRNA） 和microRNA（miRNA）， 其中，siRNA 包括ra-siRNA（重復相關小干擾RNA）、天然反義轉錄物衍生的小干擾RNA（naturalantisense transcript-derived small interfering RNA，natsiRNAs）、反式作用反義小分子干擾RNA（transactingantisense interfering RNAs，ta-siRNAs）、異染色質小分子干擾RNA（heterochromatic smallinterfering RNA，hc-siRNA）、初級siRNAs（PrimarysiRNAs） 和長的小分子干擾RNA（long smallinterfering RNAs，lsiRNA）[1-2]。其中，miRNA 是一類真核生物中廣泛存在長度約為20 ～ 24 個核苷酸（nucleotide，nt）的內源非編碼單鏈RNA，它通常是由RNA 聚合酶Ⅱ合成，其最初的轉錄產物稱初始miRNA（pri-miRNA 初級轉錄產物），初始miRNA 被DCL（Dicer-like）酶切割成雙鏈premiRNA（precursor miRNA），它通過堿基互補直接作用于靶基因mRNA，誘導靶基因的特異性切割，在轉錄和轉錄后水平調節基因的表達[3-4]。如，從馬鈴薯葉、花、根、塊莖的sRNA 表達圖譜發現sRNA 進化率低并且較為保守，具有組織特異表達的特性[5]，研究還發現：在龍眼花發育過程中，花芽中的sRNA 富集表達，通過信號途徑參與花器官發育轉錄調控網絡的特異調控[6]；在鵝掌楸葉和花中，miR157a-SPL 和miR160a-ARF 等miRNA差異表達[7]。因此，sRNA 在植物營養生長、生殖生長[8-9]、油脂形成[10] 等發育過程和模式中發揮重要作用[11]。

高通量測序由于可以同時測序幾十萬到幾百萬條RNA 分子，讀取較短序列，檢測全面、靈敏度高，因而加快了sRNA 的發現和功能研究[12]。目前，已經從179 種植物中鑒定出38 186 個miRNA 和7 838 個家族[13-15]。在擬南芥中miRNA參與胚形成并且抑制靶基因的表達，調控胚成熟[16]。如，在胚缺失突變體lacking DCL1 八細胞胚期， 上調較多的是miR156 靶基因、編碼SPL10 和SPL11 的轉錄因子。在胚成熟過程中，miRNA156 誘導介導抑制zygotic SPL 轉錄， 抑制轉錄的近成熟積累[16]。Lin 等[17] 以龍眼的同步胚性培養物為材料，研究體細胞胚胎發生過程的miRNA 和靶基因，發現了272 個miRNA 和2063個靶基因，其中，181 個miRNA 和862 個靶基因得到證實。這些靶基因參與植物代謝、信號傳導和應激響應。研究還發現：20 個保守miRNA 和4個新miRNA 參與體細胞胚狀體的發育過程。在油菜Brassica napus 中發現了50 個保守miRNA，其中，一些miRNA 影響油脂的形成積累[10]。然而，對多數sRNA 的表達和功能仍然了解甚少。本研究以油茶葉片和種子為材料構建sRNA文庫，利用Illumina高通量測序技術，確定已知miRNA 和新miRNA，篩選差異表達miRNA，預測靶基因，進行靶基因GO功能富集分析和KEGG通路富集分析，揭示miRNA 與光合作用、碳固定和油酯代謝過程中對應靶基因間可能的調控關系，為利用miRNA 技術進行油茶遺傳改良提供科學依據。

1 材料與方法

1.1 RNA提取和構建sRNA文庫

采集油茶‘橫沖89’的嫩葉和成熟種子，立即液氮冷凍，-70 ℃貯藏。利用CTAB 裂解法提取總RNA。以總RNA 為起始樣品，利用SmallRNA Sample Pre Kit，將Small RNA 加上接頭反轉錄合成cDNA，PCR 擴增，分離回收DNA，構建cDNA文庫。利用Qubit2.0 和Agilent 2100對文庫的insert size 進行定量檢測和完整性分析，采用HiSeq/MiSeq2000測序獲得原始序列。

1.2 測序數據質量評估和sRNA的結構特征、種類和數量、公共序列及特有序列分析

通過Phred數值，轉化獲得堿基測序錯誤率（e），每個堿基測序錯誤率通過公式Qphred=-10log10（e）轉化得到，以Phred 分值及對應的Q20、Q30評估數據質量。去除原始讀長（raw reads）中低質量、含N、5′ 接頭污染、無3′ 接頭序列、polyA/T/G/C的讀長，得到凈讀長（clean reads）。在18～40nt 的長度區間篩選sRNA，分析sRNA 的長度分布、種類、數量以及公共及特有序列等。

1.3 已知miRNA 分析、新miRNA預測以及miRNA家族分析

用bowtie 把sRNA 定位到拼接的油茶轉錄組序列上，將mapped 到參考序列上的sRNA 序列與miRBase 序列比對，比對分析已知miRNA，分析比對上的總sRNA 種類、數量等[18-19]。分析已知miRNA 首位堿基的偏好性。利用miREvo[15]和mirdeep2[20]miRNA 軟件，通過miRNA 的發夾結構，分析比對上miRNA 前體和成熟體，預測新miRNA（novel miRNA），對miRNA 進行家族分析，探索miRNA 家族在其他物種中的存在情況。

1.4 miRNA表達、差異分析和聚類分析

從miRNA 表達水平分析中得到的read count數據，采用TPM對read count 數據進行歸一化處理，歸一化表達量= Mapped read count/Totalreads×1 000 000。利用DEGseq（2010）R package進行差異分析。用q value 調整P-value，以q value ＜0.01 及|log2（fold change）| ＞ 1作為默認閾值進行顯著差異性分析。將差異miRNA 的TPM 值用于聚類分析。

1.5 miRNA靶基因預測、候選靶基因GO 功能富集分析

和KEGG 通路富集分析通過psRobot_tarin psRobot 預測miRNA 靶基因。利用GO（Gene Ontology，http：//www.geneontology.org/）數據庫對差異表達miRNA 靶基因進行功能富集分析，將差異表達miRNA 的候選靶基因按照生物學過程、細胞組分、分子功能進行富集條目（term）分析。通過KEGG（KyotoEncyclopedia of Genes and Genomes）通路顯著性富集，分析差異表達 miRNA 候選靶基因參與的主要代謝途徑和信號轉導途徑。

1.6 參與光合作用、油酯代謝miRNAs 的表達分析和靶基因篩選

根據miRNA 的差異表達、候選靶基因的KEGG pathway 富集分析，篩選油茶光合作用、脂肪酸油脂代謝通路的miRNA 靶基因，預測差異表達miRNA 與靶基因可能的調控關系。

2 結果與分析

2.1 總RNA 提取與sRNA 文庫制備

提取油茶幼葉和成熟種子的總RNA，用DNase Ⅰ除去殘留的基因組DNA（圖1）。RNA濃度A260/A280 分別為2.119 和2.093。總RNA 濃度達到了建庫要求，Agilent 2100 Bioanalyzer（AgilentTechnologies，USA）測定RNA 完整性，結果表明：RIN（RNA integrity number） 為8.2 和8.9， 符合建庫標準。分離富集葉片和種子的sRNA，構建sRNA 文庫。

2.2 測序數據質量評估和sRNA 的結構特征、種類和數量

通過測序獲得油茶葉sRNA 的原始讀長為15 659 389 條，測序錯誤檢驗的錯誤率為0.01%，Q20、Q30 分別達到99.38% 和97.69%，堿基GC 含量為49.15%；種子sRNA 的原始讀長為14 380 199條，測序檢驗錯誤率為0.01%，Q20 和Q30 分別達到99.34% 和97.59%，堿基GC 含量為49.72%（表1）。

在葉和種子sRNA 中，無5′ 接頭或無插入片段讀長分別為261 297 條（1.67%）和405 709 條（2.82%）；低質量讀長分別為14 239 條（0.09%）和8 125 條（00.6%）、含N% ＞ 10% 的讀長分別為39條和22條；含5′ 接頭的讀長分別為11942 條（0.08%）和9，253 條（0.06%），含有polyA/T/G/C的讀長分別為67733條（0.43%）和48174條（0.34%），除去上述讀長，獲得葉和種子sRNA 的凈讀長分別為15304139條和13908934條（圖2）。

在篩選獲得油茶葉和種子的sRNA 中，18 ～ 40 nt 的凈讀長數量分別為14 256 189 條和10 822 313條，屬于6 138179和4980518 種sRNA。與其他植物相似，油茶sRNA 的長度分布在18 ～ 40 nt，其中，24-nt sRNA 最多，在葉和種子中分別為4678323條（32.82%）和3 516 310條（32.49%）；其次是21-nt sRNA，分別為954 749條（6.7%） 和1142946 條（10.56%）；22-ntsRNA 居于第3，分別為952210條（6.68%）和941611（6.12%）條；23-nt 的居于第4，分別為872 524 條（8.7%）和749 356 條（6.92%）（圖3）。油茶葉和種子sRNA 公共序列的種類為773647種（7.48%），特有序列的種類分別為4206871種（40.67%） 和5364532種（51.86%）；公共序列的數量為1489104條（59.38%）。葉和種子sRNA 特有序列的數量分別為4418438條（23.0%）和5768160條（17.62%）。

2.3 已知miRNA 分析、新miRNA預測以及miRNA家族分析

用bowtie 將sRNA 定位到trinity 拼接的轉錄組上，結果發現：葉和種子中的sRNA 分別為14 256 189 條和10 822 313 條，比對到參照序列的sRNA 分別為5 855 009 條（41.07%）和4 122 070條（38.09%）。其中，與參考序列同向的分別為933 734（6.55%） 條和770 311 條（7.12%）， 與參考序列反向的分別為4 921 275（34.52%）條和3 351 759 條（30.97%）。比對分析已知miRNA，結果發現： 比對上的總sRNA（Total sRNA）為120 102 條，在葉和種子中分別為54 035 條和66 067 條；比對上的sRNA 為518 種，在葉和種子中分別245 種和273 種；比對上的miRNA 前體為118種，其中，在葉和種子中分別為112 種和109 種；比對上的miRNA 成熟體為78 個，其中，在葉和種子中分別為74 個和72 個。其中，ath-miR399d、osa-miR399e、ath-miR399f、gma-miR169v、osamiR1863a、osa-miR1863b 的成熟體僅出現在葉中，而ath-miR171a、ath-miR5658、gma-miR172b-5p、osa- miR1863b 的成熟體僅出現在種子中。

從圖4可以發現：18-23nt 的已知miRNA 首位堿基多偏好堿基U（uridine）， 其中，21-nt 已知miRNA 的首位堿基為U 的sRNA 數量最多，分別為45 579 條和54 137 條；22-nt 已知miRNA的次之，分別為2510條和3118條；20-nt 已知miRNA 的為第3位，分別為1707條和2775條。

利用miRDeep2軟件進行新miRNA（novelmiRNA）預測，結果發現：比對上新miRNA 前體的sRNA 種類為1178種，總sRNA 為20653個。其中，在葉中比對上的sRNA 種類為596種，總sRNA 為10330個；在種子中為582種，比對上前體的總sRNA 為10323個。分析發現了54個新miRNA，長度為19-24nt，以24-nt miRNA 為主。在葉片中比對上的新miRNA 成熟體出現了9 411次，在種子中的出現了9 504 次。此外，發現33個新miRNA 含有miRNA*（miRNA 的隨從鏈），在葉和種子中miRNA* 分別出現了443 次和464次（表2）。

2.4 miRNA表達、差異分析

和層次聚類分析統計分析miRNA 的讀長數（read count），TPM 歸一化處理這些讀長數，篩選出132 個差異表達的miRNA，其中，在葉中為128 個，在種子為126 個，122 個為兩者共表達，葉和種子表達的sRNA 讀長129 568 條和154 492 條。其中，ath-miR166a 對應讀長最多，在葉和種子中分別出現了26 594 次和29 821 次。DEGseq 差異分析發現：在葉和種子中有26 個顯著差異表達的miRNA（q value ＜ 0.01，|log2（fold change）| ＞ 1），顯著上調和下調表達的miRNAs 均為13 個。其中，顯著上調表達的miRNA 包含ath-miR169b、ath-miR399d、ath-miR399f、gma-miR169v、osamiR1863a、osa-miR399e 6 個已知miRNA 和novel82、novel 32、novel 54、novel 49、novel 53、novel 41、novel 83 6 個新miRNA。其中，athmiR169b、novel 82、ath-miR399d、ath-miR399f、gma-miR169v、osa-miR1863a、osa-miR399e 的log2 倍數差異表達值大于3。顯著下調表達的miRNAs 包括ath-miR172e、gma-miR172k、gmamiR396b-3p、gma-miR172c、gma- miR172d、ath-miR172c、ath-miR399b、osa-miR172c、gmamiR172b-5p、ath-miR171a、ath-miR5658、osamiR1863b12 個已知miRNA 和novel 72。其中，gma-miR172b-5p、ath-miR171a、ath-miR5658、osa-miR1863b 的log2 倍數差異表達值小于-3。

由圖5 聚類分析可知：在葉和種子miRNA 表達模式較為相似。差異表達的miRNA 被分為2 族，一族是在葉和種子中的表達量均較低，包括athmiR5658、osa-miR1863b、ath-miR171a 等。另一族是在葉和種子中的表達量較高，包括ath-miR399b、ath-miR172c、gma-miR369b-p、novel 53、novel 32、novel 54 等。另外，novel 49、novel 41、novel 82 在僅葉中表達量高；osa-miR172c、gma-miR172c 和gma-miR172d 等在僅種子中表達量高。

2.5 miRNA靶基因預測、候選靶基因GO功能富集分析和KEGG通路富集分析

對差異表達的miRNA 進行靶基因預測，得到67個差異表達miRNA 和它們對應的1346個靶基因。26 個顯著差異表達的miRNA 對應的靶基因為644個。結果還發現：一個miRNA 調控多個靶基因，如：ath-miR5658調控532個靶基因；同時，一個靶基因受多個miRNA 調控，如調控胚胎發育與細胞增殖的真核翻譯起始因子3類-L- 亞基基因（comp69386_c0） 受到ath-miR156g、athmiR156j、osa-miR156k 3 個miRNA 的調控。

根據差異miRNA 候選靶基因在GO 數據庫中的分布特征，計算顯著富集的GO條目（term）。從圖6可以看出，葉和種子差異表達miRNA 的候選靶基因在生物學過程、細胞組分、分子功能中的51 個條目中富集。其中，生物學過程中的生物調控過程、細胞過程、代謝過程等23 個條目富集較多；細胞組分中的細胞、細胞部分、細胞器、細胞器等17個條目富集較多；分子功能過程的結合、催化活性等11個條目顯著富集。KEGG 通路富集分析結果表明：候選靶基因在萜類主鏈生物合成、甘油酯代謝、氨酰tRNA 生物合成、RIG-I樣受體、核黃素代謝、角質軟木脂和蠟質生物合成等條目中顯著富集。

2.6 參與光合作用、油酯代謝miRNAs的表達分析和靶基因篩選

從表3可以發現：ath-miR5658 在葉片中沒有表達，而在種子中表達量為6.47。它對應靶基因的靶基因分別是光合碳固定、光合作用途徑中的類NADP 依賴型蘋果酶（comp89613_c0）、果糖二磷酸醛縮酶1（comp46464_c0）、光系統I 反應中心亞基ⅣB（comp71063_c0），甘油磷酸酯代謝中的甘油-3- 磷酸酰基轉移酶1（comp77968_c0），脂肪酸延長途徑中的3- 酮酰基-CoA 合成酶4（comp80917_c1），甘油脂代謝和甘油磷酸酯代謝中的1- 酰基-sn- 甘油-3- 磷酸酰基轉移酶5comp83553_c0 等。而gma-miR169v 在葉片中的表達量為7.72，而在種子中沒有表達。它的靶基因是甘油脂代謝中的單半乳糖酰二酰基甘油合酶2（comp90506_c0）。

3 討論

3.1 sRNA的數量、長度、種類鑒定和特征分析

sRNA 是植物發育和種子形成的重要調節因子。通過長度分布能夠鑒定sRNA 的種類，如，21～24 nt 的sRNA 主要是miRNA，多數植物的siRNA（small interfering RNAs） 為24 nt [2，17，26]。研究發現：油茶葉和種子中sRNA 最多的是24-ntsRNA，其次是21-nt sRNA。這與油茶抗炭疽病[27]略有差別，與Brassica napus[10]、Arabidopsis thaliana、Oryza sativa[28] 的sRNA 相似，第3、4 大分布群體是22-nt 和23-nt sRNA，在葉種子中分別占6.68%和6.12%；在種子中占8.7% 和6.92%（圖3）。

3.2 miRNA 的堿基偏好特征分析

植物miRNA 的5′ 端首個堿基偏好U 而抵抗G。可以通過miRNA 的U 堿基偏好性分析，獲得典型性miRNA 序列比例。從圖4 可以發現，長度在18-23 nt 的已知miRNA 首位堿基偏好U，這與中國甜柿[29] 中的miRNA 有一定的相似性，因此，增強了研究鑒定結果的可信度。

3.3 已知miRNAs家族分析、新miRNA 鑒定和靶基因分析

本研究從油茶葉和種子sRNA 中鑒定出22 個已知miRNAs 家族，105 個miRNAs 家族成員。其中，多數已知miRNA 與其他物種保守miRNA 具有高度相似性[27]。在鑒定出的miRNA 家族中，以多家族成員形式存在多個物種，這與相關研究結果[27] 相一致。本研究初步發現和證實了“一個miRNA 可能存在多個靶基因，而多個miRNA 可能均調控同一靶基因”的現象[27]。其中，67 個表達差異miRNA 調控1 346 個靶基因。在油茶中發現54 個油茶新miRNA，其長度以24 nt 為主，新miRNA 的數量和長度與楊利利[27] 對油茶的研究結果相一致，但與中國甜柿新miRNA 的長度分布（以21 nt 為主）略有差異[29]。篩選出的miRNAs 及其候選靶基因還有待進一步利用實時熒光定量PCR和特殊莖環qRT-PCR 進行驗證。由于缺少油茶的基因組序列，本研究無法進行莖環結構準確預測，導致所得的miRNA 可能是其他類型的sRNA。此外，由于缺少基因組注釋，可能會導致靶基因注釋不全。但是，隨著油茶的核基因組、葉綠體基因組、線粒體基因組等基因組組裝和測序解析研究的深入[30]，sRNA 種類的預測，特別miRNA 會更加準確。

3.4 調控光合作用和油酯代謝miRNA靶基因的表達分析

從表3 可以發現：ath-miR5658 在葉片中沒有表達，而在種子大量表達。而gma-miR169v 在葉片中大量表達，而在種子中沒有表達。前期的轉錄表達研究結果[19] 顯示：類NADP依賴型蘋果酶（comp89613_c0）、果糖二磷酸醛縮酶1（comp46464_c0）、光系統I反應中心亞基ⅣB（comp71063_c0）、單半乳糖基二酰基甘油合酶2（comp90506_c0）等基因在油茶葉中的轉錄表達較種子中均上調表達。這表明：ath-miR5658 沒有顯著抑制這些靶基因的轉錄表達，對油茶葉片光合產物碳固定、光合作用的影響較小。而靶基因甘油磷脂代謝1- 酰基-sn- 甘油-3- 磷酸酰基轉移酶5（comp83553_c0）、甘油-3- 磷酸酰基轉移酶1（comp77968_c0）、3- 酮酰基-CoA 合成酶4（comp80917_c1）在葉和種子中的轉錄表達沒有明顯差異，ath-miR5658 可能調控了這些靶基因在種子中的轉錄表達。據此我們推測：ath-miR5658既可能調控葉片中光合作用基因的表達[20]，但也會對種子中的甘油磷酸酯代謝、脂肪酸延長和甘油酯代謝產生影響。而gma-miR169v 盡管在葉片中表達，但是，其靶基因半乳糖基二酰基甘油合酶2（comp90506_c0）在葉片中的轉錄表達仍然遠高于種子中的表達。由此可見，它的調控作用局限在一定范圍內。

在前期油茶轉錄組表達分析研究[19] 的基礎上，本研究進行miRNA 靶基因篩選， 初步確定了mRNA 和靶基因的調控關系。下一步研究將是利用基因組測序組裝序列，stem loop qRT-PCR 和靶轉錄表達方法，一是驗證ath-miR5658 及其靶基因類NADP 依賴型蘋果酶（comp89613_c0）、果糖二磷酸醛縮酶1（comp46464_c0）的轉錄表達的相關性，從而確認miRNA 對靶基因的調控程度。二是在油茶果實發育和茶油形成不同階段，檢測gma-miR169v 及其靶基因的表達模式，串聯分析揭示miRNA-mRNA 調控關系。探討油脂合成是否直接或間接地受到多種miRNA 家族的調控。

4 結論

本研究構建了油茶葉和種子的sRNA 文庫，鑒定出54個新miRNA，22 個已知miRNA 家族，包括105 個miRNA 家族成員；67個差異表達miRNA， 調控1346個靶基因， 其中，26個miRNAs 為顯著差異表達，調控644 個靶基因。初步證實了“一個mi RNA 調控多個靶基因，以及多個miRNA 調控一個靶基因的現象”；發現了athmiR5658可能調控碳固定代謝、光合作用、甘油磷酯代謝、甘油酯代謝、脂肪酸延長等途徑中靶基因的表達，gma-miR169v 可能調控甘油酯代謝中靶基因的表達。