Decoding the transcriptome of calcified atherosclerotic plaque at single-cell resolution | Communications Biology – Nature.com

Tissue source and processing

Paired sections of tissue, including both artery and plaque, were recovered from the atherosclerotic core (AC) and proximally adjacent (PA) region of three patients with asymptomatic type VII calcified plaques who underwent carotid endarterectomy (Fig. S1a, TableS1). Due to the rich cellular composition of carotid artery and plethora of debris in plaque (i.e., lipid, fibrinogen, etc.), dissociation and generation of single-cell suspensions amenable to single-cell RNA sequencing were difficult. After tissue collection, enzymatic digestion, RBC lysis, and filtration were the initial steps required to generate single cells (see Methods and Fig. S1b). However, despite efficient enzymatic dissociation and significant filtering of our sample, we were still challenged by abundant plaque debris, which ultimately resulted in poor single-cell capture rates. In order to overcome this issue without isolating specific cell types through cell-marker antibody labeling, we devised a strategy to label all cells in the sample with a far-red excitation-emission live/dead cell nuclear stain (DRAQ5). All cells in the sample were stained, with debris being left unstained by the dye. Previous studies have used nuclear staining in library preparation and sequencing experiments to discriminate single versus doublet cells during cell sorting without adverse effects for downstream applications such as single-cell and bulk RNA sequencing17,18,19,20. Subsequently, DRAQ5+ cells were manually gated and sorted from the remainder of the debris using FACS. Cells isolated from the entire filtered sample represented <1% of the total particles in the sample (Figs. S2aS2f). Viability of remaining cells was assessed and was always >80% using this technique for cell separation (see Methods). The cells were then processed for single-cell sequencing.

The analytical approach in this manuscript is depicted in Fig.1a. Generation of single cells from three patient-matched AC and PA samples (batched per patient on a single NextSeq flow cell) yielded 51,981 cells total, with an average of ~13,000 AC cells/patient and ~5000PA cells/patient. Cell number disparities are due to the difference in size of the AC vs PA tissue itself. Given the abundance of AC versus PA cells, down-sampling was performed to equalize the contribution of each sample and condition to the unsupervised discovery of cell types and to mitigate bias due to class imbalance. UMAP-based clustering (see Methods) of this down-sampled dataset reveals 15 distinct cell partitions (Fig. S1c, d), representing 17,100 cells total. In order to assign partitions to major cell types we examined genes expressed in >80% of cells per partition and at a mean expression count >2. A dotplot representing three marker genes selected for each partition is presented in (Fig. S1e). A comparison of VSMC marker genes used in our study with those in the literature15 is provided (Fig. S1f). Cell-type labels assigned to these 15 initial partitions based on these marker genes include: T-lymphocytes (2 partitions), macrophages, VSMCs (2 partitions), ECs (2 partitions), B-lymphocytes, natural killer T-cells, B1-lymphocytes, mast cells, lymphoid progenitors, plasmacytoid dendritic cells, and an unidentified partition (TableS2). Following doublet filtering using a marker-gene exclusion method (see Supplemental Methods), removal of partitions with too few cells for differential gene expression analysis (mast cells, lymphoid progenitors, plasmacytoid dendritic cells, and the unidentified partition), and merging of partitions assigned to the same cell-type, we assessed differential gene expression between AC and PA regions across the 6 remaining major cell types: macrophages, ECs, VSMCs, NKT cells, T- and B-lymphocytes (Fig.1bd, Fig. S4, Supplementary Data16). We performed a number of independent partitioning experiments using various algorithmic variations to confirm the reproducibility of these partitions and cell-label assignments (see Supplemental Methods).

a Schematic diagram of analytical steps from tissue dissociation to key driver analysis. b, c UMAP visualization of 6 major cell types following doublet removal via gene exclusion criteria (see Supplemental Methods), separated by anatomic location (b), and by cell type (c). d Dotplot depicting cell-type marker genes, resulting in the identification of macrophages, ECs, VSMCs, NKT cells, T- and B-Lymphocytes. Dot size depicts the fraction of cells expressing a gene. Dot color depicts the degree of expression of each gene. n=3 for PA and AC groups.

