Variations in genome-wide gene expression in identical twins - a study of primary osteoblast-like culture from female twins discordant for osteoporosis.

BACKGROUND: Monozygotic twin pairs who are genetically identical would be potentially useful in gene expression study for specific traits as cases and controls, because there would be much less gene expression variation within pairs compared to two unrelated individuals. However the twin pair has to be discordant for the particular trait or phenotype excluding those resulting from known confounders. Such discordant monozygotic twin pairs are rare and very few studies have explored the potential usefulness of this approach.
RESULTS: We studied genome-wide gene expression in primary osteoblast-like culture from marrow aspirates obtained from three pairs of monozygotic twins. We used the latest Affymetrix microchip contains probe sets for more than 20,000 genes. Two pairs were discordant for bone mineral density at the hip by more than one standard deviation, and the third pair was unrelated concordant and used as control. Only 1.5% on average of genes showed variation in expression within pairs as compared to 5% between pairs or over 15% from the literature. Importantly we identified several groups of genes showing variations within the discordant pairs and not within the concordant pair such as chondroitin beta 1,4 N-acetylgalactosaminyltransferase, inhibin beta A, interleukin 1 beta and colony stimulating factor 1 macrophage. These genes are known to have potential roles in bone physiology relating to bone density, osteoporosis and osteoarthritis.
CONCLUSION: Using the example of osteoblast-like cells in our monozygotic discordant twins for osteoporosis, we identified genes showing differential expression. Although without further experiment, we cannot confirm or conclude these are genes definitely related to bone physiology, we believe we have shown the potential and cost-effectiveness of further gene expression studies in discordant monozygotic twin pairs. A replication study for confirmation is essential.

Background

Studies of gene expression in humans have been problematic due to wide individual variation between subjects and random noise [1]. Characteristically screening studies of genome-wide expression in cases and controls have produced in excess of 500 genes of interest to study further. These numbers are unfeasible to study in humans and better efficient methods are sought to narrow down the number of gene targets. One potential solution is using pairs of identical twins, discordant for specific complex traits or diseases, which has yet to be properly explored.

As an example of a complex trait we chose osteoporosis, a common condition resulting in fractures and affecting one in four women in developed countries [2]. One of the main risk factors for fracture is reduced bone mineral density (BMD), which is easily measured by Dual Energy X-ray Absorbtiometry (DXA). Twin studies have consistently found bone density to be strongly heritable with most estimates in the range of 70–80% [3], whilst only a small amount of variance was due to unique environmental influences.

There are only a handful cellular gene expression studies reported in monozygotic (MZ) twins [4-8] and just two of them using microchips for expression study [6,7]. One study has looked at genome-wide gene expression changes in lympoblastoid cell lines derived from unrelated individuals and MZ twins [8]. The main aim of this study was to reveal the normal expression variation of 5184 genes in unrelated individuals, and MZ twins were used to compare the variability of 5 selected highly variable genes, as expected MZ twins showed much less variation [8]. Recently, Munshi et al reported genome-wide gene expression differences between one pair of MZ discordant for multiple myeloma in peripheral blood mononuclear cells using the Affymetrix hu95av2 arrays, which covers 12,600 genes.

The aims of our study were therefore to explore the feasibility and value of using the identical twin model to uncover altered expression of key genes in the tissue of relevance to a common complex trait. We examined and compared the mRNA expression in primary osteoblast-like cells of three pairs of monozygotic female twins, two of which were discordant for bone density, and one pair concordant for bone density acted as a control.

Results

From 200 MZ pairs screened with DXA, ten were identified as discordant and 7 excluded for environmental confounders. Two MZ pairs were consented for study; their bone mineral densities (BMD) at the hip were different by at least 1SD with no obvious environmental cause. These 2 MZ twin pairs (2A/2B and 3A/3B with hip BMD T-score differed by 1.42 and 1.40 SDs respectively) were the case pairs, whereas the concordant pair (1A/1B differed by only 0.14 SD) was used as the control pair. For case pair 3A/3B, T-score also differed by 2.40 for spine BMD. For convenience of comparison, 1A, 2A and 3A have lower BMD. All three pairs of MZ twins were of the same age (all born in 1930), post-menopausal, and had similar body mass index (BMI) (Table 1).

