Spatial metatranscriptomics resolves host-bacteria-fungi interactomes

Abstract

The interactions of microorganisms among themselves and with their multicellular host take place at the microscale, forming complex networks and spatial patterns. Existing technology does not allow the simultaneous investigation of spatial interactions between a host and the multitude of its colonizing microorganisms, which limits our understanding of host-microorganism interactions within a plant or animal tissue. Here we present spatial metatranscriptomics (SmT), a sequencing-based approach that leverages 16S/18S/ITS/poly-d(T) multimodal arrays for simultaneous host transcriptome- and microbiome-wide characterization of tissues at 55-µm resolution. We showcase SmT in outdoor-grown Arabidopsis thaliana leaves as a model system, and find tissue-scale bacterial and fungal hotspots. By network analysis, we study inter- and intrakingdom spatial interactions among microorganisms, as well as the host response to microbial hotspots. SmT provides an approach for answering fundamental questions on host-microbiome interplay.

Fig. 1. Overview of the method.

a, SmT uses capture arrays on glass slides. Each capture array contains 4,992 spots that are 55 µm in diameter and 100 µm from center to center. Cells are permeabilized to release RNA molecules that hybridize to the barcoded capture probes in the spots. The captured molecules are then processed into a sequencing library. b, Capture probes consist of a sequencing adapter, a spatial barcode, a UMI and a capture moiety. Polyadenylated mRNAs are captured with poly-d(T) probes that comprise 10% of all the capture probes. Ribosomal RNAs from fungi are captured with P-ITS7 and P-ITS1 probes targeting the 18S rRNA and ITS regions, respectively. Ribosomal RNAs from bacteria and archaea are captured with P479, P799, P902 and P1205 probes targeting bacterial 16S rRNA. Bacterial and archaeal probes and fungal probes each comprise 45% of the capture probes. c, A bioinformatic workflow designed to assign the reads to host or microbial modalities. First, low-quality reads are filtered out, the remaining reads are mapped against the A. thaliana TAIR10 reference genome, and spatial barcodes are demultiplexed. Second, mapped A. thaliana reads are filtered based on their UMI and compiled to obtain a gene-count matrix. Third, the reads not mapping to A. thaliana are mapped to a universal database to remove those that are not clearly of microbial origin. The remaining microbial reads are classified with LCA based on their identity and UMIs, and unique taxa are counted to generate separate unique taxa-count matrices for fungi and for bacteria and archaea.

Fig. 2. SmT resolves the microbial profile and host transcriptome at microscopic resolution.

a, A Toluidine blue-stained bright-field image of a 14-µm thick longitudinal A. thaliana leaf section (left) and the fluorescent cDNA footprint (right) of the same section from the tissue optimization experiment. b, Left, fluorescence image of an intact A. thaliana leaf syringe-infiltrated (yellow square) with mCherry-tagged Pst DC3000 bacteria. Middle, a 14-µm thick longitudinal section from the same leaf was analyzed using a 50% poly-d(T), 25% P799 and 25% P902 array, revealing the spatial capture of Pst DC3000 16S rRNA molecules. Right, spatial distribution of PR1 gene expression in the same leaf section. Scale bars: 1 mm. c, Pearson correlation coefficient and the corresponding two-tailed significance test of bacterial 16S rRNA, eukaryotic 18S rRNA/ITS and A. thaliana molecules captured in leaf 1 with a multimodal array containing 10% poly-d(T), 45% 16S rRNA and 45% 18S rRNA/ITS probes and with 100% 16S rRNA or 18S rRNA/ITS and poly-d(T) arrays. In all correlations, P = 0. d, Bray–Curtis similarity for bacterial and fungal taxa captured on different arrays, organized by hierarchical clustering. e, Experimental validation of SmT by amplicon sequencing. f,g, Numbers of bacterial and archaeal taxa detected using the two methods in a representative sample of four leaves from two plants are compared qualitatively using a Venn diagram (f) and quantitatively using NMDS (g). NMDS, non-metric multidimensional scaling.

Fig. 3. Microbial interactions are driven by spatial organization.

a, Bacterial and fungal profiles for each of the sections of four leaves (‘L’) from two plants (‘P’). ‘Other’ denotes binned bacterial and archaeal genera and fungal genera having ≤1% abundance. b, Numbers of unique microbial taxa per capture spot. c, Significant hot- and cold-spots for bacteria and fungi in a representative leaf section. NS, not significant. Scale bars: 500 µm. d, Percentages of shared and unique hotspots among bacteria and fungi across 13 different leaf sections. Of note is the variance in shared regions across the sections. The sections were taken from two leaves from two different plants (four leaves in total); ‘P’ denotes plant and ‘L’ denotes leaf. e, The proportion of interkingdom (bacteria–fungi) interactions as a function of the proportion of shared interkingdom hotspots and the number of reads. SRCC ρ = 0.72, P = 0.0058 (two-tailed test).

Fig. 4. Host response is associated with microbial colonization pattern.

a, UMAP clustering of host gene expression. b, Projection of UMAP clusters mesophyll 1 and vascular on a representative leaf section. c, Mesophyll and vascular cell-type proportions projected on a representative leaf tissue section. d, Overrepresented GO terms for microbial-associated genes (n denotes the number of genes labeled with the indicated GO term). e, Spatial distributions of significant bacterial and fungal hotspots together with hotspots for the expression of the defense-related genes CA1 (AT3G01500), LURP1 (AT2G14560) and ACD6 (AT4G14400). NS, not significant. Scale bars: 500 µm.