上官燕, 李嬌愛, 王 哲, 孫億敬, 劉祉辛, 孫旭武

(上海師范大學 生命與環(huán)境科學學院 植物種質資源開發(fā)協(xié)同創(chuàng)新中心,上海 200234)

0 引 言

microRNA(miRNA)是21～24個核苷酸的RNA非蛋白質編碼基因的產(chǎn)物[1].在植物和動物的基因組中含有大量miRNA編碼基因,它們首先被轉錄成較長的primary microRNAs(pri-miRNA),隨后在RNA酶Ⅲ樣DICER酶的作用下從pri-miRNA中釋放出長度大約為20～22 nt的miRNA.成熟的miRNAs進一步被整合到ARGONAUTE(AGO)效應蛋白中以形成一個RNA-induced silencing complex(RISC),RISC的主要調節(jié)目標是轉錄后的信使RNA(mRNA)[1].miRNA介導的基因表達控制參與調節(jié)了植物生長發(fā)育的諸多方面,包括根和葉的發(fā)育、激素應答、發(fā)育轉變和應激反應[1].

在模式植物擬南芥(Arabidopsisthaliana)中,4種DICER-like(DCL)蛋白之一的DCL1主要參與miRNA的產(chǎn)生[2-4].蛋白DAWDLE(DDL),TOUGH(TGH),HYPONASTIC LEAVES 1(HYL1)和SERRATE(SE)與DCL1相互作用來精確和有效地調節(jié)miRNA的產(chǎn)生[5-12].在miRNA產(chǎn)生期間,長的pri-miRNA轉錄物被進一步切割形成precursor miRNA(pre-miRNA),并從其中釋放出miRNA雙鏈體.然后,miRNA通過HUA增強子1(HEN1)被甲基化,并且主要與擬南芥中10種AGO蛋白之一的AGO1結合[11,13-14].HEAT SHOCK PROTEIN 90(HSP90)和cyclophilin 40蛋白SQUINT(SQN)與AGO1相互作用并確保有效的miRNA加載[15-17].AGO1去除miRNA雙鏈體的passenger strand(或miRNA *),允許miRNA guide strand識別靶RNA內的互補序列區(qū)域[1,18].AGO1與靶序列的結合導致miRNA結合位點內的mRNA切割,或在內質網(wǎng)上進行的翻譯被抑制.

除了miRNA,植物還富含另一種類型的small interfering RNA(siRNA),其與miRNA的結構、基因組和功能均相似.siRNA源自轉基因的模板,以及內源性重復序列或轉座子[19-21].在模板、重復序列及轉座子方面,miRNA與來自轉基因的siRNA的一個關鍵區(qū)別是miRNA的靶基因不是產(chǎn)生miRNA的基因,而siRNA是針對生成它們的序列的.最近,一種新的與miRNA相似的siRNA,被命名為trans-acting siRNA(ta-siRNA)[22-23].ta-siRNAs起源于產(chǎn)生非編碼模板的基因座,它們是miRNA的靶基因[24].由miRNA介導的轉錄本的剪切招募一個依賴RNA的RNA聚合酶 (RdRP),將剪切的轉錄本作為模板,以形成長的雙鏈RNA(dsRNAs),它們是多個ta-siRNAs的來源.miRNA介導的pre-RNA的剪切對ta-siRNAs的形成至關重要[25].

SPETH等[26]利用酵母雙雜交篩選的策略,鑒定到了一個與SE互作的蛋白receptor of activated C kinase 1 (RACK1).進一步分析表明RACK1通過和SE互作,作為AGO1復合物的一部分,參與調節(jié)一些pri-miRNA的積累和加工.rack1突變體導致一些miRNA的靶基因的miRNA的過量富集,并改變了植物對脫落酸(abscisic acid,ABA)的響應和葉序的發(fā)育[26].RACK1存在于幾乎所有的真核生物中,具有7個WD40-b-propellers,其可以與多種蛋白質相互作用[27-29].RACK1本身不具有酶活性,而是充當?shù)鞍踪|之間相互作用的橋梁或與其他蛋白質競爭結合位點,從而抑制特定的相互作用.RACK1通常在蛋白質組學方法中被鑒定,表明RACK1是許多復合物的一部分,其中,它控制復合物的形成和穩(wěn)定性,或者為調節(jié)子創(chuàng)建對接位點[30].直接與RACK1相互作用的幾種類型的蛋白已經(jīng)得到了描述[27].RACK1可以作為真核生物40S核糖體亞基的核心成分,能直接調節(jié)翻譯來響應各種刺激[27,29].在細胞核中,RACK1作為銜接子連接激酶與其底物,并調節(jié)轉錄[31-34].