Table 1

Bone mineral density (BMD) and T score of the 3 pairs of MZ twins.

Sample	Year of Birth	BMI*	Hip BMD	T Score	Spine BMD	T Score
Concordant pair
1A	1930	22.5	0.6545	-2.67	0.8313	-1.96
1B	1930	23.1	0.6373	-2.81	0.8861	-1.46
Discordant pairs
2A	1930	24.5	0.6992	-2.30	0.8844	-1.48
2B	1930	22.3	0.8692	-0.88	0.9358	-1.01
3A	1930	24.6	0.6473	-2.73	0.8236	-2.03
3B	1930	22.8	0.8145	-1.33	1.0881	+0.37

* BMI: body mass index

Primary osteoblast-like cell culture was established from bone marrow precursor cells. The cell differentiation was assessed by the expression of osteogenic marker such as alkaline phosphatase (ALP) and stromal cell marker STRO1 demonstrated by flow cytometric analysis (data not shown), also more than half of the colonies were stained positive with ALP. No sub-culture or passage has been carried out.

In each of the 6 microchips, the mRNA in the samples hybridized to about half of the probe sets. From the comparison results with MAS5 using a nominal 2 fold cut-off, there were 417 and 283 probe sets showing at least a 2 fold decrease or increase in expression respectively within the concordant control pair, 1A/1B (Table 2). The respectively decrease or increase probe sets for the discordant case pairs were 90 and 45 for 2A/2B pair and 133 and 59 for 3A/3B (Table 2). In total, there were 700, 135 and 192 probe sets changes in 1A/1B, 2A/2B and 3A/3B respectively, or 3.1%, 0.6% and 0.9% of the total probe sets. Of particular note is that one probe set showed increased expression in both the discordant case pairs but not the concordant pair – this was chondroitin beta 1,4 N-acetylgalactosaminyltransferase (ChGn, NM_018371), a gene which regulates an important reaction in cartilage synthesis.

Table 2

Numbers of probe sets showed either decrease of increase in gene expression comparing within and between the MZ pairs. 2A/2B and 3A/3B are the discordant pair for BMD, and 1A/1B is the concordant pair.

Comparison	Changes		No. of probe (% of total)
		Twin	1B	2B	3B
Within pair	Decrease*	1A	417 (1.8%)
		2A		90 (0.4%)
		3A			133 (0.6%)
		Twin	1B	2B	3B
	Increase*	1A	283 (1.3%)
		2A		45 (0.2%)
		3A			59 (0.3%)
	Total	Twin	1B	2B	3B
		1A	700 (3.1%)
		2A		135 (0.6%)
		3A			192 (0.9%)
Between pair		Twin	2A	2B
	Total	1A	791 (3.6%)
		3A		1251 (5.9%)
		3B	1150 (5.6%)	1114 (5.2%)

* Decrease or increase is referred to 1B, 2B and 3B who had higher BMD than 1A, 2A and 3A respectively.

The probe sets in Table 2 were submitted to the NetAffx Centre for gene ontology (GO) analysis. Results are shown in Table 3. The patterns of GO molecular functions for decrease or increase in the three pairs are very similar with about 80% of the probe sets categorized under 3 GO terms: 5488 binding activity, 3824 enzyme activity and 4871 signal transducer activity. Hierarchy maps showed nodes branched out from these 3 main GO terms also had similar patterns for the three pairs (data not shown). However there were nodes with several probe sets grouped especially to the 2A/2B and 3A/3B case pairs but not in 1A/1B control pair. These nodes were 5125 cytokine activity, 8083 growth factor activity, 5194 cell adhesion molecular activity and 3700 transcription factor activity (Table 4). These molecular functions which might involve bone physiology included cytokine/growth factors, transcription factors and cell adhesion molecules.

Table 3

Numbers of probe sets showing decrease or increase of the three most frequent gene ontology (GO) molecular function in each pair's comparison. 2A/2B and 3A/3B are the discordant pair for BMD, and 1A/1B is the concordant pair.

		GO molecular function*	No of probe sets
