Gene Expression Analysis Suggests Bone Development-Related Genes GDF5 and DIO2 Are Involved in the Development of Kashin-Beck Disease in Children Rather than Adults

Objective To investigate the differences in gene expression between children and adults with Kashin-Beck disease (KBD). Methods 12 children with KBD and 12 healthy children were selected and divided into 4 KBD vs. control pairs matched according to age and gender, with each pair having 3 KBD children and 3 healthy children. Additionally, 15 adults with KBD and 15 healthy adults were selected and divided into 5 KBD vs. control pairs matched according to age and gender, with each pair having 3 KBD adults and 3 healthy adults. Total RNA was isolated from peripheral blood mononuclear cells (PBMCs) respectively. A total of 367 target genes were selected based on previous genome-wide gene expression profile analysis. Expression levels of the 367 genes were evaluated by customized oligonucleotide microarray and the differentially expressed genes were identified. Quantitative real-time reverse transcription polymerase chain reaction (qRT-PCR) was conducted to validate the microarray data. Results A total of 95 (25.9%) genes in KBD children and 158 (43.1%) genes in KBD adults were found to exhibit more than two-fold change in gene expression level relative to healthy controls. By comparing differentially expressed genes identified in KBD children to those of KBD adults, 42 genes were found to be differentially expressed only in KBD children. And 105 genes were found to be differentially expressed only in KBD adults. Further, 16 differentially expressed genes common to both KBD children and adults were found to be asynchronously expressed in KBD children compared to KBD adults. Conclusion Significant differences in gene expression pattern were identified between KBD children and KBD adults, indicating different molecular mechanisms underlying cartilage lesions of KBD children and KBD adults. In addition, bone development-related genes GDF5 (expression ratio = 2.14±0.02) and DIO2 (expression ratio = 0.11±0.05) may contribute to the development of KBD in children rather than in adults.


Introduction
Kashin-Beck disease (KBD) is an endemic and chronic osteochondropathy with unknown etiology. The disease mostly occurs in children between the ages of 3 and 13 in a diagonal beltlike area ranging from Northeast to Southwest China, with no observed difference in occurrence between genders. Few new KBD patients are observed among adolescents and adults [1]. It has been reported that more than 2.5 million people in China suffer from KBD and about 30 million people are at risk. The prevalence of KBD in children aged from 7 to 13 years reached 50% in serious KBD areas [2].
KBD is characterized by focal chondronecrosis in the mature chondrocytes of articular cartilage and the growth plate cartilage, which can result in impaired endochondral ossification during childhood [3]. Moreover, the premature closure of the epiphyseal plate in KBD children leads to irregular and impaired skeletal development. Enlarged and deformed joints as well as shortened long bones resulting from irreversible bone dysplasia occur during growth. Secondary osteoarthritis becomes the primary clinical manifestation of adult KBD patients.
Compared to adult KBD patients, abnormal radiological signs in the metaphyseal area and retarded skeletal development are more representative in child KBD patients [4]. Therefore, child KBD samples should provide more useful information regarding the mechanism of primary cartilage lesions and enable development of effective prevention and treatment measures for KBD. However, due to the difficulty in collecting child KBD samples and recent low KBD incidence rate, most previous KBD studies have focused on adult KBD, which provides limited information regarding the mechanism of KBD cartilage lesions at an early stage.
To identify the genes involved in the development of KBD at an early stage, the expression levels of 367 genes in KBD children and adults were evaluated and the differentially expressed genes of each group were identified. 367 KBD genes were selected from results of previous genome-wide gene expression profile analysis in PBMCs and articular cartilage from KBD patients [5][6][7]. Microarray analysis was conducted to evaluate the expression levels of those genes in peripheral blood mononuclear cells (PBMCs) of KBD children and adults respectively. Quantitative real-time reverse transcription polymerase chain reaction (qRT-PCR) was performed to validate the accuracy of the microarray data. In this study, new evidence is presented suggesting that bone development-related genes GDF5 (expression ratio = 2.1460.02) and DIO2 (expression ratio = 0.1160.05) are involved in the development of early-stage cartilage lesions in KBD children. This is believed to be the first comparative gene expression analysis between children and adults with KBD.

