Strain-specific tropism and transcriptional responses of enterovirus D68 infection in human spinal cord organoids

The mechanisms by which Enterovirus D68 (EV-D68) infection leads to acute flaccid myelitis (AFM), a severe neurological condition characterized by sudden muscle weakness and paralysis, remain poorly understood. To investigate the cellular tropism and infection dynamics of EV-D68, we profiled naive and EV-D68-infected human spinal cord organoids (hSCOs) derived from induced pluripotent stem cells (iPSCs) using single-cell RNA sequencing (scRNA-seq). Examining the cellular composition of healthy hSCOs, we found that hSCOs include diverse cell types, including neurons, astrocytes, oligodendrocyte progenitor cells (OPCs), and multipotent glial progenitor cells (mGPCs). Upon infection with 2 EV-D68 strains, US/IL/14-18952 (a B2 strain) and US/MA/18-23089 (a B3 strain), we observed distinct viral tropism and host transcriptional responses. Notably, US/IL/14-18952 showed a significant preference for neurons, while US/MA/18-23089 exhibited higher rates of infection in cycling astrocytes and OPCs. These findings provide novel insights into the host cell tropism of EV-D68 in the spinal cord, offering insight into the potential mechanisms underlying AFM pathogenesis. Understanding the dynamics of infection at single-cell resolution will inform future therapeutic strategies aimed at mitigating the neurological impact of enterovirus infections.

Introduction

Acute flaccid myelitis (AFM) is a polio-like illness characterized by muscle weakness and paralysis, primarily affecting children (Murphy et al., 2021; CDC, 2025). Increased cases of AFM suspected to be due to enterovirus infection were first recorded in 2014 (Aliabadi et al., 2016; Messacar et al., 2016, 2018) and many of such cases have been associated with Enterovirus D68 (Taxonomy: Enterovirus deconjuncti), or EV-D68 (McKay et al., 2018; Mishra et al., 2019; Schubert et al., 2019; Vogt et al., 2022). Like most enteroviruses, the vast majority of EV-D68 infections are asymptomatic or mild, making it difficult to estimate the total number of cases. Flu-like and respiratory symptoms occurring in a small proportion of individuals, perhaps between 30% and 40% (Teoh et al., 2025), potentially associated with preexisting asthma or reactive airway disease (Clopper et al., 2025). Neurologic illness, including AFM, occurs in a much smaller proportion, perhaps around 0.1% of cases, although the prevalence is difficult to estimate due to the number of asymptomatic cases, the infrequency of identified pathogen, and the temporal variability of AFM cases, perhaps due to strain differences year to year or other factors (Freeman and Messacar, 2025).

EV-D68 and AFM had coinciding biennial outbreaks from 2014 to 2018 (Shah et al., 2021), with the 2018 AFM outbreak associated with nearly twice as many confirmed AFM cases compared to 2014 (CDC, 2025). Another spike of cases was expected in 2020, but transmission was likely impeded by isolation policies during the SARS-CoV-2 pandemic (Olsen et al., 2021). EV-D68 had an additional outbreak in 2022, but AFM cases did not increase as expected (CDC, 2025; New Vaccine Surveillance Network Collaborators and Hall, 2022).

The mechanism by which enterovirus infection contributes to the development of AFM is unknown. Previous studies have suggested both direct damage to spinal cord neurons after viral infection, and subsequent cytotoxic T-cell responses to infected neurons both contribute to disease (Vogt et al., 2022; Woods Acevedo et al., 2024). One fundamental question relevant to EV-D68 pathogenesis is the cell types that contribute to virus replication and production in the CNS. Studies in multiple model systems have demonstrated that EV-D68 can target and replicate in neurons (Wang et al., 2004; Feng et al., 2016; Hixon et al., 2017; Brown et al., 2018; Rosenfeld et al., 2019). Astrocyte infection has also been identified during EV-D68 infection of murine brain slice cultures and primary human astrocytes in vitro (Rosenfeld et al., 2019; Liu et al., 2025).

We have previously shown that contemporary strains of EV-D68 can replicate in an induced pluripotent stem cell (iPSC)-derived human spinal cord organoid (hSCO) model, which provides a valuable human-derived, multicellular model in which to explore EV-D68 pathogenesis (Aguglia et al., 2023). Analysis of marker gene expression suggests hSCOs comprise multiple cell types, including neurons and glial cells, but the specific cell types present, and which contribute to enterovirus infection in the hSCO model are unknown (Aguglia et al., 2023).

To better define the cell types infected by EV-D68, and to identify potential strain-specific differences in cell tropism and pathogenesis, we infected hSCOs with two contemporary strains of EV-D68 associated with AFM, US/IL/14-18952 (18952) and US/MA/18-23089 (23089). These strains are genetically distinct and represent circulating strains from the 2014 and 2018 outbreak years (Brown et al., 2014; Hadfield et al., 2018; Sagulenko et al., 2018) (Figure 1). Across the coding sequence, the two isolates share ~92% nucleotide and 99% amino acid identity, consistent with their classification as distinct subclades within the EV-D68 B lineage (Supplementary Table 1). These strains were selected to reflect clade-level diversity among EV-D68 isolates associated with AFM outbreaks, with 18952 representing clade B2 and 23089 representing clade B3. Their use also aligns with the availability of well-characterized viral stocks, and were previously used in the development of the hSCO infection model (Aguglia et al., 2023). Notably, in these previous studies, representatives of early isolates, such as Fermon (1962), and clade B1 viruses did not exhibit robust infection in the hSCO model. While not intended to represent the full diversity of circulating strains, these two representatives provide a useful comparison of potential strain-specific tropism and host responses.

Figure 1

Figure 1. Next strain Phylodynamics of enterovirus D68 clade B. (A) phylogenetic tree of publicly available EV-D68 Clade B sequences based on VP1 region (N = 1,311 sequences; collection dates: April 2008 to December of 2024), highlighting the two strains used in this study USA/IL/2014-18952 (18952) and USA/MA/2018-23089 (23089). (B) Temporal distribution of EV-D68 Clade B sequences, shown as the proportion of sequences per year.

Using single-cell RNA sequencing (scRNAseq), we captured host and viral transcripts within different cell types and subtypes. Our analysis revealed the complex cellular composition of hSCOs, which includes neuronal and glial cell lineages. Analyzing viral transcript abundance within these cell populations demonstrated that the EV-D68 strains exhibit markedly different tropisms. Although EV-D68 18952 (B2) was primarily associated with neuronal infection, 23089 (B3) exhibited a preference for cycling astrocytes. Subsequent analysis of host cell transcriptional responses in both infected and bystander cell populations revealed further differences between these strains. Together, these findings clarify the shifting cell tropism of EV-D68 strains and provide insight into the mechanism of AFM pathogenesis in a highly relevant human model system.

Results

24-day-old spinal cord organoids comprise diverse cell lineages

Previous characterization of marker gene and protein expression in hSCOs suggested the presence of several fully-differentiated cell types, including neuronal lineages, such as motor neurons and interneurons, and roof plate-like structures and neuro epithelium (Ogura et al., 2018; Aguglia et al., 2023). However, such analyses of gene expression may not identify minor cell types or developmental intermediates of specific cell lineages, and do not provide quantitative measures of relative cell abundance. Therefore, to better characterize the diversity and abundance of individual cell types in these organoids, we performed scRNAseq. We generated single-cell suspensions by dissociating pools of 12 hSCOs differentiated for 24 days, capturing between 8,000 and 10,000 cells per sample for sequencing (Figure 2A).