1A/1B	Decrease**	All	417
		5488	159
		3824	91
		4871	85
	Increase**	All	283
		5488	111
		3824	69
		4871	45
2A/2B	Decrease	All	90
		5488	35
		4871	20
		3824	17
	Increase	All	45
		5488	18
		3824	6
		4871	6
3A/3B	Decrease	All	133
		5488	53
		3824	27
		4871	15
	Increase	All	59
		5488	23
		3824	11
		4871	10

* Gene ontology molecular function terms are: 5488 – binding activity; 3824 – enzyme activity and 4871 – signal transducer activity. ** Decrease or increase is referred to 1B, 2B and 3B who had higher BMD than 1A, 2A and 3A respectively.

Table 4

List of genes grouped under gene ontology molecular functions that showed changes in the discordant monozygotic (MZ) twin pairs but not the concordant pair.

Probe Set ID	Gene Name	Fold changed	Gene Symbol
5125 Cytokine activity in 2A/2B decrease
204926_at	inhibin, beta A (activin A, activin AB alpha polypeptide)	2.1	INHBA
205067_at	interleukin 1, beta	4.3	IL1B
205476_at*	chemokine (C-C motif) ligand 20	3.0	CCL20
205767_at	Epiregulin	3.3	EREG
206569_at	interleukin 24	3.3	IL24
210557_x_at	colony stimulating factor 1 (macrophage)	2.1	CSF1
211839_s_at	colony stimulating factor 1 (macrophage)	2.8	CSF1
221577_x_at*	prostate differentiation factor	2.3	PLAB
39402_at	interleukin 1, beta	4.0	IL1B
8083 Growth factor activity in 2A/2B decrease
204926_at	inhibin, beta A (activin A, activin AB alpha polypeptide)	2.1	INHBA
205067_at	interleukin 1, beta	4.3	IL1B
205767_at	Epiregulin	3.3	EREG
210557_x_at	colony stimulating factor 1 (macrophage)	2.1	CSF1
211839_s_at	colony stimulating factor 1 (macrophage)	2.8	CSF1
221577_x_at*	prostate differentiation factor	2.3	PLAB
39402_at	interleukin 1, beta	4.0	IL1B
5194 Cell adhesion molecular activity in 2A/2B increase
202877_s_at	complement component 1, q subcomponent, receptor 1	2.0	C1QR1
203477_at	collagen, type XV, alpha 1	2.0	COL15A1
204677_at**	cadherin 5, type 2, VE-cadherin (vascular epithelium)	5.7	CDH5
208983_s_at**	platelet/endothelial cell adhesion molecule (CD31 antigen)	3.5	PECAM1
212713_at	microfibrillar-associated protein 4	3.7	MFAP4
3700 Transcription factor activity in 3A/3B decrease
207034_s_at	GLI-Kruppel family member GLI2	2.0	GLI2
207147_at	distal-less homeo box 2	2.6	DLX2
208003_s_at	nuclear factor of activated T-cells 5, tonicity-responsive	3.3	NFAT5
208859_s_at	alpha thalassemia/mental retardation syndrome X-linked (RAD54 homolog, S. cerevisiae)	2.8	ATRX
209088_s_at	ubinuclein 1	2.1	UBN1
212420_at	E74-like factor 1 (ets domain transcription factor)	2.1	ELF1
215698_at	retinoblastoma binding protein 2	2.1	RBBP2
3700 Transcription factor activity in 3A/3B increase
205991_s_at	paired mesoderm homeo box 1	2.0	PMX1
207001_x_at	delta sleep inducing peptide, immunoreactor	2.0	DSIPI
207558_s_at	paired-like homeodomain transcription factor 2	3.7	PITX2
208502_s_at	paired-like homeodomain transcription factor 1	8.6	PITX1
210135_s_at	short stature homeobox 2	7.5	SHOX2
221950_at	empty spiracles homolog 2 (Drosophila)	4.3	EMX2
5194 Cell adhesion activity in 3A/3B increase
204677_at	cadherin 5, type 2, VE-cadherin (vascular epithelium)	3.5	CDH5
208983_s_at	platelet/endothelial cell adhesion molecule (CD31 antigen)	2.6	PECAM1