GWAS results have highlighted biological processes in the vessel wall as key drivers of coronary artery disease (CAD)21. Our prior work has demonstrated the vascular wall to be involved in the most impactful common genetic risk factor for CAD22. Our results here also demonstrate extensive differential expression in these cell types across anatomic locations compared to the remaining cell types. Therefore, we chose to focus our efforts on dissecting expression alterations in VSMCs and ECs in order to illuminate pathogenic genomic signatures within these cell-types. As above, each cell type is compared across anatomic location (Fig.2a, e), and the top differentially expressed genes are shown (Fig.3b, f), revealing interesting spatial and expression magnitude differences between AC and PA cells.

a, e UMAP visualization of VSMCs (a) and ECs (e), separated by anatomic location. b, f Volcano plots of the top differentially expressed genes in VSMCs (b) and ECs (f). Dotted lines represented q-value 0.5 and <0.5 corresponding to PA and AC cells, respectively. c, d UMAP visualization of the top 4 upregulated genes in AC VSMCs (c), and PA VSMCs (d). Gray-colored cells indicate 0 expression of designated gene, while color bar gradient indicates lowest (black) to highest (yellow) gene expression level. g, h UMAP visualization of the top 4 upregulated genes in AC ECs (g), and PA ECs (h). Color scheme is similar to the above-described parameters. VSMCs=3674 cells; ECs=2764 cells. n=3 for PA and AC groups.

a, b Normalized enrichment score (NES) ranking of all significant PA and AC Hallmarks generated from GSEA analysis of differentially expressed genes for VSMCs (a) and ECs (b) (FDR q-value<0.05). c Fully clustered on/off heatmap visualization of overlap between leading edge EMT hallmark genes generated by GSEA. Heatmaps are downsampled and represent 448 cells from each cell type and anatomic location (1792 total cells). A dotplot corresponding to gene expression levels for each cell type in the heatmap is included. Dot size depicts the fraction of cells expressing a gene. Dot color depicts the degree of expression of each gene. d Volcano plot of differentially expressed genes between the two groups of VSMCs in (c). Dotted lines represented q-value<0.01 and normalized effect >0.5 and <0.5. e, f Gene co-expression networks generated from VSMC Module 13 (d) and EC Module 1 (e) representing the EMT hallmark from GSEA analysis. Genes are separated by anatomic location (red=AC genes, cyan=PA genes), differential expression (darker shade=higher DE, gray=non-significantly DEGs), correlation with other connected genes (green line=positive correlation, orange line=negative correlation) and strength of correlation (connecting line thickness). Significantly DEGs (q<0.05) with high connectivity scores (>0.3) are denoted by a box instead of a circle. n=3 for PA and AC groups.

VSMCs generate three subclusters in the UMAP plot. A large fraction of PA VSMCs form a PA-specific VSMC subcluster. In contrast, AC VSMCs form 2 separate clusters both of which are intermingled with PA VSMCs. This suggests VSMCs occupy three major cell states, including one completely distinct PA subtype, and two that are predominantly AC VSMCs (Fig.2a). The top four upregulated genes in the AC are sparsely expressed and include SPP1, SFRP5, IBSP, and CRTAC1 (Fig.2b, c), while APOD, PLA2G2A, C3, and MFAP5 are upregulated in many PA VSMCs (Fig.2b, d).

The spatial clustering of upregulated genes in AC VSMCs suggests the presence of separate subpopulations of matrix-secreting VSMCs involved with ECM remodeling (Fig.2c). SPP1 (osteopontin) is a secreted glycoprotein involved in bone remodeling23 and has been implicated in atherosclerosis for inhibiting vascular calcification and inflammation in the plaque milieu24. IBSP (bone sialoprotein) is a significant component of bone, cartilage, and other mineralized tissues25. CRTAC1 is a marker to distinguish chondrocytes from osteoblasts and other mesenchymal stem cells26,27. These findings suggest the presence of cartilaginous and osseous matrix-secreting VSMCs in the AC region. SFRP5, an adipokine that is a direct WNT antagonist, reduces the secretion of inflammatory factors28 and is thought to exert favorable effects on the development of atherosclerosis29. The high expression of SFRP5 in the AC suggests a deceleration of these inflammatory processes in the core of the plaque, and an overall shift in the AC to calcification and matrix remodeling.