Figure 2

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) (10⁵ PFU/group) for 48 h at 33 °C, then dissociated mechanically and enzymatically. Cells were partitioned using the 10× Next GEM procedure and resulting libraries were sequenced using short-read sequencing. Flowchart created with Bio Render. (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.

Following scRNAseq, we used gene expression profiles from individual cells to perform cell type identification. We based our assignments on curated cell types from a previously published study of the developing human spinal cord (Andersen et al., 2023) through “label transfer” (Stuart et al., 2019) (Figures 2B,C). This analysis revealed that hSCOs exhibit a diverse cellular composition, encompassing neuronal lineages, along with astrocytes, oligodendrocytes and their progenitors (OPCs), glial populations such as mid plate cells and multipotent glial progenitor cells (mGPCs), and vascular leptomeningeal cells (VLMCs), which contribute to blood–brain barrier integrity.

Astrocytes and cycling astrocytes were identified as the majority cell type in the 24-old day hSCOs (54.3%). We confirmed their cell identity by assessing expression of known astrocyte-specific markers, including SOX9, FGFR3 (Pringle et al., 2003; Stolt et al., 2003) and TOP2A (Lee and Berger, 2019) (Supplementary Figure S1). Neurons comprised the next most common fully-differentiated cell type (13.9%), which we confirmed by examining expression of MAP2 and TUBB3 (Menezes and Luskin, 1994; Dehmelt and Halpain, 2005) (Supplementary Figure S1). In addition to differentiated cell types, mGPCs also made up a considerable proportion of the hSCO composition (20.1%) (Figures 2B,C; Supplementary Figure S2A). These proportions were consistent across replicate pools of hSCOs (Supplementary Figure S2B) of 24-days-old hSCOs. Fewer than 10% of cells were not classified as a specific cell type due to low predicted cell type score (Figure 2, “Unclassified”). These were not considered for subsequent analysis.

Focused reanalysis of the mGPC subset revealed five distinct clusters which may represent intermediate states in the differentiating organoid (Supplementary Figure S3A). Consistent with this interpretation, we have found the proportion of mGPCs are reduced in hSCOs in later development days (data not shown). mGPCs in clusters 0 and 4 display features of neural progenitors and early neuronal differentiation, with Cluster 0 exhibiting expression of genes linked to mature neuronal identity, while Cluster 4 shows expression patterns indicative of active neurogenesis. Among the expression differences, we identified DCX and NEUROG1 in Cluster 0, and NKX1-1 and GATA2 in Cluster 4, reflecting their roles in early neural development (Schubert et al., 1995; Gleeson et al., 1999; Ma et al., 1999; Nardelli et al., 1999). Cells in Cluster 0 showed high expression of NEUROD4 and ELAVL3 (Okano and Darnell, 1997; Miyata et al., 1999; Mulligan and Bicknell, 2023), markers of mature neurons. In contrast, Cluster 4 displays a broader range of functions, including neurotransmitter synthesis (GAD2 and SLC32A1) (Erlander and Tobin, 1991; McIntire et al., 1997), and cell adhesion and signaling (GPR83 and CNTNAP5) (Müller et al., 2013; Gomes et al., 2016; Chatterjee et al., 2019), indicating a more diverse cell population (Supplementary Figure S3B). Clusters 1, 2, and 3 are likely astrocyte progenitor cells, with diverse gene expression profiles related to stress response, lipid metabolism, cell signaling, extracellular matrix remodeling, DNA repair, cell cycle regulation, neural development, and cell polarity. Cluster 2 exhibited markers of proliferation, more consistent with cycling astrocytes.

Two contemporary EV-D68 strains show distinct tropism in hSCOs

Upon infection with two distinct strains of EV-D68 [US/IL/14-18952 (B2) and US/MA/18-23089 (B3)], we observed similar patterns of broad cellular susceptibility across the three main cell types identified within the hSCOs (Figure 3A). Overall, we identified 144 viral RNA (vRNA)-positive cells out of 8,698 total cells in the 18952 (B2)-infected organoids, representing 1.65% of the captured cell population. For strain 23089 (B3), we observed 205 vRNA-positive cells out of 15,805 total cells, corresponding to 1.3% of the captured cells. The small proportion of infected cells is similar to our prior microscopic observations of hSCO infection with these strains (Aguglia et al., 2023). Although the proportion of infected cells is similar, the two strains exhibited distinct preferences for specific cell types. EV-D68-18952 (B2)-infected cells were significantly enriched in neurons (based on permutation tests), whereas EV-D68-23089 (B3)-infected cells were significantly enriched among cycling astrocytes and oligodendrocyte progenitor cells (OPC Oligo.) (Figure 3B). Notably, these enrichments were also consistent with the number of viral RNA reads (vRNA) originating from these cell types despite significant heterogeneity in vRNA reads per individual cell (Figure 3C). Together these observations suggest marked strain differences in host cell preference in the spinal cord, which may contribute to differential pathogenesis although the consequences are not clear from this observation alone.

Figure 3

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.

Strain-specific enrichment analysis reveals distinct responses to EV-D68 infection

To address the distinct cellular and molecular responses triggered by the two strains of EV-D68 [18952 (B2) and 23089 (B3)], we performed a comprehensive enrichment pathway analysis on infected spinal cord organoids to identify strain-specific transcriptional programs. Due to the low proportion of vRNA-positive cells, we combined all cell types to compare gene expression differences based on the infecting strain. While this approach does not resolve cell type-specific changes, it allows for a more statistically robust comparison of strain-specific differences in host response. The normalized enrichment scores (NES) of significantly enriched pathways (adjusted p < 0.05) revealed three major functional clusters (Figure 4A), capturing pathways where genes are primarily up-regulated (Cluster 1) or down-regulated (Cluster 3) in response to infection with either strain, or pathways primarily upregulated in response to 23089 (B3) infection relative to 18952 (B2)-infected (Cluster 2). Notably, our analysis did not identify any pathways enriched only in the context of 18952 (B2) infection.

Figure 4

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. Mid plate 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).

The most striking difference between the two strains is in Cluster 2, which corresponded to proliferative and biosynthetic transcriptional programs uniquely enriched in 23089 (B3)-infected cells. This included pathways such as DNA replication, chromosome organization, cell cycle, and ribosome biogenesis (Figure 4A). High NES values were observed for DNA-templated DNA replication (NES = 2.07), DNA repair (NES = 1.56), and cell cycle process (NES = 1.52). Despite 18952 (B2)-infected cells showing slightly higher expression of cell cycle markers per cell (Figure 4B), the enrichment in 23089 (B3) is likely driven by the larger number of infected cycling astrocytes, increasing the overall representation of these proliferative pathways. 23089 (B3)-infected cells showed modest upregulation of RNA processing and splicing pathways across all infected cell types (e.g., SNRNP70, HNRNPA1 and SF3B1) (Figure 4B). These findings suggest that 23089 (B3) induces a transcriptionally-active and biosynthetically engaged state, possibly reflecting its preferential infection of cycling astrocytes.

Clusters 1 and 3 captured similarities in response to the two strains. Pathways identified in Cluster 1 included neuronal development and communication—such as synapse assembly, axon development, neurotransmitter secretion, G protein-coupled receptor signaling, and neurogenesis—which were positively enriched in both 18952 (B2)- and 23089 (B3)-infected cells, consistent with the neurotropic nature of EV-D68. However, strain 18952 (B2) showed higher enrichment for synaptic signaling (NES = 1.87 vs. 1.13) and neurotransmitter transport (NES = 1.95 vs. 1.13), suggesting a more robust neuronal response (Figure 4A). Cluster 3 encompassed pathways related to translation, metabolism, and developmental signaling, with 18952 (B2) showing stronger negative enrichment, including cytoplasmic translation (NES = −2.22), AMP metabolic process (NES = −2.26), and ubiquitin ligase regulation (NES = −2.30). While genes encoding ribosomal proteins (e.g., RPL10, RPLP1 and RPS6) appeared downregulated in neurons from both infected and mock conditions (Figure 4B), broader suppression of metabolic homeostasis pathways was more pronounced in the infected cells. Together, these clusters suggest that both strains impact neuronal function and basal cellular processes, but with greater intensity in 18952 (B2)-infected cells.

These transcriptional differences align with the distinct cellular tropism observed between the 2 EV-D68 strains. In organoids infected with 18952 (B2), neurons represented the major significantly enriched infected cell population, which is consistent with the prominent upregulation of synaptic signaling, axonal development, and neurotransmitter-related pathways. The robust enrichment of these neuronal processes suggests a direct viral impact on neurons, potentially contributing to neuronal dysfunction or degeneration. In contrast, 23089 (B3) predominantly infected cycling astrocytes and, to a lesser extent, oligodendrocyte progenitor cells (OPCs). This cell-type specificity is mirrored by the enrichment of DNA replication, cell cycle progression, and RNA processing pathways in 23089 (B3)-infected cells—hallmarks of transcriptionally active, proliferative glial populations. Thus, the pathway signatures not only reflect divergent viral-host interactions but also underscore how strain-specific cellular targeting shapes the overall transcriptional landscape of infected organoids.

Distinct transcriptional impacts of EV-D68 strains across bystander cell types

Differential gene expression (DEG) analysis was performed on bystander cells stratified by cell type. A summary of the number of significant DEGs identified per cell type is presented in Figure 5A. This panel reflects DEGs filtered by both adjusted p-value <0.05 and absolute log₂ fold-change greater than 0.25, ensuring that the displayed genes represent biologically relevant transcriptional shifts. Significant DEG counts were observed for neurons, astrocytes, and cycling astrocytes, whereas OPCs did not show any DEGs that met these criteria. Accordingly, they were analyzed similarly to the vRNA-positive cells.

Figure 5

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 log₂FC > 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 log₂FC = ±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.