* Also show decrease in 1A/1B pair but the changes were less than 2 folds ** Also show increase in 3A/3B pair

We have also performed between pair comparisons: 1A/2A, 2A/3B, 2B/3A and 2B/3B. As expected the number of probe sets showing at least 2 fold changes were higher than the within pair comparisons. The numbers of changes for 1A/2A, 2A/3B, 2B/3A and 2B/3B were 791, 1150, 1251 and 1114 respectively or 3.6, 5.6, 5.9 and 5.2% (Table 2). GO analysis of these probe sets, the same 3 GO terms (5488 binding activity, 3824 enzyme activity and 4871 signal transducer activity) came out as the commonest molecular functions annotated. The GO molecular functions, which found with the discordant case pairs but not the concordant control pair (Table 4), were present in the hierarchy maps in all 4 between pair comparisons (data not shown).

Discussion

We have performed a genome-wide gene expression study on 6 RNA samples extracted from primary osteoblast-like cultures prepared from 3 pairs of MZ twins ascertained for BMD. The gene expression is expected to be very similar if not identical within pairs because MZ twins should for practical purposes be 100% genetically identical and their tissue type should have the same expression profile. Indeed in the 2A/2B and 3A/3B case pairs, despite being discordant for BMD at the hip, the percentage of probe sets showing alteration more than 2 folds were only 0.6 and 0.9% respectively (Table 2), in other words the gene expression was over 99% similar. This less than 1% difference in expression in the osteoblast-like cells is likely to be due to the discordance in BMD. However paradoxically the control pair 1A/1B who were concordant for BMD showed more variation – 3.1%, or about 97% similar. There is little data in the literature to compare, but a very recent publication [8] found 3.2% of genes showed expression changes of more than 2 folds for a single MZ twin pair discordant for multiple myeloma. The authors used Affymetrix hu95av2, which tested smaller number of genes than the Affymetrix human U133A chip we used, and the experimental conditions might be different as well. Nevertheless it showed that variation could be up to 3%. Another report using lymphoblastoid cells from 35 unrelated individuals [7], who are parents from the Centre d'Etude du Polymorphism Humain (CEPH), showed 15.7% (814 out of 5184) genes were significantly variable in expression. The 3.1% gene variation of the 1A/1B control pair is much lower than this unrelated individuals' figure of 15.7%. The 1A/1B pair is concordant for BMD, but they could be discordant for other phenotypes or environmental factors not recorded by us, or this might reflect normal variation in both in vivo gene expression and microarray assay.

The majority (about 80%) of changes within the three pairs fall into 3 GO molecular terms (Table 3) – 5488 binding activity, 3824 enzyme activity and 4871 signal transducer activity, and in fact the hierarchy maps are similar in terms of the branches of nodes. It has been suggested that allelic differences might underlie gene expression variations [9]. But as MZ pairs should not have any allelic differences, hence it is possible that the genes of these 3 GO terms could express differently independent of polymorphisms of the genetic sequences, and possibly regulated by external, environmental or epigenetic factors. It will be of interest to find out whether the gene expression pattern in other tissue types is similar or not. If as likely this were a universal pattern for different types of tissues, it would suggest those differentially expressed genes are possibly regulated by external non-heritable genetic factors.

Exploring the hierarchy maps in more detail, there are particular nodes in the discordant case pairs with several probe sets grouped under them (Table 4), and these nodes are either absent from the control pair, or contain only one or two probe sets. As the two case pairs are discordant for bone density, it is of interest that these genes are all related to bone physiology. The first group of genes belong to cytokine and growth factors showing decrease in 2A/2B, they are the inhibin beta A (also known as activin A, activin AB alpha polypeptide) (INHBA) interleukin 1 beta (IL1B) and colony stimulating factor 1 macrophage (CSF1). For IL1B and CSF1, two probe sets for each gene showed consistent decrease strongly suggesting a genuine change was there. In fact, allelic variations of the IL1B gene have been reported to be associated with BMD and osteoporosis [10] and with osteoarthitis [11,12] a common joint disease which has associated sub-chondral bone changes and increased bone turnover and bone mass [13]. Mutations in the CSF1 gene can cause severe osteoporosis in rats [14-16]. It has also shown that INHBA may play crucial roles in the transformation of osteoprogenitor cells to osteoblasts during bone and bone marrow regeneration [17].

