摘要：

解析調控復雜性狀的機制對于推進作物改良至關重要。在此，我們提出了一個全面的水稻（Oryza sativa）調控組圖譜，涵蓋了來自三個代表性品種的23種不同組織的染色質可及性。我們的研究揭示了117,176個獨特的開放染色質區域（OCRs），占水稻基因組的約15%，這一比例顯著高于之前在植物中的報告。通過整合匹配組織的RNA-seq數據，我們自信地預測了59,075個OCR到基因的聯系，其中增強子構成了這些關聯的69.54%，包括許多已知的增強子到基因的聯系。利用這一資源，我們重新評估了全基因組關聯研究的結果，并發現了一個之前未知的功能：OsbZIP06在種子萌發中的作用，我們隨后通過實驗驗證了這一發現。我們優化了深度學習模型以解碼調控語法，實現了對組織特異性染色質可及性的穩健建模。這種方法使得我們能夠從基因組序列中預測跨品種的調控動態，揭示了順式調控分歧和品種間形態差異的遺傳基礎。總的來說，我們的研究為水稻功能基因組學和精準分子育種建立了一個基礎性資源，為調控復雜性狀的機制提供了寶貴的見解。

引言

水稻（Oryza sativa）不僅是世界上最重要的農作物之一，也是研究植物生長和發育的杰出模式物種。在過去的二十年中，人們付出了巨大的努力來理解水稻中重要農藝性狀的遺傳基礎。全基因組關聯研究（GWAS）在這一探索中發揮了關鍵作用，幫助將遺傳變異與表型多樣性聯系起來。這些研究鑒定出大量有望用于性狀改良的候選基因。然而，盡管取得了這些進展，我們對水稻中復雜性狀的調控機制的理解仍然不完整。

基因調控網絡（GRNs）主要由順式調控DNA序列決定，例如啟動子和增強子，這些序列被特定的轉錄因子（TFs）結合。解析這些調控序列中的調控密碼，并將調控序列與目標基因聯系起來，對于重新構建用于作物改良和性狀優化的基因調控網絡至關重要。然而，目前在水稻中對調控組（包含基因組中所有調控元件）的分析仍然受到限制。這些努力往往集中在特定組織上，忽視了整個發育階段和組織的全面景觀。同樣，建立調控區域與其目標基因之間聯系的努力在水稻中也受到限制。

與此同時，許多與水稻農藝性狀相關的功能性遺傳變異位于非編碼調控區域（例如qSH1、DROT1和FZP），這使得它們的解釋變得具有挑戰性，并強調了對調控序列進行系統剖析的必要性。鑒于不同的性狀在不同的發育階段和組織中表現出來，目前由于缺乏跨越各種組織和生長階段的全面表觀基因組圖譜，水稻中非編碼調控變異的系統注釋受到了阻礙。

為了彌合這些差距，我們利用UMI-ATAC-seq方法（一種改良的ATAC-seq協議，由我們實驗室開發）系統地繪制了三種代表性水稻品種整個生命周期中各種組織的染色質可及性圖譜。通過對145個ATAC-seq數據集的分析，我們總共獲得了117,176個獨特的開放染色質區域（OCRs），占水稻基因組的約15%。通過整合匹配組織的RNA-seq數據，我們根據基因表達和鄰近染色質可及性在組織間的相關性，預測了OCRs的潛在目標基因。通過轉錄因子足跡分析，我們推斷了組織或階段特異性的調控網絡，并通過比較秈稻和粳稻亞種之間的調控圖譜，識別了品種多態性/性狀相關的OCRs。值得注意的是，我們的分析揭示了GWAS相關變異在組織特異性OCRs中的偏好，利用這一OCR圖譜，我們能夠識別出209個復雜農藝性狀與非編碼調控變異之間的因果關聯。利用優化的深度學習模型，我們通過對組織特異性染色質可及性的建模和基于序列的跨品種預測，解碼了調控語法。這種建模方法揭示了導致順式調控分歧的關鍵遺傳變化。總體而言，這些數據不僅為植物研究界提供了一個基石資源，還為精準分子育種提供了寶貴的調控變異。

結果

1.繪制水稻染色質可及性的參考圖譜