The transcriptional response in bystanders varied notably according to both cell type and viral strain. Among the different cell types, neurons demonstrated a short list of DEGs, while astrocytes and especially cycling astrocytes showed a more pronounced transcriptional response. Neurons and astrocytes displayed a greater number of DEGs in response to the 23089 (B3) strain, suggesting a more pronounced impact of this viral strain on these cell populations. In contrast, cycling astrocytes exhibited more DEGs associated with the 18952 (B2) strain. However, the magnitude of transcriptional changes differed between strains: DEGs associated with 18952 (B2) tended to have more modest fold changes, whereas 23089 (B3) DEGs spanned a broader range of fold changes (Figure 5D). These patterns indicate that both the extent and intensity of the bystander transcriptional response differ depending on cell type and viral strain.

Bystander neurons displayed a comparatively limited transcriptional response, with fewer differentially expressed genes than astrocytes or cycling astrocytes. In the 18952 (B2) condition, upregulated genes included SLC5A7, essential for acetylcholine synthesis and synaptic transmission (Okuda et al., 2000), which is consistent with the cholinergic nature of spinal motor neurons. ECEL1, implicated in motor neuron axon development Kiryu-Seo et al. (2019), was also elevated. Other transcripts, such as HOXA4 and CCBE1, related to positional identity and extracellular matrix organization Bos et al. (2011), Philippidou and Dasen (2013), suggest a role in preserving neuronal identity and structure.

In contrast, 23089 (B3)-exposed bystander neurons upregulated stress-responsive and protein quality control genes. These included PTPN11 (Leahy et al., 2024), CHORDC1, and PRKCSH, associated with MAPK signaling, protein folding, and ER stress (Ferretti et al., 2010; Cressey et al., 2024). Increases in MT-ND6 and PSMC4 expression, involved in mitochondrial respiration and proteasomal degradation (Bai and Attardi, 1998; Zavodszky et al., 2021), may reflect compensatory responses to maintain homeostasis (Figure 5B). These trends suggest that while 18952 (B2) maintains neuronal identity, 23089 (B3) induces a mild stress-adaptive transcriptional state.

Astrocytes exhibited a broader transcriptional response than neurons, with distinct gene expression patterns between strains. In the 18952 (B2) condition, upregulated genes included CHST9 and PDGFA, involved in extracellular matrix remodeling and glial signaling, alongside HRH3 and SLC2A3, linked to histaminergic modulation and glucose transport. Additional genes such as BCL3, PFKFB4, and MT1X suggested mild activation of redox regulation and metabolic homeostasis. Together, these changes reflect a modestly reactive but structurally supportive astrocyte state. In contrast, astrocytes exposed to the 23089 (B3)-strain upregulated a distinct set of genes, including HSPA1A, CNTNAP2, PCDH19, and ROBO2, which are involved in synaptic regulation, axon-glia interaction, and cellular stress responses. The transcriptional profile also included ERBB4, LINGO2, and TOP2A, further implicating altered signaling and adhesion programs. These differences suggest that 23089 (B3)-exposed astrocytes enter a more transcriptionally active state, potentially affecting communication with the neuronal environment. The overall transcriptional divergence reinforces the strain-specific modulation of astrocyte identity in the bystander context (Figure 5C).

Cycling astrocytes displayed a unique transcriptional profile distinct from both neurons and non-cycling astrocytes, marked by broader identity and functional shifts between strains. In the 18952 (B2) condition, upregulated genes included SERPINB9, DPAGT1, and LOX, along with histone-associated transcripts such as H2AC13, H4C3, and HIST1H1B. These changes support a transcriptional program favoring immune regulation, extracellular matrix organization, and cell cycling. 23089 (B3)-exposed cycling astrocytes showed a marked shift in identity, with upregulation of neuronal and stress-related genes including STMN2, NEFM, DCX, and CRABP1. Additional activation of HSPA1A, JUN, and ERBB4 suggested elevated stress signaling and cytoskeletal remodeling, consistent with transcriptional dysregulation. The emergence of neuronal lineage markers in this glial progenitor population may reflect stress-induced unexpected activation of developmental genes not typically active in astrocyte progenitors (Sardi et al., 2006).

Notably, this mixed transcriptional phenotype in cells predicted to be cycling astrocytes may be linked to the fact that these cells are also the primary targets of infection by the 23089 (B3) strain. Even in bystander cells, signaling cues from nearby infected cells could perturb their transcriptional program and disrupt normal progenitor dynamics. Compared to non-cycling astrocytes, which maintained more canonical glial features, the 23089 (B3)-exposed cycling astrocytes exhibited broader deviations in transcription. These distinctions highlight the unique vulnerability of proliferative astrocyte populations to strain-specific viral influence (Figure 5D).

In OPC Oligodendrocytes, pathway enrichment analysis revealed distinct transcriptional programs between strains, even in the bystander population. Due to the relatively low number of OPC Oligo in this group, we adopted the same strategy used for infected cell analysis, performing enrichment pathways analysis to enable broader biological interpretation since no gene by itself stands out from our DEG analysis. OPCs exposed to the 18952 (B2)-strain showed upregulation of pathways associated with lipid metabolism and extracellular matrix organization, including foam cell differentiation, regulation of cholesterol efflux, and aminoglycan biosynthetic process. These processes are important for membrane dynamics and maintenance of oligodendrocyte identity, suggesting that in the 18952 (B2) condition they may retain a more metabolically stable and structurally supportive state.

In contrast, OPCs exposed to the 23089 (B3)-strain were enriched for pathways related to cell cycle regulation (mitotic spindle elongation, cytokinesis), microtubule organization, and developmental signaling, including glial cell fate specification. These transcriptional signatures are consistent with a disturbed cellular state, which could potentially reflect direct infection-related stress or dysregulation, as we observe for cycling astrocytes. This is in line with the fact that OPCs, alongside cycling astrocytes, represent a significant infected population in the 23089 (B3) condition. The enrichment of mitotic and developmental signaling programs may indicate that viral presence disrupts the normal proliferative or lineage-committed state of these glial progenitors. Together, these findings suggest that OPCs respond to 23089 (B3) exposure with transcriptional shifts that may compromise their homeostatic roles (Figure 5E).

Finally, we noted a lack of altered expression of innate immune genes, particularly interferon stimulated genes (ISGs), in our organoid model. Further analysis revealed we could detect expression of genes known to be regulated by Type I, II, and III interferon in both infected and uninfected organoids suggesting the lack of induction was not due to limited detection in our sequencing (Supplementary Figure S4). Although, our dataset captures a single time point (48 hpi), this lack of ISG induction is consistent with previous studies highlighting the effectiveness of viral innate immune suppression in EV-D68 infection, through protease cleavage of innate immune signaling factors and the action of other viral proteins (Kang et al., 2023; Li et al., 2024).

Discussion

Our analysis of 24-day old organoids showed a diverse cell composition, with both neuronal and glial lineages, consistent with previous IF marker analysis (Ogura et al., 2018). While we identified neurons, mGPCs, OPCs, VLMCs, and more, the most common cell type in the hSCO at this age was cycling astrocytes and astrocytes (54.3%). Given the relatively nascent nature of the hSCO model, this majority may be due to astrocytes’ vital role in developing and maintaining neuronal functions throughout spinal cord development (Allen and Barres, 2009; Chung et al., 2013; Allen, 2014; Allen and Eroglu, 2017). mGPCs (20.1%) also have the capacity to differentiate into both astrocytes and oligodendrocytes so cellular composition in more mature hSCO may shift (Trevino et al., 2021). Although this cell population has not yet committed to a differential pathway, they are still more specialized and committed than a broad progenitor cell. 24-day old organoids allow for the growth and development of CNS cell types without compromising on overall cell viability (Aguglia et al., 2023). The abundance of CNS cell types present and interacting in hSCO allow for us to better understand EV-D68 tropism and dynamics in the complex human spinal cord.

To understand which CNS cell-types are infected by EV-D68, we performed scRNA-seq on hSCO infected with contemporary EV-D68 strains US/IL/14-18952 (B2) and US/MA/18-23089 (B3). We found that while both strains had the capacity to infect neuronal and glial cell lineages, 18952 (B2) preferentially infected neurons while 23089 (B3) preferentially infected cycling astrocytes and to a lesser extent, oligodendrocyte precursor cells. These differences extended to bystander cells—i.e., uninfected cells within infected hSCOs—where we observed more DEGs in 18952 (B2)-bystander cycling astrocytes, and even larger fold changes in 23089 (B3)-bystander neurons and astrocytes. These findings indicate that glial cell populations such as astrocytes and oligodendrocytes play an important role in EV-D68 associated AFM pathogenesis, varying amongst different strains and clade classification. Previous studies have implicated spinal cord neurons in EV-D68-induced paralysis, but our results suggest additional cell types are infected and transcriptionally altered, potentially playing unidentified roles in EV-D68 pathogenesis (Feng et al., 2016; Hixon et al., 2017; Brown et al., 2018; Rosenfeld et al., 2019; Woods Acevedo et al., 2024).