Ethics Statement
This study was approved by the Institutional Review Board (IRB) of Xi'an Jiaotong University. All adult subjects or the respective guardians of child subjects gave their informed written consent by signing a document that had been carefully reviewed by the IRB of Xi'an Jiaotong University.

Sample Collection
All child KBD subjects were randomly chosen from the KBDendemic area of Tongde County and Guide County in Qinghai Province of China, while adult KBD subjects were chosen from another KBD-endemic area of Yongshou County in Shaanxi Province of China. Healthy subjects came from the same KBDendemic areas with the KBD case group. All subjects were Chinese Han people. KBD patients were diagnosed strictly according to national diagnostic criteria of Kashin-Beck disease in China [WS/T 207-2010]. Any subjects who had a history of other bone or joint diseases were excluded from this study. Child subjects were divided into six pairs, and adult subjects were divided into seven pairs, each consisting of three KBD patients (1 male and 2 females) and three healthy individuals (1 male and 2 females) matched according to age and gender. Microarray analysis was performed for four of the six child pairs, and qRT-PCR was performed for the remaining two pairs. Likewise, microarray analysis was performed for five of the seven adult pairs, and qRT-PCR was performed for the remaining two pairs. The basic characteristics of the study subjects are shown in Table 1.

RNA Extraction
Three milliliters of peripheral blood from every child and adult subject was collected and stored at 280uC. Blood samples were thawed at room temperature for 2 hours before RNA extraction. Total RNA was extracted with the Paxgene Blood RNA Kit Table 1. Basic characteristics of the study subjects. (PreAnalytix; Qiagen) following the manufacturer's instructions. The integrity of RNA was validated using 1% agarose gel electrophoresis. Extracted RNA was stored at 280uC until cDNA synthesis.

Selection of Target Genes
In this analysis, 367 genes were selected as target genes from 3 additional databases of previous microarray analysis carried out by other researchers [5][6][7]. In those analyses, a target gene with $2 or #0.5 fold-change in gene expression level and P#0.05 was regarded as differentially expressed. Overlapping differentially expressed genes and some differentially expressed genes whose function may be not related to cartilage were excluded (detailed in Figure 1).

Microarray Hybridization
Extracted total RNA was reverse-transcribed into complementary DNA (cDNA), and then transcribed into cRNA and labeled with CyDye with the Amino Allyl MessageAmp aRNA Kit (Ambion) according to the manufacturer's instructions. For each group, 0.5 mg of labeled cRNA was purified separately and then mixed together with hybridization buffer (Agilent In Situ Hybridization Plus kit) before being applied to microarrays. The customized primer array from the National Engineering Research Center for Miniaturized Detection System in China, which contains 367 oligonucleotide probes representing 367 human genes, was used to perform microarray hybridization following the recommended protocol. This is a two-color microarray which uses a two-channel labeling system. The control cRNA was labeled with Cy5 while the KBD cRNA was labeled with Cy3. The array design information has been loaded into the repository of ArrayExpress with the accession number A-MEXP-2399. Hybridization signals were recorded by an Agilent G52565BA scanner, and analyzed using Feature Extraction 9.3 (Agilent Technologies) and Spotfire 8.0 (Spotfire Inc., Cambridge, MA, USA) software. The fluorescent spots that failed to pass the quality control procedure were not used for further analysis. Linear and LOWESS normalization was conducted to eliminate possible dye-related bias in the microarray data. The gene expression data obtained from Spotfire 8.0 was imported into Excel spreadsheets (Microsoft Corp., Redmond, WA, U.S.A.) for downstream data analysis. The dataset of this study has been deposited in ArrayExpress under the accession number E-MTAB-2494.