為了生成水稻（Oryza sativa）中全面的染色質可及性景觀，我們利用了一種改進的ATAC-seq協議（UMI-ATAC-seq，該技術通過引入獨特分子標識符來提高常規ATAC-seq技術的定量和足跡分析的準確性），對水稻整個生命周期中的23種組織/器官進行了染色質可及性分析。這些代表性組織包括愈傷組織、胚根、胚芽、葉片、葉鞘、根、頂端分生組織（AM1/AM2）、休眠芽（DBuds）、莖頂端分生組織（SAM1/SAM2/SAM3）、穗頸節（PNN）、莖、幼穗（Panicle1/Panicle2/Panicle3/Panicle4）、穎殼、內稃、雌蕊、雄蕊和種皮（Seed1/Seed2/Seed3）。實驗在三種代表性水稻品種中進行，分別是日本晴（NIP；粳稻亞種）、明恢63（MH63；秈稻亞種II型）和珍汕97（ZS97；秈稻亞種I型），每個實驗至少包含兩個生物學重復（見圖1a和補充數據1）。總共生成了145個全基因組染色質可及性數據集，平均測序深度約為3070萬條reads。我們應用了ENCODE標準來建立分析流程（見方法部分）。與存儲在ChIP-Hub數據庫中的已發表植物ATAC-seq數據集相比，我們的數據展現出顯著更高的信噪比（見補充圖1）。通過對三種栽培品種的相應參考基因組進行數據分析，我們平均識別出每個實驗有40,676個（范圍從28,991到49,737）可重復的開放染色質區域（OCRs，其不可重復發現率[IRDR]小于0.05）（見圖1b）。正如預期的那樣，所有實驗中識別出的OCRs主要位于轉錄起始位點（TSS）的近端上游區域或遠端基因間區域（見圖1b、f，補充圖2和補充數據2），分別類似于啟動子或增強子。值得注意的是，來自基因內區域的OCRs所占比例相對較小（約為15.7%），而這些OCRs大多來自內含子區域（見補充圖2b）。這些觀察結果表明，絕大多數OCRs源自水稻基因組的非編碼區域。

我們估計大約15%的水稻基因組可以被注釋為OCRs，并且在每個品種中都觀察到了一致的模式（見圖1c），并且這種估計似乎在水稻中已經達到了飽和（見圖1d）。OCRs包含多個轉錄因子結合位點，并且負責調控目標基因的表達。我們收集了56種不同轉錄因子的公開可用的ChIP-seq數據（見補充數據3），并從ChIP-Hub數據庫中預測了458種轉錄因子的DNA基序，并表明OCRs顯著富集了轉錄因子結合位點（見圖1e）。此外，我們發現與側翼基因組區域相比，OCRs在進化上受到高度限制（見圖1g），支持了之前的研究結果，即保守的非編碼序列（CNSs）可以預測植物中的OCRs。

接下來，我們評估了不同品種和組織之間染色質可及性的整體相似性和差異性。我們基于從同一參考基因組（即日本晴）中鑒定出的合并開放染色質區域（OCRs，共117,176個）對所有數據集進行了量化，并利用t分布隨機鄰域嵌入（t-SNE）技術對其全局模式進行了可視化。t-SNE結果的第1和第2維度通常反映了秈稻（MH63和ZS97）和粳稻（NIP）亞種之間的差異，而第2和第3維度主要區分了不同組織類型的獨特簇（圖1h）。例如，NIP的營養組織和生殖組織的染色質可及性模式被分成了不同的簇，而幼穗和愈傷組織的模式則表現出相似性，無論它們來自何種品種。我們進一步根據Jensen-Shannon散度（JSD）指數計算了每個OCR的組織特異性。顯然，遠端OCR的特異性得分顯著高于近端OCR（圖1i和補充圖3a、b），這與之前的發現一致。

總之，水稻中全面的染色質可及性景觀為作物功能基因組學研究提供了一個有價值的資源。

