https://doi.org/10.1038/s41586-024-08383-z
تاريخ الاستلام: 26 ديسمبر 2023
تم القبول: 8 نوفمبر 2024
نُشر على الإنترنت: 8 يناير 2025
الوصول المفتوح
تحقق من التحديثات

التفاعلات بين المضيفين والميكروبات غير متجانسة لأسباب متعددة. تتكون أنسجة المضيف متعددة الخلايا من أنواع خلايا متنوعة لديها قدرات متميزة للاستجابة للميكروبات، ويمكن أن تشغل الميكروبات موائل موزعة بشكل غير متجانس في المضيف. علاوة على ذلك، قد تحدث التفاعلات الفردية بين الخلايا بشكل غير متزامن. وبالتالي، يمكن أن تتواجد حالات خلوية متنوعة في نسيج واحد. يمكن أن تخفي هذه اللامتجانسة المبادئ الأساسية للتفاعلات الخلوية عندما يتم تحليل المضيفين والميكروبات على مستوى النسيج.

تمت دراسة تفاعلات النبات والجراثيم لفهم الآليات الجزيئية التي تكمن وراء مناعة المضيف وفتك الجراثيم. لقد كشفت تسلسلات RNA أحادية الخلية (scRNA-seq) للبروتوبلاستات الورقية عن استجابات غير متجانسة للنباتات المصابة بالعوامل الممرضة البكتيرية والفطرية الضارة. “. ومع ذلك، فإن فهمنا لتنوع حالات الخلايا محدود إلى حد كبير بنوع خلية محدد (الميزوفيل) ) مصابة بالعوامل الممرضة الفتاكة المثبطة للمناعة. علاوة على ذلك، فإن الديناميات الزمانية والمكانية للاستجابات المناعية للنباتات غير واضحة بسبب الطبيعة ذات الإنتاجية المنخفضة للاختبارات المبلغ عنها بواسطة الجينات المعدلة وراثيًا. . علاوة على ذلك، لدينا فهم محدود للآليات التنظيمية الجينية التي تكمن وراء تنوع حالات الخلايا، مثل الارتباط المحدد لعوامل النسخ (TFs) بالعناصر التنظيمية الجينية الخاصة بحالة الخلية (CREs). يُستخدم اختبار النواة المفردة للوصول إلى الكروماتين القابل للوصول بواسطة الترانسبوز (snATAC-seq) غالبًا لتحديد العناصر التنظيمية الجينية المحتملة في أنواع وحالات الخلايا الفردية. لكن لم يتم تطبيقها لدراسة استجابات المناعة في النباتات. تمثل هذه الفجوات في المعرفة عقبات كبيرة لفهم كيفية عمل نظام المناعة في النباتات ككيان جماعي من تجمعات الخلايا ذات الوظائف المتميزة. .

في هذه الدراسة، هدفنا إلى سد هذه الفجوات من خلال تحليلات متعددة الأوميات على مستوى النواة الواحدة (snMultiome) والترانسكريبتوميات المكانية لنبات مضيف مصاب بعوامل مرضية بكتيرية ضارة أو غير ضارة في تجربة زمنية. على وجه التحديد، نستخدم تسلسل RNA على مستوى النواة الواحدة (snRNA-seq)، وتسلسل ATAC على مستوى النواة الواحدة (snATAC-seq) والتهجين الفلوري المتعدد الخطأ القوي في الموقع (MERFISH). والأمراض البكتيرية الثلاثة التالية: Pseudomonas syringae pv. tomato DC3000 (المشار إليه فيما بعد بـ DC3000)، DC3000 AvrRpt2 و DC3000 AvrRpm1 (المشار إليهما فيما بعد بـ AvrRpt2 و AvrRpm1، على التوالي). DC3000 هو ممرض ضار يمكنه قمع مناعة النبات من خلال المؤثرات والسموم. بينما تعتبر AvrRpt2 و AvrRpm1 مسببات أمراض غير ضارة تحمل عوامل مؤثرة يتم التعرف عليها بشكل غير مباشر بواسطة مستقبلات النباتات التي تحتوي على نطاقات ربط النوكليوتيدات وكرات ليوسين الغنية لبدء المناعة المستحثة بالعوامل المؤثرة (ETI). يظهر أطلسنا متعدد الأبعاد للخلايا المفردة حالات خلوية غير موصوفة سابقًا، بما في ذلك حالة خلوية نادرة تظهر في مركز مناطق أنسجة ETI (المعينة بخلايا PRIMER). يتم وصف حالة خلوية أخرى تحيط بخلايا PRIMER (المعينة بخلايا المتفرج)، ونقدم رؤى حول وظائف وآليات تنظيم الجينات لهذه الحالات الخلوية.

قمنا بإنشاء أطلس زمني مفصل لـ A. thaliana من أوراق مصابة بـ DC3000 أو AvrRpt2 أو AvrRpm1 (الشكل 1a). تم تلقيح العوامل الممرضة في الأوراق بجرعة منخفضة (الكثافة الضوئية عند 600 نانومتر )، حيث من المتوقع أن يكون فقط مجموعة فرعية من خلايا النبات
مختبر بيولوجيا النباتات، معهد سالك للدراسات البيولوجية، لا جولا، كاليفورنيا، الولايات المتحدة الأمريكية. مختبر التحليل الجينومي، معهد سالك للدراسات البيولوجية، لا جولا، كاليفورنيا، الولايات المتحدة الأمريكية. معهد هوارد هيوز الطبي، معهد سالك للدراسات البيولوجية، لا جولا، كاليفورنيا، الولايات المتحدة الأمريكية. برنامج المعلوماتية الحيوية وبيولوجيا النظم، جامعة كاليفورنيا، سان دييغو، لا جولا، كاليفورنيا، الولايات المتحدة الأمريكية. مختبر علم أمراض النبات، كلية الدراسات العليا للزراعة، جامعة كيوتو، كيوتو، اليابان. العنوان الحالي: مختبر ساينسبري، جامعة شرق أنجليا، نورويتش، المملكة المتحدة. العنوان الحالي: معهد آرك، بالو ألتو، كاليفورنيا، الولايات المتحدة الأمريكية. البريد الإلكتروني: ecker@salk.edu

تعبير زائف للكتل الفرعية. يتم عرض الجينات المتعلقة بالدفاع التي تم توصيفها بشكل جيد. انظر الشكل 3c في البيانات الموسعة لتحليل شامل. تشير الأشرطة العلوية إلى نوع الخلية والكتلة الرئيسية التي تم اشتقاق كل كتل فرعية منها. يتطابق نظام الألوان للكتل الرئيسية مع . e، تقسيم فرعي للعناقيد الرئيسية 6 (جينات مرتبطة بالدفاع) تظهر تعبيرًا محددًا لكل تقسيم فرعي. f، مخطط تقسيم فرعي للعناقيد 3 و7 و11 (خلايا الميزوفيل النشطة مناعياً). g، تعبير الجينات المشاركة في خطوات مختلفة من تخليق وإفراز المستقلبات الثانوية المشتقة من التريبتوفان موضحة في h.h. مخطط مبسط لتخليق وإفراز المستقلبات الثانوية المشتقة من التريبتوفان.
تواجه خلايا مسببات الأمراض. تم أخذ عينات من الأوراق المصابة في أربع نقاط زمنية مختلفة و 24 ساعة كتحكم، قمنا بإعداد نباتات تم إصابتها بشكل وهمي بالماء وأخذنا عينات بعد 9 ساعات. تم إجراء نسختين مستقلتين لشرط الإصابة بـ AvrRpt2 لمدة 9 ساعات وللشرط الوهمي. قمنا بتطوير بروتوكول لعزل النوى بسرعة بحيث يتم تقليل التغيرات في النسخ الجيني والإيبيجينومي أثناء إعداد العينة (الطرق). اجتاز ما مجموعه 65,061 خلية من 15 عينة مراقبة الجودة لكل من تحليلات بيانات RNA-seq وATAC-seq (الشكل البياني الممتد 1a-d). حددت بيانات snATAC-seq لدينا مناطق الكروماتين الأكثر وصولاً (ACRs) مقارنةً ببيانات ATAC-seq الكمية المبلغ عنها سابقًا من الأوراق المنشطة مناعياً. علاوة على ذلك، أظهرت مجموعات البيانات تداخلات كبيرة، مما يشير إلى أن snATAC-seq يمكن أن يلتقط كل من ACRs المعروفة وغير المعروفة (الشكل 1e من البيانات الموسعة). تم الكشف عن قمم ATAC-seq المرتبطة بجين housekeeping (ACTIN2) وجين مرتبط بالمناعة (ICS1) بشكل متسق عبر النسخ المكررة، مما قدم دعماً للازدواجية العالية لبيانات snATAC-seq (الشكل 1f من البيانات الموسعة). أظهرت نسختان مستقلتان إعادة برمجة نسخية متسقة ناجمة عن عدوى AvrRpt2 (الشكل 1g، h من البيانات الموسعة)، مما أكد مرة أخرى على تكرار بياناتنا.

تم إجراء تقليل الأبعاد والتجميع باستخدام بيانات تسلسل RNA الأحادي. كانت بعض المجموعات غنية بعدوى محددة.
الشروط وفي نقاط زمنية محددة (الشكل البياني الموسع 1i)، مما يشير إلى أن تحليل التجميع التقط حالات خلوية متميزة ناتجة عن عدوى مسببات الأمراض. حددنا الجينات المعبر عنها بشكل خاص في التجمعات الفردية (تم عرض العلامات الرئيسية في الشكل البياني الموسع 2a)، مما أوضح بشكل أكبر هوية كل تجمع (تم عرض التوصيفات الرئيسية لأنواع الخلايا في الشكل 1b). تم أيضًا توقع أنواع الخلايا بناءً على بيانات ATAC-seq. على سبيل المثال، كان هناك ACR محدد للتجمعات (ذروة في الكروموسوم 2 في الموقع 11172821-11173529) للتجمعات 0، 12، 19، 21 و29 مرتبطًا بـ جين علامة للجلد (الشكل 2b، c من البيانات الموسعة؛ الأسماء الكاملة للجينات المميزة في هذه الدراسة موضحة في الطرق). بشكل عام، صنفت بيانات snMultiome أنواع الخلايا وحالات أوراق A. thaliana المصابة بعامل ممرض.

تحليل إثراء علم الأحياء الجيني (GO) لجينات العلامة لكل مجموعة حددت تجمعات خلايا الميزوفيل (المجموعات 3 و 7 و 11) والجلد (المجموعة 12) الغنية بجينات مرتبطة بالدفاع، بما في ذلك تلك المشاركة في مسار هرمون الدفاع حمض الساليسيليك (SA). (الشكل 1c والشكل الإضافي 2d). كانت هذه المجموعات ممثلة بشكل جيد في ظروف ETI (عدوى AvrRpt2 وAvrRpm1) (الشكل الإضافي 1i)، مما يشير إلى أن هذه الخلايا كانت تستجيب للعوامل الممرضة المحفزة للمناعة. وقد تم دعم هذه النتيجة أيضًا من خلال التعبير القوي عن ICS1، وهو جين رئيسي لتحفيز SA بواسطة العوامل الممرضة.
التمثيل الحيوي (الشكل البياني الموسع 2e). لاحظنا زيادة في الوصول إلى مناطق الكروماتين upstream من موضع ICS1 بعد العدوى بسلالات ETI (الشكل البياني الموسع 1f)، وهو نتيجة تتماشى مع دراسة سابقة باستخدام ATAC-seq الشامل. معًا، قامت بياناتنا بتوثيق تغييرات غير متجانسة ومنسقة في تعبير جينات الدفاع والوصول إلى الكروماتين خلال عدوى العوامل الممرضة.

على الرغم من أن الكتل الرئيسية التقطت خلايا نشطة مناعياً في الميزوفيل والجلد (الشكل 1c)، إلا أن الكتل لأنواع خلايا أخرى، مثل الأوعية الدموية، احتوت على خلايا نشطة وغير نشطة مناعياً. من المحتمل أن تكون هذه النتيجة ناتجة عن توقيعات تطويرية قوية (على سبيل المثال، يتم التعبير عن جينات علامات الأوعية الدموية بمستويات عالية بغض النظر عن تنشيط المناعة). لالتقاط استجابات مناعية محددة لنوع الخلية، قمنا بإجراء جولة ثانية من التجميع لكل كتلة رئيسية، مما أسفر عن 429 تحت كتلة بأنماط تعبير جيني متنوعة (الشكل 1b (يمين) والشكل الممتد 3a). أظهرت تحليلات الجينات المختارة المتعلقة بالمناعة تعبيراً جينياً محدداً لنوع الخلية وتنوعاً في التعبير بين أنواع الخلايا (أي، تنوع حالة الخلية) (الشكل 1d). كشفت عملية التجميع الفرعي لبيانات snRNA-seq لدينا عن استجابات مناعية معقدة في أنسجة الأوراق لم يتم التقاطها بواسطة RNAseq الكلي (الشكل الممتد 3h). علاوة على ذلك، في بعض الحالات، بدت الجينات وكأنها معبّرة بشكل عالٍ على مستوى النسخ الكلي ولكنها معبّرة بشكل محدد على مستوى تحت الكتلة (الشكل الممتد 3i، j)، وهو اكتشاف يبرز المزيد من قيمة تحليلات الخلايا الفردية.

أظهر تحليل أكثر تفصيلاً لمجموعات فرعية محددة تعبيراً قوياً عن علامات خلايا رفيقة اللحاء (PCCs) (الشكل التوضيحي الممتد 3d) في المجموعة الرئيسية 6. يمكن فصل هذه المجموعة إلى 12 مجموعة فرعية، بعضها كان غنياً بالجينات المرتبطة بالمناعة (الشكل 1e). كانت الجينات المهمة للمقاومة المكتسبة نظامياً (SAR)، مثل ALD1 و FMO1، معبراً عنها بشكل خاص في المجموعة الفرعية 8 من PCC (الشكل 1e). تشير هذه النتيجة إلى أن مجموعة فرعية من PCCs تساهم في إرسال إشارات بعيدة المدى إلى الأوراق النظامية. ومن الجدير بالذكر أننا لم نلاحظ تعبيراً قوياً عن ALD1 أو FMO1 في مجموعات الأوعية الأخرى (الشكل التوضيحي الممتد 3e). كان ILL6 من بين الجينات العلامة الغنية في المجموعة الفرعية 8 من PCC (الشكل 1e والشكل التوضيحي الممتد 3f)، وهذا الجين مرتبط بـ SAR. تشير هذه النتائج معًا إلى أن هذه المجموعة الخلوية قد تلعب دورًا في SAR، وقد يكون ILL6 منظمًا محددًا لـ SAR في PCC. لقد حددنا جينات إضافية غنية بشكل خاص في تحت مجموعة PCC 8 ولكن ليس في المجموعات الرئيسية الأخرى (الشكل 3g من البيانات الموسعة). نظرًا لأن SAR يتطلب إشارة متحركة تنتقل من الأوراق المصابة محليًا إلى الأوراق النظامية عبر الأوعية الدموية، فمن الممكن أن تكون مجموعة PCC المحددة والجينات المؤشر لها لها دور محدد في SAR.

كمثال آخر على الاستجابات المناعية المحددة لسكان الخلايا، قمنا بتحليل مسارات الأيض الثانوي المتعلقة بالدفاع المشتقة من التربتوفان في خلايا الميزوفيل النشطة مناعياً (المجموعات 3 و7 و11) (الشكل 1f). تم تحديد تعبير مميز للجينات المشاركة في خطوات أو مسارات مختلفة من تخليق أو إفراز المستقلبات الدفاعية الكاماليكسين والإندول جلوكوسينولات (الشكل 1g وh). تشير هذه النتيجة إلى وجود تقسيم في تفعيل مسارات الدفاع التي يمكن أن تتنافس على الموارد (التربتوفان). معاً، تبرز هذه الأمثلة تنوع حالات خلايا المناعة في النباتات، مما يؤكد أهمية فهم مسارات الإشارات المناعية والشبكات على مستوى حالة الخلية.

مكنتنا بيانات snMultiome من مقارنة تعبير mRNA و ACRs بشكل مباشر لتحديد CREs المحتملة في أنواع خلايا مختلفة، وظروف العدوى ونقاط زمنية مختلفة (الشكل 4a من البيانات الموسعة). حددنا ما مجموعه 29,002 زوج من ACR-جين المرتبطين بشكل كبير (أو المرتبطين) ضمن 500 كيلوبايت من كل جين عبر الخلايا في كل حالة عدوى. كانت معظم الروابط ضمن مسافة قصيرة من مواقع الجينات (<400 نقطة أساسية؛ مناطق المحفز المحتملة)، بينما كانت أخرى أكثر بعدًا (مناطق المعزز المحتملة).
(البيانات الموسعة الأشكال 4ب و5أ). قمنا بتلخيص بيانات ارتباط القمة بالجين لكل جين باستخدام قيم معامل الارتباط بيرسون الأقصى (درجة ارتباط القمة بالجين) (البيانات الموسعة الشكل 4ج). كانت الجينات التي أظهرت روابط في جميع الظروف (العنقود 8؛ البيانات الموسعة الشكل 4ج) غنية بجينات علامات نوع الخلية مثل FDH وBCA2 وMAM1 (البيانات الموسعة الشكل 4ج، د). كانت الجينات التي أظهرت روابط بشكل خاص في ظروف تنشيط ETI (العنقود 4؛ البيانات الموسعة الشكل 4ج) غنية بجينات متعلقة بالمناعة (البيانات الموسعة الأشكال 4د و5ب). كان CBP60G، وهو منظم نسخي للمناعة، لديه عدة ACRs حيث كانت الوصولية مرتبطة بشكل كبير بتعبير mRNA الخاص به (البيانات الموسعة الشكل 4هـ). كانت الجينات التي أظهرت روابط بشكل خاص في عدوى DC3000 (العنقود 2؛ البيانات الموسعة الشكل 4ج) غنية بجينات متعلقة بحمض الجاسمون (JA) (البيانات الموسعة الشكل 5ب). هذه النتيجة تتماشى مع قدرة DC3000 على تنشيط مسار JA في النباتات باستخدام السم coronatine والعوامل المؤثرة لقمع مناعة النبات. تشير هذه النتائج إلى أن إعادة برمجة الوصول إلى الكروماتين والتعبير الجيني بشكل منسق ومحدد لمجموعات الخلايا هي سمة رئيسية من سمات تطور الأوراق والمناعة واستغلالها من قبل مسببات الأمراض الضارة.

على الرغم من أننا حددنا العديد من ACRs المرتبطة ارتباطًا وثيقًا بجينات الدفاع، ( 304 من 627 ) من جينات الدفاع (جينات العلامة لمجموعات النشاط المناعي و 12) لم تظهر روابط مهمة مع ACRs. غالبًا ما كانت هذه الجينات الدفاعية غير المرتبطة تحتوي على كروماتين مفتوح بشكل دائم (الشكل البياني الممتد 4f) على الرغم من التغيرات المحددة في التعبير الجيني (الشكل البياني الممتد 4f، مخطط الكمان). تشير هذه النتيجة إلى وجود طبقة إضافية من تنظيم الجينات، على سبيل المثال، تعبير TFs العلوية. كانت أنماط مختلفة غنية في ACRs على بعد 2 كيلو بايت من الجينات الدفاعية المرتبطة وغير المرتبطة (الشكل البياني الممتد 4g، h)، مما يشير إلى أن بعض TFs الدفاعية تعمل معًا مع إعادة برمجة الكروماتين بينما لا تفعل ذلك الأخرى.

لتحديد آليات تنظيم الجينات المحددة لنوع الخلية والحالة، قمنا بإجراء تحليل إثراء الأنماط للـ ACRs المرتبطة بالمجموعات الفردية. على سبيل المثال، مقارنة مجموعة الميزوفيل النشطة مناعياً (المجموعة 3) ومجموعة الميزوفيل غير النشطة مناعياً (المجموعة 1) كشفت عن إثراء الأنماط للعديد من عوامل النسخ المعروفة بمشاركتها في المناعة، بما في ذلك WRKYs وCAMTAs (الشكل 2a). تسلط هذه النتيجة الضوء على فائدة هذه الاستراتيجية. قمنا بتوسيع هذا التحليل ليشمل ACRs العلامة لجميع المجموعات وحددنا إثراء الأنماط المحددة لكل مجموعة (الشكل 2b). تشير هذه النتائج إلى أن أنواع الخلايا والحالات المختلفة تستخدم كل من تنظيم الجينات المشترك والمتميز من خلال ارتباط TFDNA. علاوة على ذلك، يمكن أن تسرع هذه الطريقة من تحديد عوامل النسخ ذات الوظيفة المحددة لنوع الخلية أو الحالة.

لتحليل وصول الأنماط على مستوى الخلية الواحدة – وليس على مستوى التجمع – قمنا بحساب درجة إثراء الأنماط لكل خلية باستخدام ChromVAR. (الشكل 2ج)، الذي كشف عن وصول غير متجانس إلى 465 نمطًا بين المجموعات وداخلها (البيانات الموسعة الشكل 6أ). ثم تساءلنا عما إذا كان هناك ارتباط بين درجة إثراء نمط TF وعبارة TF، وأظهرت النتائج أيضًا تعبيرًا غير متجانس (البيانات الموسعة الشكل 6ب). حددنا TFs التي أظهرت أعلى ارتباطات بين درجات إثراء الأنماط وعبارات mRNA الخاصة بها في أنواع خلايا مختلفة (الشكل 2د والبيانات الموسعة الشكل 6ج). على سبيل المثال، تداخل تعبير mRNA لـ WRKY46 والوصول إلى مواقع ارتباط WRKY46 بشكل رئيسي في خلايا الميزوفيل النشطة مناعياً وخلايا البشرة عبر نقاط زمنية مختلفة (الشكل 2هـ). تشير هذه النتيجة إلى أن نظام WRKY46 له دور رئيسي خلال الاستجابات المناعية في هذه الأنواع من الخلايا. من خلال البحث عن الجينات المرتبطة بـ ACRs التي تحتوي على مواقع ارتباط WRKY46 (الشكل 2و)، حددنا جينات مستهدفة محتملة لـ WRKY46، بما في ذلك العديد من الجينات المعروفة المتعلقة بالدفاع (الشكل 2ز). تتماشى هذه النتيجة مع الوظيفة المعروفة لـ WRKY46 في مسار SA والدفاع ضد . سيرينغاي . أخيرًا، قمنا بإجراء توقع للجينات المستهدفة لجميع عوامل النسخ الموضحة في الشكل 2d ووجدنا أن العديد من عوامل النسخ تستهدف جينات مرتبطة بدفاعات النبات (الشكل 2h). كانت عوامل النسخ التي تنتمي إلى نفس العائلة تميل إلى إظهار جينات مستهدفة متداخلة (الشكل 2h). بالإضافة إلى WRKYs، الأنماط

الشكل 3 | الترانسكتوميات المكانية الزمنية في الأوراق المصابة بالعوامل الممرضة. أ، مخطط لتجارب MERFISH. لكل حالة. انظر الطرق للإجراءات التفصيلية. ب، مخططات ثنائية الأبعاد للنسخ الجينية لـ 16 جينًا مختارًا (الجدول التكميلي 1) تم اكتشافها باستخدام MERFISH في كل عينة. تم توفير المخططات التي تحتوي على نسخ لجميع الجينات الـ 500 في الشكل البياني الممتد 7ب. ج، يسار، مجال رؤية تمثيلي يظهر إشارة غامضة لنوى ملونة بـ DAPI. في المنتصف، تقسيم قائم على النسخ في نفس مجال الرؤية (الطرق). تم تلوين النسخ بناءً على الخلايا المعينة. يمين، مراكز الخلايا المكتشفة باستخدام طريقة تقسيم قائمة على النسخ (نقاط حمراء) و
نتيجة الفشل في تقسيم الصبغة DAPI (المنطقة الصفراء). لوحظت أنماط مشابهة عبر مجالات الرؤية والعينات. تم تقديم تحليل كمي منهجي في الشكل 7f من البيانات الموسعة. د، مخططات الكمان تظهر عدد النسخ لكل خلية (يسار) والجينات الفريدة لكل خلية (وسط) المكتشفة في كل عينة MERFISH وحجم المنطقة لكل خلية (يمين). هـ، تمثيلات UMAP للخلايا في كل عينة بناءً على تعبير 500 جين تم اكتشافه باستخدام MERFISH. تم دمج جميع عينات MERFISH، وتم تلوين الخلايا بناءً على مجموعات Leiden الجديدة. و، رسم خرائطي مكاني لمجموعات Leiden في كل عينة باستخدام نفس نظام الألوان كما في هـ. مقياس الرسم، (ج) أو 1 مم (ب، ف).

الشكل التوضيحي الممتد 7d). بعد التقسيم، تم تخصيص النسخ إلى الخلايا، مما أنتج مصفوفة خلية مقابل جين مع وجود إحداثيات مكانية لكل خلية. بشكل عام، اكتشفنا وسطي 161 نسخة و79 جينًا لكل خلية من إجمالي 121,998 خلية من 5 مقاطع ورقية (1 عدوى وهمية و4 نقاط زمنية بعد عدوى AvrRpt2) (الشكل 3d). قمنا بإجراء دمج وتجمع جديد لجميع عينات MERFISH، مما حدد 14 مجموعة (الشكل 3e والشكل التوضيحي الممتد 7g). تم رسم هذه المجموعات مكانيًا في عينات فردية، مما كشف عن تجمعات خلوية منظمة مكانيًا (الشكل 3e). على سبيل المثال، تم التقاط مجموعة MERFISH جديدة خلايا الأوعية الدموية في مقاطع الأنسجة (الشكل التوضيحي الممتد 7h). توفر خصائص التعبير الجيني على مستوى الخلية وكمية MERFISH الفرصة لدمج بيانات MERFISH الموزعة مكانيًا وبيانات snMultiome الغنية بالمعلومات الجزيئية.