Conversely, the upregulated genes in PA VSMCs are more ubiquitously expressed by VSMCs in a PA-specific region of the UMAP plot (Fig.2d). C3 is highly differentially expressed in many PA cells (Fig.2d). Complement activation has long been appreciated for its role in atherosclerosis30, with maturation of plaque shown to be dependent, in part, on C3 opsonization for macrophage recruitment and stimulation of antibody responses31. Its predominance in our PA samples suggests complement activation in atherosclerosis is anatomically driven by VSMCs located adjacent to areas of maximal plaque build-up. PLA2G2A (phospholipase A2 group IIA), also selectively expressed by this group of cells, is pro-atherogenic, modulates LDL oxidation and cellular oxidative stress, and promotes inflammatory cytokine secretion32, further facilitating the inflammatory properties of this group of VSMCs. Full differential expression results for VSMCs are provided (Supplementary Data5).

Overall, we identify increased calcification and ECM remodeling by VSMCs in the AC versus pro-inflammatory signaling by VSMCs in the PA. These differences in biological processes are strongly supported further in the systems analyses below.

In contrast to VSMCs, for ECs we observe a more complete separation of cells into two distinct subgroups (Fig.2e). PA ECs significantly outnumber the AC ECs (2316 vs 448 cells, respectively), possibly due to intimal erosion and loss of endothelial cell layer integrity during advanced disease5,33,34,35 resulting in fewer captured ECs in the AC. Cellular transdifferentiation may also cause a subpopulation of ECs to lose common EC marker expression, resulting in lower numbers of ECs identified in AC compared to the PA counterpart. Histologic assessment of AC plaque collected from our patients supports the assertion of endothelial layer attenuation as the principal reason for lower AC EC capture (Fig. S3b, c). In contrast to VSMCs, there is a skew toward higher magnitude expression changes in AC ECs vs PA ECs. The top four upregulated genes are ITLN1, DKK2, F5, and FN1 in the AC and IL6, MLPH, HLA-DQA1, and ACKR1 in PA ECs (Fig.2g, h).

The upregulated genes in AC ECs again suggest a synthetic profile. ITLN (omentin) is an adipokine enhancing insulin-sensitivity in adipocytes36. Interestingly, circulating plasma omentin levels were shown to negatively correlate with carotid intima-media thickness37, inhibit TNF-induced vascular inflammation in human ECs38, and promote revascularization39, suggesting an anti-inflammatory and intimal repair role in AC ECs. DKK2 further indicates intimal repair as it stimulates angiogenesis in ECs40. The significant upregulation of FN1 (fibronectin) in this group further suggests active ECM remodeling and may serve as a marker for mesenchymal cells and EMT-related processes41.

Similarly to PA VSMCs, the upregulated genes in PA ECs suggest an overall inflammatory profile. Central players in inflammation and antigen presentation are upregulated specifically in PA ECs (Fig.2h). IL6, a key inflammatory cytokine associated with plaque42, is the most upregulated gene. Furthermore, ACKR1, highly upregulated in many PA ECs, binds and internalizes numerous chemokines and facilitates their presentation on the cell surface in order to boost leukocyte recruitment and augment inflammation43. Antigen presentation on ECs via HLA-DQA1 (MHC class II molecule) may support activation and exhaustion of CD4+ T-cells44,45 as previously described. Full differential expression results for ECs are provided (Supplementary Data6).

Overall, we identify two main EC subtypes: synthetic ECs in the AC that appear to participate in intimal repair, revascularization, and ECM modulation, and inflammatory ECs in the PA region that likely facilitate inflammation via antigen/chemokine presentation and recruitment of immune cells, including CD4+ T-cells. These differences in biological processes are strongly supported further in the systems analyses below.

In order to explore the anatomic differences for these cell types further, gene set enrichment analysis (GSEA) was used to asses hallmark processes most significantly altered in VSMCs and ECs (Fig.3a, b). Epithelial to mesenchymal transition (EMT), oxidative phosphorylation, and myogenesis gene upregulation were strongly enriched in both AC VSMCs and ECs, collectively suggesting an increase in cellular metabolic activity and proliferation. In contrast, a distinctly inflammatory profile was seen in PA VSMCs and ECs, with IFN gamma/alpha responses and TNFa signaling via NFkB dominating the enriched processes in these groups of cells. Because EMT and TNFa signaling were both shared and strongly enriched processes in the two cell types, the gene signatures associated with these hallmarks were further scrutinized through generation of heatmaps consisting of leading-edge differentially expressed genes from each hallmark process (EMTFig.3c, TNFa signaling via NFkBFig. S5a).