Analysis of Microarray Data
To identify differentially expressed genes, expression ratio was calculated for each gene through dividing the fluorescent value of KBD patient by that of control subject. Genes with expression ratios #0.5 or $2.0 were regarded as differentially expressed in this study. P values were calculated using the Agilent P value log ratio algorithm [5], with Erf(x) calculated as: Erf(x) represents twice the integral of the Gaussian distribution, with a mean value of 0, variance of 0.5, and xdev denotes the deviation of the log ratio from 0. This calculation provides the statistical significance of the log ratio for each feature (i.e., transcript level) between red and green channels. Only P values, 0.05 were regarded as statistically significant.

Quantitative Real-Time Reverse Transcription PCR
To validate the microarray results of the adult group, six significant differentially expressed genes were selected for parallel qRT-PCR analysis, including Gpx4 (glutathione peroxidase 4, NM_001039847), CYB5A (cytochrome b5 type A, NM_001914), MGAT3 [mannosyl (beta-1,4-)-glycoprotein beta-1,4-N-acetylglucosaminyltransferase, NM_002409], COL1A1 (collagen, type I, alpha 1, NM_000088), HAPLN1 (hyaluronan and proteoglycan link protein 1, NM_001884) and HBB (hemoglobin, beta, NM_000518). In the validation experiment of the child group, the selected genes are Gpx4, COL1A1, HAPLN1, MGAT3, and HBB. Total RNA was isolated in the same manner as for oligonucleotide array analysis. Extracted total RNA was then converted into cDNA using superscript II reverse transcriptase (Invitrogen, Carlsbad, CA) and random primers. The ABI7500 Figure 1. Flow chart of the procedure used to select target genes. This study is based on three additional gene microarray analyses [5][6][7]. doi:10.1371/journal.pone.0103618.g001 Real-Time PCR system (Applied Biosystems, Foster City, CA, USA) was used to perform qRT-PCR experiment according to the manufacturer's instructions. All primer and probe sets were purchased from TaqMan Gene Expression Assays (Applied Biosystems). GAPDH (Glyceraldehyde-3-phosphate dehydrogenase, NM_002046) was simultaneously assayed as an endogenous invariant control for data normalization. The expression values of the selected genes were normalized to the expression level of GAPDH. Paired t tests were performed to determine significance levels of expression differences for the selected genes between KBD and healthy controls.

Clinical and Radiological Characteristics of Representative Children and Adults with KBD
Children and adults with KBD are usually at different stages of KBD and thus experience different typical pathological changes ( Figure 2). Flexion of the terminal parts of the finger (Figure 2A) typically emerges earlier than enlarged finger joints. The primary radiographic changes seen in KBD children are metaphyseal lesions in phalanges, including a cone-shaped and blurred margin of the epiphysis, and abnormal closure of the epiphyseal line from the center ( Figure 2C). The premature closure of the metaphyseal plate in KBD children results in irregular and impaired development of bones, such as enlarged and crooked joints, as well as shortened fingers ( Figure 2B). As seen in radiograph of KBD adults ( Figure 2D), extensive articular cartilage damages appear in advanced cases of KBD. The metaphysis has also completely closed and become enlarged, and the margin shows irregular marginal sclerosis. In addition, irregular bone ends, narrowed joint space, and marginal blurred carpal bones are also characteristics for adult KBD patients.