Astrocytes participate in both innate and adaptive immune system, such as regulating the release of cytokines and chemokines, thereby leading to antigen presentation that lead to recruitment of helper T-cells in pathogen or damage-affected areas of the CNS (Dorf et al., 2000; Dong and Benveniste, 2001; McKimmie and Graham, 2010; Sofroniew and Vinters, 2010; Choi et al., 2014). Several of the differentially expressed genes identified in bystander astrocytes and cycling astrocytes, particularly in the 23089 (B3) condition, are associated with stress signaling pathways. These changes, combined with the disruption of glial identity in 23089 (B3)-exposed cycling astrocytes, may create an environment more permissive to immune cell recruitment. In contrast, the transcriptional profile in 18952 (B2)-exposed astrocytes suggested a more limited activation of broad immune-modulatory programs, although cycling astrocytes did express immune regulatory markers such as CD70 and SERPINB9, potentially reflecting a more localized or cell-type stress response rather than a widespread inflammatory activation. These findings raise the possibility that differential modulation of innate immune signaling by each strain could further shape their pathogenic outcomes in vivo. Since we do not know yet to what degree the immune response mediates EV-D68-associated AFM, further studies are needed to assess astrocyte function and subsequent immune response during EV-D68 infection. Additionally, because oligodendrocytes form myelin to aid neuron conductivity and communication, its functionality upon EV-D68 infection should also be assessed (Raine, 1984; Simons and Nave, 2015).

We also identified strain-specific differences in cellular tropism within hSCO, indicating that there has been a change in EV-D68 viral infection dynamics from 2014 to 2018. EV-D68-18952 (B2) strain primarily infected neurons, while 23089 (B3) preferentially infected cycling astrocytes and, to a lesser extent, OPCs. Together, these results suggest that the two strains may induce distinct forms of cellular vulnerability—18952 (B2) through direct neuronal targeting and 23089 (B3) through glial destabilization. While further studies will be needed to determine the long-term impact of these responses, the transcriptional divergence observed across cell types points to fundamentally different modes of pathogenesis. These differences may be broadly attributed to clade specific pathogenic differences between 18952 (B2) and 23089 (B3) or be pathogenic variations between two specific isolates. While both viruses are expected to cause the same clinical paralysis phenotype, the differences in tropism may have caused varying mechanisms of pathogenesis and we are not able to assess clinical differences between these two isolates.

Additional strain specific differences are reflected in the transcriptional impacts in bystander cells. Astrocytes and cycling astrocyte bystander cells had significantly more transcriptional changes compared to neuronal bystanders, with 23089 (B3) having more impact on astrocyte and neuron bystanders and 18952 (B2) having more of an impact on cycling astrocyte bystanders. Many of the upregulated genes in 23089 (B3) bystander cells were cellular stress indicators, primarily found in astrocytes. These cellular responses may also indicate a downstream loss in communication or homeostasis in neurons, as astrocytes are integral to maintaining neuronal function and stability. 23089 (B3) had a larger fold transcriptional shift in bystander cells compared to 18952 (B2). Although these are only two isolates, these differences may point to clade- or outbreak-specific variation in EV-D68 pathogenesis. The higher number of AFM cases in 2018 compared to 2014 could reflect such viral differences, but may also stem from multifactorial causes, including co-circulation of other enteroviruses (e.g., EV-A71) or improved case detection and reporting (Morens et al., 2019).

The OPC bystander population did not have any significant shifts, 18952 (B2) exposed OPCs were more stable and 23089 (B3) exposed OPCs had evidence of a stress response that interfered with glial cell differentiation. EV-D68 may be interfering with oligodendrocyte differentiation and maturation in the hSCO and thus prevent a fully functioning oligodendrocyte population. Considering the hSCO population of oligodendrocytes were precursor cells and still differentiating at the time of EV-D68 infection, cellular response may differ in mature oligodendrocytes. To better understand oligodendrocyte response to EV-D68 infection, studies in more aged hSCO will be necessary.

While hSCO scRNAseq has allowed us to further understand EV-D68 infection dynamics in the spinal cord, it does have its limitations. Since hSCO cells are relatively immature and represent a developing human spinal cord, EV-D68 infection dynamics may differ in hSCO to that of a child. Despite this limitation, its multicellular complexity and physiological relevance for AFM has and can let us learn more about EV-D68 pathogenesis in the CNS. This limitation in hSCO also lends itself to be an advantage when looking at enteroviruses that target neonates in order to better understand their tropism and pathogenesis.

Another limitation in our study is the relatively low number of infected cells identified for both EV-D68 strains, with vRNA-positive cells comprising only ~1.3–1.6% of the total population. Despite leveraging an aggregated, bulk analysis approach to enhance statistical power, this low infection frequency constrains the resolution of cell type-specific comparisons, particularly for rarer populations such as OPCs. As a result, strain-specific conclusions should be interpreted with caution in these cases. Moreover, the small number of infected cells may reflect an under representation of cell types that are more fragile or sensitive to dissociation, such as neurons or oligodendrocytes, further limiting subtype-level analysis.

Despite these limitations, our study provides valuable insights into the cellular and molecular landscape of EV-D68 infection in human spinal cord organoids, revealing distinct cell-type tropism for two contemporary strains. By integrating pathway enrichment analysis with single-cell transcriptomic profiling, we demonstrate that EV-D68-18952 (B2) exhibits a pronounced preference for neurons, driving synaptic signaling and axonal development pathways, while EV-D68-23089 (B3) predominantly targets cycling astrocytes, triggering transcriptional programs associated with cell cycle progression and RNA processing. These findings represent a step forward in understanding the strain-specific interactions of EV-D68 with neural populations, which may have implications for viral spread and neuro pathogenesis. While spinal cord organoids represent an early developmental CNS state, their cellular diversity and structural organization enable the study of viral tropism and host responses relevant to pediatric neurotropic infections. Furthermore, as the brain and spinal cord are comprised of unique cells which may respond to infection in specific manners, our use of spinal cord organoids rather than cerebroids as a model provides a physiologically relevant system for a paralytogenic virus to dissect host–virus interactions at single-cell resolution, underscoring the utility of this platform for studying neurotropic viruses.

These findings reveal distinct transcriptional programs associated with infection, but they are not yet supported by phenotypic validation. Future studies will be necessary to assess whether these transcriptional shifts correspond to functional changes in viral replication, cell viability, or neural function.

Our analysis focused on an early time point (48 hpi) to capture initial transcriptional responses to viral entry and replication. However, it is possible that key host responses—particularly delayed innate immune signaling or neuronal damage—may emerge at later stages. Follow-up studies incorporating extended time points will be important to assess the full temporal dynamics of EV-D68 in the CNS.

Although we did not perform in vivo validation in this study, our single-cell transcriptomic findings align with emerging clinical and experimental data on AFM pathogenesis. Notably, the neuron-enriched infection profile of EV-D68-18952 (B2) mirrors reports of motor neuron involvement in mouse models and human cases. In particular, a case report by Vogt et al. (2022) identified EV-D68 in anterior horn cells of a child with AFM but found only sparse detection of viral RNA and protein, suggesting that even low-level infection may contribute to significant neurological damage (Vogt et al., 2022). This aligns with prior mouse model studies demonstrating neurotropic infection and paralysis upon EV-D68 inoculation (Hixon et al., 2017; Brown et al., 2018; Rosenfeld et al., 2019), as well as findings that immune cell infiltration can drive neuropathology even when infection levels are low (Woods Acevedo et al., 2024). In human AFM cases, direct detection of EV-D68 in CSF or spinal cord tissue remains rare, yet serological studies provide strong evidence of prior CNS exposure, with enterovirus-specific antibodies detected in the cerebrospinal fluid (Mishra et al., 2019; Schubert et al., 2019). These observations reinforce the idea that rare but targeted infections events—especially in vulnerable neuronal populations—can have disproportionate pathological effects. Thus, our hSCO model, despite low infection frequencies, captures a clinically relevant aspect of EV-D68 biology. This highlights the value of single-cell approaches in modeling strain-specific tropism and generating testable hypotheses for downstream validation in more complex in vivo systems.

Materials and methods

Viruses and cells

EV-D68 US/IL/14-18952 (CDC) and US/MA/18-23089 (B3) (CDC) strains were propagated using HeLa cells incubated at 33 °C and 5% CO₂ and purified using sucrose-cushion as previously described (Morosky et al., 2016). These stocks were previously sequenced to confirm their identity with VP1 primers.

HeLa 7b (ATCC, CCl-2) cells were maintained in MEM medium (Thermo Fisher, 11095-072), supplemented to contain 5% FBS (Phenonix Scientific, PS-100), 1% penicillin/streptomycin (Corning, 30-002-Cl), and 1% NeAA (Corning, 25-025-Cl). Cells were grown at 37 °C and 5% CO2.