While overlapping at the hallmark level, separation of cells by cell type as well as anatomic location in the EMT hallmark heatmap suggests the overlapping processes are mediated by distinct gene sets in each cell type. Moreover, analysis of EMT hallmark genes further supports the presence of 2 cellular subtypes of AC VSMCs as they appear to cluster into two distinct groups of cells with dichotomous expression of contractile (MYL9, TPM2, TAGLN, FLNA) versus synthetic/EMT (POSTN, LUM, FBLN2, DCN, PCOLCE2, MGP, COL3A1) gene signatures (Fig.3c, d). These results indicate a group of VSMCs in the AC may perform the contractile functions of the blood vessel wall, while the other group of VSMCs may be involved with CTD and ECM remodeling. Furthermore, cells with an ACTA2+Thy1 gene signature in Fig.3c may be, in part, plaque-stabilizing myofibroblasts (orange line), indicating that these contractile cells may also have a large role in ECM remodeling.

In contrast to distinct subclustering of cells by EMT-related genes, there appears to be a common gradient of genes involved in inflammation and response to inflammation expressed throughout the atherosclerotic tissue, with higher levels of TNF-related inflammatory genes expressed in PA VSMCs and ECs compared to AC cells, indicating a predominance of inflammatory processes occurring in the PA region overall (Fig. S5a). Collectively these genes (EIF1, FOS, JUN, JUNB, ZFP36, PNRC1, KLF2, IER2, CEBPD, NFKBIA, GADD45B, EGR1, PPP1R15A, and SOCS3), in addition to IL6 expression in PA ECs, appear to coordinate the inflammatory response pathways in plaque and its adjacent structures. All cell types analyzed thus far are coordinated along this gradient of inflammation.

To further dissect VSMC and EC anatomical gene expression differences in order to identify candidate key genes driving the significant hallmark processes, we reconstructed gene co-expression networks using a partial correlation-based approach (see Methods), defined modules by clustering, and overlaid differential expression analysis results on these modules to identify those enriched in genes differentially expressed between AC and PA tissues.

Using this strategy, 31 and 39 distinct gene network modules were generated in our VSMC and EC datasets, respectively (see Supplemental Methods, Supplementary Data7, 8). Of these, 8 modules in VSMCs, and 5 modules in ECs were enriched with differentially expressed genes (p-value<0.05, Fishers exact test, see Methods). Furthermore, differentially expressed EMT-related hallmark genes overlapped significantly and specifically with a single VSMC and EC module. Differentially expressed TNFa signaling via NFkB-related hallmark genes also overlapped significantly with one VSMCs and EC module (p-value<0.05, Fishers exact test). No other hallmark processes overlapped with generated network modules.

The EMT gene signature generated from GSEA analysis of network modules and the robust upregulation of genes found in matrix-secreting cells in this cohort suggests the possibility of CTD occurring and/or completing in the atherosclerotic core. Therefore, in order to further characterize genes which may stimulate CTD in AC VSMCs and ECs we examined gene co-expression networks in conjunction with differential expression data from the modules enriched with EMT hallmark genes. In VSMCs we identified 9 genes (SPP1, IBSP, POSTN, MMP11, COL15A1, FN1, COL4A1, SMOC1, TIMP1) whose expression was significantly upregulated in AC cells and with strong network connectivity (see Methods). Among these genes we identify POSTN, SPP1, and IBSP as possible key gene drivers of CTD processes in AC VSMCs due to their strong central connectivity and high degree of differential expression in the network module (Fig.3e). POSTN (periostin) is expressed by osteoblasts and other connective tissue cell types involved with ECM maturation46 and stabilization during EMT in non-cardiac lineages47,48. POSTN, SPP1, and IBSP are highly interconnected in our network and likely serve as drivers of CTD by modulating other correlated genes such as TIMP1, VCAN, TPST2, SMOC1, MMP11, FN1, and COL4A1 (Fig.3e), all genes which are involved with cellular differentiation49 and extracellular matrix remodeling50,51.