قمنا بدمج بيانات MERFISH و snMultiome (خمسة ظروف تتطابق مع عينات MERFISH؛ الشكل 4a) باستخدام 500 جين مشترك، وتم نقل تسميات الكتل المحددة بواسطة snMultiome إلى خلايا MERFISH (الشكل 4a وb والطرق). استنادًا إلى هذا الدمج البياني، قمنا برسم خرائط مكانية للكتل المحددة في بيانات snMultiome. تم رسم الأنواع الرئيسية من الخلايا المحددة بواسطة snMultiome بنجاح على المناطق المتوقعة في عينة نسيج MERFISH (الشكل 4c)، مما يدل على نجاح دمج البيانات. هذا الدمج بين بيانات MERFISH و snMultiome مكننا من استكشاف التوزيع المكاني لمجموعات الخلايا المحددة في تحليل snMultiome.

باستخدام مجموعة البيانات المتكاملة هذه، قمنا بتقدير المعلومات الخاصة بالترانسكريبتوم والوصول إلى الكروماتين مكانيًا (الشكل البياني الممتد 8a-d). تم تقدير تعبير ICS1 (غير المدرج في لوحة MERFISH) بدقة، مما توقع التعبير المكاني الحقيقي لـ ICS1 (استنادًا إلى smFISH) (الشكل البياني الممتد 8a,b)، وهو نتيجة أكدت دقة تقدير البيانات. كما قمنا بتقدير درجات نشاط ATAC مكانيًا (الشكل البياني الممتد 5f) لـ ICS1 و ALD1، وأظهرت هذه أنماطًا متسقة مع تعبير mRNA (الشكل البياني الممتد 8c). تم توقع نشاط نمط HSFB2b ليكون مرتفعًا في المناطق النشطة مناعيًا في ورقة عند 24 ساعة بعد التلقيح (h.p.i.) (الشكل البياني الممتد 8d)، وهو ما يتماشى مع نمط تعبير mRNA لـ تم التحقق منه بواسطة MERFISH (الشكل الإضافي 8d). بشكل عام، تشير هذه النتائج إلى أن تقدير البيانات المكانية للترانسكريبتوم والإبيجينوم كان دقيقًا، مما مكن من تحليل آليات تنظيم الجينات بدقة خلوية ومكانية. لتسهيل استكشاف بياناتنا، قمنا بتقدير جميع 25,299 ترانسكريبت تم اكتشافه في تحليل snMultiome و465 درجة إثراء نمطية على الأنسجة الخمسة المستخدمة في تجارب MERFISH وجعلنا البيانات متاحة على متصفح البيانات (https:// plantpathogenatlas.salk.edu).

لفهم الديناميات الزمنية للاستجابات المناعية، قمنا بتطبيق تحليل الزمن الزائف على بيانات snMultiome الخاصة بنا. يقوم تحليل الزمن الزائف بمحاذاة الخلايا كمسار استنادًا إلى تعبيرها الجيني، وهو ما يُستخدم عادة لاستنتاج المسارات التطورية. . اقترحنا

الشكل 4 | دمج بيانات snMultiome وبيانات النسخ الجزيئي المكاني. أ، دمج بيانات snMultiome (المصابة بشكل وهمي والمصابة بـ AvrRpt2؛ خمس عينات في المجموع) وبيانات MERFISH الموضحة في UMAP. ب، UMAP المدمج. اليسار، يتم تلوين النوى والخلايا بناءً على تسميات المجموعات، التي تم نقلها من snRNA-seq إلى MERFISH (الطرق). اليمين، نفس UMAP المدمج ملون حسب أنواع الاختبارات. ج، اليسار، UMAP لبيانات snMultiome ملون على
أساس أنواع الخلايا الرئيسية: البشرة (أرجواني)، الميزوفيل (أخضر) والأنسجة الوعائية (أصفر). على اليمين، رسم خرائطي مكاني لخلايا snMultiome ملونة بناءً على أنواع الخلايا الرئيسية. تم استخدام عينة AvrRpt2 بعد 9 ساعات من الإصابة. مقياس الرسم، 1 مم. د، قيم الزمن الزائف المحسوبة لخلايا الميزوفيل في بيانات snRNA-seq. هـ، مخططات UMAP تظهر قيم الزمن الزائف في خلايا من كل نقطة زمنية. و، رسم خرائطي مكاني لقيم الزمن الزائف بناءً على دمج البيانات ونقل العلامات.
من خلال تطبيق تحليل الزمن الزائف على نوع خلية تطورية واحد مع حالات مناعية متنوعة، تمكنا من نمذجة الديناميات الزمنية للاستجابة المناعية والعدوى المتنوعة في ورقة مصابة بشكل أفضل. أظهرت حسابات درجات الزمن الزائف لخلايا الميزوفيل (الشكل 4d والطرق) أن توزيع درجات الزمن الزائف المتوقعة في كل عينة كان متسقًا مع ما هو متوقع من نقطة زمنية حقيقية (الشكل 4e). أي أن الخلايا ذات درجات الزمن الزائف المنخفضة كانت متزايدة في النقاط الزمنية المبكرة، بينما ظهرت الخلايا ذات درجات الزمن الزائف العالية في النقاط الزمنية اللاحقة. تشير هذه النتيجة إلى أن الديناميات الزمنية للاستجابات المناعية تم نمذجتها بنجاح. كانت هناك خلايا تتراوح درجات الزمن الزائف فيها بشكل واسع تتواجد في 9 و 24 ساعة بعد الإصابة، مما يشير إلى أن الاستجابات المناعية الخلوية غير متزامنة في الأوراق المصابة.

لفهم التوزيع المكاني للحالات المناعية غير المتجانسة، قمنا برسم خرائط زمنية زائفة مكانيًا باستخدام نقل العلامات من بيانات snMultiome إلى بيانات MERFISH (الشكل 4f). كشفت الخرائط المكانية لدرجات الزمن الزائف عن مناطق نشطة مناعيًا كانت موزعة في أنسجة الأوراق المصابة بالعوامل الممرضة بطريقة محدودة (الشكل 4f). توسعت هذه المناطق النشطة مناعيًا مع مرور الوقت وبدت وكأنها تندمج عند 24 ساعة بعد الإصابة (h.p.i) (الشكل 4f). من الجدير بالذكر أن الحالات المناعية تتغير أيضًا ديناميكيًا مع مرور الوقت في كل منطقة نشطة مناعيًا، حيث تكون الخلايا النشطة مناعيًا الأقدم (درجات زمن زائف أعلى) محاطة بخلايا نشطة مناعيًا أصغر سنًا (درجات زمن زائف أقل) (الشكل 4f). تشير هذه النتائج إلى أن أقدم الخلايا النشطة مناعيًا هي على الأرجح خلايا نباتية كانت على اتصال مباشر مع خلايا العوامل الممرضة في نقاط زمنية مبكرة وأن الاستجابات المناعية تنتشر إلى الخلايا المحيطة من خلال التواصل بين الخلايا مع مرور الوقت.

سعينا للتحقيق فيما إذا كان يمكن تفسير التعبير الجيني المناعي المحدود مكانيًا من خلال توزيع الخلايا البكتيرية. وقد اكتشفت تقنية smFISH التي تستهدف الجينات الميتاجينية البكتيرية (19 جينًا عالي التعبير) مستعمرات بكتيرية بعد 24 ساعة من الإصابة، والتي قمنا بتراكبها مع الخريطة المكانية لدرجات الزمن الزائف (الشكل 8g من البيانات الموسعة). كتحكم،
قمنا بإجراء تجربة MERFISH أخرى باستخدام ورقة من A. thaliana مصابة بالعامل الممرض المثبط للمناعة DC3000 بعد 24 ساعة من الإصابة. تداخل توزيع السلالة المحفزة للمناعة AvrRpt2 مع مناطق الأنسجة ذات درجات الزمن الزائفة العالية (مُعززة للمناعة) على عكس سلالة DC3000 المثبطة للمناعة. أكدنا هذه الملاحظة من خلال تحليل كمي لخلايا النبات المجاورة لمستعمرات البكتيريا الفردية. مجتمعة، تشير هذه النتائج إلى أن المناطق النشطة مناعياً التي حددها تحليل الزمن الزائف تتفاعل مع العامل الممرض المحفز للاستجابة المناعية؛ كما قمنا بالتقاط احتمالية تثبيط المناعة بواسطة العامل الممرض الفتاك DC3000.

قمنا بتحديد الجينات التي تغير تعبيرها بشكل ملحوظ عبر مسار الزمن الزائف في بيانات MERFISH باستخدام استنتاج التعبير التفاضلي المحدد لنوع الخلية (C-SIDE) بشكل منهجي. (الشكل 5أ والطرق). كانت BON3 وALD1 وFMO1 من بين الجينات التي تم تحديدها والتي كانت تعبيرها أعلى نحو مركز المناطق النشطة مناعياً (الشكل 5أ). تم تحفيز BON3 في مناطق محدودة للغاية في الأنسجة بعد الإصابة بـ AvrRpt2 (عند 9 ساعات بعد الإصابة) (الشكل 5ب)، وهو ما يختلف عن ALD1 وFMO1 (الشكل 5ب-د والشكل التمديدي 9أ). أكدت تحليلنا الفرعي لبيانات snMultiome النمط الذي تم ملاحظته مع نتائج MERFISH (الشكل 5هـ)، مما يشير إلى أننا حددنا حالتين متميزتين من خلايا المناعة. كان تعبير BON3 متزايداً في الخلايا ذات أعلى درجات الزمن الزائف (أي، أقدم خلايا نشطة مناعياً) (الشكل 4هـ، و4ف والشكل التمديدي 9ب)، مما يشير إلى أن هذه الخلايا هي المستجيبة المبكرة لغزو العوامل الممرضة. لذلك، نحدد هذه الحالة الخلوية كخلايا PRIMER، التي قد تنتشر منها الاستجابات المناعية. تُسمى الخلايا المحيطة بخلايا PRIMER بالخلايا المتفرجة.

أظهرت خلايا PRIMER وخلايا المتفرجين توقيعات نسخية وإبيجينية مميزة. كانت خلايا PRIMER غنية في أنماط CAMTA ونمط GT-3A، في حين كانت خلايا المتفرجين غنية في WRKY.

الشكل 5 | تحديد وتوصيف وظيفي لمجموعات خلايا المناعة. أ، الجينات التي تم تحليلها بواسطة MERFISH والتي تظهر تغييرات ملحوظة في التعبير عبر مسار الزمن الزائف (الشكل 4f) في عينة AvrRpt2 بعد 9 ساعات من الإصابة. تشير قيم التغير الإيجابي إلى تعبير أعلى في الخلايا ذات قيم الزمن الزائف العالية. أعلى الجينات المعبر عنها بشكل مختلف (DEGs؛ ) تم وسمها باللون الأحمر. ب-د، التعبير المكاني لـ ALD1 و FMO1 و BON3 بعد 9 ساعات من الإصابة. (ب) وصور مقطوعة لموقعين مختلفين في ( و د). هـ، تعبير ALD1 و FMO1 و BON3 في تحت المجموعات النشطة مناعياً من بيانات تسلسل RNA المفرد. تعبير mRNA (f)، صور MERFISH (g) ونشاط الأنماط (h) لـ GT-3A في خلايا الميزوفيل النشطة مناعياً في بيانات snMultiome. يتم أيضاً عرض تعبير BON3 في g. i، الأنماط الغنية المرتبطة بالجينات الغنية في خلايا PRIMER أو خلايا المتفرج. j، مخطط فين يوضح تداخل الجينات المعبر عنها بشكل كبير في خلايا PRIMER أو خلايا المتفرج مع الجينات التي تم إظهار أنها مثبطة بواسطة CAMTA3 (المرجع 39). k، تسلسل RNA الكمي للنباتات من النوع البري (WT) ونباتات GT-3A-ox.
تمت معالجته بالماء (M) أو DC3000 (DC) بعد 24 ساعة من الإصابة. DEGs (معدلة ، -تغيير الطي المعدل ) بين WT و GT-3A-ox بعد العدوى موضحة. يتم عرض مصطلحات GO الغنية في مجموعات الجينات المحددة مع التعديل القيم. ثلاث تكرارات لكل حالة. نمو DC3000 (I) أو AvrRpt2 (m) في نباتات Col-0 (WT) مقارنة بنباتات GT-3A-ox (I) ونباتات gt3a-KO (m) في الأوقات المحددة. و أوراق من 3 و 4 تجارب مستقلة، على التوالي. RLU، وحدة الضوء النسبية. k-m، تم حقن مسببات الأمراض باستخدام الحقن في جرعة تعليق من معدل تم حساب القيم باستخدام اختبار الطالب ذو الطرفين -اختبار مع تصحيح بنجاميني-هوشبرغ. النتائج موضحة كرسوم بيانية صندوقية، حيث تمثل الصناديق النسب المئوية من 25 إلى 75، والخط المركزي يشير إلى الوسيط، وتمتد الشعيرات إلى القيم الدنيا والقصوى ضمن 1.5 مرة من النطاق الربعي. النموذج المقترح للدور المحتمل وتنظيم حالات خلايا المناعة. شريط القياس، (ج، د)، (غ) أو 1 مم (ب).
الموتيفات (الشكل 5i والشكل الإضافي 9i). كانت الجينات التي تم إظهار أنها مثبطة بواسطة CAMTA3 مفرطة التمثيل بشكل كبير في خلايا المتفرج مقارنة بخلايا PRIMER (الشكل 5j؛ معدل الاكتشاف الخاطئ (FDR) ، مع تصحيح اختبار هايبرجيومايتريك باستخدام طريقة بنجاميني-هوشبرغ)، وهو نتيجة تدعم الدور المثبط للنسخ لـ CAMTA3 في خلايا PRIMER.

لفهم وظيفة خلايا PRIMER وآليات تنظيم الجينات الخاصة بها، قمنا بالتحقيق في مجموعة فرعية من snMultiome حيث تم إثراء (المجموعة الفرعية 4). على الرغم من أن العديد من جينات علامات خلايا PRIMER لم تكن موصوفة سابقًا، وجدنا عدة جينات معروفة، بما في ذلك WRKY8 و LSD1، كعلامات لخلايا PRIMER (الشكل 9c من البيانات الموسعة). ومن الجدير بالذكر أن BON3 و WRKY8 و LSD1 جميعها من المنظمات السلبية للمناعة. على الرغم من أن الآليات الدقيقة لكيفية قمعها للاستجابات المناعية لا تزال غامضة.

تحليل جينات العلامة الإضافية لمجموعة خلايا PRIMER حدد GT-3A (الذي يشفر عامل نسخ مرتبط بالحمض النووي ثلاثي الحلزونات) وظيفتها في مناعة الأوراق لا تزال غير موصوفة (الشكل 5f). أظهرت GT-3A نمط تعبير مكاني مشابه لذلك الخاص بـ BON3 (MERFISH؛ الشكل 5g)، وكان نمط GT-3A متاحًا بشكل كبير في خلايا PRIMER (snATAC-seq؛ الشكل 5h)، مما يوحي بأن له وظيفة تنظيم جيني في هذه المجموعة الخلوية. لفهم دور GT-3A، قمنا بإنشاء نباتات A. thaliana معدلة وراثيًا تعبر عن هذا الجين بشكل مفرط (GT-3A-ox). كشفت تحليلات RNA-seq الشاملة عن ضعف في تحفيز الجينات المعنية في مسار SA في نباتات GT-3A-ox المصابة بـ DC3000 (الشكل 5k). علاوة على ذلك، تم الحصول على سطرين مستقلين من

كانت خلايا GT-3A-ox أكثر عرضة للإصابة بـ DC3000 (الشكل 5ل)، مما يشير إلى أن هذا العامل النسخي المحدد لحالة الخلية غير المعروف سابقًا يمكن أن ينظم المناعة سلبًا. كما اختبرنا النباتات التي تم فيها حذف GT-3A (gt3a-KO)، وكانت هذه الطفرة أكثر عرضة للإصابة بـ AvrRpt2 (الشكل 5م). تشير هذه النتيجة إلى أن وظيفة GT-3A في خلايا PRIMER ضرورية للدفاع الأمثل ضد مسببات الأمراض غير الفتاكة. تحليل snRNA-seq لـ أظهرت نباتات -KO جينات قد تكون منظمة بواسطة GT-3A، إما بشكل مباشر أو غير مباشر. من بين هذه الجينات، كان تعبير PUB36 في خلايا PRIMER متأثراً بشكل كبير في -KO النباتات في كل من 9 و 24 ساعة بعد الإصابة. من الجدير بالذكر، يمتلك نمط ربط GT-3A في المنطقة العلوية (الشكل 9d من البيانات الموسعة). تم التعبير عن الجينات المرتبطة بـ SA، بما في ذلك ALD1، بمستويات أقل في خلايا المتفرجين في نباتات gt3a-KO مقارنة بالنباتات البرية (الشكل 9e، f من البيانات الموسعة). تشير هذه النتيجة إلى أن تحفيز GT-3A في خلايا PRIMER مهم للتحفيز الصحيح لجينات الدفاع في الخلايا المحيطة. باختصار، من خلال دمج بيانات snMultiome الزمنية والبيانات النصية المكانية، حددنا حالات خلايا مناعية غير معروفة سابقًا وعامل نسخ محدد لسكان الخلايا ينظم مناعة النبات.

تشمل مناعة النباتات شبكة متعددة الخلايا حيث تقوم الخلايا الفردية بتفسير إشارات المدخلات من خلال شبكاتها الجزيئية المميزة وتتواصل مع خلايا أخرى. إن تعريفنا الجزيئي المكاني الزمني
كشف أطلس الأوراق المصابة بالعوامل الممرضة عن حالات خلوية متنوعة بتفاصيل النسخ الجيني والحمض النووي الريبوزي. توفر هذه الموارد أيضًا وسيلة للتحقيق في حالات الخلايا الفردية التي تم إخفاؤها في التحليلات التقليدية للأنسجة الكلية أو المقطعة ومن خلال التصوير الحي لعدد محدود من خطوط التقارير. على سبيل المثال، حددنا مجموعة فرعية من PCC بحالة مميزة تتسم بتحفيز جينات SAR (الشكل 1e) ومجموعات فرعية من الميزوفيل التي تنشط فروعًا مختلفة من مسارات الأيض الدفاعية المشتقة من التربتوفان (الشكل 1g).

وصفنا أيضًا حالة خلوية نادرة تقع في مركز النقاط الساخنة النشطة مناعياً، والتي أطلقنا عليها اسم حالة الخلية PRIMER (الشكل 5). بالإضافة إلى رسم خرائط للجينات التي تم تحديدها سابقًا (BON3، WRKY8، LSD1 وCAMTA3) والتي لها وظائف مثبطة للمناعة شائعة، حددت تحليلاتنا المتكاملة snMultiome وMERFISH جين علامة خلوية أخرى لـ PRIMER، GT-3A، الذي يشفر عامل نسخ، وأظهرت أنه يساهم في مناعة النبات ضد عدوى مسببات الأمراض (الشكل 5). كشفت المقارنات بين خلايا PRIMER والخلايا المحيطة بها (خلايا المتفرج) عن مناظر نسخية وإيبيجينية مميزة (الشكل 5). من الممكن أن تكون هناك حالات خلوية متخصصة إضافية في خلايا PRIMER وخلايا المتفرج. يتطلب فهم أعمق لحالات خلايا المناعة تطوير وتطبيق منهجيات جديدة، مثل تقنيات التصوير التي تصور كل من بروتينات الفاعل الممرض والعديد من جينات النبات في الوقت نفسه في ثلاثة أبعاد بدقة خلوية فردية. (انظر القسم ‘قيود هذه الدراسة’ في المعلومات التكميلية).

أظهرت الدراسات السابقة أن GT-3A ينظم سلبًا دفاع النبات ضد عدوى الديدان الخيطية في الجذر. ، مما يشير إلى أن GT-3A قد يلعب دورًا في الدفاع ضد أنواع مختلفة من مسببات الأمراض التي تصيب أنسجة أخرى. دعمًا لهذه الفرضية، كانت نباتات GT-3A-ox أكثر عرضة لمسبب الأمراض الفطرية Colletotrichum higginsianum (الشكل 9j من البيانات الموسعة).
لقد أظهرت الاكتشافات الحديثة دور مستقبلات مجال ربط النوكليوتيدات وتكرارات الغلوتامين الغنيّة بالليوسين كقنوات كالسيوم أو NADases. . ومع ذلك، فإن الآليات الدقيقة التي تكمن وراء تنشيط ETI وكيفية قمع هذه العملية لنمو مسببات الأمراض ليست مفهومة تمامًا. وقد تم اقتراح أن المقاومة المكتسبة الموضعية (LAR) قد تكون مهمة لـ استجابة LAR هي استجابة دفاعية قوية في الخلايا المحيطة بتلك التي تعرضت لعوامل مسببة للأمراض. كشفنا عن انتشار استجابات ETI من خلايا PRIMER (الشكل 5e) والتقطنا استجابات LAR المحتملة بمعلومات جزيئية مفصلة. من الجدير بالذكر أن ALD1 و FMO1، اللذان يشفران مكونات SAR التقليدية، تم التعبير عنهما في الخلايا المجاورة ولكن ليس في خلايا PRIMER (الشكل 5b-e)، مما يشير إلى أن مسار SAR له دور محدد في استجابات LAR. هذه الفكرة معقولة نظرًا لقدرة SAR على الإشارة على مسافات طويلة. نقترح أن خلايا PRIMER قد تمر بموت خلوي فرط الحساسية، مما يرسل بعد ذلك إشارات إلى الخلايا المجاورة التي تنشط الاستجابات المناعية، بما في ذلك مسار SAR. حيث يمكن أن يؤدي تنشيط مسار SA إلى قمع موت الخلايا. قد يساهم GT-3A في موت الخلايا في خلايا PRIMER من خلال قمع إشارات SA (الشكل 5n). يمكن أن تكشف التحليلات الإضافية لبياناتنا عن أدوار حالات الخلايا المناعية المختلفة وكيف تتواصل هذه الخلايا مع الخلايا المحيطة لتحقيق دفاع ناجح.

بالإضافة إلى تحديد مجموعات خلايا المناعة غير المعروفة سابقًا، توقعت بيانات snMultiome لدينا العديد من وحدات TF-ACR-gene المحتملة. لقد مكنتنا التكامل الناجح بين بيانات snMultiome وMEFISH من رسم خرائط مكانية للتعبير الجيني وحالات الكروماتين (الشكل التوضيحي الممتد 8a-d). يمكن استخدام هذه النتائج لاكتشاف جينات متعلقة بالمناعة لم يتم تصنيفها سابقًا، كما يمكن أيضًا تحديد عناصر CRE المتعلقة بالدفاع من خلال تحليل الأنماط de novo.

أخيرًا، قمنا بإنشاء قاعدة بيانات (https://plantpathogenatlas.salk.edu) لتسهيل استكشاف تجمعات الخلايا غير الموصوفة سابقًا المرتبطة بالمرض والمقاومة مع معلومات مكانية وزمنية وآليات تنظيمية محتملة. يمكن استخدام قاعدة بياناتنا لتوليد الفرضيات واختبارها، وستحفز اكتشافات جديدة للآليات الجزيئية التي تكمن وراء تفاعلات النباتات والميكروبات بدقة عالية.

أي طرق، مراجع إضافية، ملخصات تقارير Nature Portfolio، بيانات المصدر، بيانات موسعة، معلومات تكميلية، شكر وتقدير، معلومات مراجعة الأقران؛ تفاصيل مساهمات المؤلفين والمصالح المتنافسة؛ وبيانات توفر البيانات والرموز متاحة علىhttps://doi.org/10.1038/s41586-024-08383-z.

ملاحظة الناشر: تظل شركة سبرينغر ناتشر محايدة فيما يتعلق بالمطالبات القضائية في الخرائط المنشورة والانتماءات المؤسسية.

الوصول المفتوح هذه المقالة مرخصة بموجب رخصة المشاع الإبداعي النسب 4.0 الدولية، التي تسمح بالاستخدام والمشاركة والتكيف والتوزيع وإعادة الإنتاج بأي وسيلة أو صيغة، طالما أنك تعطي الائتمان المناسب للمؤلفين الأصليين والمصدر، وتوفر رابطًا لرخصة المشاع الإبداعي، وتوضح إذا ما تم إجراء تغييرات. الصور أو المواد الأخرى من طرف ثالث في هذه المقالة مشمولة في رخصة المشاع الإبداعي الخاصة بالمقالة، ما لم يُشار إلى خلاف ذلك في سطر الائتمان للمواد. إذا لم تكن المادة مشمولة في رخصة المشاع الإبداعي الخاصة بالمقالة وكان استخدامك المقصود غير مسموح به بموجب اللوائح القانونية أو يتجاوز الاستخدام المسموح به، فستحتاج إلى الحصول على إذن مباشرة من صاحب حقوق الطبع والنشر. لعرض نسخة من هذه الرخصة، قم بزيارةhttp://creativecommons.org/licenses/by/4.0/.
(ج) المؤلف(ون) 2025

تم استخدام الكواشف والمجموعات التالية: بيرسولفات الأمونيوم (سيغما، 09913)؛ محلول BSA (سيغما، A1595)؛ مجموعات Chromium Next GEM Single Cell Multiome ATAC + Gene Expression (10x Genomics، PN-1000283)؛ مصافي خلايا Corning Falcon (كورنينغ، 08-771-2)؛ ثيوثريتال (ثيرمو، R0861)؛ EDTA، pH 8.0 خالي من RNase (إنفيتروجين، AM9260G)؛ KCl (2 M) خالي من RNase (إنفيتروجين، AM9640G)؛ مصافي MACS SmartStrainers (ميلتيني بيوتيك، 130-098-458)؛ مجموعات تصوير الجينات MERSCOPE 500 (فيزجن، 10400006)؛ لوحة جينات MERSCOPE 500 (فيزجن، 10400003)؛ مجموعات تحضير العينات MERSCOPE (فيزجن، 10400012)؛ -تترا ميثيل إيثيلين ديا مين (سيغما، T7024-25ML)؛ NaCl (5 M) خالي من RNase (إنفيتروجن، AM9759)؛ NP40 (IGEPAL CA-630) (سيغما، I8896)؛ بارا فورمالدهيد (سيغما، F8775)؛ PBS ( ) pH 7.4 خالي من RNase (ثيرمو فيشر، AM9625)؛ خليط مثبط للبروتيناز (سيغما، P9599)؛ مثبط RNase Protector (سيغما، 3335402001)؛ مركب Scigen Tissue-Plus OCT (فيشر، 23-730-571)؛ سبيرميدين (سيغما، S2626)؛ سبيرمين (سيغما، 85590)؛ تريتون-X (سيغما، 93443)؛ ومخزن UltraPure 1 M Tris-HCI pH 7.5 (إنفيتروجين، 15567027).

تم تسليط الضوء على الجينات التالية في هذه الدراسة: ALD1 (بروتين الاستجابة الدفاعية الشبيه بـ AGD2 1، At2g13810)؛ BCA2 (الكربونات الأنزيمية بيتا 2؛ At5g1474O)؛ BON3 (بونزاي 3؛ At1gO886O)؛ CBP6Og (بروتين ربط الكالمودولين الشبيه بـ 6O-G؛ At5g2692O)؛ FDH (فيديليهاد؛ At2g26250)؛ FMO1 (مونوأوكسيجيناز معتمد على الفلافين 1؛ At1g19250)؛ GT-3A (عامل ربط الحمض النووي ثلاثي الحلزونات (عامل GT) 3A؛ At5gO138O)؛ ICS1 (سينثاز الإيزوشوريسمات 1؛ At1g7471O)؛ LSD1 (مرض محاكاة الآفات 1؛ At4g2O38O)؛ HSFB2b (عامل النسخ الناتج عن الصدمة الحرارية B2b؛ At4g1166O)؛ ILL6 (جين مقاوم لحمض الإندول-الليوسين (ILR) الشبيه بـ 6؛ At1g4435O)؛ MAM1 (سينثاز ميثيلثيوالكيلماليت 1؛ At5g23010)؛ WRKY8 (At5g46350)؛ و WRKY46 (At2g46400).

تم زراعة A. thaliana Col-0 في غرفة عند مع فترة خفيفة و الرطوبة النسبية لـ أيام. تم زراعة سلالات البكتيريا في وسط سائل كينغ ب مع المضادات الحيوية (ريفامبيسين وتتراسيكلين) في تم وصف ثلاثة سلالات بكتيرية – P. syringae pv. tomato DC3000 مع متجه فارغ (pLAFR3)، avrRpt2 (pLAFR3) و avrRpm1 (pLAFR3) – سابقًا. تم جمع البكتيريا عن طريق الطرد المركزي وإعادة تعليقها في ماء معقم إلى 0.001 (تقريباً وحدة تشكيل المستعمرات ). في المجموع، تم حقن 20 ورقة من A. thaliana (4 أوراق متوسعة بالكامل لكل نبات) باستخدام محاقن بدون إبرة. تم إجراء حقن المحقنة في زاوية أو زاويتين من كل ورقة، وتم توزيع المعلق البكتيري في جميع أنحاء الورقة بأكملها. اخترنا أربعة نقاط زمنية مختلفة ( و 24 ساعة)، تمثل المراحل المبكرة من العدوى وعندما تم ملاحظة إعادة برمجة النسخ الديناميكية في دراسة سابقة استخدمت تسلسل RNA الشامل. لكل سلالة، تم أخذ أربع نقاط زمنية في نفس الوقت من اليوم عن طريق تسريب البكتيريا في أوقات مختلفة لتقليل تأثير الإيقاعات اليومية. تم جمع 20 ورقة مصابة باستخدام الملاقط ومعالجتها على الفور لاستخراج النوى. بالنسبة للحالة الضابطة، تم جمع الأوراق التي تم تسريبها بالماء بعد 9 ساعات.

كان pWAT206 هدية من A. Takeda. تم بناء البلازميدات التالية باستخدام مجموعات تجميع الحمض النووي HiFi (New England Biolabs) وتم التحقق منها بواسطة تسلسل سانجر. تم تضخيم ثلاث قطع PCR من pBICRMsG. باستخدام أزواج البرايمر (mEGFP-Nter_1_F بالإضافة إلى mEGFP_1_R؛ mEGFP_2_F بالإضافة إلى mEGFP_2_R؛ وmEGFP_3_F بالإضافة إلى mEGFPNter_3_R) ثم تم تجميعها في pBICAscII المقطوع بواسطة Stul/Ascl. . الناتج
تم استخدام البلازميد كقالب لردين من تفاعل البلمرة المتسلسل (PCR) باستخدام أزواج من البرايمرات mEGFP_Nter_1_F بالإضافة إلى mEGFP_4R و mEGFP_4F بالإضافة إلى mEGFP_Nter_3R. تم تجميع هذه المنتجات من تفاعل البلمرة المتسلسل في pBICAscII المقطوع بواسطة Stul/Ascl، مما أدى إلى إنتاج pBIC_mEGFP_Nter. تم تضخيم قطعة من تفاعل البلمرة المتسلسل من pBIC_mEGFP_Nter باستخدام mEGFP_Cter_F1 بالإضافة إلى mEGFP_Cter_R، ثم تم استخدامها كقالب لتفاعل البلمرة المتسلسل باستخدام mEGFP_Cter_F2 بالإضافة إلى mEGFP_Cter_R. تم تجميع قطعة تفاعل البلمرة المتسلسل الناتجة في pBIC_mEGFP_Nter المقطوع بواسطة Stul/AscI للحصول على pBIC_mEGFP_Cter. تم تضخيم تسلسل الترميز لـ GT-3A (At5g1380) بواسطة تفاعل البلمرة المتسلسل باستخدام زوج البرايمرات GT3a_mEGFP_F بالإضافة إلى GT3a_mEGFP_R وتم تجميعه في pBIC_mEGFP_Cter المقطوع بواسطة Stul للحصول على pBIC_GT3a_mEGFP. تم تضخيم التسلسل الذي يشفر GT-3A المدمج مع mEGFP من pBIC_GT3a_mEGFP بواسطة تفاعل البلمرة المتسلسل باستخدام زوج من البرايمرات ( وتم تجميعه بعد ذلك في pWAT206 المهضوم بواسطة Xhol/Xbal للحصول على pWAT206_GT3a_mEGFP. تم تحويل نباتات A. thaliana Col-0 باستخدام Agrobacterium tumefaciens GV3101 (pMP90، pSoup) مع pWAT206_GT3a_mEGFP كما هو موصوف سابقًا. تم استخدام سطرين مستقلين مقاومين لباستا لإجراء اختبارات نمو البكتيريا باستخدام سلالة بيو لومينسنت من P. syringae pv. tomato DC300043 ولتحليل تطور الآفات باستخدام C. higginsianum MAFF 305635. تم زراعة النباتات في غرفة مناخية بدرجة حرارة الرطوبة النسبية وشدة الإضاءة 6000 لوكس (حوالي ) لمدة 10 ساعات. تعليقات بكتيرية في تم حقن الماء المعقم في أوراق نباتات عمرها أسبوع. تم قياس نمو البكتيريا كضوء متألق باستخدام جهاز قياس اللمعان GloMax Navigator Microplate (Promega) كما تم وصفه سابقًا. لتحليل تطور الآفات، تم تلقيح أوراق النباتات التي تتراوح أعمارها بين 4-5 أسابيع. تعليق كونيدي لـ C.higginsianum كونيديا لكل مل). تم الاحتفاظ بالنباتات الملقحة تحت رطوبة عالية في غرفة مناخية بدرجة حرارة وشدة الإضاءة 6000 لوكس لمدة 10 ساعات. تم قياس حجم الآفة بعد 6 أيام من التلقيح باستخدام ImageJ.

تم إنشاء نباتات gt3a-KO من خلال تحرير الجينوم باستخدام CRISPR-Cpf1. تم توفير البلازميد pWAT235-Cas9-HF-cc بواسطة A. Takeda. تم تضخيم تسلسل Cpf1 من pY004 (Addgene، 69976) بواسطة PCR باستخدام زوج من البرايمرات (Spel-NLS-Cpf1-F و BamHI-NLS-Cpf1-R؛ الجدول التكميلية 2)، تلاها هضم باستخدام Spel و BamHI. تم ربط شظية الحمض النووي المهضوم مع pWAT235-Cas9-HF-cc، الذي تم هضمه باستخدام Spel و BamHI، مما أدى إلى pWAT235-Cpf1-cc. تم تضخيم شظيتين من PCR من pWAT235-Cpf1-cc باستخدام أزواج البرايمرات Gb-ccdB-F و Gb-ccdB-R، و Gb-U626-F و Gb-U626-R. تم تجميع هذه الشظايا من PCR و pWAT235-Cas9-HF-cc المهضوم بواسطة Bsal في بلازميد واحد باستخدام مجموعة تجميع الحمض النووي HiFi، مما أدى إلى pWAT235-Cpf1-U626-cc-polyT. تم تزاوج وربط أوليغونوكليوتيدين، GT3a_gRNA2_F و GT3a_gRNA2_R، مع pWAT235-Cpf1-U626-cc-polyT المهضوم بواسطة Bsal، مما أدى إلى pWAT235-Cpf1-U626-GT3a-gRNA2-polyT. تم تحويل نباتات A. thaliana Col-O باستخدام A. tumefaciens GV3101 (pMP90، pSoup) مع pWAT235-Cpf1-U626-cc-polyT كما هو موصوف سابقًا. تم اختيار خط ترانسجيني يحتوي على كودون توقف مبكر بسبب حذف مكون من 5 قواعد في جين GT-3A واستخدامه طوال هذه الدراسة.

تم حقن أوراق النباتات من النوع البري وواحدة من الخطوط المعدلة وراثيًا التي تعبر عن GT-3A-mEGFP بالماء أو DC3000 عند 0.001 وتم جمعها بعد 24 ساعة من التسلل. تم استخراج RNA الكلي باستخدام مادة TRI (سيغما). تم استخدام ميكروغرام واحد من RNA الكلي لتحضير المكتبة باستخدام طريقة BrAD-seq لإنشاء مكتبة محددة الاتجاه. مكتبات التعبير الجيني الرقمية . تم تسلسل المكتبات على منصة DNBSEQ-G400 في BGI، والتي أنتجت قراءات بطول 100 قاعدة من الطرف. نظرًا لأن جودة القراءات العكسية كانت ضعيفة بسبب البوليمر تم استخدام القراءات الأمامية فقط للتحليل. تم تقليم أول 8 قواعد والمهايئات وتصنيف الجودة.
تم التنفيذ باستخدام fastp (الإصدار 0.19.7) مع المعلمات -x-f8-q30-b50. تم تعيين القراءات المصفاة والمفلترة من حيث الجودة إلى جينوم الأرابيدوبسيس (TAIR10) باستخدام STAR (الإصدار 2.6.1b) مع المعلمات الافتراضية وتم تحويلها إلى عدد لكل جين لكل مكتبة باستخدام featureCounts (الإصدار 1.6.0) تم إجراء التحليل الإحصائي لبيانات RNA-seq في بيئة R (الإصدار 4.1.3). نظرًا لأن طريقة BrAD-seq تتضمن إثراء البوليمر (A)، تم استبعاد الجينات الميتوكوندرية والكلوروبلاستية. تم استبعاد الجينات التي تحتوي على متوسط عدد قراءات أقل من عشرة لكل مكتبة من التحليل. تم إخضاع بيانات العد الناتجة لتطبيع TMM باستخدام الدالة calcNormFactors في حزمة edgeR، تلاها تحويل لوغاريتمي باستخدام الدالة voomWithQualityWeights في حزمة limma. تم ملاءمة نموذج خطي لكل جين باستخدام الدالة ImFit في حزمة limma. من أجل تقليص التباين في حساب القيم، تم استخدام دالة eBayes في حزمة limma. النتيجة تم تصحيح القيم بعد ذلك لاختبار الفرضيات المتعددة عن طريق حساب ستوري. القيم باستخدام الدالة qvalue في حزمة qvalue. لاستخراج الجينات التي تظهر تغييرات ملحوظة في التعبير، قيم القطع من و -تم تطبيق التغير النسبي المحول | > 1.

محلول تنقية النواة الطازج (NPB؛ 15 مللي مولار تريس pH 7.5، 2 مللي مولار EDTA، تم تحضير (سبيرميدين، 0.2 مللي مول من سبيرمين، 1:100 من BSA و 1:100 من خليط مثبطات البروتياز) قبل التجربة وتم تبريده على الثلج. تم تنفيذ جميع الإجراءات اللاحقة على الثلج أو عند تم تقطيع عشرين ورقة NPB بارد مع 1:500 من Protector RNase IN باستخدام شفرة حلاقة على الثلج لمدة 5 دقائق لتحرير النوى ثم تم حضنه في 20 مل من NPB. تم تصفية المستخلص الخام للنوى بشكل متسلسل من خلال و مصفاة خلايا مرشحات خلايا كورنينغ فالكون، كورنينغ، 08-771-2؛ ، مصافي MACS الذكية، ميلتيني بيولوجي، 130-098-458). تم إضافة Triton-X و NP40 إلى المستخلص بتركيز نهائي من كل منها، وتم حضن المستخلص في لمدة 5 دقائق مع دوران. تم طرد التعليق في جهاز الطرد المركزي عند 50 جرام لمدة 3 دقائق في جهاز طرد مركزي ذو دوار متأرجح لجمع الحطام غير النووي وتم استرداد السائل العلوي. تم جمع النوى عن طريق الطرد المركزي عند 500 جرام لمدة 5 دقائق في جهاز طرد مركزي ذو دوار متأرجح. عندما كان الحبيبة خضراء، تم إعادة تعليق الحبيبة في 20 مل من NPBd (NPB مع 0.1% Triton-X و0.1% NP40) مع 1:1,000 Protector RNase IN عن طريق السحب، تلاها الطرد المركزي عند 500 جرام لمدة 5 دقائق. عندما كانت الحبيبة شبه شفافة، تم تخطي غسل NPBd. ثم تم غسل الحبيبة عن طريق إعادة تعليقها في 20 مل من NPB مع 1:1,000 Protector RNase IN والطرد المركزي عند 500 جرام لمدة 5 دقائق في جهاز طرد مركزي ذو دوار متأرجح. تم إعادة تعليق الحبيبة في من عازلة النوى (10x Genomics، PN-2000207) مع 1:40 من مثبط RNase IN و 1 مللي مول من ثيوثريتال. تم طرد تعليق النوى في جهاز الطرد المركزي عند 50 جرام لمدة 3 دقائق في جهاز طرد مركزي دوار لجمع الحطام غير النووي وتم استرداد السائل العلوي. تم تكرار هذه الخطوة مرة أخرى. تم عد النوى الناتجة يدويًا باستخدام جهاز عد الخلايا. تم جمع النوى عن طريق الطرد المركزي عند 500 جرام لمدة 5 دقائق في جهاز طرد مركزي دوار وتم إزالة السائل العلوي، مما ترك حوالي للمخزن. تم عد النوى مرة أخرى، وتم استخدام ما يصل إلى 16,000 نواة للخطوات التالية. ومع ذلك، في معظم العينات، لم نقم بتحميل الحد الأقصى من النوى التي يمكن أن تم قبول مجموعة الجينوميات لتجنب خطر انسداد الجهاز، مما قد يؤدي إلى أعداد متغيرة من النوى المستعادة بين العينات (الشكل البياني الممتد 1ب). تم بناء مكتبات scRNA-seq و ATAC-seq وفقًا لتعليمات الشركة المصنعة (10x Genomics، CG000338). تم تسلسل مكتبات scRNA-seq باستخدام جهاز Illumina NovaSeq 6000 في وضع الفهرسة المزدوجة مع عشرة دورات لمؤشرات i7 و i5. كما تم تسلسل مكتبات snATAC-seq باستخدام جهاز Illumina NovaSeq 6000 في وضع الفهرسة المزدوجة مع 8 و 24 دورة لمؤشر i7 ومؤشر i5، على التوالي.

معالجة البيانات الخام. تم معالجة بيانات التسلسل للحصول على عدد ميزات الخلايا الفردية عن طريق تشغيل cellranger (الإصدار 6.0.1) و cellranger-arc.
(الإصدار 2.0.0) لبيانات snRNA-seq وبيانات snATAC-seq، على التوالي. بالنسبة لـ snRNA-seq، تم استخدام خيار -include-introns لمحاذاة القراءات إلى النسخة النووية لـ A. thaliana التي تم بناؤها باستخدام جينوم TAIR10 ونسخة Araport 11. تم إزالة جينوم البلاستيد من الجينوم المرجعي لتحليل بيانات كل من snRNA-seq وsnATAC-seq. تم تنزيل جينوم .thaliana TAIR10 منhttps://plants.ensembl.org/، وتمت إزالة جينوم البلاستيد يدويًا. الـ تم تعديل ملف توضيح الجينات لـ .thaliana TAIR10 يدويًا عن طريق إزالة جينات البلاستيدات الخضراء واستبدال الفواصل المنقوطة في عمود ‘اسم الجين’ بشرطات لمنع الأخطاء أثناء معالجة cellranger-arc. نلاحظ أن المتوسط هو و تمت محاذاة عدد من القراءات على جينوم البلاستيد في بيانات snRNA-seq و snATAC-seq، على التوالي. لم يؤثر إزالة جينوم البلاستيد من المرجع على توزيع عدد RNA و ATAC بشكل عام. تم تحليل بيانات العد باستخدام حزم R Seurat (الإصدار 5). وسينياك .

مراقبة الجودة وتصنيف الخلايا. قبل دمج مجموعات البيانات، تم توقع الازدواجيات باستخدام DoubletFinder. وتم تصفيتها. تم إنشاء مصفوفات مراقبة الجودة لـ snATAC-seq باستخدام نسخة معدلة من دالة loadBEDandGenomeData في حزمة R Socrates تم تحديد ACRs باستخدام MACS2 (المرجع 53) مع المعلمات التالية: -g (حجم الجينوم) تحول ، حجم الإضافات ، و –qvalue ، –nomodel، –keep-dup جميعها. تم حساب نسبة القراءات التي تتطابق مع 2 كيلو بايت upstream أو 1 كيلو بايت downstream من موقع بدء النسخ (TSS). تم تصفية النوى باستخدام المعايير التالية: 200 < عدد RNA UMI < 7,000؛ عدد جينات RNA > 180؛ 200 < عدد ATAC UMI < 20,000؛ ونسبة قراءات RNA المتطابقة مع الجينوم الميتوكوندري. تم دمج كائنات Seurat من العينات الفردية باستخدام دالة Merge_Seurat_List من حزمة scCustomize.
تجميع snRNA-seq. تم إجراء تجميع snRNA-seq باستخدام حزمة R Seurat. تم تطبيع مصفوفة عدد RNA لكل خلية باستخدام SCTransform. تم إجراء تقليل الأبعاد باستخدام تحليل المكونات الرئيسية مع RunPCA. تم تقليل التباين الفني بين العينات باستخدام Harmony. استخدام المكونات الرئيسية (PCs) 1-20. تم إجراء تجميع قائم على الرسم البياني على المكونات الرئيسية 1-20 المصححة بواسطة هارموني من خلال حساب رسم بياني للجيران الأقرب المشترك باستخدام الفضاء منخفض الأبعاد للمكونات الرئيسية (مع تجميع لوفيان (الدقة ) تم تطبيقه بعد ذلك، وتم إسقاط المجموعات في مساحة مخفضة إضافيًا باستخدام UMAP (عدد الجيران و min.dist ).
استدعاء قمم snATAC-seq. تم استدعاء القمم بشكل مستقل على كل مجموعة تم تعريفها بواسطة بيانات snRNA-seq ثم تم دمجها باستخدام وظيفة CallPeaks من Signac، التي تستخدم MACS2 مع المعلمات التالية: effective.genome.size ، حجم الإضافات تحول تم قياس عدد القمم باستخدام وظيفة FeatureMatrix. بالمقارنة مع خط أنابيب استدعاء القمم الافتراضي لـ cellranger-arc (24,394 قمة)، كانت هذه الطريقة الخاصة باستدعاء القمم المحددة للمجموعة قادرة على التقاط المزيد من القمم (35,560 قمة)، وهو ما يتماشى مع تقرير سابق. .
تجميع snATAC-seq. تم إجراء تقليل الأبعاد باستخدام الفهرسة الدلالية الكامنة (LSI) . أولاً، الأعلى تم اختيار أكثر الميزات شيوعًا باستخدام دالة FindTopFeatures. ثم تم حساب تردد المصطلح عكسي تردد الوثيقة (TF-IDF) باستخدام RunTFIDF مع عامل المقياس. . تم تحليل مصفوفة TF-IDF الناتجة باستخدام تحليل القيم الفردية مع RunSVD، الذي يستخدم حزمة R المعروفة باسم irlba. تم تقليل التباين الفني بين العينات باستخدام Harmony مع مكونات LSI من 2 إلى 10. تم إجراء التجميع القائم على الرسم البياني على مكونات LSI المصححة بواسطة Harmony من 2 إلى 20 من خلال حساب رسم بياني للجيران الأقرب المشترك باستخدام الفضاء منخفض الأبعاد LSI (مع تجميع لوفين (الدقة ثم تم تطبيق )، وتمت
تم إسقاطه في مساحة مخفضة إضافيًا باستخدام UMAP (عدد الجيران = 30 و min.dist ).

تجميع RNA-ATAC المشترك. تم دمج النمطين من خلال تحليل الجوار الأقرب الموزون باستخدام FindMultiModalNeighbors من Seurat مع مكونات PCs 1-20 المصححة بواسطة Harmony لـ RNA ومكونات LSI 2-20 المصححة بواسطة Harmony لـ ATAC. ثم، SLM (الدقة ) تم تطبيقه وإسقاطه باستخدام UMAP (عدد الجيران و min.dist=0.1).

درجة نشاط جين ATAC. تم حساب درجة نشاط جين ATAC باستخدام دالة GeneActivity من Signac مع extend.upstream=400.

تحليل ارتباط القمة بالجين. تم استخدام LinkPeaks من Signac لاستدعاء ارتباط القمة بالجين المهم لكل حالة عدوى (وهمية، DC3000، AvrRpt2 وAvrRpm1). تم حساب معاملات الارتباط بيرسون المصححة للخلفية بين تعبير الجين لكل جين والوصول إلى كل قمة ضمن 500 كيلوبايت من نقطة بدء النسخ للجين. تم حساب القيمة لكل زوج من القمة والجين باستخدام اختبار أحادي الجانب -اختبار، وأزواج الجينات القمة مع ومعامل ارتباط بيرسون تم الاحتفاظ بها كروابط هامة.

تحليل إثراء الأنماط. الأنماط الموجودة في قاعدة بيانات JASPAR2020 تم استخدام Arabidopsis (رمز النوع 3702). تم تحديد القمم المحددة للمجموعة أولاً باستخدام FindMarkers من Seurat مع المعلمات الافتراضية. تم استخدام اختبار هايبرجومتري لاختبار الزيادة في تمثيل الأنماط في مجموعة القمم القابلة للوصول بشكل مختلف باستخدام FindMotifs من Signac. تم إنشاء مخططات الأنماط باستخدام MotifPlot. تم حساب درجات إثراء الأنماط (درجات انحراف الأنماط) للأنماط الفردية لعوامل النسخ في الخلايا الفردية باستخدام chromVAR. لدمج درجات إثراء الأنماط والتعبير عن mRNA، تم مطابقة أسماء الأنماط في JASPAR2020 مع أسماء الجينات. تم إزالة الأنماط التي لم يكن من الممكن ربطها بشكل فريد مع الجينات من التحليل.

تنبؤ هدف عامل النسخ. للتنبؤ بالجينات التي يتم تنظيمها بواسطة عامل النسخ، تم استخراج مناطق التحكم المعزولة (ACRs) التي تحتوي على كل نمط لعامل النسخ. ثم تم تحديد الجينات التي تتوافق تعبيرها مع هذه المناطق. اعتبرنا الجينات الموجودة ضمن 5 كيلو بايت من منطقة التحكم المعزولة مع درجة ارتباط كمرشحين مهمين. تم إجراء تحليل إثراء GO لهؤلاء المرشحين لكل TF باستخدام دالة enrichGO من clusterProfiler مع توضيح org.At.tair.db. بالنسبة للتحليل في الشكل 2h، تم تطبيق عتبة أكثر صرامة (درجة الربط > 0.2). تم إنشاء مخططات إثراء GO باستخدام ggplot2 (المرجع 57).

تحليل تحت العنقود. تم تجميع النوى التي تحمل نفس علامة تحت العنقود، و تم حساب قيم النسخ المعدلة لكل مليون (TPM). تم إزالة المجموعات الفرعية التي تحتوي على أكثر من 18,000 جين غير مكتشف من التحليل. بالنسبة للمجموعات الفرعية التي اجتازت خطوة التصفية، تم تجميع الجينات باستخدام -تجميع المتوسطات مع (تم تحديده باستخدام طريقة الكوع) (الشكل البياني الإضافي 3ب). ثم تم إجراء تحليل إثراء GO للجينات في كل مجموعة (الشكل البياني الإضافي 3ب). أظهرت ثلاث مجموعات (1 و 5 و 8) إثراء لوظيفة متعلقة بالمناعة (استجابات لـ SA)؛ تم تعريف الجينات في هذه المجموعات على أنها ‘جينات مناعية محتملة’ واستخدمت للتحليل اللاحق.

تحليل الزمن الزائف. لحساب الزمن الزائف، تم الحصول على مصفوفة الخلايا حسب الجينات للخلايا الموصوفة بأنها نسيج ورقي من بيانات snMultiome. استخدمنا الدالة scanpy.tl.dpt مع n_dcs=2. تم بناء مسارات الزمن الزائف مع استخدام كل خلية نسيج ورقي في العينات الوهمية كخلية بداية. بعد ذلك، تم حساب متوسط هذه المسارات الفردية في كل خلية لإنشاء مسار زمن زائف موحد واحد. تم حساب اتجاهات الجينات في خريطة الحرارة لجين معين من خلال ملاءمة نموذج خطي GAM باستخدام دالة pygam LinearGAM مع (0، lam=400) لملاءمة نموذج GAM لمستويات التعبير عن RNA عبر قيم الزمن الزائف المرتبة.

بين المجموعات الفرعية من خلايا الميزوفيل النشطة مناعياً في بيانات snMultiome، تم تعريف مجموعات خلايا PRIMER وbystander بناءً على تعبير الجينات المؤشرة BON3 وALD1، على التوالي. من خلال مقارنة هذه الحالات الخلوية، تم تحديد الجينات المعبر عنها بشكل مختلف (DEGs) باستخدام وظيفة FindMarkers في Seurat. تم إجراء تحليل إثراء الأنماط على ACRs ضمن 2 كيلو بايت من DEGs الغنية بحالة الخلية.

لتقييم مساهمة CAMTA3 في حالات الخلايا المختلفة، تم أولاً تحديد الجينات التي تنظمها CAMTA3 باستخدام مجموعة بيانات RNA-seq الكمية المنشورة. ، مقارنة بين النوع البري و camta3-D (طافرة سلبية سائدة لـ CAMTA3). تم تداخل الجينات المثبطة بواسطة CAMTA3 بعد ETI (الذي تم تحفيزه بواسطة عدوى AvrRpm1 أو AvrRps4) مع DEGs الخاصة بحالة الخلية المحددة أعلاه. تم إجراء اختبارات هايبرجومتريكية لتقييم أهمية التداخلات.

تم إجراء تسلسل snRNA على نباتات gt3a-KO المصابة بـ AvrRpt2 أو المعالجة بالماء (شاهد) بعد 9 و 24 ساعة من الإصابة. تم إعداد نسختين مستقلتين لكل حالة عدوى. تم تصفية البيانات باستخدام نفس المعايير المستخدمة في تحليل snMultiome دون النظر في معلومات ATAC-seq. تم دمج بيانات -KO snRNA-seq مع بيانات Col-0 snMultiome (ظروف الشاهد و AvrRpt2)، وتم إجراء التجميع من جديد بنفس الطريقة الموضحة أعلاه. بناءً على تعبير الجينات المناعية، تم تحديد تجمعات الميزوفيل النشطة مناعياً وتم تقسيمها إلى تجمعات فرعية. من هذه التجمعات الفرعية، تم تعريف خلايا PRIMER وخلايا المتفرج بناءً على تعبير الجينات المؤشرة BON3 و ALD1، على التوالي. لكل حالة خلوية، تم إجراء تحليل DEG لمقارنة نباتات Col-O و gt3a-KO لتقييم تأثير طفرة GT-3A.

تمت مقارنة بيانات snATAC-seq الخاصة بنا مع بيانات ATAC-seq الكثيفة المنشورة للبالغين أوراق الثاليانا التي تنشط المناعة المستندة إلى النمط (PTI)، أو المناعة المستندة إلى العامل (ETI) أو كليهما PTI وETI، بالإضافة إلى الأوراق غير النشطة مناعياً تم دمج جميع مناطق الربط المعزولة (ACRs) المحددة في مجموعات بيانات ATAC-seq الكبيرة ومقارنتها بمناطق الربط المعزولة المحددة في مجموعات بيانات snATAC-seq الخاصة بنا.

تصميم لوحة MERFISH. قمنا بتجميع 500 جين مستهدف تضمنت الجينات التالية: (1) علامات محددة مسبقًا لأنواع خلايا أوراق A. thaliana. ; (2) الجينات المشاركة في عمليات مختلفة مثل المناعة، ومسارات الهرمونات، والتنظيم الجيني؛ و (3) مجموعة متنوعة من عوامل النسخ التي تم تحليلها سابقًا باستخدام DAP-seq (اختبار تفاعل عامل النسخ مع الحمض النووي) تم استبعاد الجينات التي لم يكن من الممكن تصميم أكثر من 25 مسبارًا محددًا لها بناءً على برنامج تصميم المسبار من Vizgen من لوحة الجينات المستهدفة. يمكن أن تتسبب الجينات ذات التعبير العالي في ازدحام إشارات smFISH وتعيق قياس MERFISH. لتجنب تضمين الجينات ذات التعبير العالي في اللوحة، قمنا بتقييم تعبير الجينات المستهدفة باستخدام مجموعة بيانات RNA-seq العامة المتاحة لجينوم A. thaliana المصاب بـ AvrRpt2 في نفس الإعداد الذي تم فيه الدراسة الحالية عند ثماني نقاط زمنية مختلفة. . تم استخدام أعلى قيمة تعبير بين النقاط الزمنية الثمانية لكل جين. الجينات التي أظهرت قيم TPM من لم يتم تضمينها في اللجنة. كان إجمالي TPM لـ 500 جين حوالي 22,000. تم استهداف ICS1 بجولة واحدة من smFISH حيث أن هذا الجين معبر عنه بشكل كبير. تم تقديم نتيجة smFISH كصورة دون معلومات كمية. تم تصور الخلايا البكتيرية من خلال استهداف 19 جينًا معبرًا عنه بشكل كبير (استنادًا إلى بيانات RNA-seq الجماعية المنشورة سابقًا في النباتات). ) كهدف واحد. تحتوي الجدول التكميلية 1 على قائمة بالجينات المستهدفة بواسطة MERFISH. تم تصميم جميع المجسات وبناؤها بواسطة Vizgen.

تقطيع الأنسجة، التثبيت والتركيب. تم زراعة النباتات وفقًا للطرق الموضحة أعلاه. تم قطع الأوراق المطابقة للعلاجات والنقاط الزمنية المذكورة أعلاه على الفور.
تم تحضينها وتكييفها في OCT (Fisher) لمدة 5 دقائق. بعد التحضين، تم تجميد الأوراق على الفور كما هو موصوف سابقًا. تم تكييف كتل الأنسجة إلى في غرفة كريوستات مبردة مسبقًا (لايكا) لمدة ساعة. تم تقليم كتل الأنسجة حتى تم الوصول إلى الأنسجة، بعد ذلك تم فحص الأقسام بصريًا حتى تم الكشف عن منطقة الاهتمام. تم إجراء تركيب العينة والتحضير وفقًا لدليل مستخدم MERSCOPE، ولكن مع بعض التعديلات الطفيفة. باختصار، تم إذابة القسم وتركيبه على شريحة MERSCOPE بدرجة حرارة الغرفة (Vizgen، 20400001)، ووضعها في طبق بتري وتم إعادة تجميده عن طريق الحضانة في غرفة التجميد لمدة 5 دقائق. تم تنفيذ الخطوات التالية مع العينات المثبتة في طبق بتري. ثم تم خبز العينات عند لمدة 5 دقائق وتم حضنها في محلول التثبيت بي بي إس و الفورمالديهايد) لمدة 15 دقيقة في درجة حرارة الغرفة. ثم تم غسل العينات بـ PBS يحتوي على 1:500 مثبط RNase (مثبط RNase Protector، ميلليبور سيغما) لمدة 5 دقائق في درجة حرارة الغرفة في ثلاث نسخ. بعد الغسل النهائي بـ PBS، تم تجفيف العينات عن طريق الحضانة في 70% إيثانول في بين عشية وضحاها.

تجربة MERFISH. تم معالجة مقاطع الأنسجة وفقًا لبروتوكول فيزجن. بعد إزالة 70% من الإيثانول، تم تحضين العينة في محلول غسل تحضير العينة (PN20300001) لمدة دقيقة واحدة ثم تم تحضينها في محلول غسل الفوراميد (PN20300002) عند لمدة 30 دقيقة. بعد إزالة محلول غسل الفوراميد، تم تحضين العينة في مزيج لوحة جينات MERSCOPE عند لمدة 42 ساعة. بعد تفاعل التهجين، تم غسل العينة مرتين باستخدام محلول الغسيل بالفورماميد في لمدة 30 دقيقة ومرة واحدة مع محلول غسل تحضير العينة في درجة حرارة الغرفة لمدة دقيقتين. بعد خطوة الغسل، تم تضمين العينة في الهيدروجيل عن طريق الحضانة في محلول تضمين الجل (خليط تضمين الجل (PN20300004)، محلول بيرسلفات الأمونيوم و -تترا ميثيل إيثيلين ديامين) في درجة حرارة الغرفة لمدة 1.5 ساعة. ثم تم تنظيف العينة عن طريق أولاً حاضنتها في مزيج الهضم (مزيج الهضم (PN20300005) ومثبط RNase بنسبة 1:40) في درجة حرارة الغرفة لمدة ساعتين، تلاها الحضانة في محلول التنظيف (مزيج التنظيف (PN20300003) وبروتيناز K) عند لمدة 24 ساعة ثم عند لمدة 24 ساعة. تم غسل العينة المنظفة مرتين باستخدام محلول غسل تحضير العينة وتم صبغها بـ DAPI و PolyT في درجة حرارة الغرفة لمدة 15 دقيقة. ثم تم غسل العينات بمحلول غسل الفوراميد في درجة حرارة الغرفة لمدة 10 دقائق وشطفها بمحلول غسل تحضير العينة. تم تصوير العينة باستخدام جهاز MERSCOPE، وتم فك تشفير النسخ المكتشفة على جهاز MERSCOPE باستخدام كتاب رموز تم إنشاؤه بواسطة Vizgen. تم تصور النسخ باستخدام أداة التصور من Vizgen.

تقسيم ومعالجة MERFISH. تم إجراء تقسيم حدود الخلايا لكل مخرج بيانات MERSCOPE. أظهرت المجسات المستهدفة لـ DAPI والمستهدفة لـ poly(A) نجاحًا متغيرًا في صبغ النوى والسيتوبلازم، على التوالي، اعتمادًا على العينات ومناطق الأنسجة التي تم فحصها (الشكل 3ب والشكل الممتد 7هـ يظهران تقسيمًا فاشلاً وناجحًا، على التوالي). بالمقابل، كانت المناطق الكثيفة من النسخ تشير إلى مواقع النواة بشكل أكثر قوة (الشكل الممتد 7هـ). لذلك، تم استخدام طريقة تقسيم قائمة على النسخ. لكل عينة، تم إنشاء مصفوفة Numpy ثنائية الأبعاد من الأصفار، والتي نمذجة إجمالي مساحة البكسل المصورة. تم تغيير إحداثيات النسخ RNA المحددة من 0 إلى 1 في هذه المصفوفة. بعد ذلك، تم تشويش المصفوفة باستخدام دالة cv2.GaussianBlur في OpenCV مع ksize=(5,5). تم تقسيم المصفوفة الناتجة إلى مناطق بكسل . تم تحميل هذه المناطق في Cellpose ، أداة تقسيم تعتمد على التعلم العميق، وتم تدريب نموذج تقسيم مخصص عن طريق تقسيم يدوي لأجسام النواة عبر عشرة مناطق بكسل في واجهة مستخدم Cellpose.

ثم تم استخدام النموذج المخصص للتنبؤ بحدود التقسيم في جميع المناطق المتبقية، مع المعلمات diameter=22.92، flow_threshold و cell probability threshold . بعد ذلك، تم حساب العدد الإجمالي للخلايا في تجربة عن طريق جمع
عدد الخلايا الفريدة عبر جميع المناطق. تم استخدام هذا الرقم بعد ذلك لتهيئة الأمر –num-cells-init في أداة تقسيم أخرى، Baysor ، التي تأخذ في الاعتبار الاحتمالية المشتركة لتكوين النسخ وشكل الخلية للتنبؤ بحدود الخلايا. تم تشغيل Baysor باستخدام صورة Docker تم تنزيلها والمعلمات- n clusters 1، -i1، –force-2d، min-molecules-per-gene=1، min-molecules-per-cell=50، scale ، scale-std=” “، estimate-scale-from-centers=true، min-molecules-per-segment ، new-component-weight ، new-component-fraction=0.3.
لاختبار جودة طريقة تقسيمنا القائمة على النسخ، استخدمنا FOV مع صبغ DAPI ناجح (وهو نادر في عيناتنا) وأجرينا تقسيمًا قائمًا على DAPI وتقسيمًا قائمًا على النسخ (الشكل الممتد 7د). كانت النتائج من هاتين الاستراتيجيتين للتقسيم متوافقة مع بعضها البعض بشكل عام، حيث كانت الطريقة القائمة على النسخ تلتقط النسخ في السيتوبلازم بالإضافة إلى تلك الموجودة في النواة (الشكل الممتد 7د-ف). تشير هذه النتيجة إلى أن نهج التقسيم لدينا يمكنه التقاط الخلايا بشكل موثوق.
بعد تقسيم Baysor، تم إنشاء مصفوفة خلية-بجين من تعيينات خلايا النسخ. تم إزالة الخلايا التي تحتوي على أقل من 50 نسخة معينة. تم استخدام Scanpy للمعالجة اللاحقة لتجارب MERFISH الخاصة بنا. بعد تحميل مصفوفة الخلية-بجين المعنية في كائن Anndata لكل تجربة، قمنا بتخزين الإحداثيات المكانية لكل خلية تم الحصول عليها من Baysor. تم تطبيع عدد النسخ الفردية في كل خلية بواسطة إجمالي عدد النسخ لكل خلية. ثم تم تحويل مصفوفة الخلية-بجين Anndata إلى مقياس لوغاريتمي.

دمج بيانات MERFISH-snMultiome ونقل العلامات. لدمج تجارب MERFISH مع بعضها البعض، استخدمنا FindIntegrationAnchors تليها دالة IntegrateData في Seurat (v.5). لدمج بيانات تجربة MERFISH مع بيانات snMultiome الخاصة بنا، قمنا بدمج النقاط الزمنية المقابلة في كل نمط بشكل منفصل. لكل نقطة زمنية في بيانات العدوى (وهمية، 4، 6، 9 و 24 ساعة)، استخدمنا gimVI من scvi-tools لإسقاط كل من خلايا snMultiome و MERFISH في نفس الفضاء الكامن باستخدام أعداد تعبير RNA الخاصة بهم. تم تدريب gimVI على مدى 200 حقبة مع حجم 10 فضاء كامن لكل نقطة زمنية. لنقل الملاحظات المستمرة (الزمن الزائف، أعداد نشاط جينات RNA و ATAC، وغنى النمط) من بيانات snMultiome إلى بيانات MERFISH، قمنا بتعيين المتوسط العددي لأقرب 30 خلية snMultiome في الفضاء الكامن gimVI لكل خلية MERFISH. بالمثل، لنقل الملاحظات الفئوية، بما في ذلك أنواع الخلايا وتسميات المجموعات، قمنا بتعيين كل خلية MERFISH لأكثر تسمية شيوعًا في مجموعة أقرب 30 خلية snMultiome لها. لرسم التضمين المشترك لكلا النمطين (MERFISH و snMultiome) لنقطة زمنية واحدة، قمنا بتشغيل UMAP على الفضاء الكامن gimVI مع n_neighbors و min_dist . لاحظ أنه في التخطيط المكاني لأنواع الخلايا المتوقعة في snMultiome (الشكل 4ج)، على الرغم من أن المقطع كان في منتصف الورقة، تم تضمين بعض خلايا البشرة لأن المقطع لم يكن مسطحًا تمامًا.
تقدير smFISH وتحديد مستعمرات البكتيريا. تم إجراء تقدير للنسخ المسمى بواسطة smFISH باستخدام حزمة Python Big-FISH. تم إسقاط سبعة طائرات من صور smFISH الخاصة بـ MERSCOPE في بعدين باستخدام numpy.max على المحور وتم تقسيمها إلى مناطق . ثم تم استدعاء النقاط باستخدام دالة bigfish.detection.detect_spots مع العتبة ، spot_radius= و voxel_size . لتحديد مستعمرات البكتيريا، تم استدعاء bigfish.detection.detect_spots لتسمية مواقع ‘جين ميتا بكتيري’ مع العتبة ، log_kernel_size و minimum_distance=(1.456,1.456). لأخذ في الاعتبار الفلورية الذاتية من نسيج النبات، تم استخدام نفس مستدعي النقاط لاستدعاء النقاط على قنوات DAPI و ICS1. تم تجميع النقاط من جميع القنوات الثلاث، وتم استخدام DBSCAN من scikit-learn.cluster على جميع النقاط مع eps و min_samples=5. احتفظنا بجميع مجموعات DBSCAN التي تشكلت فيها نقاط DAPI و ICS1 أقل من من إجمالي نقاط المجموعة.

تمثل هذه المجموعات المتبقية مستعمرات بكتيرية محتملة. قمنا بتقييم كل مجموعة يدويًا لدمج تلك التي تشير إلى نفس المستعمرة وإزالة المجموعات التي تشير إلى الفلورية الذاتية الواضحة، على سبيل المثال، الإشارات من الثغور. تم إنشاء نوافذ من بكسل حول كل مجموعة نهائية، وتم استخدام detect_spots مع threshold=95، spot_radius=(17,17) و voxel_size=(3,3) للتقدير الدقيق للبكتيريا الفردية لكل مجموعة.

تحليل جيرة البكتيريا. لتحديد مستوى المناعة في الجوار حول كل مستعمرة بكتيرية، حددنا أقرب 1-1,000 خلية في القرب من المستعمرة. ثم تم حساب القيم الزمنية الزائفة الملساء المستنتجة لهذه الخلايا المجاورة للحصول على قيمة زمنية زائفة متوسطة واحدة لكل حجم جار. تم حساب منحنيات الخطأ باستخدام الخطأ القياسي مقسومًا على المتوسط عند كل حجم جار. تعتبر هذه القيم مؤشرات لمستوى المناعة العام للمنطقة حول كل مستعمرة.

تحليل التعبير التفاضلي المكاني. لتحديد الجينات المرتبطة مكانيًا بالمناطق الغنية بالمناعة (المتعلقة بالشكل 5أ)، قمنا أولاً بتنعيم زمننا الزائف المكاني المستنتج على أقرب 100 خلية مكانيًا للعثور على مناطق النقاط الساخنة النشطة مناعيًا. ثم استخدمنا RCTD لفك تشفير الثنائيات في بياناتنا المكانية باستخدام أنواع الخلايا في بيانات snMultiome الخاصة بنا كمرجع. أنشأنا كائن مرجعي مع وقمنا بتشغيل RCTD على بيانات MERFISH الخاصة بنا مع المعلمات gene_cutoff ، gene_cutoff_reg ، fc_cutoff ، fc_cutoff_reg و doublet_mode ‘doublet’. ثم قمنا بتشغيل C-SIDE من spacexr مع run.CSIDE.single مع cell_type_threshold لإيجاد أعلى الجينات المعبر عنها بشكل مكاني مختلف لكل نوع خلية على طول الزمن الزائف المكاني الملسّح (المتغير التفسيري) ورسمت اللوغاريتم السالب المحوّل القيمة مقابل التغير المضاعف المحول لوغاريتمياً.

معلومات إضافية حول تصميم البحث متاحة في ملخص تقارير مجموعة ناتشر المرتبط بهذه المقالة.

جميع المعلومات التي تدعم الاستنتاجات متوفرة مع الورقة. تم إيداع بيانات تسلسل الخلايا المفردة والكمية التي تم إنشاؤها في هذه الدراسة في قاعدة بيانات المعهد الوطني للمعلومات الحيوية الخاصة بالتعبير الجيني (أرقام الوصول GSE226826 و GSE248054). تتوفر الجينومات المرجعية والتعليقات التوضيحية والبيانات المعالجة بالكامل لتحليلات snMultiome من موقع Salk (http:// neomorph.salk.edu/download/Nobori_etal_merfish). بيانات MERFISH متاحة على موقع سالك (http://plantpathogenatlas. سالک.إديوتم إدراج الجينات المستهدفة باستخدام MERFISH في الجدول التكميلي 1. تتوفر معلومات حول البرايمرات المستخدمة في هذه الدراسة في الجدول التكميلي 2. تم توفير بيانات المصدر مع هذه الورقة.

الرموز لتحليل بيانات snMultiome و MERFISH متاحة على GitHub (https://github.com/tnobori/snMultiome و https://github. com/amonell/Spatial_Plant_Pathogen_Atlas).
40. Mine، أ. وآخرون. شبكة إشارات هرمونات الدفاع النباتية تمكّن إعادة برمجة النسخ بسرعة وبشدة عالية خلال المناعة المستحثة بواسطة المؤثرات. خلية النبات 30، 1199-1219 (2018).
41. كايدو، م.، تسونو، ي.، ميس، ك. وأوكونو، ت. استهداف الشبكة الإندوبلازمية لبروتين حركة فيروس الفسيفساء النخرية لزهرة البرسيم الأحمر مرتبط بتكرار RNA1 الفيروسي ولكن ليس RNA2. علم الفيروسات 395، 232-242 (2009).
42. Mine، أ. وآخرون. الأدوار المختلفة لـ Hsp 70 و Hsp 90 في تجميع مركب النسخ لفيروس RNA النباتي ذو السلسلة الموجبة. J. Virol. 86، 12091-12104 (2012).
43. زانغ، إكس.، هنريكس، ر.، لين، س.-س.، نيو، كيو.-و. وتشوا، ن.-هـ. تحويل أريابيدوبسيس ثاليانا بواسطة بكتيريا الأجروبكتيريوم باستخدام طريقة الغمر الزهري. بروتوكولات الطبيعة 1، 641-646 (2006).
44. مatsu ماتو، أ. وآخرون. أداة وسم حيوية تعتمد على ناقل Tn7 متعددة الاستخدامات للكشف الكمي والمكاني عن البكتيريا في النباتات. اتصالات النبات. 3، 100227 (2022).
45. تاونزلي، ب. ت.، كوفينغتون، م. ف.، إيتشيهاشي، ي.، زومشتاين، ك. وسينها، ن. ر. BrAD-seq: تسلسل توجيهي لمهايئ التنفس: بروتوكول تحضير مكتبة مبسط وسريع للغاية لبناء مكتبة mRNA محددة الاتجاه. فرونت. علوم النبات 6، 366 (2015).
46. تشين، س.، تشو، ي.، تشين، ي. وغو، ج. fastp: معالج مسبق FASTQ فائق السرعة شامل. المعلوماتية الحيوية 34، i884-i890 (2018).
47. دوبين، أ. وآخرون. ستار: محاذي تسلسل RNA عالمي فائق السرعة. المعلوماتية الحيوية 29، 15-21 (2013).
48. لياو، ي.، سميث، ج. ك. وشي، و. featureCounts: برنامج عام فعال لتعيين قراءات التسلسل إلى الميزات الجينومية. المعلوماتية الحيوية 30، 923-930 (2014).
49. هاو، ي. وآخرون. تحليل متكامل لبيانات الخلايا الفردية متعددة الأنماط. خلية 184، 3573-3587 (2021).
50. ستيوارت، ت.، سريفاستافا، أ.، ماداد، س.، لاريه، ك. أ. وساتيجا، ر. تحليل حالة الكروماتين على مستوى الخلية الواحدة باستخدام Signac. نات. ميثودز 18، 1333-1341 (2021).
51. مكغينيس، سي. إس.، مورو، إل. إم. وغارتنر، ز. ج. دابلت فايندر: اكتشاف الدوبلت في بيانات تسلسل RNA أحادي الخلية باستخدام الجيران الأقرب الاصطناعيين. نظام الخلية 8، 329-337 (2019).
52. ماراند، أ. ب.، تشين، ز.، جالافوتي، أ. وشميتز، ر. ج. أطلس تنظيمي في الذرة بدقة خلوية فردية. خلية 184، 3041-3055 (2021).
53. زانغ، ي. وآخرون. تحليل قائم على النموذج لبيانات ChIP-seq (MACS). جينوم بيو. 9، R137 (2008).
54. كورسونسكي، إ. وآخرون. دمج بيانات الخلايا المفردة بسرعة وحساسية ودقة مع هارموني. نات. ميثودزhttps://doi.org/10.1038/s41592-019-0619-0 (2019).
55. كوزانوفيتش، د. أ. وآخرون. تحليل الوصول إلى الكروماتين على مستوى الخلية الواحدة باستخدام الفهرسة الخلوية التوافقية. ساينس 348، 910-914 (2015).
56. فورنس، أ. وآخرون. JASPAR 2020: تحديث قاعدة البيانات المفتوحة الوصول لملفات ارتباط عوامل النسخ. أبحاث الأحماض النووية. 48، D87-D92 (2020).
57. ويكهام، هـ. ggplot2: رسومات أنيقة لتحليل البيانات (سبرينجر، 2009).
58. بروكو، سي. وآخرون. تحليل التعبير الجيني الخلوي المحدد في الأوراق والتعبير الجيني على مستوى الخلية الواحدة يكشف عن دور طبقة العمود الفقري في حماية الأشعة فوق البنفسجية. خلية النباتhttps://doi.org/10.1093/plcell/koac167 (2022).
59. أومالي، ر. س. وآخرون. ميزات السيتروما والإبيسيتروما تشكل المشهد التنظيمي للحمض النووي. خلية 165، 1280-1292 (2016).
60. نوبوري، ت. وآخرون. مشهد النسخ لجهاز بكتيري ممرض تحت مناعة النبات. وقائع الأكاديمية الوطنية للعلوم في الولايات المتحدة الأمريكية 115، E3055-E3064 (2018).
61. جيكاميلو، س. ولوندبرغ، ج. إعداد أنسجة النباتات لتمكين تحليل النسخ المكاني باستخدام المصفوفات الدقيقة المشفرة. نات. بروتوك. 13، 2425-2446 (2018).
62. باتشيتاريو، م. وسترينجر، س. Cellpose 2.0: كيفية تدريب نموذجك الخاص. نات. ميثودز 19، 1634-1641 (2022).
63. بتوكوف، ف. وآخرون. تقسيم الخلايا في النسخ الجيني المكاني القائم على التصوير. نات. بيولوجيا التكنولوجيا. 40، 345-354 (2022).
64. كابل، د. م. وآخرون. تحليل قوي لمزيج أنواع الخلايا في النسخ الجيني المكاني. نات. بيوتكنولوج.https://doi.org/10.1038/s41587-021-00830-w (2021).

الشكر والتقدير نشكر ك. تسودا وج. ووكر على تعليقاتهم حول المخطوطة. تم دعم ت.ن. من خلال زمالة طويلة الأمد من برنامج العلوم الحدودية البشرية (HFSP) (LTOO0661/2020-L). ج.ر.إ. هو باحث في معهد هوارد هيوز الطبي.

مساهمات المؤلفين: قام ت.ن. بتصميم الدراسة والتجارب بمساعدة من ج.ر.إ. قام ت.ن. بتنفيذ تجارب snMultiome و MERFISH وتحليل البيانات. قدم ت.أ.ل. مقاطع الأنسجة لتجارب MERFISH. قام أ. مونييل بتحليل بيانات MERFISH ودمج البيانات بمساعدة من ج.ز. قام ي.س. و س.س. و أ. مين بإنتاج طفرات GT-3A وأجروا اختبارات نمو مسببات الأمراض. ساعد ج.ر.ن. في التسلسل وإدارة البيانات. كتب ت.ن. المسودة الأولية للبحث. قام ت.ن. و ج.ر.إ. بتحرير البحث.

المصالح المتنافسة: يشغل J.R.E. عضوية في المجلس الاستشاري العلمي لشركة Zymo Research وشركة Cibus. يعلن المؤلفون الآخرون عدم وجود مصالح متنافسة.

معلومات إضافية النسخة الإلكترونية تحتوي على مواد إضافية متاحة فيhttps://doi.org/10.1038/s41586-024-08383-z.
يجب توجيه المراسلات والطلبات للحصول على المواد إلى جوزيف ر. إكر. تشكر مجلة نيتشر غيتا كوكير، شيوهوا تشونغ والمراجعين الآخرين المجهولين على مساهمتهم في مراجعة الأقران لهذا العمل. تقارير مراجعي الأقران متاحة.
معلومات إعادة الطبع والتصاريح متاحة علىhttp://www.nature.com/reprints.

الشكل البياني الممتد 1| مراقبة جودة بيانات snMultiome. أ، مخططات الكمان تظهر عدد النسخ لكل نواة، وعدد الجينات لكل نواة، وقراءات ATAC لكل نواة، ومناطق الكروماتين القابلة للوصول (ACRs) لكل نواة في كل عينة. ب، عدد الخلايا لكل عينة. تم تحليل ودمج نسختين مستقلتين لظروف Mock وAvrRpt2 بعد 9 ساعات من الإصابة. ج، مخططات تشتت الكثافة تظهر درجة قراءة النسبة في القمم (FRiP) في كل عينة قبل تصفية الخلايا. المحور السيني: تحولت أعماق القراءة. المحور الصادي: نسبة مواقع تكامل Tn5 في ACRs. د، مخططات تشتت الكثافة لـ تم تحويل أعماق القراءة (محور س) بواسطة نسبة مواقع تكامل Tn 5 التي تتواجد ضمن 2 كيلو بايت upstream و1 كيلو بايت downstream من مواقع بدء النسخ (TSSs) (محور ص). البيانات في كل عينة قبل تصفية الخلايا موضحة. ج-د، تم دمج نسختين من ظروف Mock وAvrRpt29. هـ، رسم بياني عمودي يوضح عدد قمم ATAC-seq التي تم تحديدها في بيانات ATAC-seq السابقة وفي بيانات snATAC-seq الحالية. القمم المشتركة والفريدة من نوعها للاختبار موضحة باللون الأزرق والأحمر، على التوالي.
f، الوصول إلى الكروماتين المجمّع حول ACTIN2 (يسار) و ICS1 (يمين). g، تحليل المكونات الرئيسية للترانسكريبتوم الزائف لكل عينة. تم وضع علامات على النسخ المستقلة من Mock و AvrRpt2 بعد 9 ساعات. h، مخططات مبعثرة تقارن الترانسكريبتومات الزائفة لعينات Mock و AvrRpt2 بعد 9 ساعات. تم عرض قيم معامل ارتباط بيرسون. i، مخططات شريطية مكدسة تظهر تمثيل مجموعات Leiden المستندة إلى تعبير الجينات في كل عينة. j، تضمين ثنائي الأبعاد لتشابه الوصول إلى الكروماتين بين النوى من جميع العينات باستخدام تقريب متعدد الأشكال الموحد والإسقاط (UMAP). يتم تلوين النوى حسب مجموعات Leiden. k، تضمينات UMAP بناءً على رسم بياني للجيران المشترك الذي يمثل كل من قياسات تعبير الجينات والوصول إلى الكروماتين. يتم تلوين النوى حسب مجموعات Leiden الجديدة بناءً على التحليل المشترك (يسار) ومجموعات Leiden المحددة فقط من خلال قياس تعبير الجينات (يمين؛ الشكل 1b).

أ

الشكل البياني الممتد 2 | تحليل التعبير الجيني والوصول إلى الكروماتين على مستوى الخلية الواحدة. أ، رسم نقطي يظهر أهم جينات العلامة في المجموعات الفردية. ب، الوصول إلى الكروماتين المجمّع حول المجموعات. ، جين علامة معروف لجلد الورقة. ج، رسم UMAP يوضح عدد ATAC-seq
على قمة بالقرب من FDH (الكروموسوم 2، الموقع 11172821-11173529) في كل نواة. د، تحليل إثراء GO لجينات العلامة لكل مجموعة. هـ، تعبير mRNA لـ ICS1 في كل نواة في كل عينة. تم عرض القيم p المعدلة من اختبار هايبرجومتريك أحادي الجانب متبوعًا بتصحيح بنجاميني-هوشبرغ.

الشكل التمديدي 3 | توصيف شامل لمجموعات الخلايا المتميزة. أ، خريطة حرارية تظهر الارتباط الثنائي بين النسخ الجينية الزائفة بين المجموعات الفرعية لمجموعات فردية. تظهر الأشرطة العلوية والجانبية تسميات المجموعات الرئيسية. ب، (يسار) مخطط سير العمل لاختيار الجينات المتعلقة بالمناعة ذات التباين العالي. أولاً، تم اختيار الجينات التي كانت غنية بشكل ملحوظ في واحدة على الأقل من المجموعات الرئيسية أو الفرعية. ثم، تم تجميع هذه الجينات بناءً على تعبير النسخ الزائفة للمجموعات الفرعية باستخدام تجميع k-mean (تم تحديد قيمة k بواسطة طريقة الكوع)، تلاها تحليل إثراء GO لكل مجموعة k-mean. أخيرًا، تم اختيار مجموعات الجينات التي تحتوي على مصطلحات GO المتعلقة بالمناعة. (يمين) أهم مصطلحات GO الغنية في 12 مجموعة k-mean. تم اختيار المجموعات 1 و2 و10 كجينات مناعية. ج، خريطة حرارية تظهر التعبير المنظم للجينات المتعلقة بالمناعة المختارة في (ب) عبر جميع المجموعات الفرعية. تشير الأشرطة العلوية إلى نوع الخلية والمجموعة الرئيسية التي تنحدر منها كل مجموعة فرعية. د، تعبير اللحاء.
علامات خلايا الرفيق SUC2 وFTIP1. e، تعبير ALD1 وFMO1. f، تعبير ILL6 عند الإصابة بـ AvrRpt2 على مدى الزمن، مما يظهر تحفيزًا محددًا في العنقود 6. g، تعبير الجينات المعبر عنها بشكل خاص في تحت العنقود 6-8. h، خريطة حرارية تظهر تعبير الجينات المتعلقة بالمناعة الموضحة في الشكل 1d في دراسة سابقة لتتابع RNA-seq الكمي، حيث تمت إصابة أوراق .thaliana بواسطة DC3000 أو AvrRpt2 أو AvrRpm1 أو تم معالجتها بالتحكم الوهمي (الماء). لا تلتقط هذه البيانات وجود مجموعات الخلايا المختلفة المحددة في الشكل 1d.i، مخططات التشتت التي تظهر العلاقة بين CRK12 و CRK41 في تحليل تسلسل RNA الكمي على مدى الزمن (الأعلى) وتسلسل RNA أحادي النواة (snRNA-seq؛ الأسفل). بالنسبة لتحليل snRNA-seq، تم استخدام بيانات التعبير الجيني الزائفة لمجموعات فرعية (موضحة في c). تعبير CRK12 و CRK41 عند الإصابة بـ AvrRpm1 كما هو موضح في UMAP. هذه الجينات المرتبطة ظاهريًا في تسلسل RNA الكمي أظهرت نمط تعبير مختلف للغاية على مستوى الخلية الفردية.

الشكل الممتد 4 | ربط التعبير الجيني والوصول إلى الكروماتين على مستوى الخلية الواحدة. أ، مخطط توضيحي لمناطق الكروماتين القابلة للوصول (ACRs) وجين. يتم “ربط” ACR وجين عندما يكون هناك ارتباط كبير بين وصول الكروماتين وتعبير mRNA عبر خلايا فردية. ب، مخطط كثافة يظهر تكرار الروابط على مسافات مختلفة من مواقع بدء النسخ (TSSs). ج، خريطة حرارية تظهر درجة الربط (معامل الارتباط بيرسون بين عدد ACR وتعبير mRNA) للجينات التي أظهرت على الأقل رابطًا واحدًا كبيرًا في واحدة على الأقل من ظروف العدوى أو الشاهد. عندما كان للجين روابط متعددة، تم عرض أعلى درجة ربط. الشريط الجانبي يظهر تجمع k-mean.
التعليق (الـ تم تحديد القيمة بواسطة طريقة الكوع). د، تعبير mRNA الذي يشفر VSP2 و UGT85A1 و MAM1. هـ، و، الوصول إلى الكروماتين المجمّع حول CBP6Og (هـ) و AT1G56660 (و). تظهر الرسوم البيانية على الجانب التعبير المجمّع لـ mRNA لكل جين. تم التعبير عن كلا الجينين بشكل مرتفع بطريقة محددة للخلايا، لكن فقط CPB6Og أظهر أنماط وصول كروماتين مرتبطة (مرتبطة). ز، ح، الأنماط العليا الغنية في مناطق المحفز (2 كيلو بايت أعلى من TSS) لجينات الدفاع (علامات الكتل النشطة مناعياً في الميزوفيل والجلد). و 12) التي مرتبطة (g) وغير مرتبطة (h).

الشكل البياني الممتد 5 | الروابط بين النسخ الجيني والوصول إلى الكروماتين. أ، مخطط مبعثر يظهر توزيع درجة الربط (معامل الارتباط بيرسون بين عدد ACR والتعبير عن mRNA) على مسافات مختلفة من TSSs. ب، تحليل إثراء GO للجينات في كل مجموعة k-mean الموضحة في (الشكل البياني الممتد 4c). يتم عرض القيم المعدلة p من اختبار هايبرجومتري أحادي الجانب تليها تصحيح بنجاميني-هوشبرغ. ج، مخطط كثافة يظهر تكرار الروابط على مسافات مختلفة من TSS لجينات علامات المجموعات والجينات غير العلامات. د، مخطط كثافة يظهر تكرار حجم ACRs للجينات المرتبطة أو غير المرتبطة.
e، رسم بياني للكثافة يظهر تكرار الجينات مع أعداد مختلفة من ACRs المرتبطة. f، مخطط تخطيطي لتحليل نشاط ATAC. لكل جين، تم تجميع قراءات ATAC المرسومة على جسم الجين أو المنطقة التي تبعد 400 نقطة أساسية للأعلى لحساب الدرجة. g، رسم بياني مبعثر يظهر العلاقة بين عدد الروابط (محور x) ودرجة ارتباط نشاط RNA-ATAC (محور y). h، رسم بياني مبعثر يظهر العلاقة بين أقصى معامل ارتباط بين كل جين وACRs المرتبطة ( -المحور) ودرجة ارتباط نشاط RNA-ATAC ( -محور).

يظهر (أ) درجات الإثراء لـ 465 نمطًا و(ب) تعبير عوامل النسخ المقابلة (TF) في كل نواة. تُظهر الأشرطة العلوية توضيح التجمع المحدد في الشكل 1ب. ج، خريطة حرارية تُظهر معامل ارتباط بيرسون بين
درجات إثراء النمط والتعبير عن mRNA لعوامل النسخ المقابلة في كل نوع من الخلايا (موضحة في الشريط العلوي). جميع عوامل النسخ موضحة. د، تحليل إثراء GO للجينات المتوقعة أن تكون مستهدفة بواسطة عوامل النسخ الموضحة في الشكل 2h.

عينة (# النصوص المكتشفة)

البيانات الموسعة الشكل 7|تحليل بيانات MERFISH. أ، صور خام مثال لMERFISH التي تلتقط نفس منطقة الأنسجة عبر ثلاث جولات تصوير. النقاط البيضاء هي إشارات مشتقة من جزيئات mRNA الفردية. لاحظنا نتائج مماثلة في جميع العينات المستقلة. ب، مخططات ثنائية الأبعاد لجميع النسخ المتماثلة التي تم اكتشافها بواسطة MERFISH في كل عينة. يتم عرض عدد النسخ المتماثلة. ج، نمط التعبير المكاني لـ ALD1 الذي تم اكتشافه بواسطة MERFISH في كل عينة. د، (يسار) مجال رؤية (FOV) يظهر إشارة تلوين نوى DAPI غير واضحة. (وسط) تقسيم قائم على النسخ في نفس مجال الرؤية (انظر الطريقة). تم تلوين النسخ المتماثلة حسب الخلايا المعينة. (يمين) مراكز الخلايا المكتشفة.
من خلال تقسيم البيانات المستندة إلى النص (النقاط الحمراء) ونتيجة تقسيم البيانات المستندة إلى صبغة DAPI الفاشلة (المنطقة الصفراء). لوحظت أنماط مماثلة عبر مجالات الرؤية والعينات. تم تقديم تحليل كمي منهجي في (f). (a,d) قضبان القياس م.ع. , انخفضت نسبة النسخ في الخلايا الم segmented بواسطة Cellpose ضمن الخلايا الم segmented بواسطة Baysor. (e) تم استخدام الخلايا المكتشفة في مجال الرؤية (FOV) الموضح في (d). (f) تم استخدام الخلايا المكتشفة في 100 FOV. g، UMAP المتكامل لجميع بيانات MERFISH. h، رسم خرائط مكاني لمجموعات MERFISH الجديدة المعلنة كأوعية دموية. b، c، h القضبان المقياس .

الشكل البياني الموسع 8 | رسم خرائط مكاني للترنسكبتوم الكامل والإيبيجينوم. a، التعبير عن mRNA لـ ALD1 (الأعلى) توقع نمط التعبير عن mRNA الذي تم قياسه باستخدام MERFISH (الأسفل). تظهر الصور المكبرة على اليمين. b، التعبير عن mRNA لـ ICS1 (الأعلى) توقع نمط التعبير عن mRNA الذي تم قياسه باستخدام smFISH (الأسفل). تظهر الصور المكبرة على اليمين. c، نشاط ATAC المتوقع (الشكل البياني الموسع 5f) لـ ICS1 (الأعلى) و ALD1 (الأسفل)، والذي أظهر أنماط متسقة مع التعبير عن mRNA (a، b). d، كانت درجات إثراء النمط المتوقع لـ HsfB2b (الأعلى) متسقة مع التعبير عن mRNA لـ تم تأكيده بواسطة MERFISH (الأسفل). e، رسم خرائط مكاني لقيم الزمن الزائف بناءً على دمج البيانات ونقل العلامات بين snRNA-seq وبيانات MERFISH لنباتات مصابة بـ DC3000.
f، smFISH لـ ICS1 في أوراق مصابة بـ AvrRpt2 أو DC3000 بعد 24 ساعة من الإصابة. g، رسم خرائط مكاني للترنسكبتات البكتيرية المكتشفة باستخدام smFISH في نباتات مصابة بـ AvrRpt2 (يسار) و DC3000 (يمين) بعد 24 ساعة من الإصابة (hpi). يتم أيضًا تصور قيم الزمن الزائف في الخلفية. تعكس حجم النقاط عدد الترنسكبتات البكتيرية المكتشفة. a-g، القضبان المقياس . h، رسم بياني يوضح عدد الترنسكبتات البكتيرية (محور x) ومتوسط قيم الزمن الزائف لخلايا النبات الخمس الأقرب ( -محور) لكل مستعمرة بكتيرية في نباتات مصابة بـ AvrRpt2 (أزرق) و DC3000 (أحمر) بعد 24 ساعة من الإصابة. i، متوسط درجات الزمن الزائف للخلايا المحيطة بكل مستعمرة بكتيرية. يظهر محور x عدد خلايا النبات الأقرب التي تم تحليلها لكل مستعمرة بكتيرية. تشير الأشرطة المظللة إلى الخطأ القياسي للمتوسط.

الشكل البياني الموسع 9 | توصيف خلايا PRIMER والخلايا المجاورة.
a، خرائط الحرارة التي تظهر التعبير عن ALD1 و FMO1 و BON3 في كل خلية بعد 9 ساعات من الإصابة. b، التعبير عن يظهر في UMAP في الأشكال 4a و 5d. يتم إثراء تعبير BON3 في الخلايا ذات أعلى درجات الزمن الزائف. c، التعبير عن WRKY8 و LSD1 في المجموعات الفرعية من المجموعات 3 و 7 و 11 في بيانات snRNA-seq (الشكل 1f). d، رسم بياني على شكل كمان يظهر التعبير عن PUB36 في مجموعة خلايا PRIMER من نباتات WT و GT-3A knockout (gt3a-KO). يتم عرض نموذج الجين لـ PUB36، قمم ATAC-seq، ونمط ربط GT-3A أدناه. e، تحليل إثراء GO للجينات التي تم تقليلها في مجموعة خلايا مجاورة من gt3a-KO مقارنة بـ WT. يتم عرض قيم p المعدلة من اختبار هايبرجومتري أحادي الجانب تليها تصحيح بنجاميني-هوشبرغ. f، رسم بياني على شكل كمان يظهر التعبير عن ALD1 في مجموعة مجاورة. d-f، تم تحليل نسختين مستقلتين من كل عينة باستخدام RNA-seq أحادي النواة. g، التعبير عن BON3 و GT-3A في خلايا في حالة وهمية (الأعلى)
أو حالة الإصابة بـ DC3000 (الأسفل). h، التعبير عن /CS1 في خلايا مصابة بـ AvrRpt2. c، g، h، تم دمج جميع النقاط الزمنية. i، نشاط النمط لـ GT-3A في خلايا الميزوفيل النشطة مناعياً. j، منطقة الإصابة الناتجة عن عدوى Colletotrichum higginsianum في Col-0 وسطرين مستقلين من زيادة التعبير عن GT-3A وطفرة pad3 بعد 6 أيام من التلقيح. تشير الحروف المختلفة إلى الأهمية الإحصائية (المعدلة ). تم حساب قيم p المعدلة باستخدام اختبار t لطلاب ذو طرفين تليه طريقة بنجاميني-هوشبرغ. , و 51 ورقة لـ WT و GT-3A-ox1 و GT-3A-ox2 و pad3، على التوالي، من ثلاث نسخ مستقلة. يتم عرض النتائج كرسوم بيانية صندوقية مع صناديق تعرض النسب المئوية من 25 إلى 75، والخط المركزي يشير إلى الوسيط، والشعيرات تمتد إلى القيم الدنيا والقصوى لا تتجاوز 1.5 نطاق ربعي. , نماذج الجينات لـ CBP6Og و SARD1 مع قمم ATAC-seq وأنماط ربط GT-3A.

المؤلف(المؤلفون) المراسلون: جوزيف ر. إكر

آخر تحديث من قبل المؤلف(المؤلفين): 7 نوفمبر 2024

تسعى Nature Portfolio لتحسين قابلية إعادة إنتاج العمل الذي ننشره. يوفر هذا النموذج هيكلًا للاتساق والشفافية في التقرير. لمزيد من المعلومات حول سياسات Nature Portfolio، انظر سياسات التحرير وقائمة مراجعة سياسة التحرير.

لجميع التحليلات الإحصائية، تأكد من أن العناصر التالية موجودة في أسطورة الشكل، أسطورة الجدول، النص الرئيسي، أو قسم الطرق.
n/a
□

حجم العينة الدقيقة لكل مجموعة/شرط تجريبي، معطاة كرقم منفصل ووحدة قياس
□

بيان حول ما إذا كانت القياسات قد تم أخذها من عينات متميزة أو ما إذا كانت نفس العينة قد تم قياسها عدة مرات
□

اختبار(اختبارات) الإحصائية المستخدمة وما إذا كانت أحادية الجانب أو ثنائية الجانب
يجب وصف الاختبارات الشائعة فقط بالاسم؛ وصف تقنيات أكثر تعقيدًا في قسم الطرق.

□ وصف لجميع المتغيرات التي تم اختبارها
□

وصف لأي افتراضات أو تصحيحات، مثل اختبارات الطبيعية والتعديل للمقارنات المتعددة
□

وصف كامل للمعلمات الإحصائية بما في ذلك الاتجاه المركزي (مثل المتوسطات) أو تقديرات أساسية أخرى (مثل معامل الانحدار) و التباين (مثل الانحراف المعياري) أو تقديرات عدم اليقين المرتبطة (مثل فترات الثقة)
□

لاختبار الفرضية الصفرية، إحصائية الاختبار (مثل ) مع فترات الثقة، أحجام التأثير، درجات الحرية و القيمة المذكورة. أعطِ القيم كقيم دقيقة كلما كان ذلك مناسبًا.

□ لمعلومات تحليل بايزي، معلومات حول اختيار الأولويات وإعدادات سلسلة ماركوف مونت كارلو

□ للتصاميم الهرمية والمعقدة، تحديد المستوى المناسب للاختبارات والتقارير الكاملة للنتائج
□ X
تقديرات أحجام التأثير (مثل حجم كوهين، ، بيرسون )، تشير إلى كيفية حسابها
تحتوي مجموعتنا على الإنترنت حول الإحصائيات لعلماء الأحياء على مقالات حول العديد من النقاط أعلاه.
تم التأكيد
حجم العينة الدقيقة لكل مجموعة/شرط تجريبي، معطاة كرقم منفصل ووحدة قياس
□
بيان حول ما إذا كانت القياسات قد تم أخذها من عينات متميزة أو ما إذا كانت نفس العينة قد تم قياسها عدة مرات

اختبار(اختبارات) الإحصائية المستخدمة وما إذا كانت أحادية الجانب أو ثنائية الجانب
وصف لجميع المتغيرات التي تم اختبارها
L


آخر تحديث من قبل المؤلف(المؤلفين): الاتساق والشفافية

يجب أن تتضمن جميع المخطوطات بيانًا حول توفر البيانات. يجب أن يوفر هذا البيان المعلومات التالية، حيثما ينطبق:

جميع المعلومات التي تدعم الاستنتاجات مقدمة مع الورقة. تم إيداع بيانات تسلسل الخلايا الفردية والكتل التي تم إنشاؤها في هذه الدراسة في قاعدة بيانات التعبير الجيني لمركز المعلومات البيولوجية الوطنية (رقم الوصول GSE226826 و GSE248054). الجينوم المرجعي، التعليق، البيانات المعالجة بالكامل لتحليلات snMultiome متاحة على neomorph.salk.edu/download/Nobori_etal_merfish. بيانات MERFISH متاحة على plantpathogenatlas.salk.edu. الجينات المستهدفة باستخدام MERFIHS مدرجة في الجدول التكميلي 1. معلومات عن البرايمرات المستخدمة في هذه الدراسة مقدمة في الجدول التكميلي 2. بيانات المصدر مقدمة مع هذه الورقة. الجينوم المرجعي، التعليق، البيانات المعالجة بالكامل لتحليلات snMultiome متاحة على neomorph.salk.edu/download/Nobori_etal_merfish.

معلومات السياسة حول الدراسات التي تشمل المشاركين في البحث البشري والجنس والنوع في البحث.

يرجى ملاحظة أنه يجب أيضًا تقديم معلومات كاملة حول الموافقة على بروتوكول الدراسة في المخطوطة.

يرجى اختيار الخيار أدناه الذي يناسب بحثك بشكل أفضل. إذا لم تكن متأكدًا، اقرأ الأقسام المناسبة قبل اتخاذ اختيارك.
علوم الحياة □ العلوم السلوكية والاجتماعية □ العلوم البيئية والتطورية والبيئية
لنسخة مرجعية من الوثيقة بجميع الأقسام، انظرnature.com/documents/nr-reporting-summary-flat.pdf

يجب على جميع الدراسات الإفصاح عن هذه النقاط حتى عندما يكون الإفصاح سلبياً.

محفظة الطبيعة | ملخص التقرير مارس 2021

نحتاج إلى معلومات من المؤلفين حول بعض أنواع المواد والأنظمة التجريبية والأساليب المستخدمة في العديد من الدراسات. هنا، يرجى الإشارة إلى ما إذا كانت كل مادة أو نظام أو طريقة مدرجة ذات صلة بدراستك. إذا لم تكن متأكدًا مما إذا كان عنصر القائمة ينطبق على بحثك، يرجى قراءة القسم المناسب قبل اختيار رد.

محفظة الطبيعة | ملخص التقرير مارس 2021

https://doi.org/10.1038/s41586-024-08383-z
Received: 26 December 2023
Accepted: 8 November 2024
Published online: 8 January 2025
Open access
Check for updates

Interactions between hosts and microorganisms are heterogeneous for multiple reasons. Multicellular host tissues are composed of diverse cell types that have distinct capacities to respond to microorganisms, and microorganisms can occupy niches heterogeneously distributed in the host. Moreover, individual interactions between cells may occur asynchronously. Thus, diverse cell states can co-exist in a tissue. Such heterogeneity can mask fundamental principles of cellular interactions when hosts and microbes are analysed at the tissue scale.

Plant-pathogen interactions have been studied to understand the molecular mechanisms that underlie host immunity and pathogen virulence . Single-cell RNA sequencing (scRNA-seq) of leaf protoplasts has revealed heterogeneous responses of plants infected by virulent bacterial and fungal pathogens . However, our understanding of cell-state diversity is largely limited to a specific cell type (the mesophyll ) infected by immunosuppressive virulent pathogens. Moreover, the spatiotemporal dynamics of plant immune responses are unclear owing to the low-throughput nature of transgenic reporter assays . Furthermore, we have a limited understanding of the gene-regulatory mechanisms that underlie cell-state diversity, such as specific binding of transcription factors (TFs) to cell-state-specific cis-regulatory elements (CREs). Single-nucleus assay for transposase-accessible chromatin followed by sequencing (snATAC-seq) is often used to identify potential CREs in individual cell types and states , but it has not been applied to study plant immune responses. These gaps in knowledge represent substantial roadblocks to understanding how the plant immune system operates as a collective entity of cell populations with distinct functions .

In this study, we aimed to fill these gaps through single-nucleus multiomic (snMultiome) and spatial transcriptomic analyses of a host plant infected by virulent or avirulent bacterial pathogens in a time-course experiment. Specifically, we use single nucleus RNA-seq (snRNA-seq), snATAC-seq and multiplexed error robust fluorescence in situ hybridization (MERFISH) and the following three bacterial pathogens: Pseudomonas syringae pv. tomato DC3000 (hereafter DC3000), DC3000 AvrRpt2 and DC3000 AvrRpm1 (hereafter AvrRpt2 and AvrRpm1, respectively). DC3000 is a virulent pathogen that can suppress plant immunity through effectors and toxins , whereas AvrRpt2 and AvrRpm1 are avirulent pathogens that carry effectors indirectly recognized by plant nucleotide-binding domain and leucine-rich repeat receptors to initiate effector-triggered immunity (ETI) . Our single-cell multidimensional atlas reveals previously uncharacterized cell states, including a rare cell state that emerges at the centre of ETI tissue regions (designated as PRIMER cells). Another cell state surrounding PRIMER cells (designated as bystander cells) is described, and we provide insights into the functions and gene-regulatory mechanisms of these cell states.

We generated a time-resolved snMultiome atlas of A. thaliana leaves infected by DC3000, AvrRpt2 or AvrRpm1 (Fig. 1a). Pathogens were inoculated into leaves at a low dose (optical density at 600 nm ), at which only a subset of plant cells are anticipated to
Plant Biology Laboratory, The Salk Institute for Biological Studies, La Jolla, CA, USA. Genomic Analysis Laboratory, The Salk Institute for Biological Studies, La Jolla, CA, USA. Howard Hughes Medical Institute, The Salk Institute for Biological Studies, La Jolla, CA, USA. Bioinformatics and Systems Biology Program, University of California, San Diego, La Jolla, CA, USA. Laboratory of Plant Pathology, Graduate School of Agriculture, Kyoto University, Kyoto, Japan. Present address: The Sainsbury Laboratory, University of East Anglia, Norwich, UK. Present address: Arc Institute, Palo Alto, CA, USA. e-mail: ecker@salk.edu

pseudobulk expression of subclusters. Well-characterized defence-related genes are shown. See Extended Data Fig. 3c for a comprehensive analysis. The top bars indicate the cell type and major cluster from which each subcluster is derived from. The colour scheme for the major clusters matches with . e, Subclustering of major cluster 6 (PCCs). Defence-related genes showing subcluster-specific expression are shown.f, Schematic of subclustering of clusters 3,7 and 11 (immune-active mesophyll cells).g, Expression of genes involved in different steps of the biosynthesis and secretion of tryptophanderived secondary metabolites shown inh.h. Simplified schematic of the biosynthesis and secretion of tryptophan-derived secondary metabolites.
encounter pathogen cells. Infected leaves were sampled at four different time points and 24 h . As a control, we prepared plants that were mock-infected with water and sampled after 9 h . Two independent replicates were made for the AvrRpt2-infection 9-h condition and the mock condition. We developed a protocol to rapidly isolate nuclei so that transcriptome and epigenome changes during sample preparation were minimized (Methods). A total of 65,061 cells from 15 samples passed quality control for both RNA-seq and ATAC-seq data analyses (Extended Data Fig. 1a-d). Our snATAC-seq data identified more accessible chromatin regions (ACRs) than previously reported bulk ATAC-seq data from immune-activated leaves . Moreover, the datasets showed large overlaps, which suggests that snATAC-seq can capture both known and unknown ACRs (Extended Data Fig. 1e). ATAC-seq peaks associated with a housekeeping gene (ACTIN2) and an immune-related gene (ICS1) were consistently detected across replicates, which provided support for the high reproducibility of the snATAC-seq data (Extended Data Fig.1f). Two independent replicates showed consistent transcriptional reprogramming caused by AvrRpt2 infection (Extended Data Fig.1g,h), which further confirmed the reproducibility of our data.

Dimensionality reduction and clustering were performed using the snRNA-seq data. Some clusters were enriched in specific infection
conditions and at specific time points (Extended Data Fig. 1i), which suggests that the clustering analysis captured distinct cell states induced by pathogen infection. We identified genes specifically expressed in individual clusters (top markers are shown in Extended Data Fig. 2a), which further clarified the identity of each cluster (major cell-type annotations are shown in Fig.1b). Cell types were also predicted on the basis of ATAC-seq data. For instance, a cluster-specific ACR (peak at chromosome 2 position 11172821-11173529) of clusters 0,12,19,21 and 29 was associated with , a marker gene for the epidermis (Extended Data Fig. 2b,c; full names of genes highlighted in this study are provided in the Methods). Overall, the snMultiome data classified cell types and states of A. thaliana leaves infected with a pathogen.

Gene ontology (GO) enrichment analysis of marker genes of each cluster identified mesophyll (clusters 3, 7 and 11) and epidermis (cluster 12) cell populations enriched with defence-related genes, including those involved in the defence hormone salicylic acid (SA) pathway (Fig. 1c and Extended Data Fig. 2d). These clusters were well represented in ETI conditions (AvrRpt2 and AvrRpm1 infection) (Extended Data Fig. 1i), which suggests that these cells were responding to the immune-activating pathogens. This finding was also supported by the strong expression of ICS1, a key gene for pathogen-induced SA
biosynthesis (Extended Data Fig. 2e). We observed increased accessibility to chromatin regions upstream of the ICS1 locus after infection by the ETI strains (Extended Data Fig. 1f), a result consistent with a previous bulk ATAC-seq study . Together, our data captured heterogeneous and coordinated changes in defence gene expression and chromatin accessibility during pathogen infection.

Although the major clusters captured immune-active cells in the mesophyll and the epidermis (Fig. 1c), clusters for other cell types, such as the vasculature, contained both immune-active and non-active cells. This result is probably due to strong developmental signatures (for example, vasculature marker genes are expressed at high levels regardless of immune activation). To capture cell-type-specific immune responses, we performed a second round of clustering for each major cluster, which resulted in 429 subclusters with diverse transcriptome patterns (Fig.1b (right) and Extended Data Fig. 3a). Analyses of selected immune-related genes showed both cell-type-specific gene expression and diverse expression in cell types (that is, cell-state diversity) (Fig. 1d). The subclustering of our snRNA-seq data revealed complex immune responses in leaf tissue that were not captured by bulk RNAseq (Extended Data Fig. 3h). Moreover, in some cases, genes seemed to be highly co-expressed at the bulk transcriptome level but specifically expressed at the subcluster level (Extended Data Fig. 3i, j), a finding that further highlights the value of single-cell analyses.

A more detailed analysis of specific subclusters revealed strong expression of markers for phloem companion cells (PCCs) (Extended Data Fig. 3d) in major cluster 6. This cluster could be separated into 12 subclusters, some of which were enriched with immune-related genes (Fig. 1e). Important genes for systemic acquired resistance (SAR), such as ALD1 and FMO1, were specifically expressed in PCC subcluster 8 (Fig. 1e). This result indicates that a subset of PCCs contribute to sending long-distance signals to systemic leaves. Notably, we did not observe strong expression of ALD1 or FMO1 in other vasculature clusters (Extended Data Fig. 3e). ILL6 was among the marker genes enriched in PCC subcluster 8 (Fig.1e and Extended Data Fig.3f), and this gene is involved in SAR . Together, these findings suggest that this cell population may have a role in SAR, and ILL6 may be a PCC-specific SAR regulator. We identified additional genes specifically enriched in PCC subcluster 8 but not in other major clusters (Extended Data Fig. 3g). As SAR requires a mobile signal that travels from locally infected leaves to systemic leaves through the vasculature, it is possible that the identified PCC population and marker genes have a specific role in SAR.

As another example of cell-population-specific immune responses, we analysed tryptophan-derived defence-related secondary metabolite pathways in immune-active mesophyll cells (clusters 3,7 and 11) (Fig.1f). Distinct expression of genes involved in different steps or pathways of the biosynthesis or secretion of the defence metabolites camalexin and indole glucosinolate was identified (Fig. 1g,h). This result indicates that there is compartmentalization in the activation of defence pathways that can potentially compete for resources (tryptophan). Together, these examples highlight the diversity of plant immune-cell states, thereby confirming the importance of understanding immune signalling pathways and networks at the cell-state level.

Our snMultiome data enabled us to directly compare mRNA expression and ACRs to identify potential CREs in different cell types, infection conditions and time points (Extended Data Fig. 4a). We identified a total of 29,002 significantly correlated (or linked) ACR-gene pairs within 500 kb of each gene across cells in each infection condition. Most links were within a short distance of gene loci (<400 bp; potential promoter regions), whereas others were more distal (potential enhancer regions)
(Extended Data Figs. 4b and 5a). We summarized peak-to-gene linkage data for each gene by using the maximum Pearson’s correlation coefficient values (peak-to-gene linkage score) (Extended Data Fig. 4c). Genes that showed links in all the conditions (cluster 8; Extended Data Fig. 4c) were enriched with cell-type marker genes such as FDH, BCA2 and MAM1 (Extended Data Fig. 4c,d). Genes that showed links specifically in ETI-activated conditions (cluster 4;Extended Data Fig. 4c) were enriched with immunity-related genes (Extended Data Figs. 4d and 5b). CBP60G, a transcriptional regulator of immunity, had multiple ACRs for which accessibility significantly correlated with its mRNA expression (Extended Data Fig. 4e). Genes that showed links specifically in DC3000 infection (cluster 2; Extended Data Fig.4c) were enriched with jasmonic acid (JA)-related genes (Extended Data Fig. 5b). This finding is consistent with the ability of DC3000 to activate the JA pathway in plants using the toxin coronatine and effectors to suppress plant immunity . These results indicate that coordinated and cell-population-specific reprogramming of chromatin accessibility and gene expression is a key feature of leaf development, immunity and exploitation by virulent pathogens.

Although we identified many ACRs that were closely associated with defence genes, ( 304 out of 627 ) of defence genes (marker genes of immune-active clusters and 12 ) did not show significant links with ACRs. Such non-linked defence genes often had constitutively opened chromatin (Extended Data Fig. 4f) despite cluster-specific changes in gene expression (Extended Data Fig. 4f, violin plot). This finding suggests that there is an additional layer of gene regulation, for example, the expression of upstream TFs. Different motifs were enriched in ACRs 2 kb upstream of linked and non-linked defence genes (Extended Data Fig. 4g,h), which suggests that some defence TFs function together with chromatin reprogramming whereas others do not.

To identify cell-type-specific and state-specific gene-regulatory mechanisms, we performed motif enrichment analysis for linked ACRs specific to individual clusters. For instance, comparing an immune-active mesophyll cluster (cluster 3) and a non-immune-active mesophyll cluster (cluster 1) identified the enrichment of motifs for many TFs known to be involved in immunity, including WRKYs and CAMTAs (Fig. 2a). This result highlights the utility of this strategy. We extended this analysis to marker ACRs for all the clusters and identified cluster-specific motif enrichment (Fig. 2b). These findings suggest that different cell types and states use both shared and distinct gene regulation through TFDNA binding. Moreover, this approach can accelerate the identification of TFs with cell-type-specific or state-specific function.

To analyse motif accessibility at the single-cell level–not the cluster level–we calculated the motif enrichment score for each cell using ChromVAR (Fig. 2c), which revealed heterogeneous accessibility to 465 motifs between and within clusters (Extended Data Fig. 6a). We then asked whether the TF motif enrichment score correlates with TF expression, and the results also showed heterogeneous expression (Extended Data Fig. 6b). We identified TFs that showed top correlations between motif enrichment scores and their mRNA expression in different cell types (Fig. 2d and Extended Data Fig. 6c). For instance, WRKY46 mRNA expression and accessibility to WRKY46-binding sites overlapped mainly in immune-active mesophyll and epidermal cells across time points (Fig. 2e). This finding suggests that the WRKY46 regulon has a key role during immune responses in these cell types. By searching for genes that are linked with ACRs containing WRKY46-binding sites (Fig. 2f), we identified potential target genes of WRKY46, including many known defence-related genes (Fig. 2g). This result is consistent with the known function of WRKY46 in the SA pathway and defence against . syringae . Finally, we performed target gene prediction for all the TFs shown in Fig. 2d and found that many TFs target genes related to plant defence (Fig. 2h). TFs that belong to the same family tended to show overlapping target genes (Fig. 2h). In addition to WRKYs, motifs

Fig. 3 | Time-resolved spatial transcriptomics in pathogen-infected leaves. a, Schematic of the MERFISH experiments. for each condition. See Methods for detailed procedures.b, Two-dimensional plots of transcripts for 16 selected genes (Supplementary Table 1) detected using MERFISH in each sample. Plots with transcripts for all 500 genes are provided in Extended Data Fig. 7b. c, Left, a representative field of view (FOV) that shows obscure signalling of DAPI-stained nuclei. Middle, transcript-based segmentation in the same FOV (Methods). Transcripts were coloured on the basis of assigned cells. Right, centroids of cells detected using the transcript-based segmentation method (red dots) and
the result of failed DAPI-staining-based segmentation (yellow region). Similar patterns were observed across FOVs and samples. A systematic quantitative analysis is provided in Extended Data Fig. 7f. d, Violin plots showing the number of transcripts per cell (left) and unique genes per cell (middle) detected in each MERFISH sample and the area size per cell (right). e, UMAP embeddings of cells in each sample based on the expression of 500 genes detected using MERFISH. All MERFISH samples are integrated, and cells are coloured on the basis of de novo Leiden clusters. f, Spatial mapping of Leiden clusters in each sample using the same colour scheme as in e. Scale bar, (c) or 1 mm (b,f).

Extended Data Fig. 7d). After segmentation, transcripts were assigned to cells, which produced a cell-by-gene matrix with each cell having its spatial coordinates. Overall, we detected a median of 161 transcripts and 79 genes per cell from a total of 121,998 cells from 5 leaf sections (1 mock infection and 4 time points after AvrRpt2 infection) (Fig. 3d). We performed integration and de novo clustering of all MERFISH samples, which identified 14 clusters (Fig. 3e and Extended Data Fig. 7g). These clusters were spatially mapped in individual samples, which revealed spatially organized cell populations (Fig. 3e). For instance, a de novo MERFISH cluster captured vasculature cells in tissue sections (Extended Data Fig. 7h). The single-cell and quantitative gene-expression properties of MERFISH provide the opportunity to integrate spatially resolved MERFISH data and molecular-information-rich snMultiome data.

We integrated MERFISH and snMultiome data (five conditions matching the MERFISH samples;Fig.4a) using the shared 500 genes, and cluster labels defined by snMultiome were transferred to the MERFISH cells (Fig. 4a,b and Methods). On the basis of this data integration, we spatially mapped clusters defined in the snMultiome data. Major cell types defined by snMultiome were successfully mapped on the expected regions in a MERFISH tissue sample (Fig. 4c), which indicated successful data integration. This integration of MERFISH and snMultiome data enabled us to explore the spatial distribution of cell populations defined in the snMultiome analysis.

Using this integrated dataset, we spatially imputed the entire transcriptome and chromatin accessibility information (Extended Data Fig. 8a-d). Imputed ICS1 (not included in the MERFISH panel) expression accurately predicted the real spatial expression of ICS1 (based on smFISH) (Extended Data Fig. 8a,b), a result that confirmed the accuracy of data imputation. We also spatially imputed ATAC activity scores (Extended Data Fig. 5f) of ICS1 and ALD1, and these showed consistent patterns with mRNA expression (Extended Data Fig. 8c). The motif activity of HSFB2b was predicted to be high in the immune-active regions of a leaf at 24 h post inoculation (h.p.i.) (Extended Data Fig. 8d), which was consistent with the mRNA expression pattern of validated by MERFISH (Extended Data Fig. 8d). Overall, these results indicate that spatial data imputation of the transcriptome and the epigenome was accurate, which enabled the analysis of gene-regulatory mechanisms at single-cell and spatial resolution. To facilitate exploration of our data, we imputed all 25,299 transcripts detected in the snMultiome analysis and 465 motif enrichment scores on the 5 tissues used in MERFISH experiments and made the data available on a data browser (https:// plantpathogenatlas.salk.edu).

To understand the temporal dynamics of immune responses, we applied pseudotime analysis to our snMultiome data. Pseudotime analysis aligns cells as a trajectory based on their gene expression, which is commonly used for inferring developmental trajectories . We proposed

Fig. 4 | Integration of snMultiome and spatial transcriptome data. a, Integration of snMultiome data (mock infected and AvrRpt2 infected; five samples in total) and the MERFISH data shown in a UMAP. b, Integrated UMAP. Left, nuclei and cells are coloured on the basis of cluster labels, which were transferred from snRNA-seq to MERFISH (Methods). Right, the same integrated UMAP coloured by assay types. c, Left, UMAP of snMultiome data coloured on
the basis of major cell types: epidermis (magenta), mesophyll (green) and vasculature (yellow). Right, spatial mapping of snMultiome cells coloured on the basis of major cell types. AvrRpt2 9 h.p.i. sample was used. Scale bar, 1 mm . d, Pseudotime values calculated for mesophyll cells in the snRNA-seq data. e, UMAP plots showing pseudotime values in cells from each time point.f, Spatial mapping of pseudotime values based on data integration and label transfer.
that by applying pseudotime analysis to a single developmental cell type with various immune states, we could better model the temporal dynamics of heterogeneous infection and immune responses in an infected leaf. Calculation of pseudotime scores for mesophyll cells (Fig. 4d and Methods) showed that the distribution of predicted pseudotime scores in each sample was consistent with what is expected from the real sampling time point (Fig. 4e). That is, cells with low pseudotime scores were enriched at early time points, whereas cells with high pseudotime scores emerged at later time points. This result indicated that the temporal dynamics of immune responses were successfully modelled. Cells with a wide range of pseudotime scores coexisted at 9 and 24 h.p.i., which suggests that cellular immune responses are asynchronous in infected leaves.

To understand the spatial distribution of heterogeneous immune states, we spatially mapped the pseudotime scores using label transfer from snMultiome to MERFISH data (Fig. 4f). The spatial mapping of the pseudotime scores revealed immune-active areas that were distributed in pathogen-infected leaf tissues in a restricted manner (Fig. 4f). These immune-active areas expanded over time and seemed to merge at 24 h.p.i. (Fig. 4f). Notably, immune states also dynamically change over time in each immune-active area, with older immune-active cells (higher pseudotime scores) being surrounded by younger immune-active cells (lower pseudotime scores) (Fig. 4f). These results indicate that the oldest immune-active cells are probably plant cells that had direct contact with pathogen cells at early time points and that immune responses spread to surrounding cells through cell-cell communication over time.

We sought to investigate whether spatially restricted immune-gene expression can be explained by the distribution of bacterial cells. smFISH targeting bacterial metagenes (19 highly expressed genes) detected bacterial colonies at 24 h.p.i., which we overlaid with the spatial map of the pseudotime scores (Extended Data Fig. 8g). As a control,
we performed another MERFISH experiment using an A. thaliana leaf infected by the immunosuppressive pathogen DC3000 at 24 h.p.i. (Extended Data Fig. 8e-g). The distribution of the immune-activating strain AvrRpt2 overlapped with tissue regions with high pseudotime scores (immune heightened) in contrast to the immunosuppressive DC3000 strain (Extended Data Fig. 8e-g). We confirmed this observation by quantitatively analysing the neighbouring plant cells of individual bacterial colonies (Extended Data Fig. 8h,i). Taken together, these results indicate that immune-active regions defined by the pseudotime analysis interact with the ETI-triggering pathogen; we also captured potential immunosuppression by the virulent DC3000 pathogen.

We systematically identified genes for which expression significantly changed across the pseudotime trajectory in the MERFISH data using cell-type-specific inference of differential expression (C-SIDE) (Fig.5a and Methods). BON3, ALD1 and FMO1 were among the genes identified for which expression was higher towards the centre of immune-active regions (Fig. 5a). BON3 was induced in highly restricted areas in the tissue after infection by AvrRpt2 (at 9 h.p.i.) (Fig. 5b), which was distinct from ALD1 and FMO1 (Fig. 5b-d and Extended Data Fig. 9a). Our subclustering analysis of snMultiome data confirmed the pattern observed with the MERFISH results (Fig. 5e), which indicates that we identified two distinct immune cell states. Expression of BON3 was enriched in cells with the highest pseudotime scores (that is, the oldest immune-active cells) (Fig. 4e,f and Extended Data Fig. 9b), which suggests that these cells are early responders to pathogen invasion. Therefore, we designate this cell state as PRIMER cells, from which immune responses might spread. Cells surrounding PRIMER cells are designated as bystander cells.

PRIMER and bystander cells showed distinct transcriptional and epigenetic signatures. PRIMER cells were enriched in CAMTA motifs and a GT-3A motif, whereas bystander cells were enriched in WRKY

Fig.5|Identification and functional characterization of immune-cell populations. a, MERFISH-analysed genes showing significant changes in expression across the pseudotime trajectory (Fig. 4f) in an AvrRpt2 9 h.p.i. sample. Positive fold change values indicate higher expression in cells with high pseudotime values. Top differentially expressed genes (DEGs; ) are labelled in red.b-d, Spatial expression of ALD1, FMO1 and BON3 at 9 h.p.i. (b) and cropped images of two different locations in ( and d). e, Expression of ALD1, FMO1 and BON3 in immune-active subclusters from snRNA-seq data. , mRNA expression (f), MERFISH images (g) and motif activity (h) of GT-3A in immune-active mesophyll cells in the snMultiome data. BON3 expression is also shown ing. i, Enriched motifs linked to genes enriched in PRIMER or bystander cells.j, Venn diagram showing overlapping of genes highly expressed in PRIMER or bystander cells with genes previously shown to be suppressed by CAMTA3 (ref. 39). k, Bulk RNA-seq of wild-type (WT) plants and GT-3A-ox plants
treated with water (M) or DC3000 (DC) at 24 h.p.i. DEGs (adjusted , -adjusted fold change ) between WT and GT-3A-ox after infection are shown. GO terms enriched in specific gene sets are shown with adjusted values. Three replicates per condition. I,m, Growth of DC3000 (I) or AvrRpt2 (m) in Col-0 (WT) plants compared with GT-3A-ox plants (I) and gt3a-KO plants (m) at indicated times. and leaves from 3 and 4 independent experiments, respectively. RLU, relative light unit. k-m, Pathogens were syringe-infiltrated at a suspension dose of . Adjusted values were calculated using two-tailed Student’s -test with Benjamini-Hochberg correction. I,m, Results are shown as box plots, in which boxes represent the 25th-75th percentiles, the centre line indicates the median, and whiskers extend to minimum and maximum values within 1.5 times the interquartile range. n, Proposed model for the potential role and regulation of immune-cell states. Scale bar, (c,d), (g) or 1 mm (b).
motifs (Fig. 5i and Extended Data Fig. 9i). Genes previously shown to be repressed by CAMTA3 were significantly overrepresented in bystander cells compared with PRIMER cells (Fig. 5j; false discovery rate (FDR) , with hypergeometric test correction using the Benjamini-Hochberg method), a result that supports the transcriptional-repressive role of CAMTA3 in PRIMER cells.

To understand the function of PRIMER cells and their gene-regulatory mechanisms, we investigated an snMultiome subcluster in which was enriched (subcluster 4). Although many PRIMER cell marker genes were previously uncharacterized, we found several known genes, including WRKY8 and LSD1, as PRIMER cell markers (Extended Data Fig. 9c). Notably, BON3, WRKY8 and LSD1 are all negative regulators of immunity , although the precise mechanisms of how they suppress immune responses remain elusive.

Further marker gene analysis of the PRIMER cell cluster identified GT-3A (which encodes a trihelix DNA-binding TF) , the function of which in leaf immunity remains uncharacterized (Fig. 5f). GT-3A showed a spatial expression pattern similar to that of BON3 (MERFISH; Fig. 5g), and the GT-3A motif was highly accessible in PRIMER cells (snATAC-seq; Fig. 5h), which implies that it has a gene-regulatory function in this cell population. To understand the role of GT-3A, we generated transgenic A. thaliana plants that ectopically overexpress this gene (GT-3A-ox). Bulk RNA-seq analyses revealed impairment in the induction of genes involved in the SA pathway of GT-3A-ox plants infected by DC3000 (Fig. 5k). Furthermore, two independent lines of

GT-3A-ox were more susceptible to DC3000 infection (Fig. 5l), which suggests that this previously unidentified cell-state-specific TF can negatively regulate immunity. We also tested plants in which GT-3A was knocked out (gt3a-KO), and this mutant was more susceptible to AvrRpt2 infection (Fig. 5m). This result indicates that GT-3A function in PRIMER cells is required for optimal defence against avirulent pathogens. snRNA-seq analysis of -KO plants revealed genes that are potentially regulated by GT-3A, either directly or indirectly. Among such genes, PUB36 expression in PRIMER cells was significantly impaired in -KO plants at both 9 and 24 h.p.i. Notably, has a GT-3A-binding motif in the upstream region (Extended Data Fig. 9d). SA-related genes, including ALD1, were expressed at lower levels in bystander cells in gt3a-KO plants than in wild-type plants (Extended Data Fig. 9e,f). This finding suggests that induction of GT-3A in PRIMER cells is important for the proper induction of defence genes in surrounding cells. In summary, through the integration of time-resolved snMultiome and spatial transcriptome data, we identified previously unknown immune-cell states and a cell-population-specific TF that regulates plant immunity.

Plant immunity comprises a multicellular network in which individual cells interpret input signals with their distinct molecular networks and communicate with other cells. Our molecularly defined spatiotemporal
atlas of pathogen-infected leaves revealed various cell states in transcriptome and epigenome detail. This resource also provides a means to investigate individual cell states that have been obscured in conventional bulk or dissected tissue analyses and by live imaging of a limited number of reporter lines. For instance, we identified a PCC subpopulation with a distinct state characterized by the induction of SAR genes (Fig.1e) and mesophyll subpopulations that activate different branches of tryptophan-derived defence metabolite pathways (Fig.1g).

We also described a rare cell state located at the nexus of immuneactive hotspots, which we termed the PRIMER cell state (Fig. 5). In addition to mapping previously characterized genes (BON3, WRKY8, LSD1 and CAMTA3) with common immunosuppressive functions to PRIMER cells, our integrative snMultiome and MERFISH analyses identified another PRIMER cell marker gene, GT-3A, which encodes a TF, and demonstrated that it contributes to plant immunity against pathogen infection (Fig. 5). Comparisons between PRIMER cells and their surrounding cells (bystander cells) revealed distinct transcriptional and epigenetic landscapes (Fig. 5). It is possible that there are additional specialized cell states in PRIMER and bystander cells. A deeper understanding of immune-cell states requires the development and application of new methodologies, such as imaging techniques that visualize both pathogen effector proteins and numerous plant genes simultaneously in three dimensions at single-cell resolution (see the section ‘Limitation of this study’ in the Supplementary Information).

Previous studies have shown that GT-3A negatively regulates plant defence against nematode infection in the root , which suggests that GT-3A may have a role in defence against different types of pathogens that infect other tissues. In support of this hypothesis, GT-3A-ox plants were more susceptible to the fungal pathogen Colletotrichum higginsianum (Extended Data Fig. 9j).
Recent discoveries have shed light on the role of nucleotide-binding domain and leucine-rich repeat receptors as calcium channels or NADases . However, the exact mechanisms that underlie ETI activation and how this process effectively suppresses pathogen growth are not fully understood. It has been proposed that localized acquired resistance (LAR) may be important for . LAR is a strong defence response in cells surrounding those that have been exposed to pathogen effectors. We revealed the spread of ETI responses from PRIMER cells (Fig.5e) and captured potential LAR responses with detailed molecular information. Notably, ALD1 and FMO1, encoding canonical SAR components, were expressed in bystander cells but not in PRIMER cells (Fig. 5b-e), which indicates that the SAR pathway has a specific role in LAR responses. This idea is plausible given the long-distance signalling capability of SAR. We propose that PRIMER cells may undergo hypersensitive cell death, which subsequently sends signals to neighbouring cells that activate immune responses, including the SAR pathway. As activation of the SA pathway can suppress cell death , GT-3A may contribute to cell death in PRIMER cells by suppressing SA signalling (Fig.5n). Further analysis of our data could uncover the roles of various immune cell states and how these cells communicate with surrounding cells to confer successful defence.

In addition to identifying previously unknown immune-cell populations, our snMultiome data predicted numerous putative TF-ACR-gene modules. The successful integration of snMultiome and MEFISH data enabled us to spatially map gene expression and chromatin states (Extended Data Fig. 8a-d). These results can be used to discover previously uncharacterized immunity-related genes, and defence-related CREs can also be identified through de novo motif analysis.

Finally, we built a database (https://plantpathogenatlas.salk.edu) to facilitate the exploration of previously uncharacterized cell populations associated with disease and resistance with spatial and temporal information and potential regulatory mechanisms. Our database can be used for hypothesis generation and testing and will catalyse new discoveries of molecular mechanisms that underlie plant-microorganism interactions at high resolution.

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-024-08383-z.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
(c) The Author(s) 2025

The following reagents and kits were used: ammonium persulfate (Sigma, 09913); BSA solution (Sigma, A1595); Chromium Next GEM Single Cell Multiome ATAC + Gene Expression kits (10x Genomics, PN-1000283); Corning Falcon cell strainers (Corning, 08-771-2); dithiothreitol (Thermo, R0861); EDTA, pH 8.0 RNase-free (Invitrogen, AM9260G); KCl (2 M) RNase-free (Invitrogen, AM9640G); MACS SmartStrainers (Milteny Biotec, 130-098-458); MERSCOPE 500 Gene Imaging kits (Vizgen, 10400006); MERSCOPE 500 Gene Panel (Vizgen, 10400003); MERSCOPE Sample Prep kits (Vizgen, 10400012); -tetramethylethylenediamine (Sigma, T7024-25ML); NaCl (5 M) RNase-free (Invitrogen, AM9759); NP40 (IGEPAL CA-630) (Sigma, I8896); paraformaldehyde (Sigma, F8775); PBS ( ) pH 7.4 RNase-free (Thermo Fisher, AM9625); protease inhibitor cocktail (Sigma, P9599); Protector RNase inhibitor (Sigma, 3335402001); Scigen Tissue-Plus OCT compound (Fisher, 23-730-571); spermidine (Sigma, S2626); spermine (Sigma, 85590); Triton-X (Sigma, 93443); and UltraPure 1 M Tris-HCI buffer pH 7.5 (Invitrogen, 15567027).

The following genes were highlighted in this study: ALD1 (AGD2-LIKE DEFENCE RESPONSE PROTEIN 1, At2g13810); BCA2 (BETA CARBONIC ANHYDRASE 2; At5g1474O); BON3 (BONZAI 3; At1gO886O); CBP6Og (CALMODULIN BINDING PROTEIN 6O-LIKE G; At5g2692O); FDH (FIDDLEHEAD;At2g26250); FMO1 (FLAVIN-DEPENDENT MONOOXYGENASE 1; At1g19250); GT-3A (TRIHELIXDNA-BINDING FACTOR (GTFACTOR) BINDING 3A; At5gO138O); ICS1 (ISOCHORISMATE SYNTHASE 1; At1g7471O); LSD1 (LESION SIMULATING DISEASE 1; At4g2O38O); HSFB2b (HEAT SHOCK TRANSCRIPTION FACTOR B2b; At4g1166O); ILL6 (IAA-LEUCINE RESISTANT (ILR)-LIKE GENE 6; At1g4435O); MAM1 (METHYLTHIOALKYLMALATESYNTHASE 1; At5g23010); WRKY8 (At5g46350); and WRKY46 (At2g46400).

A. thaliana Col-0 was grown in a chamber at with a light period and relative humidity for days. Bacterial strains were cultured in King’s B liquid medium with antibiotics (rifampicin and tetracycline) at . Three bacterial strains-P. syringae pv. tomato DC3000 with an empty vector (pLAFR3), avrRpt2 (pLAFR3) and avrRpm1 (pLAFR3)-have been previously described . Bacteria were collected by centrifugation and resuspended in sterile water to an of 0.001 (approximately c.f.u. ). In total, 20 A. thaliana leaves (4 fully expanded leaves per plant) were syringe-inoculated with bacterial suspensions using a needleless syringe. Syringe infiltration was performed on one or two corners of each leaf, and bacterial suspensions were spread throughout the entire leaf. We chose four different time points ( and 24 h ), representing early stages of infection and when dynamic transcriptional reprogramming was observed in a previous study that used bulk RNA-seq . For each strain, four time points were sampled at the same time of day by infiltrating bacteria at different times to minimize the influence of circadian rhythms. The 20 infected leaves were collected using forceps and immediately processed for extraction of nuclei. For the mock condition, water-infiltrated leaves were collected after 9 h .

pWAT206 was a gift from A. Takeda. The following plasmids were constructed using HiFi DNA assembly kits (New England Biolabs) and were verified by Sanger sequencing. Three PCR fragments were amplified from pBICRMsG using primer pairs (mEGFP-Nter_1_F plus mEGFP_1_R; mEGFP_2_F plus mEGFP_2_R; and mEGFP_3_F plus mEGFPNter_3_R) and were then assembled into Stul/Ascl-digested pBICAscII . The resulting
plasmid was used as a template for two PCR reactions using primer pairs mEGFP_Nter_1_F plus mEGFP_4R and mEGFP_4F plus mEGFP_ Nter_3R. These PCR products were assembled into Stul/Ascl-digested pBICAscII, which produced pBIC_mEGFP_Nter. A PCR fragment was amplified from pBIC_mEGFP_Nter using mEGFP_Cter_F1 plus mEGFP_ Cter_R, and then used as a template for PCR using mEGFP_Cter_F2 plus mEGFP_Cter_R. The resulting PCR fragment was assembled into Stul/ AscI-digested pBIC_mEGFP_Nter to obtain pBIC_mEGFP_Cter. The coding sequence of GT-3A (At5g1380) was amplified by PCR using the primer pair GT3a_mEGFP_F plus GT3a_mEGFP_R and was assembled into Stul-digested pBIC_mEGFP_Cter to obtain pBIC_GT3a_mEGFP. The sequence encoding carboxy-terminally mEGFP-fused GT-3A was amplified from pBIC_GT3a_mEGFP by PCR using a primer pair ( and 35 Ster_HiFi_R) and was then assembled into Xhol/Xbal-digested pWAT206 to obtain pWAT206_GT3a_mEGFP. A. thaliana Col-0 plants were transformed using Agrobacterium tumefaciens GV3101 (pMP90, pSoup) with pWAT206_GT3a_mEGFP as previously described . Two independent Basta-resistant transgenic lines were used for bacterial growth assays with a bioluminescent P. syringae pv. tomato DC300043 strain and for lesion development analysis using C. higginsianum MAFF 305635. Plants were grown in a climatic chamber with a temperature of relative humidity and light intensity of 6,000 lux (about ) for 10 h . Bacterial suspensions at sterile water were syringe-infiltrated into leaves of -week-old plants. Bacterial growth was measured as bioluminescence using a GloMax Navigator Microplate luminometer (Promega) as previously described . For lesion development analysis, leaves of 4-5-week-old plants were drop-inoculated with a conidial suspension of C.higginsianum ( conidia per ml). The inoculated plants were kept under high humidity in a climatic chamber with a temperature of and light intensity of 6,000 lux for 10 h . Lesion size was measured at 6 days after inoculation using ImageJ.

The gt3a-KO plants were created by genome editing with CRISPR-Cpf1. The plasmid pWAT235-Cas9-HF-cc was provided by A. Takeda. The Cpf1 sequence was amplified from pY004 (Addgene, 69976) by PCR using a primer pair (Spel-NLS-Cpf1-F and BamHI-NLS-Cpf1-R; Supplementary Table 2), followed by digestion with Spel and BamHI. The digested DNA fragment was ligated to pWAT235-Cas9-HF-cc, which had been digested with Spel and BamHI, resulting in pWAT235-Cpf1-cc. Two PCR fragments were amplified from pWAT235-Cpf1-cc using primer pairs Gb-ccdB-F and Gb-ccdB-R, and Gb-U626-F and Gb-U626-R. These PCR fragments and Bsal-digested pWAT235-Cas9-HF-cc were assembled into a single plasmid using a HiFi DNA assembly kit, which resulted in pWAT235-Cpf1-U626-cc-polyT. Two oligonucleotides, GT3a_gRNA2_F and GT3a_gRNA2_R, were hybridized and ligated to Bsal-digested pWAT235-Cpf1-U626-cc-polyT, which resulted in pWA T235-Cpf1-U626-GT3a-gRNA2-polyT. A. thaliana Col-O plants were transformed using A. tumefaciens GV3101 (pMP90, pSoup) with pWAT235-Cpf1-U626-cc-polyT as previously described .A transgenic line with a premature stop codon due to a 5 -bp deletion in the GT-3A gene was selected and used throughout this study.

Leaves of wild-type plants and one of the transgenic lines expressing GT-3A-mEGFP were syringe-infiltrated with water or DC3000 at an of 0.001 and were collected 24 h after infiltration. Total RNA was extracted using TRI Reagent (Sigma). One microgram of total RNA was used for library preparation using the BrAD-seq method to create strand-specific digital gene-expression libraries . The libraries were sequenced on a DNBSEQ-G400 platform at BGI, which produced 100 -bp end reads. Because the quality of reverse reads was poor due to the poly sequence, only forward reads were used for analysis. Trimming of the first 8 bases and adaptors and quality filtering were
performed using fastp (v.0.19.7) with the parameters-x-f8-q30-b50. The trimmed and quality-filtered reads were mapped to the Arabidopsis genome (TAIR10) using STAR (v.2.6.1b) with default parameters and transformed to a count per gene per library using featureCounts (v.1.6.0) . Statistical analysis of the RNA-seq data was performed in the R environment (v.4.1.3). Because the BrAD-seq method involves poly(A) enrichment, mitochondrial and chloroplast genes were excluded. Genes with mean read counts of fewer than ten per library were excluded from the analysis. The resulting count data were subjected to TMM normalization using the function calcNormFactors in the package edgeR, followed by log transformation by the function voomWithQualityWeights in the package limma. To each gene, a linear model was fit using the function ImFit in the limma package. For variance shrinkage in the calculation of values, the eBayes function in the limma package was used. The resulting values were then corrected for multiple hypothesis testing by calculating Storey’s values using the function qvalue in the package qvalue. To extract genes with significant changes in expression, cut-off values of and -transformed fold change| > 1 were applied.

Fresh nucleus purification buffer (NPB; 15 mM Tris pH 7.5, 2 mM EDTA, spermidine, 0.2 mM spermine, 1:100 BSA and 1:100 protease inhibitor cocktail) was prepared before the experiment and chilled on ice. All the subsequent procedures were performed on ice or at . Twenty leaves were chopped in cold NPB with 1:500 Protector RNase IN with a razor blade on ice for 5 min to release nuclei and then incubated in 20 ml NPB. The crude extract of nuclei was sequentially filtered through and cell strainers ( , Corning Falcon cell strainers, Corning, 08-771-2; , MACS SmartStrainers, Milteni Biotec, 130-098-458). Triton-X and NP40 were added to the extract to a final concentration of each, and the extract was incubated at for 5 min with rotation. The suspension was centrifuged at 50 g for 3 min in a swing-rotor centrifuge to pellet non-nucleus debris and the supernatant was recovered. The nuclei were pelleted by centrifugation at 500 g for 5 min in a swing-rotor centrifuge. When the pellet was green, the pellet was resuspended in 20 ml NPBd (NPB with 0.1% Triton-X and 0.1% NP40) with 1:1,000 Protector RNase IN by pipetting, followed by centrifugation at 500 g for 5 min . When the pellet was translucent, the NPBd wash was skipped. The pellet was then washed by resuspending it in 20 mlNPB with 1:1,000 Protector RNase IN and centrifuging at 500 g for 5 min in a swing-rotor centrifuge. The pellet was resuspended in of Nuclei Buffer (10x Genomics, PN-2000207) with 1:40 Protector RNase IN and 1 mM dithiothreitol. The suspension of nuclei was centrifuged at 50 g for 3 min in a swing-rotor centrifuge to pellet non-nucleus debris and the supernatant was recovered. This step was repeated one more time. The resulting nuclei were manually counted using a haemocytometer. Nuclei were pelleted by centrifugation at 500 g for 5 min in a swing-rotor centrifuge and the supernatant was removed, leaving approximately of the buffer. Nuclei were counted again, and up to 16,000 nuclei were used for subsequent steps. However, in most samples, we did not load the maximum number of nuclei that the Genomics kit accepts to avoid the risk of clogging the instrument, which can result in variable numbers of recovered nuclei among samples (Extended Data Fig. 1b). scRNA-seq and ATAC-seq libraries were constructed according to the manufacturer’s instructions (10x Genomics, CG000338). scRNA-seq libraries were sequenced using an Illumina NovaSeq 6000 in dual-index mode with ten cycles for i7 and i5 indices. snATAC-seq libraries were also sequenced using an Illumina NovaSeq 6000 in dual-index mode with 8 and 24 cycles for the i 7 index and the i 5 index, respectively.

Raw data processing. Sequence data were processed to obtain single-cellfeature counts by running cellranger(v.6.0.1) and cellranger-arc
(v.2.0.0) for snRNA-seq data and snATAC-seq data, respectively. For snRNA-seq, the -include-introns option was used to align reads to the A. thaliana nuclear transcriptome built using the TAIR10 genome and the Araport 11 transcriptome. The chloroplast genome was removed from the reference genome for the analysis of both snRNA-seq and snATAC-seq data. The . thaliana TAIR10 genome was downloaded from https://plants.ensembl.org/, and the chloroplast genome was manually removed. The . thaliana TAIR10 gene annotation file was manually modified by removing chloroplast genes and replacing semicolons in the ‘gene_name’ column with hyphens to prevent errors during cellranger-arc processing. We note that a mean of and of reads were mapped on the chloroplast genome in snRNA-seq data and snATAC-seq data, respectively. Removing the chloroplast genome from the reference did not affect the overall RNA and ATAC count distribution. Count data were analysed using the R packages Seurat (v.5) and Signac .

Quality control and cell filtering. Before integrating the datasets, doublets were predicted using DoubletFinder and filtered out. Quality control matrices for snATAC-seq were generated using a modified version of the loadBEDandGenomeData function in the R package Socrates . ACRs were identified using MACS2 (ref. 53) with the following parameters: -g (genomesize) , shift , extsize , and –qvalue , –nomodel, –keep-dup all. The fraction of reads mapping to within 2 kb upstream or 1 kb downstream of the transcription start site (TSS) was calculated. Nuclei were filtered using the following criteria: 200 < RNA UMI count < 7,000; RNA gene count > 180; 200 < ATAC UMI count < 20,000; and fraction of RNA reads mapped to mitochondrial genome . Seurat objects of individual samples were merged using the Merge_Seurat_List function of the scCustomize package.
snRNA-seq clustering. snRNA-seq clustering was performed using the R package Seurat. The cell-by-gene RNA count matrix was normalized using SCTransform. Dimension reduction was performed using principal component analysis with RunPCA. Technical variance among samples was reduced using Harmony using principal components (PCs) 1-20.Graph-based clustering was performed on the Harmony-corrected PCs 1-20 by first computing a shared nearest-neighbour graph using the PC low-dimensional space (with neighbours). Louvain clustering (resolution ) was then applied, and the clusters were projected into an additionally reduced space using UMAP (n.neighbours and min.dist ).
snATAC-seq peak calling. Peaks were called independently on each cluster defined by snRNA-seq data and then combined using the CallPeaks function of Signac, which uses MACS2 with the following parameters: effective.genome.size , extsize , shift . Peak counts were quantified using the FeatureMatrix function. Compared with the default peak calling pipeline of cellranger-arc ( 24,394 peaks), this cluster-specific peak calling approach was able to capture more peaks ( 35,560 peaks), which is consistent with a previous report .
snATAC-seq clustering. Dimensionality reduction was performed using latent semantic indexing (LSI) . First, the top most common features were selected using the FindTopFeatures function. Then, the term-frequency inverse-document-frequency (TF-IDF) was computed using RunTFIDF with scale.factor . The resulting TF-IDF matrix was decomposed with singular value decomposition with RunSVD, which uses the irlba R package. Technical variance among samples was reduced using Harmony with LSI components 2-10. Graph-based clustering was performed on the Harmony-corrected LSI components 2-20 by first computing a shared nearest-neighbour graph using the LSI low-dimensional space (with neighbours). Louvain clustering (resolution ) was then applied, and the clusters were
projected into an additionally reduced space with UMAP(n.neighbours= 30 L and min.dist ).

RNA-ATAC joint clustering. The two modalities were integrated by weighted nearest-neighbour analysis using FindMultiModalNeighbors of Seurat with Harmony-corrected PCs 1-20 for RNA and Harmony-corrected LSI components 2-20 for ATAC. Then, SLM (resolution ) was applied and projected with UMAP (n.neighbours and min.dist=0.1).

ATAC gene activity score. ATAC gene activity score was calculated using the GeneActivity function of Signac with extend.upstream=400.

Peak-to-gene linkage analysis. LinkPeaks of Signac was used to call significant peak-to-gene linkage for each infection condition (mock, DC3000, AvrRpt2 and AvrRpm1). Background-corrected Pearson’s correlation coefficients between the gene expression of each gene and the accessibility of each peak within 500 kb of the gene TSS were calculated. A value was calculated for each peak-gene pair using a one-sided -test, and peak-gene pairs with and a Pearson’s correlation coefficient of were retained as significant links.

Motif enrichment analysis. Motifs present in the JASPAR2020 database for Arabidopsis (species code 3702) were used. Cluster-specific peaks were first identified using FindMarkers of Seurat with default parameters. A hypergeometric test was used to test for the over-representation of motifs in the set of differentially accessible peaks using FindMotifs of Signac. Motif plots were generated using MotifPlot. Motif enrichment scores (motif deviation scores) of individual TF motifs in individual cells were calculated using chromVAR . For the integration of motif enrichment scores and mRNA expression, motif names inJASPAR2020 were matched with gene names. Motifs that could not be uniquely associated with genes were removed from analysis.

TF target prediction. To predict genes that are regulated by a TF, ACRs containing each TF motif were extracted. Then, genes for which expression correlated with these ACRs were identified. We considered genes within 5 kb from an ACR with a linkage score of as significant candidates. GO enrichment analysis was performed for these candidates for each TF using the enrichGO function of clusterProfiler with org.At.tair. db annotation. For the analysis in Fig. 2h, a more stringent threshold (linkage score > 0.2) was applied. GO enrichment plots were created using ggplot2 (ref. 57).

Subcluster analysis. Nuclei with the same subcluster label were aggregated, and -transformed transcripts per million (TPM) values were calculated. Subclusters with more than 18,000 undetected genes were removed from the analysis. For the subclusters that passed the filtering step, genes were clustered using -mean clustering with (determined using the elbow method) (Extended Data Fig. 3b). Then, GO enrichment analysis was performed for the genes in each cluster (Extended Data Fig. 3b). Three clusters (1, 5 and 8) showed enrichment of an immunity-related function (responses to SA); genes in these clusters were defined as ‘putative immune genes’ and used for downstream analysis.

Pseudotime analysis. To calculate pseudotime, the cell-by-gene matrix for cells annotated as mesophyll was obtained from the snMultiome data. We used the function scanpy.tl.dpt with n_dcs=2. Pseudotime trajectories were constructed with each mesophyll cell in mock samples being used as a starting cell. Then, these individual trajectories were averaged in each cell to create a single, unified pseudotime trajectory. Heat map gene trends for a given gene were calculated by fitting a linear GAM using the pygam function LinearGAM withs (0, lam=400) to fit a GAM to RNA expression levels across sorted pseudotime values.

Among the subclusters of immune-active mesophyll cells in the snMultiome data, PRIMER and bystander cell clusters were defined on the basis of expression of the marker genes BON3 and ALD1, respectively. By comparing these cell states, DEGs were identified using the FindMarkers function of Seurat. Motif enrichment analysis was performed on ACRs within 2 kb of the cell-state-enriched DEGs.

To assess the contribution of CAMTA3 in different cell states, genes regulated by CAMTA3 were first identified using a published bulk RNA-seq dataset , comparing wild-type and camta3-D (a dominant negative mutant of CAMTA3). Genes suppressed by CAMTA3 after ETI (triggered by AvrRpm1 or AvrRps4 infection) were overlapped with the cell-state DEGs defined above. Hypergeometric tests were performed to assess the significance of the overlaps.

snRNA-seq was performed on gt3a-KO plants infected by AvrRpt2 or treated with water (mock) at 9 and 24 h.p.i. Two independent replicates were prepared for each infection condition. The data were filtered using the same criteria as for the snMultiome analysis without considering ATAC-seq information. The -KO snRNA-seq data were integrated with the Col-0 snMultiome data (mock and AvrRpt2 conditions), and de novo clustering was performed in the same way as described above. On the basis of immune-gene expression, immune-active mesophyll clusters were identified and further subclustered. From these subclusters, PRIMER and bystander cells were defined on the basis of expression of the marker genes BON3 and ALD1, respectively. For each cell state, DEG analysis was performed comparing Col-O and gt3a-KO plants to assess the effect of the GT-3A mutation.

Our snATAC-seq data were compared with published bulk ATAC-seq data of mature . thaliana leaves activating pattern-triggered immunity (PTI), ETI or both PTI and ETI, as well as non-immune-active leaves . All ACRs identified in the bulk ATAC-seq datasets were combined and compared with the ACRs identified in our snATAC-seq datasets.

MERFISH panel design. We curated 500 target genes that included the following genes: (1) previously defined markers of A. thaliana leaf cell types ; (2) genes involved in various processes such as immunity, hormone pathways and epigenetic regulation; and (3) a variety of TFs previously analysed using DAP-seq (a TF-DNA interaction assay) . Genes for which more than 25 specific probes could not be designed based on probe design software from Vizgen were excluded from the target gene panel. Highly expressed genes could cause the overcrowding of smFISH signals and hinder MERFISH quantification. To avoid including highly expressed genes in the panel, we assessed target gene expression by using a publicly available bulk RNA-seq dataset of A. thaliana infected by AvrRpt2 in the same setup as the current study at eight different time points . For each gene, the highest expression value among the eight time points was used. Genes that showed TPM values of were not included in the panel. The total TPM of the 500 genes was approximately 22,000 . ICS1 was targeted with a single round of smFISH as this gene is highly expressed. The smFISH result was provided as an image without quantitative information. Bacterial cells were visualized by targeting 19 highly expressed genes (based on previously published in planta bulk RNA-seq data ) as a single target. Supplementary Table 1 has a list of the genes targeted by MERFISH. All the probes were designed and constructed by Vizgen.

Tissue sectioning, fixation and mounting. Plants were grown according to the methods described above. Leaves matching the aforementioned treatments and time points were excised and immediately
incubated and acclimated in OCT (Fisher) for 5 min. Following incubation, the leaves were immediately frozen as previously described . Tissue blocks were acclimated to in a pre-cooled cryostat chamber (Leica) for 1 h . Tissue blocks were trimmed until the tissue was reached, after which sections were visually inspected until the region of interest was exposed. Sample mounting and preparation were performed according to the MERSCOPE user guide, but with slight modifications. In brief, a section was melted and mounted onto a room-temperature MERSCOPE slide (Vizgen, 20400001), placed into a Petri dish and re-frozen by incubation in the cryostat chamber for 5 min . Subsequent steps were performed with the mounted samples in the Petri dish. The samples were then baked at for 5 min and were incubated in fixation buffer ( PBS and formaldehyde) for 15 min at room temperature. Samples were then washed with PBS containing 1:500 RNase inhibitor (Protector RNase inhibitor, Millipore Sigma) for 5 min at room temperature in triplicate. Following the final PBS wash, samples were dehydrated by incubation in 70% ethanol at overnight.

MERFISH experiment. Tissue sections were processed following Vizgen’s protocol. After removing 70% ethanol, the sample was incubated in the sample prep wash buffer (PN20300001) for 1 min and then incubated in the formamide wash buffer (PN20300002) at for 30 min . After removing the formamide wash buffer, the sample was incubated in the MERSCOPE Gene Panel mix at for 42 h . After probe hybridization, the sample was washed twice with the formamide wash buffer at for 30 min and once with the sample prep wash buffer at room temperature for 2 min . After the washing step, the sample was embedded in hydrogel by incubation in the gel embedding solution (gel embedding premix (PN20300004), ammonium persulfate solution and -tetramethylethylenediamine) at room temperature for 1.5 h . Then, the sample was cleared by first incubating it in digestion mix (digestion premix (PN20300005) and 1:40 protector RNase inhibitor) at room temperature for 2 h , followed by incubation in the clearing solution (clearing premix (PN20300003) and proteinase K) at for 24 h and then at for 24 h . The cleared sample was washed twice with the sample prep wash buffer and stained with DAPI and PolyT staining reagent at room temperature for 15 min . Samples were then washed with the formamide wash buffer at room temperature for 10 min and rinsed with the sample prep wash buffer. The sample was imaged using a MERSCOPE instrument, and detected transcripts were decoded on the MERSCOPE instrument using a codebook generated by Vizgen. Transcripts were visualized using Vizgen’s Visualizer.

MERFISH segmentation and processing. Cell-boundary segmentation was performed for each MERSCOPE data output. DAPI-targeting and poly(A)-targeting probes demonstrated variable success in staining nuclei and cytoplasm, respectively, depending on the samples and tissue regions examined (Fig. 3b and Extended Data Fig. 7e show failed and successful segmentation, respectively). By contrast, dense transcript areas marked nucleus locations more robustly (Extended Data Fig. 7e). Therefore, a transcript-based segmentation method was used. For each sample, a two-dimensional Numpy array of zeroes was generated, which modelled the total pixel area imaged. The coordinates of identified RNA transcripts were changed from 0 to 1 in this array. Next, the array was blurred using the cv2.GaussianBlur function in OpenCV with ksize=(5,5). The resulting array was chunked into pixel regions. These regions were loaded into Cellpose , a deep-learning-based segmentation tool, and a custom segmentation model was trained by manually segmenting nucleus objects across ten pixel regions in the Cellpose GUI.

The custom model was then used to predict the segmentation boundaries in all the remaining regions, with the parameters diameter=22.92, flow_threshold and cell probability threshold . Next, the total number of cells in an experiment was calculated by summing the
number of unique cells across all regions. This number was then used to initialize the –num-cells-init command in another segmentation tool, Baysor , which considers the joint likelihood of transcriptional composition and cell morphology to predict cell boundaries. Baysor was run using a downloaded Docker image and parameters-s n clusters 1, -i1, –force-2d, min-molecules-per-gene=1, min-molecules-per-cell=50, scale , scale-std=” “, estimate-scale-from-centers=true, min-molecules-per-segment , new-component-weight , new-component-fraction=0.3.
To test the quality of our transcript-based segmentation method, we used an FOV with successful DAPI staining (which was rare in our samples) and performed DAPI-based watershed segmentation and transcript-based segmentation (Extended Data Fig. 7d). Results from these two segmentation strategies agreed with each other in general, with the transcript-based approach capturing transcripts in the cytoplasm in addition to those in the nucleus (Extended Data Fig. 7d-f). This result indicated that our segmentation approach can reliably capture cells.
After Baysor segmentation, a cell-by-gene matrix was created from the transcript cell assignments. Cells with fewer than 50 assigned transcripts were removed. Scanpy was used for post-processing of our MERFISH experiments. After loading the respective cell-by-gene matrix into an Anndata object for each experiment, we stored the spatial coordinates of each cell obtained from Baysor. The individual transcript counts in each cell were normalized by the total number of transcript counts per cell. The Anndata cell-by-gene matrix was then log-scaled.

MERFISH-snMultiome data integration and label transfer. To integrate the MERFISH experiments with each other, we used FindIntegrationAnchors followed by the IntegrateData function in Seurat (v.5). To integrate the MERFISH experiment data with our snMultiome data, we integrated the corresponding time points in each modality separately. For each time point in the infection data (mock, 4, 6, 9 and 24 h), we used gimVI from scvi-tools to project both the snMultiome and MERFISH cells into the same latent space using their RNA expression counts. gimVI was trained over 200 epochs with a size 10 latent space for each time point. To transfer continuous observations (pseudotime, RNA and ATAC gene activity counts, and motif enrichment) from the snMultiome data to the MERFISH data, we assigned the numerical average of the nearest 30 snMultiome cells in the gimVI latent space for each MERFISH cell. Similarly, to transfer categorical observations, including cell types and cluster labels, we assigned each MERFISH cell the most common label in the set of its 30 nearest snMultiome cells. To plot the joint embedding of both modalities (MERFISH and snMultiome) for a single time point, we ran UMAP on the gimVI latent space with n_neighbors and min_dist . Note that in the spatial mapping of cell types predicted in snMultiome (Fig.4c), although the section was in the middle of the leaf, some epidermal cells were included because the section was not completely flat.
smFISH quantification and bacterial colony identification. Quantification of transcripts labelled by smFISH was performed using the Python package Big-FISH. Seven planes of MERSCOPE smFISH images were projected into two dimensions by using numpy.max along the axis and chunked into regions. Spots were then called using the function bigfish.detection.detect_spots with threshold , spot_radius= and voxel_size . To identify bacterial colonies, bigfish.detection.detect_spots was called to label ‘bacterial meta gene’ locations with threshold , log_kernel_size and minimum_distance=(1.456,1.456). To account for autofluorescence from the plant tissue, the same spot-caller was used to call spots on DAPI and ICS1 channels. The spots from all three channels were aggregated, and DBSCAN from scikit-learn.cluster was used on all spots with eps and min_samples=5. We kept all DBSCAN clusters for which the DAPI and ICS1 spots constituted less than of the total cluster spots.

These remaining clusters represented potential bacterial colonies. We manually evaluated each cluster to merge those marking the same colony and removed the clusters marking obvious autofluorescence, for example, signals from stomata. Windows of pixels were generated around each final cluster, and detect_spots was used with threshold=95, spot_radius=(17,17) and voxel_size=(3,3) for accurate quantification of individual bacteria per cluster.

Bacterial neighbourhood analysis. To determine the level of immunity of the neighbourhood around each bacteria colony, we identified the 1-1,000 nearest cells in proximity to the colony. Then, the smoothed imputed pseudotime values of these neighbouring cells were averaged to obtain a single mean pseudotime value for each neighbour size. Error curves were calculated using the standard error divided by the mean at each neighbourhood size. These values serve as indicators of the overall immunity level of the area around each colony.

Spatial differential expression analysis. To identify genes that are spatially associated with immune-rich areas (related to Fig. 5a), we first smoothed our imputed spatial pseudotime over the spatially nearest 100 cells to find areas of immune-active hotspots. We then used RCTD to deconvolve doublets in our spatial data using the cell types in our snMultiome data as a reference. We created a reference object with and ran RCTD on our MERFISH data with the parameters gene_cutoff ,gene_cutoff_reg ,fc_cutoff , fc_cutoff_ reg and doublet_mode ‘doublet’. We then ran C-SIDE of spacexr with run.CSIDE.single with cell_type_threshold to find the top spatially differentially expressed genes per cell type along the smoothed spatial pseudotime (explanatory variable) and plotted the negative log-transformed value against the log-transformed fold change.

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

All information supporting the conclusions are provided with the paper. The single-cell and bulk sequencing data generated in this study have been deposited in the National Center for Biotechnology Information Gene Expression Omnibus database (accession numbers GSE226826 and GSE248054). Reference genomes, annotation and fully processed data for snMultiome analyses are available from the Salk website (http:// neomorph.salk.edu/download/Nobori_etal_merfish). The MERFISH data are available from the Salk website (http://plantpathogenatlas. salk.edu). Genes targeted using MERFISH are listed in Supplementary Table 1. Information about primers used in this study is provided in Supplementary Table 2. Source data are provided with this paper.

The codes to analyse snMultiome and MERFISH data are available at GitHub (https://github.com/tnobori/snMultiome and https://github. com/amonell/Spatial_Plant_Pathogen_Atlas).
40. Mine, A. et al. The defense phytohormone signaling network enables rapid, high-amplitude transcriptional reprogramming during effector-triggered immunity. Plant Cell 30, 1199-1219 (2018).
41. Kaido, M., Tsuno, Y., Mise, K. & Okuno, T. Endoplasmic reticulum targeting of the Red clover necrotic mosaic virus movement protein is associated with the replication of viral RNA1 but not that of RNA2. Virology 395, 232-242 (2009).
42. Mine, A. et al. Differential roles of Hsp 70 and Hsp 90 in the assembly of the replicase complex of a positive-strand RNA plant virus. J. Virol. 86, 12091-12104 (2012).
43. Zhang, X., Henriques, R., Lin, S.-S., Niu, Q.-W. & Chua, N.-H. Agrobacterium-mediated transformation of Arabidopsis thaliana using the floral dip method. Nat. Protoc. 1, 641-646 (2006).
44. Matsumoto, A. et al. A versatile Tn7 transposon-based bioluminescence tagging tool for quantitative and spatial detection of bacteria in plants. Plant Commun. 3, 100227 (2022).
45. Townsley, B. T., Covington, M. F., Ichihashi, Y., Zumstein, K. & Sinha, N. R. BrAD-seq: Breath Adapter Directional sequencing: a streamlined, ultra-simple and fast library preparation protocol for strand specific mRNA library construction. Front. Plant Sci. 6, 366 (2015).
46. Chen, S., Zhou, Y., Chen, Y. & Gu, J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34, i884-i890 (2018).
47. Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15-21 (2013).
48. Liao, Y., Smyth, G. K. & Shi, W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923-930 (2014).
49. Hao, Y. et al. Integrated analysis of multimodal single-cell data. Cell 184, 3573-3587 (2021).
50. Stuart, T., Srivastava, A., Madad, S., Lareau, C. A. & Satija, R. Single-cell chromatin state analysis with Signac. Nat. Methods 18, 1333-1341 (2021).
51. McGinnis, C. S., Murrow, L. M. & Gartner, Z. J. DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. 8, 329-337 (2019).
52. Marand, A. P., Chen, Z., Gallavotti, A. & Schmitz, R. J. A cis-regulatory atlas in maize at single-cell resolution. Cell 184, 3041-3055 (2021).
53. Zhang, Y. et al. Model-based analysis of ChIP-seq (MACS). Genome Biol. 9, R137 (2008).
54. Korsunsky, I. et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat. Methods https://doi.org/10.1038/s41592-019-0619-0 (2019).
55. Cusanovich, D. A. et al. Multiplex single cell profiling of chromatin accessibility by combinatorial cellular indexing. Science 348, 910-914 (2015).
56. Fornes, O. et al. JASPAR 2020: update of the open-access database of transcription factor binding profiles. Nucleic Acids Res. 48, D87-D92 (2020).
57. Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer, 2009).
58. Procko, C. et al. Leaf cell-specific and single-cell transcriptional profiling reveals a role for the palisade layer in UV light protection. Plant Cell https://doi.org/10.1093/plcell/koac167 (2022).
59. O’Malley, R. C. et al. Cistrome and epicistrome features shape the regulatory DNA landscape. Cell 165, 1280-1292 (2016).
60. Nobori, T. et al. Transcriptome landscape of a bacterial pathogen under plant immunity. Proc. Natl Acad. Sci. USA 115, E3055-E3064 (2018).
61. Giacomello, S. & Lundeberg, J. Preparation of plant tissue to enable spatial transcriptomics profiling using barcoded microarrays. Nat. Protoc. 13, 2425-2446 (2018).
62. Pachitariu, M. & Stringer, C. Cellpose 2.0: how to train your own model. Nat. Methods 19, 1634-1641 (2022).
63. Petukhov, V. et al. Cell segmentation in imaging-based spatial transcriptomics. Nat. Biotechnol. 40, 345-354 (2022).
64. Cable, D. M. et al. Robust decomposition of cell type mixtures in spatial transcriptomics. Nat. Biotechnol. https://doi.org/10.1038/s41587-021-00830-w (2021).

Acknowledgements We thank K. Tsuda and J. Walker for comments on the manuscript. T.N. was supported by a Human Frontiers Science Program (HFSP) Long-term Fellowship (LTOO0661/2020-L). J.R.E. is an Investigator of the Howard Hughes Medical Institute.

Author contributions T.N. conceived and designed the study and experiments with guidance from J.R.E. T.N. performed snMultiome and MERFISH experiments and analysed data. T.A.L. provided tissue sections for the MERFISH experiments. A. Monell performed MERFISH data analyses and data integration with assistance from J.Z. Y.S., S.S. and A. Mine produced GT-3A mutants and performed pathogen growth assays. J.R.N. assisted with sequencing and data management. T.N. wrote the initial draft of the manuscript. T.N. and J.R.E. edited the manuscript.

Competing interests J.R.E. serves on the scientific advisory board of Zymo Research and of Cibus. The other authors declare no competing interests.

Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s41586-024-08383-z.
Correspondence and requests for materials should be addressed to Joseph R. Ecker. Peer review information Nature thanks Gitta Coaker, Xuehua Zhong and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.
Reprints and permissions information is available at http://www.nature.com/reprints.

Extended Data Fig. 1| Quality control of snMultiome data. a, Violin plots showing transcripts per nucleus, genes per nucleus, ATAC reads per nucleus, and accessible chromatin regions (ACRs) per nucleus in each sample.b, The number of cells per sample. Two independent replicates were analyzed and combined for Mock and AvrRpt2 9 hpi conditions. c, Density scatterplots showing fraction reads in peaks (FRiP) score in each sample before cell filtering. x-axis: transformed read depths. y-axis: fraction of Tn5 integration sites in ACRs. d, Density scatterplots of transformed read depths ( x -axis) by the fraction of Tn 5 integration sites mapping to within 2 kb upstream and 1 kb downstream of transcription start sites (TSSs) ( y -axis). Data in each sample before cell filtering is shown. c-d, Two replicates of Mock and AvrRpt29 h conditions were combined.e, Bar plot showing the number of ATAC-seq peaks identified in previous bulk ATAC-seq data and the present snATAC-seq data. Shared and assay unique peaks are shown in blue and red, respectively.
f, Sample-aggregated chromatin accessibility around ACTIN2 (left) and ICS1 (right).g, Principle component analysis of pseudobulk transcriptome of each sample. Independent replicates of Mock and AvrRpt2 9 h were labeled. h, Scatter plots comparing pseudobulk transcriptomes of Mock and AvrRpt29 h samples. Pearson’s correlation coefficient values were shown. i, Stacked bar plots showing the representation of gene expression-based Leiden clusters in each sample. j, Two-dimensional embedding of chromatin accessibility similarity among nuclei from all samples with uniform manifold approximation and projection (UMAP). Nuclei are colored by Leiden clusters. k, UMAP embeddings based on a joint neighbor graph that represents both gene expression and chromatin accessibility measurements. Nuclei are colored by de novo Leiden clusters based on the joint analysis (left) and Leiden clusters defined by gene expression measurement alone (right; Fig.1b).

a

Extended Data Fig. 2 | Single-cell analysis of gene expression and chromatin accessibility. a, Dot plot showing the top marker genes of individual clusters. b, Cluster-aggregated chromatin accessibility surrounding , a known marker gene for leaf epidermis. c, UMAP plot showing the ATAC-seq count
on a peak near FDH(chromosome 2, position 11172821-11173529) in each nucleus. d, GO enrichment analysis for marker genes of each cluster. e, Expression of ICS1 mRNA in each nucleus in each sample. Adjusted p-values from a one-sided hypergeometric test followed by Benjamini-Hochberg correction are shown.

Extended Data Fig. 3 | Comprehensive characterization of distinct cell populations. a, Heatmap showing pair-wise correlation of pseudobulk transcriptomes between sub-clusters of individual clusters. The top and side bars show major cluster labels.b, (Left) Schematic workflow for the selection of highly variable immune-related genes. First, genes that are significantly enriched in at least one of major or sub clusters were selected. Then, these genes were clustered based on pseudobulk expression of sub-clusters using k-mean clustering (k value was determined by the elbow method), followed by GO enrichment analysis of each k-mean cluster. Finally, gene clusters with enriched immunity-related GO terms were selected. (Right) Top GO terms enriched in 12 k -mean clusters. Clusters 1,2 , and 10 were selected as immune genes.c, Heatmap showing normalized expression of immune-related genes selected in (b) across all the sub clusters. The top bars indicate cell type and major cluster that each sub-cluster derived from. d, Expression of the phloem
companion cell markers SUC2 and FTIP1. e, Expression of ALD1 and FMO1. f, Expression of ILL6 upon infection of AvrRpt2 in time course, showing specific induction in cluster 6.g, Expression of genes specifically expressed in subcluster 6-8.h, Heatmap showing expression of immune-related genes shown in Fig. 1d in a previous time-course bulk RNA-seq study, where . thaliana leaves were infected by DC3000, AvrRpt2, or AvrRpm1 or treated with mock control (water). These data do not capture the presence of various cell populations identified in Fig.1d.i, Scatter plots showing the correlation between CRK12 and CRK41 in the time course bulk RNA-seq (top) and single-nucleus RNA-seq (snRNA-seq; bottom). For snRNA-seq analysis, pseudobulk gene expression data of the subclusters (shown in c) were used. , Expression of CRK12 and CRK41 upon infection with AvrRpm1 as shown in the UMAP. These apparently correlated genes in bulk RNA-seq showed a highly different expression pattern at the single-cell level.

Extended Data Fig. 4 | Linking gene expression and chromatin accessibility at the single-cell level. a, Schematic diagram of linked accessible chromatin regions (ACRs) and a gene. An ACR and a gene are “linked” when there is a significant correlation between chromatin accessibility and mRNA expression across individual cells.b, Density plot showing the frequency of linkages at different distances from the transcription start sites (TSSs). c, Heatmap showing the linkage score (Pearson correlation coefficient between ACR count and mRNA expression) for genes that showed at least one significant link in at least one of the infection conditions or Mock. When a gene had multiple links, the maximum linkage score was shown. The sidebar shows the k-mean cluster
annotation (the value was determined by the elbow method). d, Expression of mRNA encoding VSP2, UGT85A1, and MAM1. e,f, Cluster-aggregated chromatin accessibility surrounding CBP6Og (e) and AT1G56660 (f). Violin plots on the side show aggregated mRNA expression of each gene. Both genes were highly expressed in a cell-specific manner, but only CPB6Og showed correlated (linked) chromatin accessibility patterns.g,h, Top motifs enriched in the promoter regions ( 2 kb upstream from the TSS) of defense genes (markers of immune-active mesophyll and epidermis clusters , and 12 ) that are linked (g) and not linked (h).

Extended Data Fig. 5 | Links between transcriptome and chromatin accessibility. a, Scatter plot showing the distribution of the linkage score (Pearson correlation coefficient between ACR count and mRNA expression) at different distances from TSSs.b, GO enrichment analysis of genes in each k -mean cluster shown in (Extended Data Fig. 4c). Adjusted p-values from a one-sided hypergeometric test followed by Benjamini-Hochberg correction are shown. c, Density plot showing the frequency of linkages at different distances from TSS for cluster marker genes and non-marker genes. d, Density plot showing the frequency of the size of ACRs for linked or non-linked genes.
e, Density plot showing the frequency of genes with different numbers of linked ACRs.f, Schematic diagram of the ATAC activity analysis. For each gene, ATAC reads mapped on the gene body or the 400 bp upstream region were aggregated to calculate the score.g, Scatter plot showing the relationship between the number of linkages ( x -axis) and RNA-ATAC activity correlation score ( y -axis). h, Scatter plot showing the relationship between the maximum correlation coefficient between each gene and linked ACRs ( -axis) and RNA-ATAC activity correlation score ( -axis).

showing (a) enrichment scores of 465 motifs and (b) expression of corresponding transcription factors (TF) in each nucleus. The top bars show cluster annotation defined in Fig. 1b. c, Heatmap showing Pearson correlation coefficient between
motif enrichment scores and mRNA expression of the corresponding TFs in each cell type (shown in the top bar). All TFs are shown. d, GO enrichment analysis of genes predicted to be targeted by TFs shown in Fig. 2h.

Sample (# transcripts detected)

Extended Data Fig. 7|MERFISH data analysis. a, Example raw images of MERFISH that capture the same region of tissue across three imaging rounds. White spots are signals derived from single mRNA molecules. We observed similar results in all independent samples. b, Two-dimensional plots of all the transcripts detected by MERFISH in each sample. The number of transcripts is shown.c, Spatial expression pattern of ALD1 detected by MERFISH in each sample.d, (Left) A field of view (FOV) that shows obscure DAPI nuclei staining signal. (Middle) Transcript-based segmentation in the same FOV (see Method). Transcripts were colored by assigned cells. (Right) Centroids of cells detected
by the transcript-based segmentation (red dots) and the result of failed DAPI-staining-based segmentation (yellow region). Similar patterns were observed across FOVs and samples. A systematic quantitative analysis is provided in (f). (a,d) Scale bars m.e, , The fraction of transcripts in Cellpose-segmented cells fell within Baysor-segmented cells. (e) Cells detected in the field of view (FOV) shown in (d) were used. (f) Cells detected in 100 FOVs were used. g, Integrated UMAP of all MERFISH data.h, Spatial mapping of de novo MERFISH clusters annotated as vasculature.b,c,h Scale bars .

Extended Data Fig. 8 |Spatial mapping of whole transcriptome and epigenome. a, Imputed mRNA expression of ALD1 (Top) predicted its mRNA expression pattern measured with MERFISH (Bottom). Magnified images are shown on the right.b, Imputed mRNA expression of ICS1 (Top) predicted its mRNA expression pattern measured with smFISH (Bottom). Magnified images are shown on the right. c, Imputed ATAC activity (Extended Data Fig. 5f) of ICS1 (top) and ALD1 (bottom), which showed consistent patterns with mRNA expression (a,b).d, Imputed motif enrichment scores of HsfB2b (Top) was consistent with mRNA expression of confirmed by MERFISH (Bottom). e, Spatial mapping of pseudotime values based on data integration and label transfer between snRNA-seq and MERFISH data of DC3000-infected plants.
f, smFISH of ICS1 in leaves infected by AvrRpt2 or DC3000 at 24 hpi.g, Spatial mapping of bacterial transcripts detected with smFISH in plants infected by AvrRpt2 (left) and DC3000 (right) at 24 h post-infection (hpi). Pseudotime values are also visualized in the background. Dot size reflects the number of bacterial transcripts detected. a-g, Scale bars . h, Scatter plot showing the number of bacterial transcripts ( x -axis) and averaged pseudotime values of five nearest neighbor plant cells ( -axis) for each bacterial colony in plants infected by AvrRpt2 (blue) and DC3000 (red) at 24 hpi . i, Averaged pseudotime scores of cells surrounding each bacterial colony. The x-axis shows the number of nearest plant cells analyzed for each bacterial colony. Shaded error bands indicate standard error of the mean.

Extended Data Fig. 9 | Characterization of PRIMER and bystander cells.
a, Heatmaps showing expression of ALD1, FMO1, and BON3 in each cell at 9 hpi . b, Expression of shown in the UMAP in Figs. 4a and 5d. BON3 expression is enriched in cells with the highest Pseudotime scores. c, Expression of WRKY8 and LSD1 in the sub-clusters of clusters 3, 7, and 11 in the snRNA-seq data (Fig.1f). d, Violin plot showing expression PUB36 in the PRIMER cell cluster of WT and GT-3A knockout (gt3a-KO) plants. The gene model of PUB36, ATAC-seq peaks, and a GT-3A binding motif are shown below. e, GO enrichment analysis of genes downregulated in a bystander cell cluster of gt3a-KO compared to WT. Adjusted p-values from a one-sided hypergeometric test followed by Benjamini-Hochberg correction are shown.f, Violin plot showing expression of ALD1 in a bystander cluster. d-f, Two independent replicates of each sample were analyzed with single-nucleus RNA-seq. g, Expression of BON3 and GT-3A in cells in mock (top)
or DC3000-infection (bottom) condition. h, Expression of/CS1 in cells infected by AvrRpt2. c,g,h, All time points were combined. i, Motif activity of GT-3A in immune-active mesophyll cells.j, Lesion area created by Colletotrichum higginsianum infection in Col-0 and two independent GT-3A overexpression lines and pad3 mutant at 6 days post inoculation. Different letters indicate statistical significance (adjusted ). Adjusted p -values were calculated with two-tailed Student’s t-test followed by Benjamini-Hochberg method. , and 51 leaves for WT, GT-3A-ox1, GT-3A-ox2, and pad3, respectively, from three independent replicates. Results are shown as box plots with boxes displaying the 25th-75th percentiles, the centerline indicating the median, whiskers extending to the minimum, and maximum values no further than 1.5 interquartile range. , The gene models of CBP6Og and SARD1 with ATAC-seq peaks and GT-3A binding motifs.

Corresponding author(s): Joseph R. Ecker

Last updated by author(s): Nov 7, 2024

Nature Portfolio wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Portfolio policies, see our Editorial Policies and the Editorial Policy Checklist.

For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.
n/a
□

The exact sample size for each experimental group/condition, given as a discrete number and unit of measurement
□

A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly
□

The statistical test(s) used AND whether they are one- or two-sided
Only common tests should be described solely by name; describe more complex techniques in the Methods section.

□ A description of all covariates tested
□

A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons
□

A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals)
□

For null hypothesis testing, the test statistic (e.g. ) with confidence intervals, effect sizes, degrees of freedom and value noted Give values as exact values whenever suitable.

□ For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings

□ For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes
□ X
Estimates of effect sizes (e.g. Cohen’s , Pearson’s ), indicating how they were calculated
Our web collection on statistics for biologists contains articles on many of the points above.
Confirmed
exact sample size for each experimental group/condition, given as a discrete number and unit of measurement
□
statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly

The statistical test(s) used AND whether they are one- or two-sided
A description of all covariates tested
L


ast updated by author(s): sistency and transparency

All manuscripts must include a data availability statement. This statement should provide the following information, where applicable:

All information supporting the conclusions are provided with the paper. The single-cell and bulk sequencing data generated in this study are deposited in the National Center for Biotechnology Information Gene Expression Omnibus database (accession no. GSE226826 and GSE248054). Reference genome, annotation, fully processed data for snMultiome analyses are available at neomorph.salk.edu/download/Nobori_etal_merfish. The MERFISH data are available at plantpathogenatlas.salk.edu. Genes targeted with MERFIHS are listed in Supplementary Table 1. Information of primers used in this study is provided in Supplementary Table 2. Source data are provided with this paper. Reference genome, annotation, fully processed data for snMultiome analyses are available at neomorph.salk.edu/download/Nobori_etal_merfish.

Policy information about studies involving human research participants and Sex and Gender in Research.

Note that full information on the approval of the study protocol must also be provided in the manuscript.

Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection.
Life sciences □ Behavioural & social sciences □ Ecological, evolutionary & environmental sciences
For a reference copy of the document with all sections, see nature.com/documents/nr-reporting-summary-flat.pdf

All studies must disclose on these points even when the disclosure is negative.

nature portfolio | reporting summary March 2021

We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response.

nature portfolio | reporting summary March 2021