Human iPSC line SCTi003A (STEMCELL Technologies, 200-0511) was maintained in mTeSR™ Plus medium (STEMCELL Technologies, 100-0276), supplemented with 10 μM Y-27632 (Tocris, 1254). They were seeded and passaged in flasks coated with 150 μg/mL Cultrex (R&D Systems, 3434-005-02). The 3-DiSC hSCO were propagated and differentiated as described in Aguglia et al. (2023) for up to 24 days.

Viral titration by plaque assay

Plaque assays were conducted using HeLa cells overplayed with 1% agarose, incubated at 33 °C for 72 h (EV-D68) and 5% CO₂ and stained with crystal violet prior to enumeration.

hSCO infections and dissociations

hSCOs were infected in pools of 12 organoids each and inoculated with virus at 105 PFU/pool. After 1 h of incubation at room temperature, hSCOs were washed 3X with PBS and moved to new wells before incubation with fresh medium. No further media changes were performed for the rest of the experiment. The EV-D68 infected pool was incubated at 33 °C for 48 h post infection (hpi). After 48 hpi, the pools were dissociated to a single-cell suspension with Accumax (Invitrogen), incubating for 15 min in a water bath at 37 °C, followed by gentle mixing to dissociate. We then proceeded with the proposed 10X Genomics protocol (Cell Preparation Guide—CG00053 Rev. C). Cell viability was determined by Countess II (Invitrogen) with trypan blue staining. Viability for all samples was between 60% and 70%. Once the viable cell number was determined; we diluted cells to the recommended concentration before incubating on ice.

Single cell RNAseq cDNA library generation

All samples were calculated to achieve ~5,000–8,000 targeting cells in the single-cell preparations using Chromium Next GEM Single Cell 5′ standard kit. For the preparation of the cDNA and sequencing library generation, we followed the instructions from the user guide Chromium Next GEM Single Cell 5′ Reagent Kit v2 (Dual index). All other steps were followed to produce cDNA and subsequent Illumina sequencing library for single cell sequencing. The illumina library preparation was submitted to quality control in the Tape Station D1000 high sensitivity for size distribution and DNA concentration was measured by Qubit High Sensitivity dsDNA kit. The molar concentration of the libraries were determined and the samples were diluted for sequencing according to illumina sequencing protocol. We aimed to sequence each library to achieve ~50,000 reads per cell.

Cell ranger

Cell Ranger v7.0.0 was used (Zheng et al., 2017) to align reads against a composite genome, which encompassed both the GRCh38 human reference genome and the Enterovirus genome corresponding to the specific strain identity of the sample. Feature-barcode matrices were generated using GENCODE v44 GRCh38 gene models and the viral strain’s genome as a single ORF. Default parameters of the “cell ranger count” function were used.

Seurat

Seurat v5.0.2 was used (Hao et al., 2021). Cell Ranger gene counts were made compatible across different strains by setting the name of the gene encompassing all viral reads for each sample to “Viral-Gene”. Cells with mitochondrial gene expression >25% or detected gene counts <1,000 were removed. Gene counts of filtered cells were normalized using SC Transform, while specifying the “vars to regress” parameter to the “S. Scores” and “G2M. Scores” obtained from the “Cell Cycle Scoring” function, utilizing Seurat’s “cc. genes” cell cycle gene list. Reciprocal PCA integration analyses were then performed to generate integrated datasets for EV-D68. Dimensionality reduction of the integrated dataset was performed using the “Run PCA” function and “Run UMAP” function using the top 30 principal components. Cell cluster analysis was performed using the “Find Neighbors” function using the top 30 principal components and the “Find Clusters” function using a “resolution” of 0.1. Cell type annotations were predicted using cell label transfer with Seurat’s Find Transfer Anchors, Transfer Data, and Add Meta Data functions using a previously annotated spinal cord scRN Aseq dataset (Andersen et al., 2023) as our reference, considering “true” cell types above 0.5 prediction score max. The cells that were below this threshold were assigned as NA or Unclassified.

Assignment of infection status

Viral read percentages were calculated as the fraction of total UMIs per cell mapping to the viral gene (i.e., “Viral-Gene”) and were computed using the “Percentage Feature Set” Seurat function. To classify cells as “Infected” or “Non Infected”, a Poisson test was applied to these percentages using the “Estimate Non Expressing Cells” function from the Soup X R package (Young and Behjati, 2020), using an FDR of 0.05. This function accounts for ambient RNA contamination in the sample, which was estimated using Soup X’s “Auto Est Count” function with “tfidfMin” set to 1.0 and “Soup Quantile” set to 0.9.

Permutation-based enrichment analysis

To assess whether specific cell types were significantly enriched in infected (vRNA+ cells) versus non-infected (Bystander cells) conditions, we performed a permutation test on cell count distributions across infection states. For each viral strain, we computed contingency tables comparing observed cell type frequencies across infection status (FDR < 0.05). To establish a null distribution, we randomly permuted cell type labels 10,000 times while preserving the infection status labels and recomputed the contingency tables for each iteration. Median values and 95% confidence intervals (2.5th and 97.5th percentiles) were derived from the permuted distributions. Observed counts exceeding the upper confidence interval were considered significantly enriched (p value <0.05). Enrichment scores were calculated as the ratio of observed to permuted median counts for each cell type and condition. Analyses were conducted in R using the data, Seurat, and Seurat Object packages.

Differential gene expression analysis

Differential gene expression (DEG) analysis was conducted using the Find Markers function in Seurat v5. Two complementary approaches were applied. First, a pseudo-bulk-style analysis was performed by grouping cells according to their combined infection status and viral strain (e.g., “18952_Infected” vs. “18952_NotInfected”), enabling comparisons across aggregated conditions. This approach was used because the number of infected cells was insufficient for robust stratification by cell type. Second, a cell-type-specific DEG analysis was performed using only bystander cells, focusing on neurons, astrocytes, cycling astrocytes, and oligodendrocyte lineage cells (OPC Oligo.) to compare transcriptional responses between EV-D68 strains 18952 and 23089. All analyses were conducted using the RNA assay, with normalization via Normalize Data. DEG identification used a minimum expression threshold of 10% (min.pct = 0.1), no log fold-change cutoff (logfc. threshold = 0), and a minimum of 50 cells per group. Fold changes were calculated as log₂-transformed values, and significance was assessed using Bonferroni-adjusted p-values. Genes with adjusted p value <0.05 and | log₂FC | > 0.25 were considered significantly differentially expressed. Visualizations, including volcano plots and DEG count bar plots by cell type and direction, were generated using ggplot2, ggrepel, and patchwork.

Pathway enrichment analysis

Gene set enrichment analysis (GSEA) was performed using the fgsea package to identify biological pathways enriched in infected conditions or specific cell types, even in the absence of significantly differentially expressed genes. For certain DEG comparisons, particularly those involving low-abundance populations (e.g., infected cells or bystander OPC Oligo.), no genes met the adjusted p-value threshold for significance. In these cases, the full ranked DEG lists based on average log₂ fold change were used as input for GSEA, enabling the detection of coordinated pathway-level shifts. Ranked gene lists were generated from pairwise comparisons of infection conditions (e.g., 18952-Infected vs. Mock) and bystander cell types (e.g., 23089 vs. 18952 within OPC Oligo.), and enrichment was computed using 100 million permutations for robust estimation. Gene sets were sourced from the Gene Ontology Biological Process (GO-BP) category via msigdbr. Normalized enrichment scores (NES) were computed for each condition, and pathways with adjusted p-values <0.05 were considered significantly enriched. Heatmaps of NES values across conditions were generated using ComplexHeatmap, and selected pathway comparisons were visualized with ggplot2 and ggrepel.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://www.ncbi.nlm.nih.gov/geo/, GSE292051.

Ethics statement

Ethical approval was not required for the studies on humans in accordance with the local legislation and institutional requirements because only commercially available established cell lines were used.

Author contributions

ND: Writing – review & editing, Methodology, Investigation, Writing – original draft. SM: Writing – review & editing, Writing – original draft. CM: Formal analysis, Writing – review & editing, Data curation. TE: Resources, Writing – review & editing, Formal analysis, Data curation. PD: Investigation, Writing – review & editing, Methodology, Formal analysis, Supervision, Project administration, Funding acquisition, Conceptualization, Visualization, Writing – original draft. MF: Investigation, Supervision, Writing – review & editing, Funding acquisition, Writing – original draft, Resources, Methodology, Project administration, Conceptualization.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This work was supported by the Intramural Research Program of the National Institute of Allergy and Infectious Diseases at the National Institutes of Health—Project Number 1ZIAAI001360 (PTD). MCF receives support from National Institutes of Health K08AI171177. This research was supported in part by the Intramural Research Program of the National Institutes of Health (NIH). The contributions of the NIH author(s) were made as part of their official duties as NIH federal employees, are in compliance with agency policy requirements, and are considered Works of the United States Government. However, the findings and conclusions presented in this paper are those of the authors and do not necessarily reflect the views of the NIH or the U.S. Department of Health and Human Services.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Gen AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2025.1698639/full#supplementary-material

