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