The second group of genes belong to transcription factors showing both increase and decrease in 3A/3B (Table 4). They are mainly homeo-box genes. The next group of genes are cell adhesion molecules such as the collagen type XV alpha 1 (COL15A1), cadherin 5 type 2 VE-cadherin (vascular epithelium) (CDH5) and platelet/endothelial cell adhesion molecule (CD31 antigen) (PECAM1) showing increase in both the case pairs. It has been shown that several homeo box genes together with different cell adhesion genes initiate skeletal development [18]. Furthermore, collagen type XV like other collagens is an important component for skeletal muscle and bone [19]. Additionally we also noted that expression of the chondroitin beta 1,4 N-acetylgalactosaminyltransferase (ChGn) gene was increased in both case pairs but not the control pair. ChGn is the key enzyme in chondronitin sulphate biosynthesis [20,21] and therefore relates directly to cartilage formation, which is important in osteoarthritis, a related disease to osteoporosis.

All the probe sets and genes mentioned above showed more than 2 fold changes in expression between at least one of the case pairs, but not in the control pair. Whilst failure to replicate expression variations between the case pairs might reflect spurious results or false positives or negatives, it is also possible that there is genetic heterogeneity for BMD discordance in the two discordant case pairs. Indeed the 3A/3B was discordant for both the hip and spine (Table 1), whilst 2A/2B was discordant for the hip. If this were the case we might observe different genes whose expressions altered, and finding explanatory genes using this method could be more difficult than anticipated. Nevertheless the genes we highlighted are documented to have effects on bone physiology and could be influencing different BMD within the MZ twins.

Apart from comparing gene expression within pairs for the microarray results, we have performed between-pair comparisons. This is equivalent to comparing between unrelated individuals. As expected there were more gene expression changes than when comparing within pairs, on average 5% for between pairs as compared to 1.5% for within pairs (Table 2). However the overall patterns were consistent with the same 3 GO molecular terms showing most changes. From the hierarchy maps, 5125 cytokine activity, 8083 growth factor activity and 5194 cell adhesion molecular activity, the three GO molecular terms that present in discordant BMD case pairs but not the concordant control pair, were present in all four comparisons between the MZ pairs, basically unrelated individuals. This is not a surprise finding as cytokine, growth factor and cell adhesion molecules are associated and functions in a great variety of metabolic pathways, and would show changes in gene expression between unrelated individuals [7]. This further supports that the use of MZ twin pairs could select genes related to phenotypic differences efficiently.

Due to the limited amount of RNA samples obtained from primary culture, we could not perform the microarray experiments in duplicate or triplicate and therefore statistical and clustering analyses, which might give a better indication of the expression changes. For the same reason, we could not carry out further RNA analysis such as Northern blot and real time polymerase chain reaction to confirm our results. However we have performed the microarray experiments not for only one pair of MZ twins, as most previous reports did, but in 3 pairs of the same age with one as the control pair. Hopefully our results could be generalizable to cases of similar bone density differences. A replication study to confirm our results is essential.

There is little previous data on the use of MZ twins in gene discovery. One report in the literature using T cells from one MZ twin pair discordant for type I diabetes demonstrated that the interleukin-4 excretion was impaired in the diabetic cotwin and correlated to difference in gene expression [6,22]. In another report, one disease gene mutation has been pinpointed with the use of one MZ twin pair – this was one of the rare forms of cleft palate (van der Woude Syndrome) [23]. Our study is more detailed than these studies, and examined more genes in a common phenotypic trait.

We have identified 10 hip BMD discordant pairs out of 200. Seven of them with identifiable environmental confounders were excluded. These 7 pairs could have been studied in the same way looking for genes that were either influenced by environmental factors or resulted from the discordance in BMD. Although not the primary aim of the present study, it would be a useful idea for further studies looking for possible causative genes. Twin registries around the world have collected tens of thousands twin pairs with detailed phenotypes recorded, it would cost no extra to find discordant MZ twin pairs to study despite their relative rarity.

