Shaping modern human skull through epigenetic, transcriptional and post-transcriptional regulation of the RUNX2 master bone gene | Scientific Reports…

Study of RUNX2 locus in modern and ancient humans

The sequences of the RUNX2 locus (Chr6:45,318,00045,670,000 hg38, for a total of 352.000 bases) have been aligned in AMH and ancient species (Neandertal and Denisovan) genomes to map the changes occurred during recent evolution. We have identified 459 and 470 changes (including nucleotide substitutions, indels, and annotated human SNPs) acquired in AMH compared with Neandertal and Denisovan, respectively (Supplementary File 1 and Fig.2). Our analysis confirmed that no change occurred within RUNX2 exons. Selected changes map within the P1 and P2 promoter sequences, and in both proximal and distal 3UTRs (Fig.2). In addition, to confirm which of these nucleotide changes have occurred during the recent human evolution, we also analysed the RUNX2 sequences of non-human Primates [chimpanzee (Pan troglodytes), gorilla (Gorilla gorilla), Sumatran orangutan (Pongo abelii), gibbon (Nomascus leucogenys) and rhesus macaque (Macaca mulatta)] focusing our attention on changes mapping in promoters and 3UTRs (Supplementary File 1).

Divergence of RUNX2 locus in AHM compared with ancient hominins. The pie charts display the extent of genomic divergence in terms of number of sequence changes and corresponding locations found in AMH compared with the Neandertal (a) and the Denisovan (b) species identified using the Integrative Genomics Viewer (IGV) tool. Most of the sequence changes are found in intronic regions, whereas few substitutions occur in the regulatory regions outside the open reading frame: Promoter 1, Promoter 2, proximal 3UTR and distal 3UTR. The Neandertal skull in the upper left panel of this figure is exhibited at the Natural History Museum of London (personal picture, credits to W. Lattanzi); the AMHs and Denisovans skulls were created with BioRender.com (note that a reliable reconstruction of the Denisovan skull is not available).

The P1 sequence in AMH differs for as little as one nucleotide from Neandertal and Denisovan genes (Supplementary File 1 and Fig.2). In the P2 sequence, we found three changes in AMH compared with both ancient hominins, plus an additional change differing exclusively between AMH and Neandertal (Supplementary File 1 and Fig.2). Interestingly, three substitutions mapping in the P2 promoter were not present in any of the tested Primates species, which retained the same nucleotides observed in the ancient hominins species. These changes seem indeed to be specific to AMH (Supplementary File 1).

The proximal 3UTR sequence of AMH differs in 4 and 2 nucleotides, compared to Neandertal and Denisovan species, respectively (Supplementary File 1 and Fig.2). Our analysis revealed that one change in the proximal 3UTR identified in Denisovan RUNX2 sequence is also conserved in all non-human Primates analysed (Supplementary File 1). The sequence of the distal 3UTR region of AMH, instead, differs from both Neandertal and Denisovan due to two nucleotide substitutions (Supplementary File 1 and Fig.2). All the tested Primates genome sequences differ from both modern and ancient humans in one of these changes (Supplementary File 1).

RUNX2 locus also includes two additional genes, recently classified as lncRNAs antisense to RUNX2, namely, AL096865.1 and RUNX2 antisense RNA 1 (RUNX2-AS1) (Fig.1). AL096865.1 partially overlaps with the RUNX2 P2 region and contains 3 sequence changes that occurred in the AMH genome compared with both Neandertal and Denisovan (i.e. the same changes already described for P2 sequence; Supplementary File 1). RUNX2-AS1 encodes a transcript (ENST00000563807.1) with two exons, and maps on the reverse strand of the last intron of the RUNX2 ENST00000576263.5 coding isoform (see Fig.1). Our results showed that 6 changes map within the RUNX2-AS1 gene, 2 of which fall in the first exon (Supplementary File 1). Three of the substitutions mapping in the RUNX2-AS1 locus were not found in the tested Primates genomes, which shared the same nucleotides observed in the ancient hominins (Supplementary File 1).

To assess the putative functional consequence of the sequence changes observed in the RUNX2 gene regulatory regions, we have performed in silico predictive analysis.