Microarray Data Analysis
In this study, 95 (25.9%) differentially expressed genes (31 upregulated and 64 down-regulated) were identified in KBD children, and 158 (43.1%) differentially expressed genes (58 upregulated and 100 down-regulated) were identified in KBD adults. Further analysis of the differentially expressed genes in KBD children relative to KBD adults revealed that 42 of 95 differentially expressed genes were differentially expressed only in KBD children, but not in KBD adults, including 16 up-regulated genes ( Table 2) and 26 down-regulated genes (Table 3). Those 42 genes are involved in various biological processes, such as bone development, growth factors and apoptosis. In addition, 105 differentially expressed genes were differentially expressed only in KBD adults, but not in KBD children, including 46 up-regulated genes and 59 down-regulated genes ( Table 4). These 105 genes cover a wild range of biological processes, such as metabolism, apoptosis, immune function and redox system.
Additionally, this study identified 53 differentially expressed genes shared by both KBD children and adults. Sixteen of the 53 genes were found to be expressed asynchronously between KBD children and adults, including 9 genes found to be up-regulated in KBD children but down-regulated in KBD adults, and 7 genes found to be up-regulated in KBD adults but down-regulated in KBD children (detailed in Appendix S1).

qRT-PCR Validation of Microarray Data
The results of qRT-PCR experiment were consistent with microarray analysis in both KBD children and adults (Figure 3). According to qRT-PCR results, the expression levels of Gpx4, MGAT3 were higher in children with KBD than healthy children (expression ratio$2), whereas expression levels of COL1A1, HAPLN1, and HBB were lower in children with KBD than healthy children (expression ratio#0.5). Expression levels of Gpx4, MGAT3, and HBB were higher in adults with KBD than healthy adults (expression ratio$2), whereas expression levels of COL1A1, CYB5A, and HAPLN1 were lower in adults with KBD than healthy adults (expression ratio#0.5).

Discussion
Gene expression comparative analysis was conducted using customized microarray in KBD children and adults. Significant differences in gene expression pattern were observed between children and adults with KBD. Forty-two genes (16 up-regulated and 26 down-regulated) were found to be differentially expressed in KBD children, but not in KBD adults. Among these genes, bone development-related genes, such as GDF5 (expression ratio = 2.1460.02) and DIO2 (expression ratio = 0.1160.05), are of particular importance, and may provide valuable information for further investigation of KBD development at an early stage.
Endochondral bone formation is responsible for most longitudinal bone growth. During this process, chondrocytes in the growth plate undergo a process of maturation and terminal differentiation, which is precisely regulated by a series of signaling molecules. Any disturbances in this complicated regulation system could lead to shortened, or even deformed bones. In the pathology of KBD cartilage, chondronecrosis appears in the mature chondrocytes of the growth plate and articular cartilage [3], which can result in impaired endochondral ossification and skeletal development.
GDF5 is a member of the bone morphogenetic protein (BMP) family. It plays an important role in the development of the skeletal system [8]. A missense mutation of this gene results in chondrodysplasia Grebe-type (CGT), which is an autosomal recessive disorder exhibiting severe limb shortening and dysmorphogenesis [9]. Additionally, a functional polymorphism site of GDF5 is associated with susceptibility to osteoarthritis [10]. Overexpression of a homologue of GDF5 in chickens can increase the length and size of skeletal elements [11]. In humans, it has been shown that over-expression of GDF5 and its receptors is associated with dedifferentiation of human articular chondrocytes, which also can be seen in KBD [12,13]. The polymorphism of haplotype TGC on GDF5 is found to be significantly associated with KBD, confirming the potential role of GDF5 in the development of KBD [14]. Based on this evidence, it appears that over-expression of GDF5 may play a potential role in the pathology of KBD, although the mechanism needs further evaluation.
The DIO2 gene encodes the deiodinase type 2 (D2), which removes an iodine atom from tetraiodothyronine (T4) to catalyze its conversion into triiodothyronine (T3), the active form of thyroid hormone. Thyroid hormone is a systemic factor that potently regulates skeletal maturation in the growth plate [15]. Two essential trace elements, selenium and iodine, are very important for synthesis and function of DIO2 which is a selenoprotein. Selenium and iodine deficiency are two important environmental factors for KBD [16,17]. Previous studies have shown that the average level of blood T3 is lower in preschool children from a KBD-endemic country [18]. Additional studies have found that hypothyroidism occurs more frequently in KBD children and adolescents than control subjects [16]. These studies indicate that the interaction of low selenium and iodine levels, thyroid hormone, and DIO2 may contribute to cartilage lesions in KBD children.
Additionally, 105 genes were identified to be differentially expressed in KBD adults, but not KBD children. Among them, 46 genes were up-regulated and 59 genes were down-regulated. We noticed that different members of glutathione peroxidase family were involved in adult KBD and child KBD. For example, GPX1, GPX2, GPX3 were abnormally expressed only in KBD adults (expression ratio = 2.2860.67, 2.3260.65, 2.4260.48, respective-  [19] and members of the glutathione peroxidase family are responsible for antioxidant function of organism. More glutathione peroxidase members are involved in adult KBD. It may be because oxidative is more severe and extensive in KBD adult. Some genes related to immune system function were also found to be differentially expressed in KBD adults rather than children, such as HLA-C, CCL4 and EBI2 (expression ratio = 0.3660.08, 0.3960.09, 2.4160.31, respectively). HLA-C encodes the HLA class I heavy chain paralogue. Class I molecules play an important role in the immune system through presenting peptides derived from endoplasmic reticulum lumen. The protein encoded by CCL4 is a mitogen-inducible monokine, which has chemokinetic and inflammatory functions. EBI2 is expressed in B-lymphocyte cell lines but its function is not clear yet. In adult KBD patients, secondary osteoarthritis is the primary clinical manifestation. These genes may be related to the development of secondary osteoarthritis of KBD.
Between children and adults with KBD, 16 common differentially expressed genes were found to be expressed asynchronously (9 up-regulated in KBD children but down-regulated in KBD adults, and 7 up-regulated in KBD adults but down-regulated in KBD children). To study the different roles these genes may play in different stages of KBD and how their expression levels change during development of KBD may be critical to understand KBD and warrants further research. For example, APAF1, which encodes a cytoplasmic protein that initiates apoptosis, was found to be up-regulated (expression ratio = 2.9060.30) in KBD children but down-regulated (expression ratio = 0.3460.04) in KBD adults. The protein product of this gene forms an oligomeric apoptosome in the presence of cytochrome c and dATP, which binds and cleaves caspase 9 preproprotein to make it release its mature, activated form. Activated caspase 9 stimulates the subsequent caspase cascade that commits the cell to undergo apoptosis [20]. Abnormal chondrocyte apoptosis has been proven to appear in the articular cartilage of both children and adults with KBD [21,22]. Amounts of apoptosis related factors, such as BCL2 and CASP6, among others, as well as related pathways, have been shown to be expressed abnormally in articular cartilage of adult KBD patients [7,[22][23][24]. Additionally, APAF1 was found to be up-regulated in articular cartilage of KBD adults [7]. APAF1 thus appears to be critical for the abnormal apoptosis seen in articular cartilage of KBD patients, although it exerts different biological effects at different stages of KBD.
In summary, the gene expression levels of children and adults with KBD were compared, and significant differences in gene expression pattern were found between KBD children and adults, which may indicate the molecular mechanism underlying cartilage lesions. In addition, the results of this study revealed several genes that are only abnormally expressed in KBD children or KBD adults, as well as some common differentially expressed genes that are expressed asynchronously between KBD children and adults.
The study results suggest that bone development-related genes GDF5 and DIO2 may play a potential role in the development of KBD in children rather than in adults and may provide a basis for further study of the molecular mechanism of KBD, particularly in children.

Supporting Information
Appendix S1 List of common differentially expressed genes in two groups. Although these differentially expressed genes were shared by two groups, they expressed asynchronously in KBD children and KBD adults. This table includes 7 genes down-regulated in KBD children but up-regulated in KBD adults and 9 genes down-regulated in KBD adults but up-regulated in KBD children. (DOC)