在動物中研究發(fā)現(xiàn),RACK1在miRNA途徑中起作用.秀麗隱桿線蟲(Caenorhabditiselegans)的AGO蛋白ALG-1直接與核糖體相關的RACK1結合,意味著秀麗隱桿線蟲RACK 1可能是被翻譯的mRNA上的ALG-1的錨點[35].在另一項使用肝腫瘤細胞系的研究中發(fā)現(xiàn),RACK1影響miRNA效應復合物中的加載,并且參與調節(jié)KH-type splicing regulatory protein(KSRP)的定位,KSRP參于促進一小部分pre-miRNA的成熟[36].雖然在秀麗隱桿線蟲中發(fā)現(xiàn)RACK1的敲低會導致更高的miRNA水平,而在人類細胞中的結果則是相反的[35-36].這些看似矛盾的結果和發(fā)現(xiàn)暗示,在許多不同的細胞過程中,RACK1在miRNA途徑中更具體的功能可能被其多效性作用所遮蔽.為了進一步探索RACK1在植物miRNA途徑中的潛在功能,本文作者對美國國家生物技術信息中心(National Center for Biotechnology Information,NCBI)上共享的有關rack1突變體的small RNA-sequence的數(shù)據(jù)進行了系統(tǒng)的生物信息學分析,發(fā)現(xiàn)rack1突變體顯著影響著27,24,及21 nt的small RNA的富集.分析差異small RNA以及其靶基因發(fā)現(xiàn),rack1突變體主要影響了細胞核組分、葉綠體以及線粒體組分相關的基因的表達.

1 材料與方法

從NCBI網(wǎng)站下載small RNA-sequence數(shù)據(jù)包GSE40579(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE40579),根據(jù)規(guī)范的生物信息學分析流程進行數(shù)據(jù)分析.采用Illumina Hiseq 系統(tǒng)測序所得的數(shù)據(jù),稱為 raw reads 或 raw data.使用cutadapt軟件去掉接頭序列,并按長度對序列進行過濾,去掉長度小于15 bp及大于41 bp的序列.使用Fastx_Toolkit軟件,對序列進行Q20質控,保留Q20達到80%及以上的序列.使用NGS QC Toolkit過濾掉含有N堿基的reads,最終得到高質量的clean reads并用于后續(xù)分析.使用Fastx_Toolkit (version 0.0.13)對clean reads中的重復序列進行去冗余,得到uniq reads.使用自行編寫的程序對clean reads進行長度分布統(tǒng)計,對uniq reads的重復序列進行統(tǒng)計.

Clean reads 將和指定物種基因組進行比對,并統(tǒng)計比對率.由于小RNA切膠范圍為 18～35 bp,miRNA長度分布主要集中在21～22 bp,piRNA主要集中在30 bp左右,transcription initiation RNA(tiRNA)以及tRNA-derived fragments(TRFs)主要集中在30～35 bp.Small RNA測序reads中包含:miRNA、Piwi-interacting RNA(piRNA)、tRNA、rRNA、snoRNA等,得到clean reads后,使用bowtie或者blast軟件與上述small RNA 的數(shù)據(jù)比對,進行small RNA的分類注釋.標準流程中會對鑒定出來的miRNA進行差異篩選、靶基因預測,以及功能富集分析.通過分類注釋,沒有注釋的reads即unannotation 序列,將會使用mirdeep2 軟件進行novel miRNA 預測.

2 結 果

2.1 rack1突變體的small RNA sequence數(shù)據(jù)的統(tǒng)計分析

分析rack1突變體和野生型 (WT)的clean reads的長度分布,結果如圖1所示.由圖1可知,rack1突變體中的不同的small RNA的分布發(fā)生了改變.

圖1 clean reads 的長度分布統(tǒng)計.(a) rack1突變體樣本第一組;(b) rack1突變體樣本第二組;(c) WT樣本第一組;(d) WT樣本第二組