SUPPLEMENTARY FIGURE S1 | UMAP ‘FeaturePlot’ showing Log-normalized expression of key marker genes of astrocytes (SOX9 and FGFR3), cycling astrocytes (TOP2A) and neurons (MAP2 and TUBB3).

SUPPLEMENTARY FIGURE S2 | Cell type characterization in hSCOs. (A) UMAP facets of all cell types in 24-day old hSCO. (B) Donut plots showing the relative frequency of cell types in hSCOs infected with each EV-D68 strain and mock-infected controls.

SUPPLEMENTARY FIGURE S3 | (A) UMAP showing the cellular transcriptional phenotypes of mGPC clusters in 24-days old hSCO. (B) Heatmap of top 20 marker genes that distinguish each cluster. The heatmap displays scaled expression values (z-scores) of the top 20 marker genes per cluster, calculated using Seurat’s ScaleData function. Values are centered and scaled per gene across all cells, such that 0 represents the mean expression and ±2 corresponds to approximately two standard deviations above or below the mean. This highlights relative over- or under-expression patterns across clusters.

SUPPLEMENTARY FIGURE S4 | Scatter plots comparing of the proportion of cells expressing ISGs induced by Type I, II, and III interferons, and Volcano plots comparing the fold change in expression and the significance, as −log₁₀(adjusted p-value), for comparisons of (A) Mock vs. 18952 Bystander cells, (B) Mock vs. 20892 Bystander cells, (C) 18952 Bystander cells vs. 20892 Bystander cells.

References

Aguglia, G., Coyne, C. B., Dermody, T. S., Williams, J. V., and Freeman, M. C. (2023). Contemporary enterovirus-D68 isolates infect human spinal cord organoids. MBio 14:e0105823. doi: 10.1128/mbio.01058-23

PubMed Abstract | Crossref Full Text | Google Scholar

Aliabadi, N., Messacar, K., Pastula, D. M., Robinson, C. C., Leshem, E., Sejvar, J. J., et al. (2016). Enterovirus D68 infection in children with acute flaccid myelitis, Colorado, USA, 2014. Emerg. Infect. Dis. 22, 1387–1394. doi: 10.3201/eid2208.151949

PubMed Abstract | Crossref Full Text | Google Scholar

Allen, N. J. (2014). Astrocyte regulation of synaptic behavior. Annu. Rev. Cell Dev. Biol. 30, 439–463. doi: 10.1146/annurev-cellbio-100913-013053

PubMed Abstract | Crossref Full Text | Google Scholar

Allen, N. J., and Barres, B. A. (2009). Neuroscience: glia—more than just brain glue. Nature 457, 675–677. doi: 10.1038/457675a

PubMed Abstract | Crossref Full Text | Google Scholar

Allen, N. J., and Eroglu, C. (2017). Cell biology of astrocyte-synapse interactions. Neuron 96, 697–708. doi: 10.1016/j.neuron.2017.09.056

PubMed Abstract | Crossref Full Text | Google Scholar

Andersen, J., Thom, N., Shadrach, J. L., Chen, X., Onesto, M. M., Amin, N. D., et al. (2023). Single-cell transcriptomic landscape of the developing human spinal cord. Nat. Neurosci. 26, 902–914. doi: 10.1038/s41593-023-01311-w

PubMed Abstract | Crossref Full Text | Google Scholar

Bai, Y., and Attardi, G. (1998). The mtDNA-encoded ND6 subunit of mitochondrial NADH dehydrogenase is essential for the assembly of the membrane arm and the respiratory function of the enzyme. EMBO J. 17, 4848–4858. doi: 10.1093/emboj/17.16.4848

PubMed Abstract | Crossref Full Text | Google Scholar

Bos, F. L., Caunt, M., Peterson-Maduro, J., Planas-Paz, L., Kowalski, J., Karpanen, T., et al. (2011). CCBE1 is essential for mammalian lymphatic vascular development and enhances the lymphangiogenic effect of vascular endothelial growth factor-C in vivo. Circ. Res. 109, 486–491. doi: 10.1161/CIRCRESAHA.111.250738

PubMed Abstract | Crossref Full Text | Google Scholar

Brown, D. M., Hixon, A. M., Oldfield, L. M., Zhang, Y., Novotny, M., Wang, W., et al. (2018). Contemporary circulating enterovirus D68 strains have acquired the capacity for viral entry and replication in human neuronal cells. mBio 9, e01954–e01918. doi: 10.1128/mBio.01954-18

PubMed Abstract | Crossref Full Text | Google Scholar

Brown, B. A., Nix, W. A., Sheth, M., Frace, M., and Oberste, M. S. (2014). Seven strains of enterovirus D68 detected in the United States during the 2014 severe respiratory disease outbreak. Genome Announc. 2:e01201. doi: 10.1128/genomeA.01201-14

PubMed Abstract | Crossref Full Text | Google Scholar

CDC (2025). AFM cases and outbreaks. Acute flaccid myelitis (AFM). Available online at: https://www.cdc.gov/acute-flaccid-myelitis/cases/index.html (Accessed March 5, 2025).

Google Scholar

Chatterjee, M., Schild, D., and Teunissen, C. E. (2019). Contacting in the central nervous system: role in health and disease. Neural Regen. Res. 14, 206–216. doi: 10.4103/1673-5374.244776

PubMed Abstract | Crossref Full Text | Google Scholar

Choi, S. S., Lee, H. J., Lim, I., Satoh, J.-I., and Kim, S. U. (2014). Human astrocytes: secretive profiles of cytokines and chemokines. PLoS One 9:e92325. doi: 10.1371/journal.pone.0092325

PubMed Abstract | Crossref Full Text | Google Scholar

Chung, W.-S., Clarke, L. E., Wang, G. X., Stafford, B. K., Sher, A., Chakraborty, C., et al. (2013). Astrocytes mediate synapse elimination through MEGF10 and MERTK pathways. Nature 504, 394–400. doi: 10.1038/nature12776

PubMed Abstract | Crossref Full Text | Google Scholar

Clopper, B. R., Lopez, A. S., Goldstein, L. A., Ng, T. F. F., Toepfer, A. P., Staat, M. A., et al. (2025). Enterovirus D68-associated respiratory illness in children. JAMA Netw. Open 8:e259131. doi: 10.1001/jamanetworkopen.2025.9131

PubMed Abstract | Crossref Full Text | Google Scholar

Cressey, R., Han, M. T. T., Khaodee, W., Xiyuan, G., and Qing, Y. (2024). Navigating PRKCSH’S impact on cancer: from N-linked glycosylation to death pathway and anti-tumor immunity. Front. Oncol. 14:1378694. doi: 10.3389/fonc.2024.1378694

PubMed Abstract | Crossref Full Text | Google Scholar

Dehmelt, L., and Halpain, S. (2005). The MAP2/tau family of microtubule-associated proteins. Genome Biol. 6:204. doi: 10.1186/gb-2004-6-1-204

PubMed Abstract | Crossref Full Text | Google Scholar

Dong, Y., and Benveniste, E. N. (2001). Immune function of astrocytes. Glia 36, 180–190. doi: 10.1002/glia.1107

PubMed Abstract | Crossref Full Text | Google Scholar

Dorf, M. E., Berman, M. A., Tanabe, S., Heesen, M., and Luo, Y. (2000). Astrocytes express functional chemokine receptors. J. Neuro 111, 109–121. doi: 10.1016/s0165-5728(00)00371-4

PubMed Abstract | Crossref Full Text | Google Scholar

Erlander, M. G., and Tobin, A. J. (1991). The structural and functional heterogeneity of glutamic acid decarboxylase: a review. Neurochem. Res. 16, 215–226. doi: 10.1007/BF00966084

PubMed Abstract | Crossref Full Text | Google Scholar

Feng, M., Guo, S., Fan, S., Zeng, X., Zhang, Y., Liao, Y., et al. (2016). The preferential infection of astrocytes by enterovirus 71 plays a key role in the viral neurogenic pathogenesis. Front. Cell. Infect. Microbiol. 6:192. doi: 10.3389/fcimb.2016.00192

PubMed Abstract | Crossref Full Text | Google Scholar

Ferretti, R., Palumbo, V., Di Savino, A., Velasco, S., Sbroggiò, M., Sportoletti, P., et al. (2010). Morgana/chp-1, a ROCK inhibitor involved in centrosome duplication and tumorigenesis. Dev. Cell 18, 486–495. doi: 10.1016/j.devcel.2009.12.020

Crossref Full Text | Google Scholar