圖1 | 水稻開放染色質景觀的特征分析。
a. 在水稻的三個品種（日本晴、明恢63和珍汕97）中，對各種組織進行了整個生命周期的ATAC-seq和RNA-seq實驗。樣本收集的詳細描述見補充數據1。整個圖中使用了一致的組織顏色編碼。
b. 柱狀圖顯示了在三個水稻品種的每種組織中鑒定出的可重復開放染色質區域（OCRs）的數量。根據OCR峰頂到最近轉錄起始位點（TSS）的距離，將OCRs分為三類：遠端（>1 kb）、近端（<=1 kb）和基因內。沒有SAM3（日本晴）、Seed1（珍汕97）和莖（珍汕97）組織的數據。
c. 在本研究中，被注釋為開放染色質區域（OCRs）的水稻基因組的比例。
d. 每種組織中獨特的OCRs的累積數量，通過排除與OCR超集重疊的OCRs來計算。
e. 密度圖顯示了在日本晴（NIP）中，圍繞OCRs的轉錄因子結合位點（TFBSs）的富集情況。TFBSs是通過56種不同轉錄因子的ChIP-seq數據集（左側）或458種轉錄因子的DNA基序（右側）預測的，這些數據來自ChIP-Hub數據庫。兩側的側翼區域為1 kb。
f. 在三個水稻品種中，OCR峰頂到其最近TSS的距離分布。包括了已發表的水稻（日本晴）開放染色質數據進行比較。基于該分布，使用1 kb的截止值（虛線）來區分近端和遠端調控OCRs。
g. 圍繞日本晴OCRs的保守性PhastCons評分的分布。
h. t-SNE圖顯示了對不同樣本的染色質可及性進行無監督聚類分析的結果。每個點代表一個重復樣本。顏色編碼與a相同。
i. 盒須圖顯示了基因內（n=14239）、近端（n=29524）和遠端（n=57153）OCRs的組織特異性評分的分布（左側），或每個組織的中位數評分。P1=4.01e-39，P2=2.13e-96，P3=1.23e-95。所有P值均通過近端和遠端OCRs的特異性差異的雙側Mann–Whitney U檢驗計算得出。組織顏色編碼與a相同。盒須圖顯示了中位數（水平線）、第二至第三四分位數（盒）和Tukey式須（盒外）。源數據以源數據文件形式提供。

---

2.將開放染色質區域與目標基因聯系起來

為了破譯這些開放染色質區域（OCRs）可能調控哪些基因，我們在每種水稻品種中為研究的組織生成了匹配的RNA-seq數據集（見補充圖3c和補充數據4）。我們采用了一種策略，基于所有樣本中OCR可及性與基因表達之間的相關性分析來預測OCR到基因的聯系（見圖2a；詳細方法見方法部分）。基因可以通過染色質相互作用被多個OCRs（包括啟動子和增強子）調控，這些相互作用被認為發生在拓撲相關域（TADs）內。根據Hi-C數據，水稻基因組中TADs的大小估計在35千堿基對（kb）到45 kb之間，因此我們將分析范圍限制在40 kb內（即從轉錄起始位點（TSS）上游20 kb到下游20 kb），以預測OCRs的目標基因。使用絕對皮爾遜相關系數|R|≥0.4和P < 0.05作為閾值，我們總共獲得了59,075個OCRs（n = 38,437，占所有OCRs的32.8%）和基因（n = 18,781，占注釋基因的48.1%）之間的獨特聯系（見補充圖4a、b和補充數據5）。正如預期的那樣，OCR到基因的聯系更傾向于發生在近端OCRs中，因此基因表達與染色質可及性之間的相關性在近端聯系中更高（見補充圖4c–f）。

OCRs中的遺傳變異可以通過表達數量性狀位點（eQTL）影響基因表達水平。 我們將本研究中鑒定的OCR到基因的聯系與已發表的水稻eQTL數據進行了共定位分析，發現OCR到基因的聯系與eQTL基因對之間存在顯著的重疊（卡方檢驗，P < 1.55e-06）（見補充圖4g）。正如預期的那樣，與eQTL共定位的OCR到基因的聯系的相關系數顯著高于那些沒有共定位的聯系（Wilcoxon檢驗，P = 4.11e-38；見補充圖4h）。我們鑒定出許多已知的調控變異，這些變異影響與農藝性狀相關的基因的表達。例如，qSH1的一個遠端調控區域（約12 kb上游）中的變異調節了其表達動態，從而改變了水稻的落粒性9。相應地，在各種組織中，尤其是頂端分生組織（SAM）中基因表達增加的情況下，這個增強子的可及性與qSH1的表達之間存在正相關（R=0.47，P < 0.013）（見補充圖4i，l）。同樣，OsLG1與上游調控區域緊密相連，這些區域與與穗形性狀相關的強QTL共定位27（見補充圖4j，l）。IPA1的增強子活性與基因表達之間顯示出顯著的正相關（R = 0.84，P < 2.95e-8），在與產量相關的組織中表達增加（見補充圖4k，l），這證實了IPA1在塑造水稻理想株型（IPA）以及提高谷物產量方面的重要作用28。綜上所述，預測的OCR到基因的聯系為水稻農藝性狀發育提供了調控方面的見解，并突出了可用于基因組編輯的重要基因的可靶向OCRs。

