Evolutionary insights into primate skeletal gene regulation using a comparative cell culture model

Abstract

The evolution of complex skeletal traits in primates was likely influenced by both genetic and environmental factors. Because skeletal tissues are notoriously challenging to study using functional genomic approaches, they remain poorly characterized even in humans, let alone across multiple species. The challenges involved in obtaining functional genomic data from the skeleton, combined with the difficulty of obtaining such tissues from nonhuman apes, motivated us to consider an alternative in vitro system with which to comparatively study gene regulation in skeletal cell types. Specifically, we differentiated six human (Homo sapiens) and six chimpanzee (Pan troglodytes) induced pluripotent stem cell lines (iPSCs) into mesenchymal stem cells (MSCs) and subsequently into osteogenic cells (bone cells). We validated differentiation using standard methods and collected single-cell RNA sequencing data from over 100,000 cells across multiple samples and replicates at each stage of differentiation. While most genes that we examined display conserved patterns of expression across species, hundreds of genes are differentially expressed (DE) between humans and chimpanzees within and across stages of osteogenic differentiation. Some of these interspecific DE genes show functional enrichments relevant in skeletal tissue trait development. Moreover, topic modeling indicates that interspecific gene programs become more pronounced as cells mature. Overall, we propose that this in vitro model can be used to identify interspecific regulatory differences that may have contributed to skeletal trait differences between species.

Fig 1. Comparative skeletal cell culture model.

Schematic of the differentiation protocol and the stages at which single-cell RNA-seq data were collected (A), along with descriptions of the human and chimpanzee cell lines used (B), a diagram of the overall study design (C), and cell images from one human cell line and one chimpanzee cell line at each stage of differentiation (D). Pluripotent cells (Time 0) and mesenchymal cells (Time 1) are phase contrast images at 4X magnification, and osteogenic cells (Time 2) are stained with Alizarin Red and zoomed out to display the entire cell culture well. Silhouette images were adapted from http://phylopic.org/ and courtesy of T. Michael Keesey and Tony Hisgett (http://creativecommons.org/licenses/by/3.0/).

Fig 2. scRNA-seq data recovery is similar across species and reproducible across replicates.

UMAP dimensional reduction plot of scRNA-seq data with cells labeled by the stage of differentiation at which they were collected (A), along with a bar plot depicting the number of chimpanzee and human cells collected at each stage of differentiation (B), violin plots displaying the distribution of UMI counts per cell (C) and gene counts per cell (D), and the correlation of average gene expression patterns between technical replicates (human and chimpanzee) collected in pluripotent cells (Time 0) (E), mesenchymal cells (Time 1) (F), and osteogenic cells (Time 2) (G). Enrichment of GO functional categories among marker genes for cells collected at each stage of differentiation (H-J). The top 5 GO functions identified in biological processes (BP), cell components (CC), and molecular functions (MF) are displayed along with the adjusted p-value (p-adjust), the number of marker genes overlapping a GO function (Count), and the ratio of overlapping to non-overlapping marker genes for a given GO function (GeneRatio).

Fig 3. Interspecific DE across three stages of osteogenic differentiation.

Bar plot showing the number of interspecific DE genes identified for each stage of differentiation using standard methods (A). Correlation motifs based on the probability of differential expression between species for each stage of differentiation (B) and correlation motifs based on the probability of differential expression across stages of differentiations for each species (C) with the number of genes assigned to each motif shown in the bar plot on the right and the posterior probability that a gene is DE shown by the shading of each box. Box plots of the log2 transformed gene expression variance values for cells collected at each stage of differentiation for each species (*** p<0.001) (D). Enrichment of external DE gene sets among Cormotif interspecific DE genes identified for each stage of differentiation for validation (E) and functional interpretation (F) with the p-value (p.value), the number of DE genes overlapping an external gene set (DE.Interest), and the ratio of overlapping to non-overlapping DE genes for a given external gene set (GeneRatio) denoted.

Fig 4. Interspecific DE across five substages of osteogenesis.

UMAP dimensional reduction plot of scRNA-seq data with cells labeled by the stage of osteogenesis to which they were assigned (A), along with a simplified schematic of the osteogenic ad hoc assignment method. Bar plot showing the number of interspecific DE genes identified for each stage of osteogenesis using standard methods (C). Correlation motifs based on the probability of differential expression between species for each osteogenic ad hoc assignment with the number of genes assigned to each motif shown in the bar plot on the right and the posterior probability that a gene is DE shown by the shading of each box (D). Enrichment of external DE gene sets among Cormotif interspecific DE genes identified for each stage of osteogenesis for validation (E) and functional interpretation (F) with the p-value (p.value), the number of DE genes overlapping an external gene set (DE.Interest), and the ratio of overlapping to non-overlapping DE genes for a given external gene set (GeneRatio) denoted.

Fig 5. Examining cell types at different resolutions using discrete and continuous perspectives.

UMAP dimensional reduction plot of scRNA-seq data with cells labeled by the unsupervised cluster (resolutions of 0.05, 0.25, and 0.50) to which they were assigned (A), and structure plots showing the results of topic modeling at k = 3, k = 4, k = 5, k = 6, and k = 7 with each row representing the gene expression profile from one cell, each colored bar representing a topic, and the grade of membership in each topic depicted by the length of the bar along the x-axis (B). In the structure plots, cells are grouped by their species of origin and collection time point. In both sets of plots, the key notes the top GO category enrichment of marker genes for a given cluster or topic.