Freeman, M. C., and Messacar, K. (2025). Enterovirus and Par echovirus neurologic infections in children: clinical presentations and neuropath genesis. J. Pediatric Infect. Dis. Soc. 14:piae069. doi: 10.1093/jpids/piae069

PubMed Abstract | Crossref Full Text | Google Scholar

Gleeson, J. G., Lin, P. T., Flanagan, L. A., and Walsh, C. A. (1999). Doublecortin is a microtubule-associated protein and is expressed widely by migrating neurons. Neuron 23, 257–271. doi: 10.1016/S0896-6273(00)80778-3

PubMed Abstract | Crossref Full Text | Google Scholar

Gomes, I., Bobeck, E. N., Margolis, E. B., Gupta, A., Sierra, S., Fakira, A. K., et al. (2016). Identification of GPR83 as the receptor for the neuroendocrine peptide PEN. Sci. Signal. 9:ra43. doi: 10.1126/scisignal.aad0694

PubMed Abstract | Crossref Full Text | Google Scholar

Hadfield, J., Megill, C., Bell, S. M., Huddleston, J., Potter, B., Callender, C., et al. (2018). Next strain: real-time tracking of pathogen evolution. Bioinformatics 34, 4121–4123. doi: 10.1093/bioinformatics/bty407

PubMed Abstract | Crossref Full Text | Google Scholar

Hao, Y., Hao, S., Andersen-Nissen, E., Mauck, W. M. 3rd, Zheng, S., Butler, A., et al. (2021). Integrated analysis of multimodal single-cell data. Cell 184, 3573–3587.e29. doi: 10.1016/j.cell.2021.04.048

PubMed Abstract | Crossref Full Text | Google Scholar

Hixon, A. M., Yu, G., Leser, J. S., Yagi, S., Clarke, P., Chiu, C. Y., et al. (2017). A mouse model of paralytic myelitis caused by enterovirus D68. PLoS Pathog. 13:e1006199. doi: 10.1371/journal.ppat.1006199

PubMed Abstract | Crossref Full Text | Google Scholar

Kang, J., Huang, M., Li, J., Zhang, K., Zhu, C., Liu, S., et al. (2023). Enterovirus D68 VP3 targets the interferon regulatory factor 7 to inhibit type I interferon response. Micro. Spectr. 11:e0413822. doi: 10.1128/spectrum.04138-22

PubMed Abstract | Crossref Full Text | Google Scholar

Kiryu-Seo, S., Nagata, K., Saido, T. C., and Kiyama, H. (2019). New insights of a neuronal peptidase DINE/ECEL1: nerve development, nerve regeneration and neurogenic pathogenesis. Neurochem. Res. 44, 1279–1288. doi: 10.1007/s11064-018-2665-x

PubMed Abstract | Crossref Full Text | Google Scholar

Leahy, S. N., Vita, D. J., and Broadie, K. (2024). PTPN11/corkscrew activates local presynaptic mapk signaling to regulate synapsing, synaptic vesicle pools, and neurotransmission strength, with a dual requirement in neurons and glia. J. Neurosci. 44:e1077232024. doi: 10.1523/JNEUROSCI.1077-23.2024

PubMed Abstract | Crossref Full Text | Google Scholar

Lee, J. H., and Berger, J. M. (2019). Cell cycle-dependent control and roles of DNA topoisomerase II. Genes Basel 10:859. doi: 10.3390/genes10110859

PubMed Abstract | Crossref Full Text | Google Scholar

Li, X., Guo, H., Yang, J., Liu, X., Li, H., Yang, W., et al. (2024). Enterovirus D68 3C protease antagonizes type I interferon signaling by cleaving signal transducer and activator of transcription 1. J. Virol. 98:e0199423. doi: 10.1128/jvi.01994-23

PubMed Abstract | Crossref Full Text | Google Scholar

Liu, X., Li, H., Li, Z., Gao, D., Zhou, J., Ni, F., et al. (2025). MFSD6 is an entry receptor for respiratory enterovirus D68. Cell Host Microbe 33, 267–278.e4. doi: 10.1016/j.chom.2024.12.015

PubMed Abstract | Crossref Full Text | Google Scholar

Ma, Q., Fode, C., Guillemot, F., and Anderson, D. J. (1999). Neurogenin1 and neurogenin2 control two distinct waves of neurogenesis in developing dorsal root ganglia. Genes Dev. 13, 1717–1728. doi: 10.1101/gad.13.13.1717

PubMed Abstract | Crossref Full Text | Google Scholar

McIntire, S. L., Reimer, R. J., Schuske, K., Edwards, R. H., and Jorgensen, E. M. (1997). Identification and characterization of the vesicular GABA transporter. Nature 389, 870–876. doi: 10.1038/39908

PubMed Abstract | Crossref Full Text | Google Scholar

McKay, S. L., Lee, A. D., Lopez, A. S., Nix, W. A., Dooling, K. L., Keaton, A. A., et al. (2018). Increase in acute flaccid myelitis - United States, 2018. MMWR Morb. Mortal Wkly. Rep. 67, 1273–1275. doi: 10.15585/mmwr.mm6745e1

PubMed Abstract | Crossref Full Text | Google Scholar

McKimmie, C. S., and Graham, G. J. (2010). Astrocytes modulate the chemokine network in a pathogen-specific manner. Biochem. Biophys. Res. Commun. 394, 1006–1011. doi: 10.1016/j.bbrc.2010.03.111

PubMed Abstract | Crossref Full Text | Google Scholar

Menezes, J. R., and Luskin, M. B. (1994). Expression of neuron-specific tubulin defines a novel population in the proliferative layers of the developing telencephalon. J. Neurosci. 14, 5399–5416. doi: 10.1523/JNEUROSCI.14-09-05399.1994

PubMed Abstract | Crossref Full Text | Google Scholar

Messacar, K., Abzug, M. J., and Dominguez, S. R. (2016). 2014 outbreak of enterovirus D68 in North America. J. Med. Virol. 88, 739–745. doi: 10.1002/jmv.24410

PubMed Abstract | Crossref Full Text | Google Scholar

Messacar, K., Asturias, E. J., Hixon, A. M., Van Leer-Buter, C., Niesters, H. G. M., Tyler, K. L., et al. (2018). Enterovirus D68 and acute flaccid myelitis—evaluating the evidence for causality. Lancet Infect. Dis. 18, e239–e247. doi: 10.1016/S1473-3099(18)30094-X

PubMed Abstract | Crossref Full Text | Google Scholar

Mishra, N., Ng, T. F. F., Marine, R. L., Jain, K., Ng, J., Thakkar, R., et al. (2019). Antibodies to enteroviruses in cerebrospinal fluid of patients with acute flaccid myelitis. mBio 10, e01903–e01919. doi: 10.1128/mBio.01903-19

PubMed Abstract | Crossref Full Text | Google Scholar

Miyata, T., Maeda, T., and Lee, J. E. (1999). NeuroD is required for differentiation of the granule cells in the cerebellum and hippocampus. Genes Dev. 13, 1647–1652. doi: 10.1101/gad.13.13.1647

PubMed Abstract | Crossref Full Text | Google Scholar

Morens, D. M., Folkers, G. K., and Fauci, A. S. (2019). Acute flaccid myelitis: something old and something new. mBio 10, e00521–e00519. doi: 10.1128/mBio.00521-19

PubMed Abstract | Crossref Full Text | Google Scholar

Morosky, S., Lennemann, N. J., and Coyne, C. B. (2016). BPIFB6 regulates secretory pathway trafficking and enterovirus replication. J. Virol. 90, 5098–5107. doi: 10.1128/JVI.00170-16

PubMed Abstract | Crossref Full Text | Google Scholar

Müller, T. D., Müller, A., Yi, C.-X., Habegger, K. M., Meyer, C. W., Gaylinn, B. D., et al. (2013). The orphan receptor Gpr83 regulates systemic energy metabolism via ghrelin-dependent and ghrelin-independent mechanisms. Nat. Commun. 4:1968. doi: 10.1038/ncomms2968

PubMed Abstract | Crossref Full Text | Google Scholar

Mulligan, M. R., and Bicknell, L. S. (2023). The molecular genetics of naval in brain development and disease. Eur. J. Hum. Genet. 31, 1209–1217. doi: 10.1038/s41431-023-01456-z

PubMed Abstract | Crossref Full Text | Google Scholar

Murphy, O. C., Messacar, K., Benson, L., Bove, R., Carpenter, J. L., and Crawford, T. (2021). Acute flaccid myelitis: cause, diagnosis, and management. Lancet 397, 334–346. doi: 10.1016/S0140-6736(20)32723-9

PubMed Abstract | Crossref Full Text | Google Scholar