Clean reads 的長度分布有助于比較不同處理樣本的小 RNA 情況,一般情況下樣本中的長度分布峰值為21～25 bp,動物中的峰值有時為31～33 bp,可能為piRNA或tRNA.將clean reads依次和Rfam數(shù)據(jù)庫、cDNA 序列、物種重復序列庫、miRBase 數(shù)據(jù)庫等進行對比注釋.首先將clean reads與參考基因組進行比對,統(tǒng)計比對上的參考基因組上的reads數(shù)比例,以及未比對上的參考基因組的比例.使用blastn 軟件,對clean reads序列同Rfam (version 10.0) 數(shù)據(jù)庫進行比對,提取 E-value 小于等于0.01 的結果,注釋出rRNA、snRNA、snoRNA和tRNA等序列.盡可能地發(fā)現(xiàn)并去除其中可能的rRNA、scRNA、snoRNA、snRNA和tRNA,并進行small RNA的去除.最終將這些注釋到Rfam數(shù)據(jù)庫中的序列進行過濾去除,不用于后續(xù)的已知miRNA比對以及新的miRNA預測.使用bowtie軟件,將過濾后的Rfam數(shù)據(jù)庫reads同轉錄本序列進行無錯配比對,能完全匹配上轉錄本序列,并且長度大于26 bp的序列,認為可能是mRNA的降解片段序列,長度小于15 bp的序列不用于后續(xù)miRNA的預測,過濾去除這些序列.最終提取出比對結果中序列長度為15～26 bp的序列,將之與無法比對上轉錄本的序列進行合并,用于后續(xù)的已知miRNA比對,及新的miRNA預測分析.重復上述步驟,使用Bowtie軟件將過濾后的reads與miRBase中該miRNA的成熟體序列進行無錯配比對,比對上的序列被認為是已知的miRNA,并作為miRNA的表達量統(tǒng)計,進行后續(xù)差異分析.

miRNA由前體發(fā)育為成熟體的過程是由Dicer酶切完成的,一般來說,酶切位點的特異性使得miRNA成熟體序列首位堿基對于堿基U具有很強的偏向性.其他的位點也有統(tǒng)計,如:2～4 號位一般缺少U,第10號位偏向于A(第10號位一般是miRNA對其靶基因發(fā)生剪切作用時的剪切位點).對rack1突變體和野生型樣本中的已知miRNA堿基偏好性進行了統(tǒng)計分析(圖2),結果表明rack1突變體和野生型樣本之間的miRNA堿基偏好性沒有顯著差異.通過上述各數(shù)據(jù)庫的比對注釋,最終匯總出各樣本中的small RNA分布情況(圖3).

圖2 各堿基偏好性統(tǒng)計.(a) rack1突變體樣本第一組;(b) rack1突變體樣本第二組;(c) WT樣本第一組;(d) WT樣本第二組(Position 代表在smallRNA中的堿基位置,順序從5′～3′)

圖3 Reads 分類注釋餅狀圖.(a) rack1突變體樣本第一組;(b) rack1突變體樣本第二組; (c) WT樣本第一組;(d) WT樣本第二組

2.2 新miRNA的預測

將small RNA測序數(shù)據(jù)和不同的數(shù)據(jù)庫進行了比對,去除tRNA、snRNA、rRNA 等small RNA,將過濾后的序列和miRBase數(shù)據(jù)庫比對,統(tǒng)計已知miRNA.使用Mirdeep2軟件對未注釋上的序列進行新的miRNA預測.采用RNAfold軟件,對能夠比對上基因組的序列進行二級結構預測,對能夠形成miRNA發(fā)夾前體的序列,歸類為可能的新 miRNA序列.提取出預測的miRNA的mature序列和star序列,同時進行novel miRNA的定量分析.對預測得到的miRNA的結構信息、序列信息,以及樣本歸屬情況進行生物信息學分析,包括:miRNA hairpin發(fā)卡結構(二級結構)的預測結果,clean miRNA reads比對參考基因組的位置,錯配情況,樣本歸屬情況等.定位預測已知miRNA在hairpin 中的具體定位(比對結果)以及每條已知miRNA前體與成熟體的read比對情況及錯配統(tǒng)計.預測到了17個novel miRNA.圖4為所關注的葉綠體的chloroplast_35062以及在rack1突變體中顯著上調的2_12894的結構和錯配信息.其中,sample 列表示 sample 縮寫標簽.reads 列表示比對在該處的 read 數(shù)目,mm 列表示錯配堿基個數(shù)(錯配堿基在比對圖中用大寫字符表示).在比對區(qū)域內參考序列中彩色堿基區(qū)域為測序片段所覆蓋到的區(qū)域.

2.3 rack1突變體的miRNA的表達水平分析