Conclusion

In conclusion, with the latest expression chip of over 20,000 genes, it is possible to examine genome-wide gene expression. Using the example of osteoblast-like cells in our MZ twins discordant for BMD at hip, we obtained promising results – indicating that whilst 99% of expressed genes are similar, several identified genes showing different expressions within pairs have potential roles in bone physiology. Although without further experiment, we cannot confirm or conclude these are the genes definitely related to BMD, we believe we have shown the potential and cost-effectiveness of further gene expression studies using discordant MZ twins.

Methods

Twin pairs were selected from a panel of 400 unselected adult female MZ twins for discordance of BMD at the hip using a Hologic QDR 4500 densitometer. Pairs were excluded if there was an obvious environmental cause for the discordance such as smoking, anorexia, bulimia, amenorrhea or chronic diseases like colitis, surgical menopause, use of bone drugs and alcoholism. At least one of the twin pair had to be close to the World Health Organization definition of osteoporosis – namely being 2.5 standard deviations (SD) below the peak bone mass of an average 30 year-old (called T-score).

Following ethical approval from the St Thomas' Hospital Ethics committee with verbal and written consent, bone marrow aspirates from the iliac crest were obtained from three pairs of female MZ twins. The bone marrow sample was added to 25 ml of heparinised serum-free Dulbecco's Modified Eagle Medium (DMEM). Twenty ml of Ficoll-Hypaque were layered on top of the bone marrow suspension and centrifuged at 600 g for 30 minutes. Mononuclear cells at the interface were harvested and washed. The cells were counted and plated out under basal conditions in DMEM containing 15% foetal calf serum and 100 μM ascorbate-2-phosphate. The cells were seeded at the same density, 104 cells per cm2. After 28 days, cell proliferation and the total number of colonies were determined. The colony forming units (CFU's) were stained for alkaline phosphatase (ALP) and the percentage of ALP positive CFU's were derived. All experiments were done in quadruplicate. Flow cytometric analysis (FACS) was carried out to demonstrate the pattern of expression of the stromal cell marker STRO1 and ALP.

RNA was also extracted from the primary cell culture at confluence after 21 – 28 days using the Trizol method (Invitrogen, Paisley UK) according to the manufacturer's protocol. The six RNA samples were sent to the Medical Research Council (MRC) Geneservice of Human Genome Mapping Project – Resource Center at Cambridge, UK for microarray assay using Affymetrix Human U133A microchip. Experimental information can be found at the MRC Geneservice RNA Services web page [24]. Results within the same pair were compared using the Microarray Suite 5.0 (MAS5) software. The probe sets that showed decrease or increase between the paired samples with a signal ratio of equal or larger than 1.0 (that means the decrease or increase change is equal to or larger than 2 folds) were submitted to the Netaffx Analysis Centre of the Affymetrix Website [25] for analysis by the Gene Ontology Miming Tool [26]. There are 22,283 probe sets on each Human U133A microchip of which 13,615 are annotated with molecular function for the gene ontology.

List of abbreviations

ALP = alkaline phosphatase

BMD = bone mineral density

CDH5 = cadherin 5 type 2 VE-cadherin

CEPH = Centre d'Etude du Polymorphism Humain

CFU's = colony forming units

ChGn = chondroitin beta 1,4 N-acetylgalactosaminyltransferase

COL15A1 = collagen type XV alpha 1

CSF1 = colony stimulating factor 1 macrophage

DMEM = Dulbecco's Modified Eagle Medium

DXA = dual energy X-ray absorbtiometry

FACS = flow cytometric analysis

GO = gene ontology

IL1B = interleukin 1 beta

INHBA = inhibin beta A

MAS5 = Microarray Suite 5.0

MZ = monozygotic

PECAM1 = platelet/endothelial cell adhesion molecule (CD31 antigen)

Authors' contributions

YTM carried out the gene ontology analysis and drafted the manuscript. GH carried out the cell culture and RNA preparation. JNB performed the FACS analysis and phenotypic characterisation of the cultures. TDS identified the subjects for study, performed the tissue biopsy and co-wrote the manuscript. All authors read and approved the final manuscript.