Strain-Specific Tropism and Transcriptional Responses of Enterovirus D68 Infection in Human Spinal Cord Organoids

Figure 1.. Nextstrain Phylodynamics of Enterovirus D68 Clade B. (A) Phylogenetic tree of publicly available EV-D68 Clade B sequences based on VP1 region (N = 1311 sequences; collection dates: April 2008 to December of 2024), highlighting the two strains used in this study - US/IL/14–18952 (18952) and US/MA/18–23089 (23089). (B) Temporal distribution of EV-D68 Clade B sequences, shown as the proportion of sequences per year.

Figure 2.. Cell type distribution in 24-day-old hSCOs. (A) Design of scRNAseq experiment. Pooled hSCOs (n=12) were infected with EV-D68 strains, US/IL/14–18952 (18952) or US/MA/18–23089 (23089) (105 PFU/group) for 48 hours at 33°C, then dissociated mechanically and enzymatically. Cells were partitioned using the 10x NextGEM procedure and resulting libraries were sequenced using short-read sequencing. Flowchart created with BioRender. (B) Uniform Manifold Approximation Projection (UMAP) of individual cells identified from all the samples collected in this study colored by assigned cell type (41,098 cells). (C) Donut plot showing the relative abundance of cell types captured across all experiments.

Figure 3.. Infection profiling in hSCO infected with EV-D68. (A) UMAP of cellular transcriptional phenotypes in 24-day-old organoids infected with EV-D68 strains 18952 and 23089. Cells positive for viral RNA are shown colored by viral RNA content. (B) Bar plots displaying the proportion of vRNA+ cells vs. Bystander cells in various cell types per each condition, highlighting significant differences with an asterisk. (C) Stacked bar charts showing viral read counts in each cell type. The number of total stacked cells are indicated on the top of each bar.

Figure 4.. Strain-specific enrichment of host pathways in infected cells. (A) Heatmap of normalized enrichment scores (NES) from gene set enrichment analysis (GSEA) comparing two aggregated cell groups: 18952-Infected and 23089-Infected. Pathways shown represent all significantly enriched Gene Ontology Biological Process (GO-BP) terms (adjusted p < 0.05) including both positively and negatively enriched gene sets. Negative NES values (blue) reflect pathways more enriched in Mock cells relative to the infected group. (B) Dot plots showing the expression of key genes involved in neuronal signaling and communication, proliferation and transcriptional activity, and suppressed translation and metabolism processes across distinct cell types in hSCOs infected with EV-D68 strains 18952 and 23089. Dot size represents the percentage of cells expressing each gene, while color intensity reflects the average expression level, scaled (z-score) across all cells based on log-normalized values. Midplate and VLMCs were excluded from the Mock condition, as they were either absent from the infected groups or detected in only one of them (e.g., VLMCs: 1 cell in EV-D68–18952).

Figure 5.. Strain-specific transcriptional responses in bystander cells. (A) Bar plot showing the number of significantly differentially expressed genes (DEGs) between EV-D68 strains 18952 and 23089 across bystander neurons, astrocytes, cycling astrocytes, and OPCs. DEGs were defined as adjusted p-value < 0.05 and log2FC > 0.25. B–D). Volcano plots displaying DEGs between 23089- and 18952-exposed bystander cells in neurons (B), astrocytes (C), and cycling astrocytes (D). Genes upregulated in strain 23089 are shown in green, and those upregulated in strain 18952 are in orange; non-significant genes are shown in gray. Dashed lines indicate significance thresholds (adjusted p = 0.05 and log2FC = ±0.25). (E) Bar plot showing normalized enrichment scores (NES) from GSEA of OPC Oligo. bystander cells, comparing transcriptional profiles between strains 23089 and 18952. Despite the absence of individual DEGs meeting significance, GSEA revealed pathway-level differences. Bars are colored by the direction of enrichment: green for pathways enriched in 23089 and orange for pathways enriched in 18952.