Comparative analysis of osteoblast gene expression profiles and Runx2 genomic occupancy of mouse and human osteoblasts in vitro

Con en s lis s a ailable a ScienceDi ec Gene jou nal homepage: www.else ie .com/loca e/gene Resea ch pape Compa a i e analysis o os eoblas gene exp ession p ofiles and Runx2 genomic occupancy o mouse and human os eoblas s in i o Ka i Ta kkonen a , Reija Hie a b , Ville Ky ölä b , Ma i Nyk e b,c , Riku Ki i an a a,d,⁎ a Ins i u e o Biomedicine, Uni e si y o Tu ku, Tu ku, Finland b GeneVia Technologies, Tampe e, Finland c Compu a ional Biology, Ins i u e o Biosciences and Medical Technology (BioMediTech), Uni e si y o Tampe e, Tampe e, Finland d Depa men o Endoc inology, Di ision o Medicine, Tu ku Uni e si y Hospi al, Tu ku, Finland ARTICLE INFO Keywo ds: Os eoblas Mesenchymal s em cell MSC Runx2 ChIP ChIP-Seq T ansc ip ome sequencing Mic oa ay Exp ession p ofiling ABSTRACT Fas p og ess o he nex gene a ion sequencing (NGS) echnology has allowed global ansc ip ional p ofiling and genome-wide mapping o ansc ip ion ac o binding si es in a ious cellula con ex s. Howe e , limi ed numbe o eplica es and high amoun o da a p ocessing may weaken he significance o he findings. Compa a i e analyses o independen da a se s acqui ed in he diffe en labo a o ies would g ea ly inc ease he alidi y o he da a. Runx2 is he key ansc ip ion ac o egula ing os eoblas diffe en ia ion and bone o ma ion. We pe o med a compa a i e analysis o h ee published Runx2 da a se s o ch oma in immunop e- cipi a ion ollowed by deep sequencing (ChIP-seq) analysis in os eoblas s om mouse and human o igin. Mo eo e , we assessed he simila i y o he co esponding ansc ip ion da a o hese s udies a ailable online. The ChIP-seq da a analysis confi med gene al ea u es o Runx2 binding, including loca ion a genic s in e genic egions and abundan Runx2 binding on p omo e s o he highly exp essed genes. We also ound high equency o Runx2 DNA binding wi hou a consensus Runx2 mo i a he binding si e. Impo an ly, mouse and human Runx2 showed mode a ely simila binding pa e ns in e ms o peak-associa ed closes genes and hei associa ed genomic on ology (GO) pa hways. Acco dingly, he gene exp ession p ofiles we e highly simila and os eoblas ic pheno ype was p ominen in he diffe en ia ed s age in bo h species. In conclusion, ChIP-seq me hod shows good ep oducibili y in he con ex o ma u e os eoblas s, and mouse and human os eoblas models esemble each o he closely in Runx2 binding and in gene exp ession p ofiles, suppo ing he use o hese models as adequa e ools in s udying os eoblas diffe en ia ion. 1. In oduc ion Fas de elopmen and easy a ailabili y o high h oughpu sequen- cing echnologies has allowed apid accumula ion o unbiased genome wide ansc ip ion da a o se e al cell and issue ypes. Use o ch oma in immunop ecipi a ion combined o NGS sequencing (ChIP- seq) in u n has esul ed in apidly accumula ing esea ch li e a u e on global ansc ip ion ac o (TF) binding si e mapping o many diffe en cell ypes (Cao e al., 2010; Handoko e al., 2011; Heinz and Glass, 2012; Lin e al., 2010; Mikkelsen e al., 2010). Despi e he as p og ess in he field, he e has been ela i ely ew s udies u ilizing ChIP-seq in cells o he os eogenic lineage, eflec ing pe haps he challenging sample ma e ial ma u e bone ma ix p oducing os eoblas s ep esen . Ne e heless, ChIP-seq analyses o Runx2 (Håkelien e al., 2014; Meye e al., 2014b; Wu e al., 2014), C/EBPβ(Meye e al., 2014b), VDR and RXR (Meye e al., 2014a) in os eoblas s ha e ecen ly been epo ed. Because o he high amoun o da a NGS-app oaches p oduce, i is e iden ha he e a e se e al pu a i e candida e genes and mechanisms o be ollowed by u u e s udies, e en ually gene a ing as amoun s o new in o ma ion on he egula ion o os eoblas ogenesis. Howe e , he choice o significan candida es o ake o wa d is challenging as he s a is ical powe emains usually low because gene ally only ew eplica e samples a e included in he indi idual expe imen s. Thus he mo e s udies a e published and compa ed o p e ious da a by bioin o ma ic ools and app oaches, he mo e consis en and eliable in o ma ion conce ning indi idual genes and pa hways can eme ge. Os eoblas s o igina e om mesenchymal s em cells (MSC) ha can gi e ise o a numbe o specialized cell ypes such as adipocy es, h p://dx.doi.o g/10.1016/j.gene.2017.05.028 Recei ed 27 Janua y 2017; Recei ed in e ised o m 10 May 2017; Accep ed 10 May 2017 ⁎ Co esponding au ho a : Ins i u e o Biomedicine, Uni e si y o Tu ku, FI-20520 Tu ku, Finland. E-mail add ess: iku.ki i an a@u u.fi(R. Ki i an a). Abb e ia ions: NGS, nex gene a ion sequencing; ChIP, ch oma in immunop ecipi a ion; ChIP-seq, ch oma in immunop ecipi a ion ollowed by deep sequencing; TF, ansc ip ion ac o ; MSC, mesenchymal s em cell; DE, diffe en ially exp essed; TSS, ansc ip ional s a si e; ALP, alkaline phospha ase; shRNA, sho hai pin RNA; GREAT, genomic egions en ichmen o anno a ions ool; ENA, eu opean nucleo ide a chi e; GO, genomic on ology; YAP, Yes-associa ed p o ein; TAZ, ansc ip ional coac i a o wi h PDZ-binding mo i Gene 626 (2017) 119–131 A ailable online 11 May 2017 0378-1119/ © 2017 Published by Else ie B.V. MARK myoblas s and chond ocy es. Nume ous ho mones, g ow h ac o s and cy okines egula e he diffe en ia ion o an MSC o a specific cell ype, d i en by a se ies o ansc ip ion ac o s ha con ol pheno ype- specific gene exp ession. In case o os eoblas s, ea ly bipo en chond o- os eogenic p ogeni o cells exp ess ansc ip ion ac o SOX9, which is hen ollowed by he exp ession o Run amily ansc ip ion ac o (Runx2) and i s downs eam a ge Os e ix (Osx) in p eos eoblas s (Long, 2012). Runx2 and Osx a e bo h indispensable o os eoblas diffe en ia ion as hei null mu an mice show o al absence o bone and os eoblas s (Ducy, 2000; Nakashima e al., 2002). Runx2 belongs o Run amily o ansc ip ion ac o s (Runx1–3) ha egula e de elopmen and diffe en ia ion o many diffe en cell lineages. Runx2 p o ein con ains a conse ed 128 amino acid Run domain, which is esponsible o he DNA binding and he e odime iza- ion wi h CBFβ, ha enhances Runx2 binding o DNA ( e iewed in Cohen, 2009). Impo an ly, Runx2 can in e ac wi h se e al p o eins including co- egula o y p o eins and ch oma in emodeling ac o s, leading o complex ole in egula ing bone specific genes and diffe - en ia ion. Con ol o Runx2 exp ession is complex, and includes epigene ic mechanisms such as miRNAs and se e al his one modi ying enzymes (Huang e al., 2009; Rojas e al., 2015; Yang e al., 2013). Mo eo e , Runx2 ac i i y is egula ed by se e al pos ansla ional modifica ions such as phospho yla ion, ace yla ion, ubiqui ina ion ( e iewed in Jonason e al., 2009) and sumoyla ion (Kim e al., 2014). Mouse MC3T3-E1 cell line is a clonal non- ans o med cell line es ablished om new bo n mouse cal a ia (Sudo e al., 1983) and is commonly used o s udying os eoblas diffe en ia ion in i o. These cells ep esen Runx2 posi i e p e-os eoblas s commi ed o os eogenic lineage. MC3T3-E1 cells ma u a e o mine alizing os eoblas s in he p esence o s anda d os eogenic medium con aining asco bic acid and Na-β-glyce ophospha e. O human o igin, he e a e ew non- ans- o med os eoblas cell lines. Thus diffe en ia ion o human cells is o en in es iga ed by using p ima y human os eoblas s, immo alized human os eoblas lines (HOBs) o immo alized mesenchymal s em cells (iMSCs). C i ical e alua ion o he cu en os eoblas cell cul u e models is impo an o u he de elopmen o eliable and adequa e in i o models, no only o high s anda d basic esea ch bu also o use in pha maceu ical and bioma e ial esea ch. Typically, pheno ypic assess- men o os eoblas ic cells include measu emen o he exp ession le el o os eoblas ic genes (Runx2, Sp7, Ocn, Opn), alkaline phospha ase (ALP) ac i i y and o ma ion o mine alized bone ma ix in diffe en ia- ion cul u es. Compa ison o diffe en models has o en been limi ed o hese ew pheno ypic p ope ies. In iguingly, h ee pape s we e published epo ing ChIP-seq ana- lyses o genome-wide Runx2 binding in he con ex o os eoblas diffe en ia ion in sp ing 2014. Two o hem, pape s om he labo a- o ies o Wesley Pike (Meye e al., 2014b) and Jane Lian (Wu e al., 2014), desc ibed Runx2 ChIP-seq analyses o mouse MC3T3-E1 in undiffe en ia ed and diffe en ia ed s a e. Pape o Håkelien and co- wo ke s in u n (Håkelien e al., 2014), desc ibed mapping o Runx2 binding si es in iMSCs (Ska n e al., 2014)diffe en ia ed o ma u e os eoblas s in i o. In gene al, he e was high deg ee o simila i y in he esul s epo ed in hese pape s, bu due o he publishing da es close o each o he , he da a se s we e no di ec ly compa ed by any o he au ho s. Ou app oach was o objec i ely e alua e and compa e he published global gene ansc ip ion and Runx2 ChIP-seq da a se s collec ed om hese h ee selec ed os eoblas s udies. Ou specific aims we e 1) o s udy diffe ences in he Runx2 genomic occupancy in mouse MC3T3-E1 os eoblas s du ing hei diffe en ia ion in da a se s p oduced by wo diffe en labo a o ies, o e alua e he ep oducibili y o ChIP- seq expe imen s in his model, 2) o examine he in e species diffe - ences in Runx2 binding pa e ns and a ge genes be ween human and mouse samples and 3) o e alua e he simila i y o gene exp ession p ofiles o undiffe en ia ed and diffe en ia ed os eoblas s om mouse and human o igin. 2. Me hods 2.1. O iginal da a Summa y o he published Runx2 ChIP-seq da a and gene exp ession da a used in he s udy a e desc ibed in Table I. 2.1.1. Da a p e-p ocessing and alignmen Raw ChIP-seq sequencing eads we e downloaded om Eu opean Nucleo ide A chi e (ENA) and he eads we e subjec ed o quali y con ol using Fas QC so wa e (And ews, 2010). Alignmen s o e e - ence genomes we e pe o med using bow ie2 ( e sion 2.1.0) (Langmead and Salzbe g, 2012), wi h eads om Meye e al. and Wu e al. samples aligned o mouse mm9 genome assembly and eads om Håkelien e al. samples aligned o human hg19 genome assembly. All samples we e subjec ed o PCR duplica ed emo al, a e which Fas QC quali y con ol was applied again. Alignmen s a is ics o he Bow ie2 alignmen s a e shown in Supplemen al Table I. Ob ained alignmen s we e inspec ed isually using In eg a i e Genomics Viewe (IGV) (Robinson e al., 2011). 2.1.2. ChIP-seq peak de ec ion and anno a ion Peaks we e de ec ed om he alignmen files by MACS 1.4.2 so wa e (Zhang e al., 2008) using de aul se ings. Fo Håkelien e al. sample, he peak en ichmen was de e mined ela i e o a con ol ChIP wi h an i-H3 an ibody in day 28 iMSC3 cells. Fo Meye e al. samples, he peak en ichmen s we e de e mined ela i e o a con ol wi h IgG an ibody in day 0 and day 15 cells. Fo Wu e al. samples, he peak en ichmen was de e mined ela i e o sonica ed inpu DNA om day 9 MC3T3-E1 cells. Peaks wi h p- alue < 10 −10 we e used in all downs eam analyses. Peaks we e anno a ed o he closes p o ein coding genes using Bed ools closes (Quinlan and Hall, 2010) ool and gene anno a ions o he peaks lis s so ed by ascending p- alue we e compa ed. Because o Table I Summa y o he expe imen s and da a p o ided by he o iginal a icles. S udy Cell model ChIP-seq ime poin s Gene exp ession da a: me hod, ime poin s ENA s udy numbe Håkelien e al.: The Regula o y Landscape o Os eogenic Diffe en ia ion iMSC#3 Immo alized human mesenchymal s em cell line 28 d RNA-seq, 0 d cells 28 d cells ERP003787 Meye e al.: The RUNX2 Cis ome in Os eoblas s: cha ac e iza ion, down- egula ion ollowing diffe en ia ion, and ela ionship o gene exp ession MC3T3-E1 Mouse p eos eoblas ic cell line 0d 15 d Mic oa ay (Mouse 385 K mic oa ay, Roche Nimblegen) 0 d cells (POB) 15 d cells (OB) SRP016885 Wu e al.: Genomic Occupancy o Runx2 wi h Global Exp ession P ofiling Iden ifies a No el Dimension o Con ol o Os eoblas ogenesis MC3T3-E1 Mouse p eos eoblas ic cell line, subclone 4 0d 9d 28 d Mic oa ay (GeneChip Mouse Gene 1.0 ST A ay e .4, Affyme ix) ShRunx2 (9d) cells ShSc (0 d and 9 d) cells SRP035343 K. Ta kkonen e al. Gene 626 (2017) 119–131 120 highe de ia ion o Meye e al. 15d samples, he eplica e wi h highes FRiP (F ac ion o Reads in Peaks) sco es and la ges amoun o high confidence peaks we e used o compa isons wi h o he samples. Fo c oss-species sample compa ison, human gene symbols we e ansla ed o o hologous mouse gene symbols. Fo mouse samples om Meye e al. and Wu e al. expe imen s, he o e lap o he MACS de ec ed peak egions was also compa ed di ec ly using R package ChIPSeeke (Yu e al., 2015), which calcula ed s a is ical significance o he peak o e lap. The genomic loca ions o Runx2 binding peaks we e examined using a peak anno a ion unc ion o HOMER (Hype geome ic Op imiza ion o Mo i EnRichmen ) so wa e (Heinz e al., 2010). Simila analysis was also pe o med using Genomic Regions En ichmen o Anno a ions Tool (GREAT) (McLean e al., 2010). The analysis was done using de aul se ings ha assign a basal egula o y domain o −5 kb and +1 kb o he ansc ip ional s a si e (TSS) and ex end i in bo h di ec ions o he nea es gene's basal egula o y domain bu no > 1000 kb dis ance. The ChIP-seq peak egions we e hen associa ed wi h he genes in whose egula o y domains hey laid. GREAT was also used o find GO anno a ions ha a e en iched among he genes nea he ChIP- seq peak egions. Fo GO en ichmen analysis, he assigned gene egula o y domain ex ended in bo h di ec ions o he midpoin be ween he gene's TSS and he nea es gene's TSS bu no > 50 kb dis ance. 2.1.3. Runx2 binding mo i scan En ichmen analysis o conse ed Runx2 binding sequence mo i s was pe o med using HOMER so wa e (Heinz e al., 2010). The analysis was done using de aul se ings and epea -masked sequence. 2.2. DNA mic oa ay da a analysis 2.2.1. P ocessing o Meye e al. mic oa ay da a Meye e al. p o ided gene exp ession da a ob ained om mouse 385 k mic oa ays (Roche NimbleGen) o p e-os eoblas s age MC3T3- E1 cells (POB) and os eoblas s age MC3T3-E1 cells (OB) a e 15 days o diffe en ia ion. The da a consis ed o no malized gene exp ession alues wi h and wi hou log2 ans o ma ion as well as log2 old- change alues and was used in his analysis as p o ided. Fo some genes he e we e mo e han one old change alue in he able p o ided by Meye e al., which we e appa en ly ob ained om mul iple p obes de ec ing he same gene. These alues we e a e aged by aking a geome ic mean o he old change and he esul ing single old change alue o each gene was log2- ans o med. 2.2.2. P ocessing o Wu e al. mic oa ay da a Wu e al. p o ided gene exp ession da a ob ained om GeneChip Mouse Gene 1.0 ST A ay e .4 (Affyme ix) as RMA no malized alues in he Gene Exp ession Omnibus (GEO) da abase (accession numbe GSE53982). The no malized, log2- ans o med gene exp ession alues we e used o u he analyses. The samples used o his analysis we e day 0 MC3T3-E1 cells in ec ed wi h con ol Sc amble-shRNA and he same cells a e 9 days o diffe en ia ion. Diffe en ial exp ession be ween he day 0 and day 9 Sc -shRNA cells was de e mined using R Bioconduc o package limma (Ri chie e al., 2015) wi h ob ained p- alues adjus ed o mul iple es ing by Benjamini-Hochbe g me hod (Benjamini and Hochbe g, 1995). Diffe en ial exp ession was defined as significan i he adjus ed p- alue was below he significance le el 0.05 and absolu e log2 old change was a leas 1. 2.2.3. P ocessing o Håkelien e al. RNA-seq da a Fo Håkelien e al. s udy, aw RNA-seq ead files o day 0 undiffe en ia ed iMSC3 cells and day 28 diffe en ia ed iMSC3 cells ( wo eplica es o each) we e downloaded om Eu opean Nucleo ide A chi e (ENA) (accession numbe ERA246830). Quali y con ol o he eads was pe o med using Fas QC so wa e. The quali y o all ead files was ound o be easonable and no p e-p ocessing was conside ed necessa y. The eads we e aligned agains human hg19 genome using Bow ie2. Bed ools mul ico ool (Quinlan and Hall, 2010) was used o coun ing he eads aligned o all exonic egions in he genome. The coun s om all exons belonging o he same gene we e hen summa - ized o ge gene-le el ead coun s, and he ob ained coun s we e FPKM (F agmen s Pe Kilobase o exon pe Million agmen s mapped) no malized. The non-no malized ead coun s we e used o diffe en ial exp ession analysis by R Bioconduc o package DESeq2 ( e sion 1.6.3) (Lo e e al., 2014). The ob ained p- alues we e adjus ed o mul iple es ing by Benjamini-Hochbe g me hod. Diffe en ial exp ession was defined as significan i he adjus ed p- alue was below he significance le el 0.05 and absolu e log2 old change was a leas 1. 2.3. Gene exp ession da a compa ison Fo compa ison o he gene exp ession p ofiles om Meye e al., Wu e al. and Håkelien e al. s udies, he gene symbols in Håkelien e al. da a we e mapped o o hologous mouse gene symbols. Log2- ans- o med exp ession alues o all samples, scaled o he same ange, we e hie a chically clus e ed based on compu ed Euclidian dis ance o he samples wi h comple e linkage me hod. Fu he compa ison o he gene exp ession p ofiles be ween he h ee s udies u ilizing diffe en pla o ms and wo diffe en species was done using R Bioconduc o package O de edLis (Yang e al., 2006) wi h unc ion ha de ec s simila i ies be ween wo o de ed gene lis s. The unc ion compa es wo- anked lis , in his case gene lis s o de ed by dec easing gene exp ession alue, and a simila i y sco e is assigned based on he numbe o o e lapping genes in he op anks. Random sco es a e compu ed by compa ing one lis o he andomly shuffled second lis , and based on he andom sco es an empi ical p- alue can be compu ed o he obse ed sco e. 5000 pe mu a ions we e used he e o es ima ing he empi ical p- alues. Pa hway en ichmen analyses o he h ee gene exp ession da a se s we e pe o med by GSEAP e anked, which uns Gene Se En ichmen Analysis (GSEA) (Sub amanian e al., 2005) agains a use -supplied, anked lis o genes. Fo he analysis, he gene exp ession lis s o all samples con aining all he genes we e o de ed based on he log2 exp ession alues in descending o de and mouse gene symbols in Meye and Wu gene lis s we e mapped o o hologous human gene symbols. En ichmen was es ed o all a ailable KEGG pa hway gene se s wi h numbe o pe mu a ions se o 1000 and en ichmen s a is ic se as classic. Fo compa ison o he pa hway en ichmen esul s, he esul s we e fil e ed using a s ingen alse disco e y a e (FDR) cu off o 0.01 and o de ed by ascending FDR. 2.4. Compa ison o gene exp ession o Runx2 occupancy a gene p omo e s The ela ionship be ween Runx2 binding o gene p omo e egions and gene exp ession was s udied using he undiffe en ia ed and diffe en ia ed mouse samples om Meye and Wu s udies. Gene exp ession lis s o de ed by descending log2 exp ession alues we e subse in o ou g oups: (1) Genes wi h exp ession le el wi hin he highes 20% in undiffe en ia ed cells, (2) genes wi h exp ession le el wi hin he lowes 20% in undiffe en ia ed cells, (3) genes wi h exp ession le el wi hin he highes 20% in diffe en ia ed cells and (4) genes wi h exp ession le el wi hin he lowes 20% in diffe en ia ed cells. The numbe o Runx2 ChIP-seq peaks on gene p omo e a eas in samples a diffe en ime poin s (days 0 and 15 om Meye s udy and days 0, 9 and 28 om Wu s udy) we e analyzed and he peak dis ibu ions be ween he ou gene g oups we e compa ed. 3. Resul s and discussion 3.1. Gene al compa ison o he expe imen al se up be ween he s udies The expe imen al se up o os eoblas diffe en ia ion cul u e was K. Ta kkonen e al. Gene 626 (2017) 119–131 121 highly simila in all h ee Runx2 ChIP-seq s udies analyzed he e, wi h some a ia ion in he model cell line, in he leng h o he cul u e and composi ion o he os eogenic medium (Table I). In he wo MC3T3-E1 s udies, diffe en clones o MC3T3-E1 cells we e used, Wu e al. using a la e isola ed apidly mine alizing subclone om he o iginal MC3T3- E1 cells (Wang e al., 1999). iMSC#3 line used in Håkelien e al. s udy was p oduced by immo alizing human bone ma ow de i ed s omal cells by e o i al ansduc ion o elome ase e e se ansc ip ase (TERT). These cells can diffe en ia e o os eoblas s and adipocy es (Ska n e al., 2014). Os eogenic medium con ained asco bic acid and Na-β-glyce ophospha e in all s udies, bu dexame hasone was used only in he cul u es o iMSCs. Os eoblas ma u a ion was e ified simila ly in all s udies by s anda d me hods o ALP, aliza in ed and Von Kossa s aining o demons a e lineage commi men and bone nodule o ma- ion, espec i ely. Fo Runx2 ChIP-seq, in bo h MC3T3-E1 s udies he fi s cell s age chosen o analysis we e he undiffe en ia ed confluen cells, and he endpoin o sample collec ion was ei he 15d (Meye e al.) o 28d (Wu e al.) a e diffe en ia ion induc ion. In addi ion, Wu e al. included a imepoin ep esen ing ma ix-p oducing cells p io o mine aliza ion (9d). In iMSCs, only he 28d diffe en ia ed ma u e os eoblas s we e analyzed o Runx2 binding. Two diffe en Runx2 an ibodies we e used in he s udies discussed he e, a abbi polyclonal an ibody o MC3T3- E1 cells, and a mouse monoclonal an ibody o iMSCs. No ably, he epi ope o bo h hese an ibodies is in he e y same egion o Runx2 p o ein. Despi e s anda dized p o ocols o some cell ypes, ch oma in immunop ecipi a ion (ChIP) is a challenging me hod equi ing ex en- si e op imiza ion o sample and an ibody-specific condi ions including e.g. c osslinking ime and empe a u e, cell lysis, ch oma in shea ing and washing condi ions. Fu he mo e, special equi emen s conce ning inpu DNA size and quali y ha e o be me in o de o p oduce good quali y sequencing lib a ies o ChIP-seq. The la e-s age os eoblas ic cells a e igh ly su ounded by he collagenous ma ix and hus a e difficul sample ma e ial o ChIP. In he s udies discussed he e, he e was some a iabili y in he ChIP p ocedu es including e.g. cell lysis and sonica ion condi ions, which migh ha e an impac on he sample quali y and esul s, bu hese aspec s a e no in he ocus in ou s udy. The me hod and expe imen al se up o ansc ip ional p ofiling o os eoblas ic cells as well as app oaches o co ela e Runx2 binding si e da a o ansc ip ome diffe ed be ween he s udies. Meye e al. pe o med mic oa ay exp ession analysis o MC3T3-E1 a undiffe en- ia ed and diffe en ia ed s a e while Håkelien e al. used RNA-seq analysis o iMSC#3 cells o undiffe en ia ed and 28d diffe en ia ed cells. Wu e al. in u n did a mic oa ay gene exp ession p ofiling o 9d diffe en ia ed MC3T3-E1 cells wi h shRNA media ed knockdown o Runx2, and compa ed he changes in gene exp ession o con ol shRNA cells. In o de o compa e gene exp ession p ofiles be ween he s udies, only he con ol sc ambled shRNA da a o Wu e al. s udy was analyzed he e, co esponding o he pa en al MC3T3-E1 cells o Meye e al. s udy. 3.2. ChIP-seq analysis 3.2.1. Quali y o he sequencing da a Raw ChIP-seq sequencing eads we e downloaded om ENA and subjec ed o quali y con ol. In summa y, Wu e al. and Håkelien e al. samples showed e y high quali y h oughou he eads. Meye e al. samples had no iceable dec ease in quali y owa ds he end o eads, especially in 15d eplica e sample 1, bu as he majo i y o bases in eads we e o e y good quali y and he median quali y sco es we e wi hin he e y good quali y ange, no u he p e-p ocessing was conside ed necessa y. 3.2.2. Compa ison o peak de ec ion esul s A e duplica e emo al he aligned eads we e used o peak de ec ion using MACS (Model-Based Analysis o ChIP-seq) (Zhang e al., 2008). Peak de ec ion esul s o all samples a e shown in Table II. In summa y, la ge amoun s o Runx2 peaks we e de ec ed in all samples e en i mo e s ingen p- alue cu offo 10 −10 was used. No iceable is he educ ion in peak coun in Meye day 15 eplica e 1 sample a e mo e s ingen cu offwas applied, in acco dance wi h he lowe quali y o eads ha was obse ed in quali y con ol. A simple quali y me ics F IP o ChIP-seq da a is also shown in Table II. The ecommenda ion o his me ics is ≥1% by ENCODE conso ium o a success ul expe imen (Land e al., 2012), and he sco es o all samples exceed his limi . Peaks wi h p- alue < 10 −10 we e used in all down- s eam analyses. In he published pape Meye e al. epo ed con ac ion o he Runx2 cis ome du ing diffe en ia ion, demons a ed by 12,674 and 6272 genomic egions wi h Runx2 binding in he ea ly and la e s a e os eoblas s, espec i ely (Meye e al., 2014b). Wu e al. in u n epo ed app oxima ely 25,000 significan ly en iched egions a day 0 and an inc ease in he numbe o Runx2 bound si es in la e imepoin s, yielding up o 80,000 Runx2 en iched egions in MC3T3-E1 cells (Wu e al., 2014). In he e-analysis pe o med he e, we did no obse e as high a ia ion in he de ec ed peak numbe s be ween he wo MC3T3- E1 s udies as epo ed p e iously, sugges ing he diffe ence in he peak numbe s be ween he o iginal s udies o be ela ed o he da a analysis a he han o a majo diffe ence in Runx2 binding abundance (Table II). In ou analysis, Wu e al. samples showed inc eased numbe o Runx2 peaks in D9 and D28 samples compa ed o D0 sample, whe eas peak numbe in Meye e al. s udy s ayed cons an be ween he imepoin s (when excluding he one D15 eplica e o lowe quali y) (Table II). In he o iginal pape s, he e was a mino diffe ence in Runx2 p o ein exp ession pa e n, he p o ein le el being cons an in Meye 's pape bu showing inc ease in Wu's pape , which is in acco dance wi h he dynamics in Runx2 peak numbe s obse ed. In iMSC#3 cells, we ound la ge numbe o Runx2 peaks (80061) and e y high FRiP alue (48,7%), sugges ing a success ul ChIP expe imen . Howe e , Håkelien e al. epo ed only 9549 peaks in he 50 kb egion o closes 5′gene end, sugges ing majo diffe ences in he analysis me hods and/o fil e ing pa ame e s used. Thus, we will no u he compa e ou iMSC analysis esul s o he published esul s, bu use ou de no o analysis esul s o he analysis o mouse MC3T3-E1 cells in o de o pe o m in e -species compa isons. 3.2.3. Compa ison o anno a ed genes close o Runx2 peaks Peaks we e anno a ed o he closes p o ein coding genes using bed ools closes ool (Quinlan and Hall, 2010). Anno a ed peak lis s we e so ed by ascending p- alue so ha he mos en iched peaks o Table II ChIP-seq peak de ec ion esul s. Meye e al. Wu e al. Håkelien e al. MACS de ec ion D0_1 D0_2 D15_1 D15_2 D0 D9 D28 D28 Peak coun Cu offp<10–5 38,211 34,296 46,249 34,296 23,667 96,368 49,425 122,325 Cu offp<10–10 26,297 22,966 9810 22,966 14,396 53,887 32,894 80,061 FRiP (F ac ion o Reads in Peaks) 17.3% 13.6% 4.3% 11.3% 4.1% 21.0% 13.0% 48.7% K. Ta kkonen e al. Gene 626 (2017) 119–131 122 each sample could be compa ed. To allow c oss-species compa isons, he human gene symbols we e ansla ed o o hologous mouse gene symbols. Venn diag ams o he unique and common genes anno a ed o he op 5000 mos en iched peaks a e shown in Fig. 1. Meye e al. and Wu e al. undiffe en ia ed mouse samples seem o be ai ly simila , wi h 70% o op 5000 peak anno a ions sha ed be ween he samples (Fig. 1A), whe eas diffe en ia ed samples (Meye e al. day 15 and Wu e al. day 28) show a smalle , 59%, o e lap (Fig. 1B). Compa ison o Meye e al. undiffe en ia ed and diffe en ia ed samples shows an o e lap o 74% in anno a ions, wi h co esponding compa isons be ween Wu samples showing 68% and 71% o e laps o day 9 and day 28 samples, espec i ely (Supplemen a y Fig. 1). In conclusion, he compa isons be ween he diffe en ia ion s ages sugges ha Runx2 cons i u i ely occupies mo e han hal o he o al Runx2 bound genes ega dless o he diffe en ia ion s age. Compa ison o peak anno a ions o human diffe en ia ed iMSC3 samples o diffe en ia ed mouse MC3T3-E1 samples om bo h Meye and Wu shows o e laps o 39% and 43% in op 5000 anno a ed genes, espec i ely (Fig. 1C,D). The o e lap pe cen ages a e, as could be expec ed, lowe han hose o wi hin-species compa isons and in line wi h p e ious findings ha ha e shown ha mos TF binding e en s a e species-specific(Wilson and Odom, 2009). Howe e , al oge he 1344 genes om op 5000 anno a ed genes we e in common o all samples in he las diffe en ia ion ime poin , sugges ing subs an ial simila i y in Runx2 peak associa ed genes be ween he species and expe imen s. This could sugges ha hese genes a e especially impo an o he Runx2- media ed egula ion o os eoblas diffe en ia ion. Fo Meye and Wu mouse samples i was possible o compa e di ec ly he o e lap o he peak egions. The esul s indica e ha peaks in undiffe en ia ed and diffe en ia ed mouse samples ha e a highly significan o e lap bo h wi hin and be ween he s udies (Table III). Fo example, when compa ing undiffe en ia ed cells o diffe en ia ed cells in each s udy, 88% (Wu e al. s udy) and 67% (Meye e al. s udy) o he peak egions a day 0 we e p esen also in diffe en ia ed cells. Mo e- o e , a day 0, 88% o peaks obse ed in Wu e al. s udy a e de ec ed also in Meye e al. day 0 sample. In he diffe en ia ed cells, co esponding o e lap was 50%, possibly ela ed o he a iable ime cou se and he diffe en MC3T3-E1 cell clones used in he wo s udies as well as he challenging sample ma e ial a he la e ime poin s. Impo an ly, hese esul s indica e ha ChIP-seq expe imen s in MC3T3-E1 cells show ela i ely high ep oducibili y in wo indi idual s udies and again ha majo i y o he Runx2 bound si es a e cons an du ing diffe en ia ion. 3.2.4. Compa ison o genomic loca ions o Runx2 peaks The genomic loca ions o Runx2 binding peaks we e examined using a peak anno a ion unc ion o HOMER so wa e. The dis ibu ions o Runx2 peaks a e displayed in Fig. 2A. In acco dance wi h he o iginal a icles, he as majo i y o Runx2 binding occu ed a in e genic and in onic egions. The g ea es a ia ion be ween he diffe en ia ion s ages was in Runx2 occupancy a p omo e s, wi h highe occupancy in Fig. 1. Venn diag ams o peak anno a ion compa isons. Gene anno a ions o he op 5000 mos en iched Runx2 ChIP-seq peaks we e used o compa isons: A–B) Undiffe en ia ed and diffe en ia ed MC3T3-E1 cells o Meye e al. and Wu e al. Samples we e compa ed, and C–D) Diffe en ia ed Håkelien e al. iMSC3 cells we e compa ed o Meye e al. and Wu e al. diffe en ia ed MC3T3-E1 cells. Table III S a is ical es o ChIP-seq peak o e lap in mouse samples. Que y sample Ta ge sample Que y peak coun Ta ge peak coun O e lapping peaks p- alue Adjus ed p- alue Meye DO_1 Meye D15_2 26,297 34,871 17,569 0 0 Wu D0 Wu D9 14,396 53,887 13,255 0 0 Wu D0 Wu D28 14,396 32,894 12,714 0 0 Meye D0_1 Wu D0 26,297 14,396 12,609 0 0 Meye D15- _2 Wu D28 34,871 32,894 16,577 0 0 K. Ta kkonen e al. Gene 626 (2017) 119–131 123 Fig. 2. Genome-wide p ofile o Runx2 occupancy. A–B) Dis ibu ion o Runx2 binding peaks ac oss mouse and human genomes we e classified in o eigh ca ego ies: Exon, In on, P omo e (−1 kb o +100 bp om ansc ip ion s a si e), TTS egion (−100 bp o +1 kb om ansc ip ion e mina ion si e), 5′UTR exon, 3′UTR exon, In e genic egion and o he egions. Peak dis ibu ions in all samples a e plo ed by peak numbe (A) and by he pe cen age o o al peaks (B). C–E) Analysis o he Runx2 peak associa ion o genomic egions using GREAT a he las ime poin in each o he h ee s udies. In he analysis, a basal gene egula o y domain o −5 kb and 1 kb o he TSS was assigned and ex ended in bo h di ec ions o he nea es gene's basal egula o y domain wi h a maximum dis ance o 1000 kb. Runx2 peak egions whe e hen associa ed wi h he genes in whose egula o y domains hey a e loca ed. In addi ion, GO Biological P ocess e m en ichmen analysis o he las ime poin samples in each s udy using GREAT is p esen ed ( igh panels). In his analysis, he assigned gene egula o y domain ex ended in bo h di ec ions o he midpoin be ween he gene's TSS and he nea es gene's TSS wi h a maximum dis ance o 50 kb. The figu e shows he mos en iched GO e ms among he genes nea he ChIP-seq peak egions. K. Ta kkonen e al. Gene 626 (2017) 119–131 124 undiffe en ia ed cells in bo h MC3T3-E1 s udies (Fig. 2A,B, Table IV). Simila analysis was pe o med using GREAT so wa e, which finds GO anno a ions ha a e en iched among he genes nea he ChIP-seq peak egions. G aphs showing he dis ance be ween peak egions and hei pu a i ely egula ed genes as well as he mos en iched biological p ocess GO e ms ob ained om GREAT analysis o diffe en ia ed cells a e shown in Fig. 2C–E. In summa y, he pa e n o Runx2 peak en ichmen was e y simila be ween he expe imen s, showing clea en ichmen o Runx2 binding in he icini y o TSS. Among en iched GO e ms in biological p ocesses ca ego y he e we e bo h common and closely ela ed e ms and some unique pa hways, eflec ing he o e lap o he anno a ed genes close o he peaks (Fig. 1). The en iched e ms such as RNA ela ed p ocesses, cell- and issue mo phology ela ed pa hways, esembled closely he published da a o each epo . The mino diffe ences in he e m en ichmen s o he o iginal a icles esul s migh be due o diffe ences in pe o ming he analysis, o example Wu e al. used ime-dependen dynamic clus e s o Runx2 peaks, whe eas in ou analysis all Runx2 peaks a ce ain ime poin we e included. 3.2.5. Runx2 binding mo i en ichmen in Runx2 peak egions ChIP-seq peaks we e scanned o conse ed TF sequence mo i s using HOMER so wa e. HOMER includes a mo i da abase ha is mos ly based on he analysis o public ChIP-Seq da a se s (Heinz e al., 2010). The ull ables o he esul s o he mo i scan a e p esen ed upon eques . In mouse samples, he mos significan ly en iched mo i s we e h ee mo i s ha ha e been expe imen ally associa ed wi h binding o RUNX- amily o ansc ip ion ac o s ((A/T/C)TGTGGTT(A/T); (G/T) (T/C)TGTGGTTT; CTGTGGTTT(G/C))), p esen ed in Table V. These mo i s we e p esen in al oge he 36–54% o Runx2 peaks in all MC3T3- E1 samples (Table VI). In he mo i scan, wo misma ches o he consensus sequence we e allowed in he analysis, leading o inclusion o high numbe o diffe en a ian mo i s as Runx2 mo i s, e e ed he ea e al oge he as Runx2 mo i . Acco dingly, significan en ich- men o a classic Runx2 binding co e mo i TGTGGT a he Runx2- bound egions was epo ed in bo h MC3T3-E1 s udies p e iously. In iMSCs on he o he hand, he known Runx2 mo i s we e ound in only 14,4% o he peaks and he mo i s anked as mos significan ly en iched we e o AP1 and ETS amily ansc ip ion ac o s. Thus, in iMSC sample, al hough ha ing he highes numbe o peaks le a e quali y fil e ing, only mino i y o hese peaks appea o con ain Runx2 binding mo i s. This may indica e ha he Håkelien e al. sample con ains mo e alse posi i es and would equi e mo e s ingen fil e ing o he de ec ed peaks, o ha he Runx2 mo i in human DNA is highly a iable and was no ecognized in ou analysis. In e es ingly, he mo i scan e ealed also o he significan ly en iched TF consensus mo i s a Runx2 peak egions. Seconda y binding mo i s may sugges o example p o ein-p o ein in e ac ions, by which Runx2 is ec ui ed o si es wi hou Runx2 binding mo i , o possible co-ope a i e binding and/o collabo a ion o Runx2 wi h o he ansc ip ion ac o s a he same egion. When compa ing he op 20 en iched mo i s in each sample, he e we e 9 common mo i s ound in all samples (Table VII). Mos o hese we e ela ed o AP-1 amily p o eins including A 3, F a1 and Jun-AP1 mo i . These consensus mo i s a e ound a high equency in he genome and a e in ol ed in a wide a ie y o cellula p ocesses. In e es ingly, ano he significan ly en iched mo i in e e y da a se was o TEAD amily ac o s and especially o TEAD4, p esen in 9–17% o Runx2 peaks in mouse and in 10% o peaks o human samples. In Håkelien e al. o iginal a icle, TEAD2 mo i en ichmen was epo ed in egions en iched o p omo e and enhance ela ed his one modifica ions H3K4me3, H3K9ac and H3K27ac a he end o diffe en ia ion, and i s knockdown was hen u he shown o cause impai ed mine aliza ion in i o. TEAD amily o ansc ip ion ac o s a e effec o s in he Hippo signaling pa hway egula ing o gan size by con olling cell p oli e a ion and apop osis (Landin-Mal e al., 2016). TEADs o m a complex wi h Yes-associa ed p o ein (YAP) o i s pa alog ansc ip ional co-ac i a o wi h PDZ binding mo i (TAZ) o ac i a e a ge gene exp ession. Bo h YAP and TAZ in u n in e ac also wi h Runx2 (Hong e al., 2005; Zaidi e al., 2004) and hey ha e been implica ed in egula ing MSC diffe en ia ion (Hong e al., 2005) and in Wn and BMP signaling (Va elas, 2014). Ve y ecen ly, s abiliza ion o TAZ-Runx2 complex was shown o play an impo an ole in p omo ing os eoblas ogenesis (Ma sumo o e al., 2016). When inspec ing Runx2 peaks wi h he Table IV Genome-wide Runx2 peak dis ibu ions by he pe cen age o o al peaks. Meye Wu Håkelien D0_1 D0_2 D15_1 D15_2 D0 D9 D28 D28 3′UTR 0.6% 0.6% 0.6% 0.5% 0.5% 1.0% 0.8% 0.8% TTS 1.4% 1.4% 1.0% 0.9% 1.5% 1.8% 1.5% 1.2% Exon 2.2% 2.2% 0.8% 0.9% 3.4% 5.5% 4.3% 2.5% In on 36.5% 36.3% 32.9% 26.2% 30.7% 40.5% 36.2% 49.6% In e genic 38.1% 37.5% 52.2% 58.6% 32.0% 35.4% 33.1% 36.3% P omo e 19.1% 19.9% 11.8% 12.1% 28.6% 13.3% 20.4% 8.0% 5′UTR 1.7% 1.9% 0.5% 0.7% 3.0% 2.3% 3.3% 1.1% O he 0.3% 0.3% 0.3% 0.2% 0.3% 0.3% 0.4% 0.5% To al 100.0% 100.0% 100.0% 100.0% 100.0% 100.0% 100.0% 100.0% Table V Runx2 mo i s en iched in HOMER mo i en ichmen analysis* o MC3T3-E1 cells. Mo i name Consensus sequence RUNX2(Run )/PCa-RUNX2-ChIP-Seq(GSE33889)/ Home A/T/C)TGTGGTT(A/T) RUNX1(Run )/Ju ka -RUNX1-ChIP-Seq(GSE29180)/ Home (G/T)(T/C)TGTGGTTT RUNX(Run )/HPC7-Runx1-ChIP-Seq(GSE22178)/ Home CTGTGGTTT(G/C) *The analysis allowed o wo misma ch nucleo ides in he consensus sequence. Table VI Pe cen ages o MACS peak sequences wi h Runx- ela ed mo i s de ec ed by HOMER analysis. Peaks wi h RUNX mo i Sample To al numbe o peaks Peaks wi hou RUNX mo i Numbe Pe cen age Meye D0_1 26,297 16,792 9505 36.1 Meye D0_2 22,962 14,304 8658 37.7 Meye D15_1 9810 4514 5296 54.0 Meye D15_2 34,871 20,813 14,058 40.3 Wu D0 14,396 9219 5177 36.0 Wu D9 53,887 31,310 22,577 41.9 Wu D28 32,894 20,088 12,806 38.9 Håkelien D28 80,061 68,603 11,458 14.3 K. Ta kkonen e al. Gene 626 (2017) 119–131 125 TEAD4 mo i mo e closely, app oxima ely hal o he TEAD4 mo i con aining peaks con ained also Runx2 mo i in mouse samples, comp ising 4–10% o all Runx2 peaks, whe eas in human iMSCs only 1,5% o o al Runx2 peaks con ained bo h mo i s (Supplemen al Table II). These esul s sugges ha TEAD4 and Runx2 in e ac ion migh ake place a leas in pa o he Runx2 occupied egions h ough a common p o ein complex (TEAD4 mo i only egions) o h ough collabo a i e binding ia pu a i e nea by DNA binding mo i s (co-p esence o TEAD4 and Runx2 mo i s). Gi en he eme ging li e a u e o he impo ance o TEAD/TAZ/YAP p o eins in os eogenesis (Ma sumo o e al., 2016; Tang e al., 2016), ole o TEAD-Runx2 in e ac ion and p o ein complex in os eoblas s would be highly in e es ing a ge o u u e s udies. 3.2.6. Runx2 occupancy a he p omo e egions Runx2 binding si es whe e inspec ed in mo e de ail a p omo e egions, which a e mo e likely o be conse ed ac oss species han he dis al enhance si es (Cheng e al., 2014). The p omo e sequences (2 kb ups eam om TSS) o all o hologous pai s o human and mouse genes we e e ie ed and anno a ed o he nea es ChIP-seq peaks in mouse and human samples. The esul s o he ex ensi e analysis a e p o ided in Supplemen al Table III. The p omo e coo dina es, nea es peak egion's coo dina es, and he dis ance o he peak om he p omo e a e shown o all genes oge he wi h in o ma ion indica ing he numbe o samples, in which Runx2 peaks we e ound. In addi ion, he numbe o consensus Runx2 mo i s ound a each peak egion is epo ed in he able. I he peak dis ance om he p omo e egion is anno a ed as 0, i means ha he peak is o e lapping wi h he p omo e . Gene exp ession da a om he s udies (desc ibed in de ail in Sec ion 3.3.) ha e been in eg a ed in o he able as well. Ano he iew o he da a is p esen ed in Supplemen al Table IV, which shows he p esence o Runx2 a gene p omo e s a diffe en ime poin s oge he wi h he ela ed gene exp ession da a, offe ing possibili y o assess he dynamics o Runx2 binding a indi idual gene p omo e s be ween he s udies. These ables wi h anno a ed gene in o ma ion can be used as a esou ce, when explo ing possible Runx2 binding in he icini y o any gene o in e es . Fo example, by fil e ing and o ganizing columns in Tables SIII and SIV, a ious g oups o genes o in e es can be selec ed based on hei Runx2 binding p ope ies and exp ession p ofiles. We nex sea ched o genes ha had Runx2 peaks nea he p omo e egions in all 7 mouse samples indica ing high confidence o cons i u- i e Runx2 binding du ing diffe en ia ion using he Supplemen al Table IV. These genes and hei p omo e s had se e al common ea u es. Fi s , mos o hese Runx2-occupied genes we e ela i ely highly exp essed bu only mino i y (11%) o hem con ained Runx2 mo i s in he peak egions. Simila ly, he co esponding human o hologues o hese genes had nea by Runx2 peaks, we e exp essed and showed low equency o Runx2 consensus mo i s. These obse a ions sugges ei he Runx2 mo i is highly a ian on hese si es, o ha Runx2 may be ec ui ed o many egions by o he mechanisms han by di ec binding o Runx2 mo i s in he DNA. Whe he Runx2 binding is equi ed o he ansc ip ion o hese genes is no clea . In e es ingly, in he s udy o Wu e al. only 159 genes esponded o Runx2 knockdown du ing diffe en ia ion, which is a less han genes occupied by Runx2. As Runx2 has also genome- o ganizing capabili ies (Lian e al., 2003), i may well be ha on hese o he genes Runx2 modifies he genomic landscape a he han di ec ly egula es gene exp ession. Ne e heless, Wu e al. epo ed mo e Runx2 peaks in he shRunx2 down egula ed genes ( hus Runx2 up egula ed genes) compa ed o shRunx2 up egula ed o non esponsi e genes. Acco dingly, when we inspec ed he abundance o Runx2 peaks and Runx2 mo i s close o he shRunx2 esponsi e/un esponsi e genes (Table VIII), he shRunx2 down egula ed genes con ained indeed mo e Runx2 peaks pe gene, and especially mo e Runx2 mo i con aining peaks pe gene (a e age 3,1/gene) han he up egula ed o non esponsi e genes (a e age 1,1/ gene), suppo ing he high equency o Runx2 mo i s oge he wi h Runx2 binding o be an s ong indica o o di ec ansc ip ional ac i a ion by Runx2. 3.3. Gene exp ession da a The gene exp ession analysis was pe o med by using diffe en pla o ms and expe imen al se ups as desc ibed in Sec ion 3.1., esul - ing in da a p o ided in diffe en o ma s in on-line eposi o ies. Fo he compa isons, all gene exp ession da a se s we e p ocessed o p oduce no malized log2- ans o med exp ession alues o he exp essed genes. 3.3.1. Diffe en ial gene exp ession analysis Meye e al. p o ided a eady- o-use able o gene exp ession p ofiles ob ained om mouse DNA mic oa ay analysis o p e-os eo- blas s age MC3T3-E1 cells (POB) and os eoblas s age MC3T3-E1 cells (OB) a e 15 days o diffe en ia ion. On he lis , he e we e 721 ansc ip s including 498 significan ly diffe en ially exp essed genes in he Meye e al. s udy be ween he diffe en ia ion s a es. Wu e al. in u n pe o med exp ession p ofiling a 0 and 9 days o diffe en ia ion o MC3T3-E1 ea ed wi h ei he con ol sc ambled shRNA o Runx2 a ge ing shRNAs leading o Runx2 down egula ion. In he cu en s udy, in o de o compa e gene exp ession p ofiles o MC3T3-E1 cells Table VII Top 20 en iched TF mo i s in all Runx2 ChIP-seq samples. Meye Wu Håkelien 0 day 15 day 0 day 9 day 28 day 28 day RUNX RUNX RUNX RUNX RUNX A 3 RUNX1 RUNX1 RUNX1 RUNX1 RUNX1 BATF RUNX2 RUNX2 RUNX2 RUNX2 RUNX2 F a1 RUNX-AML RUNX-AML RUNX-AML RUNX-AML RUNX-AML AP-1 TEAD4 TEAD4 F a1 F a1 F a1 Fosl2 TEAD TEAD BATF A 3 BATF Jun-AP1 F a1 TEAD2 A 3 CTCF A 3 Bach2 Fosl2 A 3 Fosl2 BATF Fosl2 NF-E2 Jun-AP1 F a1 Jun-AP1 Fosl2 AP-1 Bach1 BATF BATF AP-1 AP-1 Jun-AP1 N 2 TEAD2 AP-1 REST-NRSF Jun-AP1 CTCF Ma K A 3 Fosl2 TEAD4 BORIS Bach2 ERG AP-1 Jun-AP1 Bach2 Bach2 REST-NRSF ETV1 CTCF CTCF TEAD NF1 TEAD4 ETS1 Bach2 Ap4 TEAD2 TEAD4 TEAD TEAD4 NF1-hal si e BORIS Fli1 Kl 4 Fli1 Fli1 BORIS Bach2 Elk4 KLF5 BORIS GABPA Elk1 A oh1 Elk1 EKLF TEAD2 TEAD ETS Tc 12 ELF1 TEAD Ap4 Ma A Elk4 MyoG ETS Tlx Elk1 TEAD2 Table VIII Numbe o Runx2 peaks and Runx2 consensus mo i s a he Runx2 peaks a he p omo e egions o shRunx2 esponsi e genes. A e age 200 bp egion a ound peak cen e 500 bp egion a ound peak cen e shRunx2 esponse Gene coun Anno a ed peaks Numbe o peaks/ gene Peaks wi h Runx mo i (s)/gene % o peaks wi h Runx mo i (s) Peaks wi h Runx mo i (s)/gene % o peaks wi h Runx mo i (s) Down- egula ed 44 288 6.5 3.1 47 4.3 65 Up- egula ed 115 296 2.6 1.1 44 1.6 61 Non- esponsi e 20,788 58,886 2.8 1.1 40 1.6 57 K. Ta kkonen e al. Gene 626 (2017) 119–131 126 be ween Meye and Wu s udies, we analyzed only he da a o con ol shRNA exp essing cells o Wu s udy and compa ed hem o he pa en al MC3T3-E1 cells used in Meye s s udy. Wu e al. p o ided no malized signal alues and p obe IDs. Ou analysis o diffe en ially exp essed genes fil e ed by adjus ed p- alue and he 2- old change yielded 345 diffe en ially exp essed (DE) genes be ween he ime poin s. Fo Håkelien e al. s udy, aw RNA-seq ead files o day 0 undiffe en ia ed iMSC3 cells and day 28 diffe en ia ed iMSC3 cells we e a ailable and we e analyzed as desc ibed in ma e ials and me hods. The analysis esul s we e fil e ed o adjus ed p- alue and he ex en o diffe en ial exp ession wi h 2- old change. The fil e ed gene lis con ained 3039 genes, 1504 up- egula ed and 1535 down- egula ed, he numbe s being in line wi h he 3157 diffe en ially exp essed genes epo ed in Håkelien e al. s udy. Top en diffe en ially exp essed genes in each s udy a e shown in Table IX. All DE gene lis s a e p o ided as Supplemen al Tables V–VII. 3.3.2. Hea map o gene exp ession p ofiles In o de o compa e gene exp ession p ofiles o all samples in Meye e al., Wu e al. and Håkelien e al. s udies, he gene symbols in Håkelien e al. da a we e ansla ed o o hologous mouse gene symbols. Log2- ans o med exp ession alues o all samples, scaled o he same ange and hie a chically clus e ed, a e isualized in a hea map in Fig. 3. The hea map clea ly shows ha samples coming om he same s udy a e mos simila o each o he based on hei gene exp ession p ofiles, ega dless o he diffe en ia ion s age o he cells. Thus, he effec s caused by diffe ences in he cell clones o o he echnical aspec s be ween he s udies mask he effec o diffe en ia ion on he o e all gene exp ession p ofiles e en in he samples om he same species. The finding o only mode a e changes in gene exp ession p ofiles be ween he diffe en ia ion s a es migh also eflec he ac ha MC3T3-E1 cells a e al eady commi ed o he os eogenic lineage and he mos s iking changes in gene exp ession om he s em cell s age o p eos eoblas s ha e al eady passed. Ne e heless, hese esul s demons a e ha he widely dis ibu ed cell lines may exhibi clonal shi esul ing in pheno ypic/epigene ic/ ansc ip omic changes o e ime, and highligh he impo ance o c i ical e alua ion o esul s ob ained om a single cell line. 3.3.3. Simila i y o gene exp ession p ofiles In o de o assess he simila i y o he gene exp ession p ofiles be ween he samples om diffe en s udies using diffe en pla o ms and e en om diffe en species, R Bioconduc o package O de edLis (Yang e al., 2006) was u ilized o de ec ing simila i ies o wo o de ed gene lis s. In sho , wo anked lis s (acco ding o he mRNA exp ession le el) we e compa ed and a simila i y sco e was assigned based on he numbe o o e lapping genes in he op anks. Plo s o o e lapping genes in Meye e al. and Wu e al. day 0 op o bo om anked gene exp ession lis s a e shown in Fig. 4A–B, whe e he o e lap size inc eases i he gene in op anks o lis 1 is ound wi hin he op anks o lis 2. The obse ed o e lap was compa ed o he expec ed o e lap de i ed om a hype geome ic dis ibu ion. The significance o he simila i y o he lis s wi hin op 1000 genes was e y high (Fig. 4B), wi h p- alue o 0. Wi hin op 100 genes he simila i y was lowe bu s ill highly significan (Fig. 4A), wi h p- alue 0,0026. Numbe o common genes in hese op 100 and op 1000 exp essed gene lis s was 15 and 271, espec i ely. Simila compa isons o undiffe en ia ed human os eo- Table IX Top 10 significan ly diffe en ially exp essed genes in each s udy. Meye e al. 15d s 0d Wu e al. 9d s 0d shSc amble Håkelien e al. 28d s 0d Gene name Log2FC Gene name Log2FC Gene name Log2FC Cd200 6.233142 Lum 6.024314816 MMP13 11.50375 Mmp13 5.055013 Apod 4.675722274 COL10A1 10.72343 I gbl1 4.405752 F13a1 3.826338637 DPT 10.20353 Lec 1 4.374315 Ibsp 3.764573159 OMD 10.10910 Col15a1 4.164052 Mme 3.375728736 CHI3L1 9.27727 Ig 2 4.03787 Omd 3.355162249 FAM20A 9.23509 Akap12 3.937885 Pgm5 3.349156029 BRINP1 9.21489 F13a1 3.900387 Edn a 3.213000637 SERPINF1 9.19554 Bmpe −4.09268 Ppbp -3.205871844 GGT5 9.07341 Cxcl7 −4.58329 Anxa8 -3.106959652 MMP7 9.02634 Fig. 3. Hea map o log2- ans o med, no malized gene exp ession alues o Håkelien e al., Wu e al., and Meye e al. samples. Dend og am on op o he figu e shows he esul o hie a chical clus e ing o he samples using Euclidian dis ances wi h comple e linkage me hod. K. Ta kkonen e al. Gene 626 (2017) 119–131 127