Nardelli, J., Thiesson, D., Fujiwara, Y., Tsai, F. Y., and Orkin, S. H. (1999). Expression and genetic interaction of transcription factors GATA-2 and GATA-3 during development of the mouse central nervous system. Dev. Biol. 210, 305–321. doi: 10.1006/dbio.1999.9278

PubMed Abstract | Crossref Full Text | Google Scholar

New Vaccine Surveillance Network CollaboratorsHall, A. J. (2022). Increase in acute respiratory illnesses among children and adolescents associated with rhinoviruses and enteroviruses

Google Scholar

Ogura, T., Sakaguchi, H., Miyamoto, S., and Takahashi, J. (2018). Three-dimensional induction of dorsal, intermediate and ventral spinal cord tissues from human pluripotent stem cells. Development 145:dev162214. doi: 10.1242/dev.162214

PubMed Abstract | Crossref Full Text | Google Scholar

Okano, H. J., and Darnell, R. B. (1997). A hierarchy of Hu RNA binding proteins in developing and adult neurons. J. Neurosci. 17, 3024–3037. doi: 10.1523/JNEUROSCI.17-09-03024.1997

PubMed Abstract | Crossref Full Text | Google Scholar

Okuda, T., Haga, T., Kanai, Y., Endou, H., Ishihara, T., and Katsura, I. (2000). Identification and characterization of the high-affinity choline transporter. Nat. Neurosci. 3, 120–125. doi: 10.1038/72059

PubMed Abstract | Crossref Full Text | Google Scholar

Olsen, S. J., Winn, A. K., Budd, A. P., Prill, M. M., Steel, J., Midgley, C. M., et al. (2021). Changes in influenza and other respiratory virus activity during the COVID-19 pandemic—United States, 2020-2021. MMWR Morb. Mortal Wkly. Rep. 70, 1013–1019. doi: 10.15585/mmwr.mm7029a1

PubMed Abstract | Crossref Full Text | Google Scholar

Philippidou, P., and Dasen, J. S. (2013). Hox genes: choreographers in neural development, architects of circuit organization. Neuron 80, 12–34. doi: 10.1016/j.neuron.2013.09.020

PubMed Abstract | Crossref Full Text | Google Scholar

Pringle, N. P., Yu, W.-P., Howell, M., Colvin, J. S., Ornitz, D. M., and Richardson, W. D. (2003). Fgfr3 expression by astrocytes and their precursors: evidence that astrocytes and oligodendrocytes originate in distinct neuro epithelial domains. Development 130, 93–102. doi: 10.1242/dev.00184

PubMed Abstract | Crossref Full Text | Google Scholar

Raine, C. S. (1984). “Morphology of myelin and myelination” in Myelin (Boston, MA: Springer US), 1–50.

Google Scholar

Rosenfeld, A. B., Warren, A. L., and Racaniello, V. R. (2019). Neuro tropism of enterovirus D68 isolates is independent of sialic acid and is not a recently acquired phenotype. mBio 10, e02370–e02319. doi: 10.1128/mBio.02370-19

PubMed Abstract | Crossref Full Text | Google Scholar

Sagulenko, P., Puller, V., and Neher, R. A. (2018). Tree Time: maximum-likelihood phylodynamic analysis. Virus Evol. 4:vex042. doi: 10.1093/ve/vex042

PubMed Abstract | Crossref Full Text | Google Scholar

Sardi, S. P., Murtie, J., Koirala, S., Patten, B. A., and Corfas, G. (2006). Presenilin-dependent ErbB4 nuclear signaling regulates the timing of astrogenesis in the developing brain. Cell 127, 185–197. doi: 10.1016/j.cell.2006.07.037

PubMed Abstract | Crossref Full Text | Google Scholar

Schubert, F. R., Fainsod, A., Gruenbaum, Y., and Gruss, P. (1995). Expression of the novel murine home box gene Sax-1 in the developing nervous system. Mech. Dev. 51, 99–114. doi: 10.1016/0925-4773(95)00358-8

PubMed Abstract | Crossref Full Text | Google Scholar

Schubert, R. D., Hawes, I. A., Ramachandran, P. S., Ramesh, A., Crawford, E. D., and Pak, J. E. (2019). Pan-viral serology implicates enteroviruses in acute flaccid myelitis. Nat. Med. 25, 1748–1752. doi: 10.1038/s41591-019-0613-1

PubMed Abstract | Crossref Full Text | Google Scholar

Shah, M. M., Perez, A., Lively, J. Y., Avadhanula, V., Boom, J. A., and Chappell, J. (2021). Enterovirus D68-associated acute respiratory illness─new vaccine surveillance network, United States, July-November 2018-2020. MMWR Morb. Mortal Wkly. Rep. 70, 1623–1628. doi: 10.15585/mmwr.mm7047a1

PubMed Abstract | Crossref Full Text | Google Scholar

Simons, M., and Nave, K.-A. (2015). Oligodendrocytes: myelination and axonal support. Cold Spring Harb. Perspect. Biol. 8:a020479. doi: 10.1101/cshperspect.a020479

PubMed Abstract | Crossref Full Text | Google Scholar

Sofroniew, M. V., and Vinters, H. V. (2010). Astrocytes: biology and pathology. Acta Neuropathol. 119, 7–35. doi: 10.1007/s00401-009-0619-8

PubMed Abstract | Crossref Full Text | Google Scholar

Stolt, C. C., Lommes, P., Sock, E., Chaboissier, M.-C., Schedl, A., and Wegner, M. (2003). The Sox9 transcription factor determines glial fate choice in the developing spinal cord. Genes Dev. 17, 1677–1689. doi: 10.1101/gad.259003

PubMed Abstract | Crossref Full Text | Google Scholar

Stuart, T., Butler, A., Hoffman, P., Hafemeister, C., Papalexi, E., Mauck, W. M. 3rd, et al. (2019). Comprehensive integration of single-cell data. Cell 177, 1888–1902.e21. doi: 10.1016/j.cell.2019.05.031

Crossref Full Text | Google Scholar

Teoh, Z., Midgley, C., Quigley, C. M., Conrey, S. C., Cline, A. R., Ostrow, A., et al. (2025). P-1212. Enterovirus D68 (EV-D68) incidence and spectrum of illness among children 0-2 years in a U.S. community-based birth cohort in Cincinnati, OH during 2017-2020. Open Forum Infect. Dis. 12:1394. doi: 10.1093/jpids/piaf057

Crossref Full Text | Google Scholar

Trevino, A. E., Müller, F., Andersen, J., Sundaram, L., Kathiria, A., Shcherbina, A., et al. (2021). Chromatin and gene-regulatory dynamics of the developing human cerebral cortex at single-cell resolution. Cell 184, 5053–5069.e23. doi: 10.1016/j.cell.2021.07.039

PubMed Abstract | Crossref Full Text | Google Scholar

Vogt, M. R., Wright, P. F., Hickey, W. F., De Buysscher, T., Boyd, K. L., and Crowe, J. E. Jr. (2022). Enterovirus D68 in the anterior horn cells of a child with acute flaccid myelitis. N. Engl. J. Med. 386, 2059–2060. doi: 10.1056/NEJMc2118155

PubMed Abstract | Crossref Full Text | Google Scholar

Wang, Y.-F., Chou, C.-T., Lei, H.-Y., Liu, C.-C., Wang, S.-M., Yan, J.-J., et al. (2004). A mouse-adapted enterovirus 71 strain causes neurological disease in mice after oral infection. J. Virol. 78, 7916–7924. doi: 10.1128/JVI.78.15.7916-7924.2004

PubMed Abstract | Crossref Full Text | Google Scholar

Woods Acevedo, M. A., Lan, J., Maya, S., Jones, J. E., Williams, J. V., Freeman, M. C., et al. (2024). Immune cells promote paralytic disease in mice infected with enterovirus D68. bioRxiv :10.14.618341. doi: 10.1101/2024.10.14.618341

PubMed Abstract | Crossref Full Text | Google Scholar

Young, M. D., and Behjati, S. (2020). SoupX removes ambient RNA contamination from droplet-based single-cell RNA sequencing data. Gigascience 9:giaa151. doi: 10.1093/gigascience/giaa151

PubMed Abstract | Crossref Full Text | Google Scholar

Zavodszky, E., Peak-Chew, S.-Y., Juszkiewicz, S., Narvaez, A. J., and Hegde, R. S. (2021). Identification of a quality-control factor that monitors failures during proteasome assembly. Science 373, 998–1004. doi: 10.1126/science.abc6500

PubMed Abstract | Crossref Full Text | Google Scholar

Zheng, G. X. Y., Terry, J. M., Belgrader, P., Ryvkin, P., Bent, Z. W., Wilson, R., et al. (2017). Massively parallel digital transcriptional profiling of single cells. Nat. Commun. 8:14049. doi: 10.1038/ncomms14049

PubMed Abstract | Crossref Full Text | Google Scholar