In our EC network we identified 18 genes (ITLN1, FN1, OMD, S100A4, SCX, PRELP, GDF7, TMP2, SERPINE2, SLPI, HEY2, IGFBP3, FOXC2, RARRES2, PTGDS, TAGLN, LINC01235, and COL6A2) whose expression was significantly upregulated in AC cells and with strong network connectivity. Among these genes, we identify ITLN1, S100A4, and SCX as possible gene drivers of CTD in ECs associated with the AC (Fig.3f). ITLN1 (omentin) is highly upregulated in ECs associated with the atherosclerotic core, and network data indicate it is strongly correlated with genes involved with cellular proliferation and ECM modulation. ITLN is also strongly correlated to OGN (osteoglycin) which induces ectopic bone formation52, indicating that ITLN1 may modulate ECs with osteoblast-like features in the atherosclerotic core. SCX (scleraxin), a transcription factor that plays a critical role in mesoderm formation, and the development of chondrocyte lineages53, as well as regulating gene expression involved with ECM synthesis and breakdown in tenocytes54, is co-expressed with IL11RA, an interleukin receptor implicated in chondrogenesis55, as well as with a variety of genes involved with cellular development and modulation of ECM structures. Thus, SCX may modulate chondrocyte-like ECs in the AC. S100A4 is a calcium-binding protein that is highly expressed in smooth muscle cells of human coronary arteries during intimal thickening56, and silencing this gene in endothelial cells prevents endothelial tube formation and tumor angiogenesis in mice57. Co-expression with HEY2, a transcription factor involved with NOTCH signaling and critical for vascular development58, may indicate an important role in repair via re-endothelialization of plaque-denuded artery.

Next, genes critical to stimulating TNFa signaling via NFkB in PA VSMCs and ECs were evaluated. In the VSMC module we identified 14 genes (APOLD1, MT1A, ZFP36, EGR1, JUNB, FOSB, JUN, FOS, RERGL, MT1M, DNAJB1, CCNH, HSPA1B, and HSPA1A) whose expression was significantly upregulated in PA cells and with strong network connectivity. Among these genes we identify immediate-early (IE) genes ZFP36, EGR1, JUNB, FOSB, and FOS as critical response genes in this hallmark process. Importantly, the paired-sample study design in which AC and PA samples from the same patient are processed identically at the same time ensures that these IE genes preferentially upregulated in the PA region are critical for the inflammatory response and not an artifact of tissue processing stressors.

In the EC module we identified two genes (IER2 and FOS) whose expression was significantly increased in PA EC cells (Fig. S5e), and are highly correlated with other critical transcription factors in our network, including FOSB, JUNB, EGR1, and ZFP36, further supporting this group of genes importance in the TNFa signaling hallmark (Fig. S5d).

Finally, in order to identify and characterize refined subpopulations from each anatomic region, we selected the 7 VSMC and 5 EC differentially expressed modules described above and biclustered cells and genes (Fig.4a, d). The likely biological functions of these subpopulations were then inferred based on the genes differentially expressed and subsequent gene ontology enrichment analysis across these subpopulations. A continuous gene expression model, based on the fraction of AC cells per subpopulation, and subsequent gene ontology enrichment analysis was used to evaluate these cell subtype differences (Fig.4b, c, e, f).

a, d Biclustered heatmap visualization of all significant genes (q<0.05) from VSMC (a) and EC (d) modules enriched with differentially expressed genes. a 1224 VSMCs from each anatomic location (2448 cells total). Large color bar denotes PA (cyan) and AC (orange) VSMCs. Small color bar above denotes distinct cell subpopulations (blue, forest green, lime green, brown, purple, magenta, red). d 448 ECs from each anatomic location (896 cells total) in. Large color bar denotes PA (blue) and AC (red) ECs. Small color bar above denotes distinct cell subpopulations (cyan, green, magenta). A dotplot corresponding to gene expression levels for each cell subpopulation on the heatmap is included. Colored dots next to specific genes correspond to critical genes related to the designated cell subpopulation. Continuous gene expression based gene ontology enrichment analysis of biological function performed based on the fraction of AC cells per subpopulation of VSMCs (b, c) and ECs (e, f). n=3 for PA and AC groups.