代表性組織的全面染色質可及性圖譜為我們提供了揭示組織特異性調控語法的機會。我們利用JSD評分量化了OCR（開放染色質區域）的組織特異性，從而能夠根據上述預測的OCR-基因鏈接，將目標基因從管家基因（例如GAPDH29和OsGOGAT130）區分為組織特異性基因（例如OsYABBY531和OsWRKY4732）（補充圖5和補充數據6）。我們特別關注分析高度組織特異性的OCR（n = 6686，JSD > 0.08的截止值，約占所有OCR的7%），因為它們可能編碼組織特異性調控語法。這些OCR進一步被注釋為啟動子（n = 2322）或增強子（n = 4364），根據它們與轉錄起始位點（TSS）的基因組距離。通過基于OCR-基因鏈接對染色質可及性和目標基因表達進行聯合聚類分析，我們識別出20個不同的OCR簇（圖2b和補充數據7）。每個簇包含200~500個在特定組織中高度激活的OCR-基因鏈接，并且與相應組織的已知生物學特征高度一致（圖2b-d）。例如，簇5（C5）中的稃片和穎片特異性鏈接包含位于GW8基因座的啟動子-增強子相互作用，GW8是一個已知的控制水稻粒重的基因（圖2c）。相應地，GW8在雌蕊、穎片和稃片中高表達。利用C5中的基因進行基因本體（GO）富集分析顯示，“花粉-雌蕊相互作用”和“授粉”等生物學過程被過度代表（圖2d）。同樣，我們在C19中識別出許多在分生組織樣組織（包括幼穗和莖頂端分生組織）中高度且特異性可及的OCR，相關目標基因在功能上顯著富集于“生殖系統發育”、“花發育”和“莖系統發育”（圖2b、d）。值得注意的是，RFL，一個對植物結構和開花時間至關重要的調節因子，是這些目標基因之一（圖2c）。有趣的是，我們觀察到，與組成性OCR（12.3%）相比，組織特異性OCR中有更高比例（28.9%）來自遠端基因間區域。相比之下，大約85%的組成性OCR來自近端啟動子區域（圖2e）。

為了描述可能結合這些組織特異性OCR的轉錄因子（TF），我們使用了GimmeMotifs工具，這是一個多功能工具，可以通過比較多個實驗中的TF結合活性來檢測組織特異性TF結合基序。我們將分析限制在每個組織中特異性測量（SPM）評分最高的前2500個OCR。預測的調控基序在匹配的組織類型中顯示出顯著的組織特異性富集（補充圖6和補充數據8）。我們將重點縮小到每個組織類型中富集度最高的調控因子，并發現許多推斷的鏈接對應于已知的調控關系（圖3a）。例如，OsIDS1在塑造花序結構和建立花分生組織中起著重要作用，在幼穗中表現出相對較高的活性。OsbZIP72在胚芽組織中富集，已被發現通過調節脫落酸（ABA）信號傳導來調節胚芽長度，并促進種子萌發。值得注意的是，種子和雌蕊組織顯示出關鍵調控因子的共富集模式，這些調控因子參與花和種子發育，包括MADS基因家族的MFO1和MADS6342-44（圖3a）。對于每種組織類型，我們進行了系統分析，計算TF家族內調控因子的相對偏好。我們的分析揭示了不同的組織特異性TF結合模式，表明不同組織對特定調控因子有明確的偏好（圖3b）。例如，TCP轉錄因子家族在莖、雄蕊和穗頸節點（PNN）組織中表現出富集偏好。這一觀察結果與TCP基因的已知生物學功能一致，特別是它們在發育組織中調節細胞增殖的作用。

通過足跡分析（footprinting）分析時間序列的ATAC-seq數據，可以幫助識別控制發育進程和轉變的關鍵調控因子，例如先鋒因子。我們從幼穗（這是決定水稻產量的關鍵器官）的四個連續發育階段（<1mm、1–2mm、3–5mm和5–10mm；圖1a）生成了時間序列的開放染色質數據。我們努力識別那些在幼穗發育階段的富集上表現出正相關或負相關的調控基序，使用動態變化的OCR（開放染色質區域，n = 9244；圖3c、補充圖7a和補充數據9）。富集度最高的調控因子主要表現出正相關，表明它們作為轉錄激活因子的功能。相反，一部分因子表現出負相關，暗示它們具有抑制作用。在這方面，DL（編碼OsYABBY49）、OsSPL950和OsSPL1451被鑒定為水稻幼穗發育過程中的代表性正向調控因子（圖3d和補充圖7b）。然而，需要進一步的實驗數據來驗證這些轉錄因子在幼穗發育中的潛在參與。