統(tǒng)計了已知miRNA和新預測出的miRNA序列在rack1突變體和野生型植物中的表達水平,如圖5(a)所示,miRNA表達量計算采用5個統(tǒng)計量:最小值、第一四分位數(shù)(25 %)、中位數(shù)(50 %)、第三四分位數(shù)(75 %)和最大值,來粗略地檢測數(shù)據(jù)是否具有對稱性以及分布的分散程度等信息.進一步對在rack1突變體和野生型植物中的miRNA表達水平作了transcript per million(TPM)密度分布分析,如圖5(b)所示,檢測各個樣品的基因表達模式,不同顏色的曲線代表不同的樣本,曲線上點的橫坐標表示對應樣品TPM的對數(shù)值,點的縱坐標表示概率密度.對rack1突變體和野生型植物中的miRNA 表達水平進行了相關性檢驗,結果如圖5(c)所示,驗證了實驗的可靠性和樣本選擇的合理性.對rack1突變體和野生型植物中的miRNA 表達水平作了主成分分析(principal component analysis,PCA),如圖5(d)所示,WT 和RACK 1的miRNA的分布發(fā)生了較大的差異.

圖4 預測得到的新miRNA 結構圖.(a) 為miRDeep2軟件對miRNA前體的打分和比對到成熟序列、環(huán)狀結構、star序列上的reads數(shù)及其前體的二級結構預測圖(紅色為成熟序列,黃色為環(huán)狀結構,紫色為star序列);(b) 顯示比對到本條前體上的reads分布;(c) 包含了成熟序列、環(huán)狀結構、star序列的位置(紫色為miRDeep2軟件預測的star序列,亮藍色為測序reads支持的star序列)

2.4 rack1突變體的miRNA的差異表達分析

對rack1突變體和野生型植物中的miRNA表達水平進行了差異表達分析,在得到差異表達之后,對差異表達miRNA進行靶基因預測及gene ontology (GO)功能顯著性和KEGG pathway顯著性分析.隨后對差異表達miRNA進行非監(jiān)督層次聚類分析,結果如圖6所示.計算多個樣品兩兩之間的距離,構成距離矩陣,合并距離最近的兩類為一新類,計算新類與當前各類的距離,循環(huán)以上步驟直至只有一類為止.用挑選的差異miRNA的表達情況來計算樣品之間的相關性,一般來說,同一類樣品能通過聚類出現(xiàn)在同一個簇中,聚在同一個簇中的 miRNA 可能具有相似的生物學功能.

2.5 rack1突變體的差異 miRNA 靶基因預測及靶基因功能分析

植物miRNA與靶mRNA幾乎完全互補配對,而大多數(shù)動物miRNA與其靶mRNA 不是完全互補配對的.植物miRNA可以與靶mRNA的任何區(qū)域結合(主要是蛋白編碼區(qū)),通過切割靶mRNA或抑制靶mRNA翻譯實現(xiàn)對基因表達的調控,動物miRNA主要針對靶mRNA的3′非編碼區(qū)域(3′UTR),作用機制為抑制翻譯.針對差異表達分析得到的miRNA,使用Target Finder軟件進行植物的靶基因預測.對預測的靶基因進行GO富集分析,如圖7(a)所示.進一步對預測的靶基因作了KEGG pathway功能分析,結果如圖7(b)所示,每個小點對應一個通路,p值越大,表示該通路中基因個數(shù)越多.由圖7可知,RACK1調控的miRNA的靶基因主要參與調節(jié)防御反應、鹽脅迫、鈣離子代謝、葉綠體功能、葉綠素合成,以及核基因表達等.

圖5 miRNA 表達分析.(a) 樣本箱線圖;(b) 樣本密度圖 (density指概率密度);(c) 樣品間相關系數(shù)熱圖;(d) 樣本 PCA 結果

圖6 差異miRNA 熱圖

圖7 預測miRNA 靶基因的功能通路分析.(a) GO 富集 top 10 條目圖;(b) KEGG 富集 top 20 氣泡圖

3 討 論

對rack1突變體和野生型擬南芥的small RNA-sequence數(shù)據(jù)進行了重新分析,并對分析結果進行了深入的整理和挖掘.結果表明:rack1突變體主要影響了27 nt和24 nt的small RNA產(chǎn)生.鑒定到了17個novel miRNA,其中2-12814的表達在rack1突變體中發(fā)生了顯著性差異上調,其靶基因主要是一系列transposable element gene,暗示rack1突變體除了自身直接影響miRNA的表達和其靶基因的表達以外,也可能通過影響transposable element gene的表達來更為廣泛地影響它們所調節(jié)的基因的表達.進一步的分析表明:rack1影響的miR397b和miR5028的靶基因分別是HEMC和HEMB1,它們是參與葉綠體葉綠素合成途徑中上游途徑的關鍵步驟中的酶類,暗示RACK1在調控植物的葉綠體的發(fā)育和功能中扮演了關鍵的角色.