We identified four cell subpopulations of VSMCs with some overlapping features in our analysis (Fig.4a). The four subpopulations appear to form a continuum of cell states, starting with a population that consists exclusively of PA VSMCs (Fig.4a, green bar), characterized by genes involved in recruitment of inflammatory mediators, with early signs of CTD. Specifically, C3 (opsonization and macrophage recruitment; normalized effect=6.5, q=1.74e07) is highly differentially expressed in this subpopulation and likely augments PA inflammation and macrophage recruitment. This group of VSMCs also shows evidence of early migratory and CTD-like qualities given the expression of FBN1, SEMA3C, HTRA3, and C1QTNF3, (normalized effect=2.77, 3.65, 4.0, 3.58, respectively; q=6.93e41, 1.25e20, 2.53e05, 0.00012, respectively) genes that are both highly differentially expressed in this cohort and with high signal strength in our networks (Fig.4a, Supplementary Data7). FBN1 (ECM component) is strongly correlated with TGFBR3, SEMA3C, and CD248 (modulators of EMT-like processes)59,60,61. Interestingly, this group of cells co-expresses IGSF10, a marker of early osteochondroprogenitor cells62, TMEM119 (bone formation and mineralization; promotes differentiation of myeloblasts into osteoblasts)63,64, and WNT11 (bone formation)65 (Supplementary Data7).

On the other end of this continuum, we identify a subpopulation of ~70% AC cells (Fig.4a, red bar) that have elevated expression of POSTN (osteoblasts; normalized effect=2.206, q=3.60e16), CRTAC1 (chrondrocytes; normalized effect=3.22, q=3.91e26), TNFRSF11B (bone remodeling; normalized effect=0.98, q=7.31e06)66, ENG (VSMC migration; normalized effect=0.87, q=1.41e13)67, COL4A2, and COL4A1 (cell proliferation, association with CAD; normalized effect=0.98, 1.03 and q=3.17e15, 5.68e11, respectively)68,69. Collectively, the differential gene expression data and the underlying biology behind our gene co-expression networks support this group of cells as likely representing synthetic osteoblast- and chondrocyte-like VSMCs which facilitate calcification and cartilaginous matrix-secretion and reside largely in the AC.

Furthermore, gene ontology enrichment analysis provides a clear progression from muscle system processes, extracellular structure reorganization, and catabolic processes enriched in the PA to processes involved with CTD such as ossification, fat cell differentiation, and regulation of cell motility, adhesion, and cellular transdifferentiation enriched in the AC (Fig.4b, c). The shift in cell states supports a continuum of cell state changes leading to increased CTD in the atherosclerotic core.

Overall, we observe three EC subpopulations. Like VSMCs, these cells display transitory properties as they move through a continuum of cell states (Fig.4d). First, there is a group comprised near exclusively of inflammatory PA ECs that is involved in recruitment of inflammatory mediators (Fig.4d, magenta bar). This group has a greater number of cells expressing immune genes such as the cluster of HLA genes, as well as CD74 (normalized effect=1.63, q=2.07e112), a gene which forms part of the invariant chain of the MHC II complex and is a receptor for the cytokine macrophage migration inhibitory factor (MIF)70. The upregulation of MHC class II complex in this subset of PA ECs complements our previous finding of CD4+T-cell recruitment to this subpopulation of PA ECs, leading to over-activation and exhaustion via antigen-persistence.

The next group of cells is intermediate in its composition of AC (67.5%) and PA (~32.5%) ECs with a mixed gene expression profile with characteristics similar to each of the other two groups of cells (Fig.4d, green bar), likely representing dysfunctional ECs that are in transition from the inflamed subtype to the CTD subtype described below.

The final group of cells is largely comprised of ECs from the AC (96.8%) (Fig.4d, cyan bar) and is largely devoid of endothelial-marker gene EMCN71 (normalized effect=0.86, q=1.17e09). Critical EMT genes identified earlier (ITLN1, SCX, and S100A4) are predominantly expressed in this large cluster of AC ECs alongside highly correlated genes OMD, OGN, and CRTAC1, again indicating that this population of ECs likely represents the main group of transdifferentiated ECs.

Gene ontology enrichment analysis further supports this shift in EC cell state from cells primarily involved with immune response (antigen processing and presentation, adaptive immune response, etc.) to cell states predominantly involved with proliferation, migration, vascular development, and angiogenesis (Fig.4e, f).

Excerpt from:
Decoding the transcriptome of calcified atherosclerotic plaque at single-cell resolution | Communications Biology - Nature.com