總體而言，上述結果為研究組織特異性基因調控的候選關鍵調控因子提供了一個寶貴的資源。

---

3.系統定位GWAS變異在組織特異性調控DNA中的位置

全基因組關聯研究（GWAS）已經識別出許多與水稻各種農藝性狀相關的自然變異。為了系統地將GWAS相關變異與上述注釋的調控元件共定位，特別是那些來自非編碼調控區域的元件，我們從最近的全基因組關聯薈萃分析研究以及NGDC GWAS圖譜數據庫中整理了一個全面的水稻GWAS目錄。總體而言，我們收集了4831個顯著（P < 1e?5）且具有代表性（僅考慮主要SNP）的關聯，這些關聯涉及209個不同的數量性狀，這些性狀可以分為七個主要類別：形態特征、生理特征、產量性狀、種子質量、抗性、色澤和其他（圖4a和補充數據10）。簡而言之，這些GWAS SNP主要位于基因間非編碼區域（圖4b和補充圖8a），其中24.5%的SNP要么位于非編碼OCR（開放染色質區域）內（21.1%），要么與鄰近OCR中的SNP處于連鎖不平衡（LD）狀態（3.4%）（圖4c）。此外，與蛋白編碼序列相比，OCR中GWAS SNP的富集度顯著更高（圖4d），這突顯了調控變異在決定表型特征中的關鍵作用。

此外，我們的研究結果表明，含有GWAS SNP的OCR表現出更強的組織特異性（圖4e、f和補充圖8b-d）。例如，一個含有GWAS主要變異vg072467105553（C/T，GWAS P < 9.27e?8）的OCR與穗數顯著相關。該OCR在幼穗組織中表現出高度特異性的可及性，其可及性與GW7基因表達之間存在正向OCR-基因關聯（R = 0.59，P < 9.14e?5；圖4g）。在另一個例子中，GWAS主要變異vg0431427332與葉片寬度顯著相關（P < 1.58e?8），它位于一個分生組織/幼穗特異性的OCR中，正向調控NAL1基因的表達（R = 0.72，P < 1.16e?6）（圖4h）。先前的研究已經表明，NAL1不僅與葉片寬度有關，還與產量有關，并且其表達水平存在自然變異。更多經過驗證的OCR相關關聯的例子在補充圖8e中展示。

組織特異性調控變異解釋農藝性狀關聯

開放染色質區域（OCRs）內的DNA序列變異在通過改變染色質狀態和基因表達模式推動表型創新方面發揮著重要作用，而這種改變通常是組織特異性的。為了研究與農藝性狀相關的遺傳變異與組織特異性OCR之間的關系，我們以組織特異性的方式計算了OCR內遺傳變異的富集情況。結果表明，顯著的GWAS SNP經常富集在與性狀相關的組織的OCR中（圖4f和補充圖8d）。例如，與小穗性狀相關的GWAS變異在分生組織1（SAM1）、雌蕊和幼穗組織特異性的OCR中高度富集。受此觀察結果的啟發，我們使用一種稱為CHEERS的SNP富集方法，對來自不同組織的OCR中的GWAS鑒定的SNP進行了富集分析（補充圖9）。在209個整理的與GWAS相關的性狀中，約78%（209中的163個）的表型性狀在至少一個組織中顯示出GWAS SNP的富集（補充圖10和補充數據11）。

觀察到的與農藝性狀相關的變異在調控元件中的富集高度特異性于組織類型，并且這種關聯在很大程度上與我們對組織功能的當前理解一致（圖5a）。例如，在各種GWAS研究中，與株高相關的調控變異在莖相關組織中富集；而與種子相關的性狀（如種子厚度、寬度、長度、每株癟谷數和每株飽滿粒數）的遺傳關聯在種子、穎片、雌蕊和雄蕊組織特異性的OCR中高度富集（圖5a）。同時，我們發現與根長相關的變異主要富集在根組織中。具體而言，一個顯著的SNP（vg080620195758，P < 3.98e-8）位于OsHAK12基因的根特異性增強子中，該基因已被證明參與根中的鉀離子吸收（補充圖11a）。