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.


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 Xray 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][5][6][7][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 osteoblastlike 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).
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 Nacetylgalactosaminyltransferase (ChGn, NM_018371), a gene which regulates an important reaction in cartilage synthesis.
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.
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). 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 * Also show decrease in 1A/1B pair but the changes were less than 2 folds ** Also show increase in 3A/3B pair osteoporosis in rats [14][15][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 show-ing 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 cm 2 . 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.