We first focused on the analysis of nucleotide substitutions within the two promoter regions, P1 and P2. As discussed above, P1 sequence in AMH differs from that of Neandertal and Denisovan for just 1 nucleotide. Downstream analysis using different bioinformatic tools (Meme suite https://meme-suite.org/meme/, Consite http://consite.genereg.net/ and rVista 2.0 https://rvista.dcode.org/)27,28,29 did not reveal any association of this genomic substitution with putative regulatory regions.

For P2, we found four nucleotide changes, three differing in AMH compared with both ancient hominins, and 1 differing exclusively between AMH and Neandertal. Three of these substitutions also overlap with the antisense lncRNA AL096865.1 (Supplementary File 1). Motif analysis using MEME suite27 indicated that 3 substitutions map within 3 conserved CCYCCCWCCTC sequence motifs (Fig.3a, red boxes). Interestingly, analysis of transcription factor binding using Tomtom30 identified this motif as a putative binding site for the C2H2-type zinc finger protein ZNF263 (Fig.3b). We therefore applied position weigh matrix (PWM) analysis to investigate if the identified substitutions could affect binding of this transcription factor to RUNX2 P2 promoter. Analysis of position-specific scoring matrices using available software (FIMO, Find Individual Motif Occurrences, https://meme-suite.org)31 indicated that: (i) two of these substitutions (A/G in position 45421497 and T/C in position 45421552) decreased ZNF263 matrix score in AMHs from 17.102 to 15.734, and from to 6.637 to 4.959, respectively, while (ii) the inclusion of a T in Denisovan and AMHs (position 45421406) increases the matrix score from 4.204 to 4.653 (Supplementary File 2). ZNF263 was recently shown to participate in the CCCTC binding factor (CTCF)-mediated chromatin looping32 and its motif is also enriched in lncRNAs located at chromatin boundaries33. As chromatin looping was previously associated with the regulation of RUNX234, we hypothesize that the tested nucleotide changes could alter chromatin looping and therefore influence the levels of gene expression. To support this hypothesis, we first investigated if P2 promoter was indeed bound by CTCF and ZNF263. To this end, we explored available ENCODE (https://www.encodeproject.org/) genome-wide datasets from different human cell lines (osteoblasts, HEK293 and K562 cells), using the UCSC Genome browser resources (https://genome.ucsc.edu). Data shown in Fig.3c show that CTCF binding is enriched at the P2 promoter in all cell lines analysed, whereas ZNF263 binding was found in the two cell lines for which there were ENCODE data available (HEK293 and K562). Finally, analysis of the GeneHancer database revealed different loop interactions involving the RUNX2 P2 promoter and other regulatory elements in the region, such as SUPT3H (SPT3 homolog, SAGA and STAGA complex component) promoter (Fig.3c). Altogether, these analyses suggested that the genomic changes at the ZNF263 binding sites comprised within the RUNX2 P2 promoter/AL096865.1 lncRNA sequence could influence DNA looping and therefore differentially regulate RUNX2 expression in AMH as compared with ancient hominins.

Genomic variants in RUNX2 P2 promoter overlap with ZNF263 binding sites at CTCF-mediated loop anchor regions. (a) Motif analysis using MEME suite of the P2 proximal promoter region, containing a cluster of three nucleotide variants in AMHs compared with Denisovan/Neandertal species (chr6:45,421,39845,421,567), identified three CCYCCCWCCTC motifs overlapping with the three genomic variants. (b) In silico analysis using Tomtom identified ZNF263 as a transcription factor binding the CCYCCCWCCTC motif. (c) Scheme of the RUNX2 locus (Chr6:45,318,00045,670,000 GRCh38/hg38) showing the different RUNX2 isoforms and regulatory regions and available Encode data. These include: ChIP-seq tracks (fold change over control) for CTCF binding (light blue) in osteoblasts, HEK293 and K562 cell lines, ChIP-seq tracks (fold change over control) for Znf263 binding (dark blue) in HEK293 and K562, short nucleotide clinical variants ClinVar SNVs and GeneHancer regulatory elements (GH Reg Elem (DE)), including 3D chromatin loop formation.

We then focused further analyses on the study of the nucleotide substitutions within the regulatory elements at the 3 ends of RUNX2, to evaluate possible implications in post-transcriptional regulation.

Using the UCSC Genome browser resources, we observed that the proximal 3UTR belongs to a hotspot region of variants annotated in the ClinVar database (https://www.ncbi.nlm.nih.gov/clinvar/; see Fig.3c). As mentioned above, the proximal 3UTR sequence of AMH differs in six nucleotides respect to ancient hominis (in particular, four respect to Neanderthal and two nucleotides compared to Denisovan). Of these, four nucleotides specific for ancient hominis represent common variants annotated in the Homo sapiens genome (rs144321470, rs188598788, rs537488922, rs182124295), as benign mutations in CCD cohorts. Also, the variants in the distal 3UTR of Neandertal and Denisovan correspond to common variants observed in Homo sapiens (rs74697776, rs6458457) that are not annotated in ClinVar.

The 3'UTR of a mRNA can be bound by miRNAs that modify the stability and half-life of the target transcripts and/or protein translation. We thus tested the hypothesis that the nucleotide changes found in the 3UTRs of RUNX2 could influence miRNA binding. We first used different online databases and tools (miRDB http://mirdb.org/, TargetScanHuman http://www.targetscan.org/vert_72/, miRbase http://www.mirbase.org/ and PolymiRTS Database 3.0 http://compbio.uthsc.edu/miRSNP/), to select the list of miRNAs that, based on sequence complementarity, most likely bind the RUNX2 3UTRs. Then, we annotated all the miRNAs whose seed sequence fall into one of the sequence regions that diverge between AMH and ancient hominins.

Our results showed that the substitutions fixed in the proximal 3UTR of RUNX2 affected the binding of different miRNAs, indicating that during recent human evolution, the AMH RUNX2 genomic region acquired or lost miRNA binding sites, potentially involved in post-transcriptional regulation (Supplementary File 3 and Fig.4a).

In silico prediction of miRNA binding sites within RUNX2 3UTRs in AMH, Neandertal and Denisovan. The picture shows the nucleotide changes (highlighted in red) identified in the untranslated regions (UTRs) at the 3 ends of RUNX2 transcripts in Neandertal and Denisovan hominins compared with that of anatomically modern humans (AMH). The substitutions are apparently able to modify the pattern of miRNA binding both in the proximal 3UTR (a) and in the distal 3UTR (b).

Our data showed also that the distal 3UTR does not differ between Neandertal and Denisovan species, suggesting a similar miRNA-dependent post-transcriptional regulation pattern. We predicted that the sequence divergence between AMH and ancient hominins observed in this 3UTR fall within the binding sites of 16 different miRNAs in AMH (Supplementary File 3 and Fig.4b). In particular, the in silico analysis suggested that the 2 nucleotide changes that were fixed in AMH generate new binding sites for several miRNAs (Fig.4b). The first substitution in the distal 3UTR found in the ancient hominins instead should allow the binding of an additional miRNA, namely miR-5589-5p (Fig.4b).

All the genomic loci of the selected miRNAs were comparatively analysed in the aligned sequences of AMH and ancient hominins, using the IGV tool, to evaluate the evolutionary conservation of their DNA sequences across human species. Our analysis showed that nucleotide changes occurred in the AMH miR-3143, compared with the Denisovan sequence, and in the AMH miR-149-3p, compared both with Neandertal and Denisovan species (Supplementary File 4). Nonetheless, both these fixed changes fall outside the seed sequences of these two miRNAs. All the other miRNAs did not present any nucleotide change.

To delve into functional aspects, we have assessed the differential contribution of the alternative transcript isoforms of the RUNX2 gene to the osteogenic cascade activation. To this end, we have analysed the expression of the different RUNX2 transcript variants in mesenchymal stromal cells isolated from calvarial sutures (CMSC) isolated from unfused suture tissues obtained as surgical waste from surgery of nonsyndromic craniosynostosis (NCS) patients. We relied on this cellular model as a benchmark for the study of membranous ossification, previously optimized and standardized in our lab35. P1-derived transcripts, P2-derived transcripts, splice variants containing the proximal 3UTR and the distal 3UTR, were independently amplified by Real-Time PCR using isoform-specific oligonucleotide primer pairs.

CMSC were induced toward the osteogenic lineage up to 21days in order to assess the regulation of each transcript group during the in vitro differentiation process at different timepoints (i.e. 3, 6, 10, 14 and 21days). Our results showed that total RUNX2 levels undergo an initial down-regulation during the first days of commitment, followed by an increase thereafter (Fig.5a). The upregulation of the RUNX2 P1-derived isoforms and of the RUNX2 isoforms containing the proximal 3UTR was significant until 21days of osteogenic induction (Fig.5b,c). Instead, the expression of RUNX2 P2-derived isoforms and of the RUNX2 isoforms with the distal 3UTR appeared to reach a peak at 10days of in vitro differentiation and thereafter decrease or lose statistical significance (Fig.5d,e). The expression of the Sp7 transcription factor gene (osterix, OSX) and of the alkaline phosphatase (ALP) gene were evaluated to confirm that the cells were differentiating along the osteogenic lineage (Fig.5f,g).

Relative expression of RUNX2 isoforms during osteogenic induction in vitro. The graphs show the expression of total RUNX2 levels (a), RUNX2 P1-derived transcripts (b), RUNX2 splice variants containing the proximal 3UTR (c), P2-derived transcripts (d), isoforms with the distal 3UTR (e), OSX (f) and ALP (g) transcripts measured by Real Time PCR, in cells undergoing osteogenic induction up to 21days. The fold change of expression levels is expressed as relative quantity (RQ), calculated according to the Ct method, setting the day 0 (when the standard growth medium has been replaced with the osteogenic medium) value as reference. All data are expressed as mean fold changestandard deviation across replicates. *p0.05; **p0.01; ***p0.001; ****p0.0001.

The expression of the two long noncoding transcripts AL096865.1 and RUNX2-AS1 was also assessed in CMSC and correlated to the levels of RUNX2, to evaluate the possible involvement of these two lncRNAs in regulating RUNX2 expression. We focused our analysis during the early stages of the osteogenic commitment (0-to-10days of osteogenic induction), as a functionally relevant timeframe for gene expression regulation. AL096865.1 levels were significantly upregulated after 6 and 10days of differentiation (Fig.6a); also RUNX2-AS1 levels showed an increasing trend during osteogenic differentiation (Fig.6b). Furthermore, our analysis showed a significant positive correlation between the expression of AL096865.1 and that of RUNX2 total isoforms (Fig.6c), of RUNX2 P2-derived isoforms (Fig.6d) and of RUNX2 isoforms containing the proximal 3UTR (Fig.6e). A significant positive correlation could be also outlined for RUNX2-AS1 and the RUNX2 isoforms containing the proximal 3UTR (Fig.6f). These data might suggest a feasible involvement of AL096865.1 and RUNX2-AS1 in modulating RUNX2 expression.

lncRNAs expression during osteogenic induction in vitro. (a,b) The bar graphs show the expression of AL096865.1 (a) and RUNX2-AS1 (b) measured by real time PCR, in cells undergoing osteogenic induction up to 10days. The fold change of expression levels is expressed as relative quantity (RQ), calculated according to the Ct method, setting the day 0 (when the standard growth medium has been replaced with the osteogenic medium) value as reference. All data are expressed as mean fold changestandard deviation across replicates. *p0.05; **p0.01. (cf) The line graphs draw the linear regressions during the first phases of osteogenic induction (up to 10days) of the relative expression levels of AL096865.1 with total RUNX2 levels (c), RUNX2 P2-derived transcripts (d) and RUNX2 splice variants containing the proximal 3UTR (e), and of RUNX2-AS1 with the group of RUNX2 isoforms containing the proximal 3UTR (f). Pearson correlation was evaluated and the Pearson correlation coefficient (r value) and p-value (p) were reported.

We further focused our analysis on the distal 3UTR of RUNX2 (ENST00000576263.5), given the structural peculiarity of this gene, which includes an additional coding exon (see Fig.1), and taking into consideration that the nucleotide substitutions identified in the proximal 3UTR were annotated as clinical variants (though benign, to date).

In silico protein modelling, based on the I-TASSER tool, revealed that the presence of this terminal exon introduces a leucine-zipper DNA binding motif in the protein (Fig.7a), likely affecting RUNX2 function as a transcriptional regulator.

In silico modelling and miRNA profiling. (a) The images show the three-dimensional computational modelling of the last coding exon of the RUNX2 isoform containing the distal 3UTR and its predicted interaction with DNA, by I-TASSER tool (https://zhanglab.ccmb.med.umich.edu/I-TASSER/). (b) Surface plasmon resonance data showed the changes in miRNA-target binding affinities. Representative graph showing the interaction of the synthetic sequence that mimics the miR-3150a-3p with the 3UTR of anatomically modern human (AMH) at different concentration (RU: response unit). Empty circles represent the experimental traces at different analytes concentration, while the black lines are the global fit using a 1:1 kinetic model. (ch) The bar graphs show the changes in relative expression of miR-3150a-3p, miR-6825-5p, miR-4700-5p, miR-149-3p, miR-4486 and miR-6785-5p measured by Real Time PCR, in cells undergoing osteogenic induction up to 21days. The fold change of expression levels is expressed as relative quantity (RQ), calculated according to the Ct method, setting the day 0 (when the standard growth medium has been replaced with the osteogenic medium) value as reference. All data are expressed as mean fold changestandard deviation across replicates. *p0.05; **p0.01; ****p0.0001. (i,j) The line graphs show the linear regressions during the osteogenic induction (up to 21days) between the relative expression levels of miR-3150a-3p and miR-6785-5p and RUNX2 total isoforms. Pearson correlation was analysed and the relative Pearson correlation coefficient (r value) and p-value (p) were reported.

We thereafter investigated in vitro if the nucleotide changes in this 3UTR actually affect miRNA binding. To this end, we used Surface Plasmon Resonance (SPR) to determine the affinity of the molecular interactions between the 15 selected miRNAs and the target sequences on the distal 3UTR of AMH and Neandertal/Denisovan RUNX2 genes. The resulting dissociation constant values obtained for each tested miRNA are reported in Table 1. Of the 15 miRNAs tested, eight failed to show any interaction. The remaining seven miRNAs showed differential binding affinities to the species-specific RUNX2 transcripts: (i) namely, miR-4700-5p, miR-6825-5p, miR-3150a-3p, miR-6763-5p showed interaction exclusively with the AMH RUNX2 gene; (ii) miR-149-3p showed a higher affinity for the AMH RUNX2 mRNA; and (iii) miR-6785-5p and miR-4486 showed higher affinity in the interaction with the ancient species target sequence (see Table 1). A representative figure corresponding to the interaction of the RUNX2 3UTR of AMH with the synthetic sequence that mimics the miR-3150a-3p is reported in Fig.7b.

To validate these data at the functional level in the biological context, we have analysed in CMSC the levels of the miRNAs previously selected through SPR (miR-3150a-3p, miR-6825-5p, miR-4700-5p, miR-149-3p, miR-4486 and miR-6785-5p). Our results showed that all the miRNAs are expressed in CMSC and that their expression vary during osteogenic differentiation. The expression levels measured during the osteogenic induction were variable across the replicates and it was not always possible to identify a clear trend of activation/silencing during the osteogenesis process (Fig.7ch). Due to this technical limitation, we could not detect any significant correlation between the trend of expression of the selected miRNAs and that of the RUNX2 isoform with the distal 3UTR during the osteogenic induction (data not shown). However, the expression levels of miR-3150a-3p and of miR-6785-5p inversely correlated with those of RUNX2 total isoforms: while the expression of RUNX2 increased during the osteogenesis in vitro, the levels of miR-3150a-3p and miR-6785-5p decreased (Fig.7i,j). This may suggest that these two miRNAs could have gained a function as post-transcriptional regulators of RUNX2 in AMHs, tampering excessive production of this transcription factor and delaying ossification.

Go here to see the original:
Shaping modern human skull through epigenetic, transcriptional and post-transcriptional regulation of the RUNX2 master bone gene | Scientific Reports...