Identification of mRNAs Related to Tibial Cartilage Development of Yorkshire Piglets

Abstract

Cartilage dysplasia is one of the important reasons for the weakness of pig limbs and hooves. Porcine rickets with weak limbs and hooves bring huge economic losses to the pig industry. However, research on the development of pig cartilage is lacking. This study investigated the key genes and molecular mechanisms involved in cartilage development via an RNA-seq technique. Samples of proximal tibia cartilage were collected from three normal piglets with 1 day, 14 days, and 28 days of age, respectively, and then these samples were divided into two comparison groups (1-day vs. 14-day group, 14-day vs. 28-day group). Through the transcriptome analysis, 108 differentially expressed genes (DEGs), such as FORL2, were obtained from 1-day vs. 14-day comparison group, and 3602 DEGs were obtained from 14-day vs. 28-day comparison group, including SOX9, BMP6, and MMP13. The gene ontology (GO) functional and KEGG pathway enrichment revealed that many functions of DEGs were related to bone development. The pathways of DEGs from Day 1 vs. Day 14 were mainly enriched in mineral absorption, but the DEGs of Day 14 vs. Day 28 were enriched in osteoclast differentiation. Then, the expression patterns of six candidate genes were verified via qPCR. In conclusion, candidate genes affecting cartilage development in Yorkshire pigs were obtained by transcriptome analysis, and the clues showed that Day 14 to Day 28 is a more active and extensive period in cartilage developments, which played a key role in revealing the molecular mechanism of pig cartilage development basis, also compensating for vacancies in cartilage research.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In the modern large-scale pig industry, pig’s foot disease has become one of the most common and complex problems in pigs and occurs in pigs of all ages. Piglets are susceptible to sow injury or joint infection in the bed, causing foot and hoof disease, which affects the survival and growth rate of piglets during development. At present, leg and foot disease has become a major factor restricting the production, reproduction, and longevity of pigs, and it has also caused a large number of elimination and economic losses in the pig industry. Limb and foot disease is the second most important reason for the elimination of young sows [1]; Heinonen et al. have shown that sows eliminated by leg and foot disease account for 20–50% of all eliminated sows, and young sows are more susceptible to hoof rupture [2]. Therefore, further research on the pathogenesis and pathological mechanism of swine foot and hoof disease can help to improve limb and hoof problems and reduce economic losses.

It has been reported that osteochondrosis is an important cause of weakness in pigs’ hooves [3]. Bone formation originated from mesenchymal progenitor cells, then undergoes development and differentiation into mature chondrocytes, and is eventually replaced by bone. Studies have reported that cartilage cells are differentiated in the agglutination center region of mesenchymal cells, and a cartilage primordium composed of immature circular chondrocytes is produced [4]. A number of transcriptional and signaling factors are involved in this process, in which the agglutination of mesenchymal cells can express an important transcription factor involved in chondrogenesis Sox9 [5].

The most common site of rickets is the knee joint [6], followed by the elbow joint, and the ankle and hip joints are relatively few [7]. Ytrehus and other studies have found that piglets’ growth cartilage receives blood supply through the cartilage passage, but over time, the cartilage tube is converted to cartilage by softening, which is a physiological process that does not cause rickets disease. However, it was found that there is a correlation between an abnormal cartilage tube and rickets disease, which leads to local ischemic necrosis due to blood supply interruption, which eventually leads to secondary rickets [8]. Genetic factors play an important role in cartilage lesions. Studies have found that the estimated heritability of bone strength was also moderate (0.26) [9]. However, cartilage research related to pigs is very scarce. The genetic mechanism of cartilage development has yet to be explored.

Therefore, RNA sequencing of the proximal tibial cartilage of Yorkshire piglets in different developmental stages was mainly used in this study. With a transcriptome analysis strategy, the DEGs were screened and verified by studying the developmental regulation mechanism of cartilage and related signal pathways. Revealing the genetic mechanism of cartilage development could help to improve the development of limbs and hooves, reduce pig culling rate and mortality, and reduce the economic loss of farms.

2. Materials and Methods

2.1. Experiment Material

In this experiment, Yorkshire piglets were used as the main experimental animals. All disposals of the experimental animals were approved by the Scientific Ethic Committee of Huazhong Agricultural University (HZAUSW-2017-008). The selected piglets from the experiment came from Guangxi Yangxiang Co., Ltd. The experimental animals were divided into three developmental stages: Days 1, 14, and 28 of birth. The number of piglets assigned to each group is three. Samples of each age are siblings of the same parity of the same sow. The piglets used in the experiment were all females. The selected sows had no leg disease and received normal feeding, sufficient milk, and purebred Danish Yorkshire sows, and we ensured that the selected sows had the same gestation period. When growing to Day 1, Day 14, or Day 28, respectively, the proximal tibia cartilage was collected from the left hind limb and the muscles and other contaminants were removed. All tissues were collected in 30 min and quickly placed in liquid nitrogen. RNA was extracted from the proximal cartilage of the tibia and stored at −80°C for transcriptome sequencing.

2.2. RNA Sequencing

A total amount of 1.5 μg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer’s recommendations and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNaseH-). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylylation of 3′ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure was ligated to prepare for hybridization. In order to select cDNA fragments with right length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 μl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95°C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR Primers, and Index (X) Primer. At last, products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using HiSeq 4000 PE Cluster Kit (Illumia) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 4000 platform and 150 bp paired-end reads were generated.

2.3. Bioinformatics Analysis

2.3.1. Quality Control and Sequence Alignment

Sequencing errors, low-quality sequences, sequence junctions, and so forth are common in raw reads. Therefore, it is necessary to perform basic retreatment on the FASTQ sequence file obtained by RNA-seq and filter out the sequence with poor sequencing quality to obtain clean reads. Short sequences were mapped to the pig’s reference genome version Sus Scrofa 11.1 using HISAT2 software (version 2.1.0). The mapped reads were sorted by SAMtools (version 1.9). The count was calculated by the featureCounts software (version 1.6.2).

2.3.2. Differential Expression Gene

Differential expression analysis of transcriptome was performed by “DESeq2(version 1.22.2)” R packages, and DEGs were screened, using P-adj <0.01 and $\left|\text{FC}\right|>2$ as a DEG requirement.

2.3.3. GO and Pathway Enrichment

Based on the differentially expressed gene data after screening, functional annotations were performed using online metascape software (http://metascape.org/gp/index.html), pathway enrichment with KEGG on the DAVID (https://david.ncifcrf.gov/), and reference of Sus Scrofa used in KEGG.

2.4. Validation of DEGs by RT-qPCR

According to the transcriptome sequencing data and bioinformatics analysis results, six DEGs with high values of Fold Change, which are related to bone development and disease, are chosen. And the sequences of genes in pigs were obtained in Ensembl (https://asia.ensembl.org/index.html). Primer design was carried out using Primer Premier 5 (http://www.premierbiosoft.com). The primer sequences are shown in Supplementary Table 2. β-Actin was used as an internal reference, and Quantitative PCR amplification was done using TOYOBO’s SYBR Green I real-time PCR Mix, Adopt Roche’s LightCycler480 Real-Time PCR Detection System. Gene expression on Day 1 was used as a control group. Statistical analysis method on mRNA level comparisons was t-test. The relative expression level of the gene was calculated by the 2^−ΔΔCt method and calculated as follows: ΔΔC_t = (C_t Target–C_tβ-actin) Apiece-(C_t Target–C_tβ-actin) Control.

3. Results

3.1. Transcriptome Data Detection

The clean reads of the obtained proximal tibial cartilage tissue were compared with the reference genome of the pig (Supplementary Table 1), and the total number of reads (total mapped reads) and the number of reads of the single position alignment (uniquely) were obtained. The above results indicate that the quality of the sequencing data is high, which can be used for subsequent analysis. The R function was used to calculate the distance between the samples. Figure 1 shows that the 28-day-old sample had the best biological repeatability, followed by the 14-day-old sample, and there is some deviation in the data of 1-day-old piglets; however, the results showed that the sample has a good repeat effect and meets the requirements of further analysis. It can be seen that the difference between the 1-day-old sample and the 14-day-old sample was small, and the 1-day-old and 14-day-old samples were both different from the 28-day-old sample. It can also be inferred that the cartilage development level of piglets at Day 1 and Day 14 was relatively close, and the level of cartilage development of piglets at Day 28 was different.

[figure omitted; refer to PDF]

3.2. Differentially Expressed Genes (DEGs)

The DEGs of 1-day-old, 14-day-old, and 28-day-old tibia cartilage were compared and divided into a 1-day-old and 14-day-old comparison group and a 14-day-old and 28-day-old comparison group. The most stringent screening conditions were selected by observing the volcano map (Figures 2(a) and 2(b)). Trends in DEGs were observed by making heat maps (Figures S1 and S2), and descriptive statistics were performed on the DEGs (Table 1). A total of 108 DEGs were found in the comparison group of Day 1 vs. Day 14, among which 74 DEGs were upregulated and 34 DEGs were downregulated. A total of 3602 DEGs were found in the comparison group of Day 14 vs. Day 28, of which 1897 DEGs were upregulated and 1705 DEGs were downregulated. Obviously, the number of DEGs between Day 14 and Day 28 is much larger than the number of DEGs between Day 1 and Day 14. Therefore, we speculate that Day 14 to Day 28 is a significant change in the cartilage development of the piglets.

[figures omitted; refer to PDF]

Table 1

DEGs of tibial cartilage RNA at different ages of piglets.

	Upregulated genes	Downregulated genes
Day 1 vs. Day 14	74	34
Day 14 vs. Day 28	1897	1705

3.3. GO Analysis of DEGs

GO functional enrichment of DEGs from the two comparison groups is shown in Figure 3(a). GO analysis of DEGs was performed mainly in the skeletal system development, cell chemotaxis, collagen biosynthetic process, bone development, and other GO terms. As shown in Figure 3(b), after the GO analysis of DEGs, the top 20 GO terms were obtained, and the biological processes mainly included skeletal system development, ossification, myeloid cell differentiation, and other GO terms.

[figures omitted; refer to PDF]

To further understand the relationships between the terms, a subset of enriched terms have been selected and rendered as a network plot (Figures 3(c) and 3(d)), where the terms with a similarity >0.3 are connected by edges. We selected the terms with the best $P$values from each of the 20 clusters, with the constraint that there are no more than 15 terms per cluster and no more than 250 terms in total. The network is visualized using metascape, where each node represents an enriched term and was colored first by its cluster ID. We found that the GO terms enriched by the DEGs found in the 1-day-old and 14-day-old comparison group were more closely related, while the GO articles enriched in the 14-day-old and 28-day-old comparison group were relatively independent.

By enriching the GO function of the two comparison groups, we found that more DEGs were aggregated in the comparison of 14-day-old and 28-day-old samples in the GO terms associated with bone development, and there are also more functions related to bone development that are aggregated. We also found DEGs aggregation in ossification, which is an important function of cartilage development. This indicates that more genes related to cartilage development changed from Day 14 to Day 28.

3.4. Pathway Enrichment of DEGs

Pathway enrichment of DEGs mapping is shown in Figures 4(a) and 4(b). In 1-day-old and 14-day-old comparison groups, the DEGs enrichment has fewer pathways, so the DEGs with P-adj <0.05 are used for pathway enrichment, and the DEGs are mainly enriched in mineral absorption, protein digestion and absorption, and chemokine signaling pathway. In the comparison between the 14-day-old and the 28-day-old groups, the DEGs were mainly enriched in pathways such as rheumatoid arthritis, hematopoietic cell lineage, osteoclast differentiation, and chemokine signaling pathway. The pathway enrichment of DEGs found that from Day 1 to Day 14, the change is more significant for mineral absorption, and from Day 14 to Day 28, it is more significant for osteoclast differentiation. At the same time, we selected one pathway of interest in each of the two comparison groups to indicate the location of gene enrichment as shown in Figures S3 and S4. A total of four genes were enriched in the pathway of mineral absorption, in which MT-2B and ICA were upregulated on Day 14, and ATP1A2 and ENSSSCG00000027157 were downregulated on Day 14. A total of 40 genes were enriched in the osteoclast differentiation pathway; with the exception of AKT3, FOS, LOC100626904, MAPK10, and NFATC2, all other 35 genes were upregulated at 28 days.

[figures omitted; refer to PDF]

3.5. Expression Pattern Analysis

Totally, 3602 DEGs were found by comparison of 14-day-old and 28-day-old groups, far more than the 108 DEGs obtained from comparison between 1-day-old and 14-day-old groups. Thus, compared with Day 1 to Day 14, Day 14 to Day 28 may be a special period of cartilage developmental changes. Therefore, we analyzed the expression patterns of the DEGs shared by the two comparison groups. We found that there are 77 identical DEGs between the two comparison groups (Figure 5(a)). As shown in Figure 5(b), the expression trend of these 77 DEGs is obtained. We can see that the trend of expression of these genes is clearly divided into two groups. The DEGs of the first cluster were relatively at a lower expression level on Day 1 and Day 28 and at a relatively higher expression level on Day 14, and the first cluster has 53 DEGs. There were 24 DEGs aggregated in the second cluster, which was at a relatively high expression level on Days 1 and 28 and a relatively low expression level on Day 14. Then, six candidate genes which were significantly different in at least one comparison (P-adj <0.05) and related to cartilage development were selected among the two comparison groups, and the expression trends of the three periods were observed, as shown in Figure S5. It was clearly found that Day 14 was a special period of cartilage development compared to Day 1 and Day 28.

[figures omitted; refer to PDF]

3.6. Validation of DEGs by Quantitative PCR (qPCR)

Real-time PCRs were performed to verify these six DEGs. The results showed (Figure 6) that SOX9, BMP6, and FORL2 were significantly upregulated on Day 14 and significantly downregulated on Day 1 and Day 28. The transcription factor Sox9 is required for the proliferation and differentiation of mesenchymal cells into chondrocytes. And Bmp6 is a member of the bone morphogenetic protein family and plays an important role in chondrogenesis and differentiation. FOLR2 plays a role in the regulation of folate receptors, and folate receptors have a great relationship with macrophages in the synovial membrane of osteoarthritis. The expression level of LRP5 at Day 14 was higher than that of Day 1, and the expression on Day 28 was significantly lower than that on Day 14. Lrp5 is one of the receptor proteins of the Wnt signaling pathway and can affect changes in bone mass. The expression level of MMP13 on Day 1 was higher than that on Day 14, and the expression on Day 28 was significantly higher than that on Day 14. MMP13 is involved in the degradation of collagen matrix during bone and cartilage development. The expression of TGFB3 on Day 28 was significantly lower than that on other days, and the expression level on Day 14 was lower than that on Day 1. TGFβ is associated with chondrogenesis and differentiation, and deletion may not affect cartilage formation but dissipate the function of differentiation. The regression analysis results of qPCR and RNA-seq data are shown in Figure 7. The result shows a certain reliability. Comparing the RNA sequencing results with the quantitative PCR experiments (Supplementary Table 3), it was found that the results of the five genes SOX9, LRP5, BMP6, MMP13, and FOLR2 were consistent, while TGFB3 was slightly different from the sequencing results.

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

4. Discussion

Pig’s limb and foot firmness is a quantitative trait, and genetic factors are among the most important factors. In 2013, Nikkilä et al. found that the leg heritability of pigs was 0.07–0.31 [10]. At the same time, differences in breeds can also lead to differences in heritability. Huang et al. reported that the leg and heritability of Duroc, Landrace, and Yorkshire pigs were 0.23, 0.30, and 0.39, respectively [11]. Although the heritability of these reports varies, it can be shown that the occurrence of leg and foot disease is affected by genetic factors. In 2010, Laenoi et al. used molecular biology experiments to find that MGP (Matrix gla protein) is a potential calcification inhibitor of an extracellular matrix. The mRNA abundance of healthy cartilage MGP is significantly lower than that of rickets [12]. Rangkasenee et al. found that the KRT8 gene is associated with bone mineral density and content in 310 binary and 298 commercial groups, the FAF1 gene is associated with articular cartilage disease, and the PTH1R gene is significantly associated with the articular cartilage score, and they are related to leg weakness [13].

Based on the results of transcriptome analysis, 108 DEGs were obtained in the 1-day-old and 14-day-old comparison group and 3602 DEGs in the 14-day-old and 28-day-old comparison group. It can be seen from the results that the number of DEGs in the 14-day-old and 28-day-old comparison group is much larger than that in the 1-day-old and 14-day-old comparison group. At the same time, through sample repeatability analysis, we found that the difference between the 1-day-old and 14-day-old samples was small. After GO functional enrichment and KEGG pathway enrichment analysis, the comparison group on Day 1 and Day 14 lacked bone developmental biological processes, such as osteoblast differentiation and osteoclast differentiation, as well as Wnt signaling pathway and BMP signaling pathway. This indicates that there is little difference in the molecular regulation of cartilage development on the day of birth and Day 14; however, after Day 14, DEGs are enriched in skeletal system development, ossification, myeloid cell differentiation, and other GO terms, as well as the osteoclast differentiation pathways. It is indicated that the molecular regulation of cartilage development varies with development time. Through we analyzed the common differential expression patterns in the two comparison groups and the changes in the expression levels of the selected candidate genes during the three periods, we found that Day 14 is in a relatively special time point, and it can even be said to be a turning point. Therefore, based on the above results, we speculate that the process from Day 14 to Day 28 is a process of greater changes in cartilage development, which was not found in previous studies.

Current studies have found that many signal pathways are involved in bone development, such as the Wnt signaling pathway, BMP signaling pathway, TGFβ signaling pathway, and Hedgehog signaling pathway. The Wnt signaling pathway plays an important role in the development of various organs and tissues, including bone. Its discovery provides new ideas and directions for studying bone development in the future. Wnt/β-catenin (Wnt/β-catenin) is particularly important for bone biology. Low-density lipoprotein receptor-related protein (LRP5) is an important protein molecule in the Wnt signaling pathway, and DKK1 and SOST are two endogenous inhibitors in the Wnt pathway. Lrp5 under normal conditions is combined with an endogenous inhibitor. When a missense mutation occurs in Lrp5, it can destroy its own binding to the inhibitor but does not affect the binding to the agonist and thus increases the bone mass, Currently, results in humans [14] and mice [15] have been reported in the literature. A number of studies have reported that this increase in bone mass is primarily a transmission of Wnt/β-catenin signal [16, 17]. And SMO in the Hedgehog (Hh) signaling pathway was also found to be differentially expressed during Day 14 vs. Day 28. Hedgehog was originally isolated from the body of Drosophila, and its mutations can make Drosophila embryos develop into hairy cells, which is also known as the hedgehog gene. The transmission of vertebrate Hh signal is mainly through two multichannel transmembrane proteins, Smoothened (Smo) and Patched 1 (Ptch1) [18]. The Hedgehog (Hh) signal is required for the differentiation of mesenchymal stem cells into osteoblasts during embryonic development.

In our study, the DEGs in the comparison group on Days 1 and 14 were mainly enriched in the pathway of mineral absorption. Recent studies have shown that many trace elements, including Zn, Cu, and Fe, are bone-protecting elements that lack the ability to delay bone mass increase during childhood and adolescence and accelerate bone loss at old or postmenopausal age [19]. Various trace elements affect the metabolism and production of bone by regulating the activity of different enzymes in bone metabolism, regulating calcium absorption and excretion, and promoting different pathways, such as osteogenesis and ossification [20, 21]. In 1989, Hall et al. reported that 27 children with reduced blood levels of Mn had ischemic necrotizing bone disease or selective growth disorder, and this Mn deficiency caused ischemic bone disease in chickens [22]. Another study by Staley showed that congenital skeletal malformations in Holstein calves are associated with maternal Mn deficiency [23]. Fe-deficient rats can present with bone mineralization and pathological changes in trabecular bone structure [24, 25] and severe iron deficiency anemia, and rats showed reduced bone matrix formation and decreased bone mineralization [26]. The recommended intake of Mn in the United States is 1.8 to 2.3 mg/d [27]. However, some elements, such as aluminum, cadmium, and lead, are toxic, and these elements can make the quality of bone deteriorate [28]. The DEGs in the comparison group on Days 14 and 28 were mainly enriched in the pathway of osteoclast differentiation. The osteoclast is a tissue-specific macrophage polykaryon created by the differentiation of monocyte/macrophage precursors cells at or near the bone surface [29]. Once formed, bone undergoes a process termed remodeling that involves the breakdown (resorption) and build-up (synthesis) of bone; this occurs in microscale throughout the skeleton. Bone remodeling is the predominant metabolic process regulating bone structure and function during adult life, with the key participant being the osteoclast [30, 31]. Most adult skeletal diseases are due to excess osteoclastic activity, leading to an imbalance in bone remodeling which favors resorption [32]. Such diseases would include osteoporosis, periodontal disease, rheumatoid arthritis, multiple myeloma, and metastatic cancers. The study of this pathway provides a deeper understanding of how diverse physiological and pathophysiological signals exert their effects on the RANK signaling pathway to induce osteoclastogenesis, bone resorption, and skeletal remodeling and so control bone mass [31].

According to the results of DEGs, six differentially large multiples and research-related DEGs were selected, namely, SOX9, LRP5, BMP6, MMP13, TGFβ3, and FOLR2. The transcription factor Sox9 is required for the proliferation and differentiation of mesenchymal cells into chondrocytes. If Sox9 is deleted, cartilage formation is blocked [4], and premature maturation of immature chondrocytes into mast cells takes place [33]. At the same time, Sox9 is thought to directly activate chondrocyte differentiation markers [34], which can induce the expression of Sox5 and Sox6 [5] together with Sox9 to activate the chondrocyte differentiation program [35]. Lrp5, a low-density lipoprotein receptor protein, is one of the Wnt signaling pathway receptors; Lrp5 deletion or missense mutations disrupt binding to inhibitors but do not disrupt binding to agonists, leading to changes in bone mass [36]. Bmp6 is a member of the bone morphogenetic protein family and plays an important role in chondrogenesis and differentiation. When Bmps lacks receptor binding, the expression of Sox5, Sox6, and Sox9 is absent in the cartilage anterior condensation of mice and does not produce chondrocytes [37]. TGFβ is associated with chondrogenesis and differentiation, and deletion may not affect cartilage formation but dissipate the function of differentiation [38]. The role of transforming growth factor beta (TGFβ) signaling in cartilage has been controversial. Functional studies in cell culture models have shown that TGFβ can trigger cartilage formation [39]. However, blocking TGFβ signaling in limbal mesenchymal cells by targeting receptor TGFβ2 results in a certain effect on chondrocyte differentiation and joint formation but does not seem to affect the initial stage of cartilage [38]. Matrix metalloproteinase 13 (MMP13) is secreted by multiple myeloma cells and induces osteoclast active factors produced by multiple myeloma cells. MMP13 has been implicated in degradation of the collagenous matrices during the development of bone and cartilage [40]. FOLR2 plays a role in the regulation of folate receptors, and folate receptors have a great relationship with macrophages in the synovial membrane of osteoarthritis [41]. After real-time quantitative PCR in each cartilage tissue, it was consistent with the trend of sequencing results.

5. Conclusions

In this study, a transcriptome analysis of the tibial cartilage from Yorkshire pigs with different development stages was performed. Totally, 108 DEGs were found in the 1-day-old and 14-day-old comparison group, and 3360 DEGs were found in the 14-day-old and 28-day-old comparison group as well. These clues showed that the number of DEGs from 14-day-old to 28-day-old group was much larger than that from 1-day-old to 14-day-old group. It is speculated that Day 14 to Day 28 is a critical period of cartilage developmental changes. Furthermore, transcriptome analysis and real-time PCR showed that the Fold Change of SOX9, FOLR2, LRP5, BMP6, MMP13, and TGFB3 genes was significant and there are reports related to bone, which was an important candidate gene affecting cartilage phenotype. This study identified candidate genes that affect the development of tibial cartilage, which will compensate for the vacancies in cartilage development-related research. It also inferred the critical period of cartilage development and made an indelible contribution to reveal the genetic mechanism of tibial cartilage development.

Disclosure

The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Authors’ Contributions

X. D. L. and H. Y. W. conceived and designed the experiments; G. J. M., D. D. Y., and C. W. performed the experiments and sample collection; S. F. F. analyzed the data and wrote the paper; X. Y. D. and S. H. Z. revised the paper. All authors have read and approved the manuscript for publication.

References

[1] L. Engblom, N. Lundeheim, A.-M. Dalin, K. Andersson, "Sow removal in Swedish commercial herds," Livestock Science, vol. 106 no. 1, pp. 76-86, DOI: 10.1016/j.livsci.2006.07.002, 2007.

[2] M. Heinonen, J. Oravainen, T. Orro, "Lameness and fertility of sows and gilts in randomly selected loose-housed herds in Finland," Veterinary Record, vol. 159 no. 12, pp. 383-387, DOI: 10.1136/vr.159.12.383, 2006.

[3] B. Jorgensen, J. Arnbjerg, M. Aaslyng, "Pathological and radiological investigations on osteochondrosis in pigs, associated with leg weakness," Journal of Veterinary Medicine Series A, vol. 42 no. 1–10, pp. 489-504, DOI: 10.1111/j.1439-0442.1995.tb00404.x, 1995.

[4] W. Bi, J. M. Deng, Z. Zhang, R. R. Behringer, B. de Crombrugghe, "Sox9 is required for cartilage formation," Nature Genetics, vol. 22 no. 1, pp. 85-89, DOI: 10.1038/8792, 1999.

[5] H. Akiyama, M.-C. Chaboissier, J. F. Martin, A. Schedl, B. de Crombrugghe, "The transcription factor Sox9 has essential roles in successive steps of the chondrocyte differentiation pathway and is required for expression of Sox5 and Sox6," Genes & Development, vol. 16 no. 21, pp. 2813-2828, DOI: 10.1101/gad.1017802, 2002.

[6] R. C. Schenck, J. M. Goodnight, "Osteochondritis dissecans," Journal of Bone and Joint Surgery, vol. 78 no. 3, pp. 439-456, 1996.

[7] T. S. Lindholm, K. Osterman, E. Vankka, "Osteochondritis dissecans of elbow, ankle and hip: a comparison survey," Clinical Orthopaedics and Related Research, vol. 148, pp. 245-253, 1980. https://www.ncbi.nlm.nih.gov/pubmed/7379401/

[8] B. Ytrehus, C. S. Carlson, N. Lundeheim, "Vascularisation and osteochondrosis of the epiphyseal growth cartilage of the distal femur in pigs--development with age, growth rate, weight and joint shape," Bone, vol. 34 no. 3, pp. 454-465, DOI: 10.1016/j.bone.2003.07.011, 2004.

[9] A. Storskrubb, M.-L. Sevon-Aimonen, P. Uimari, "Genetic parameters for bone strength, osteochondrosis and meat percentage in Finnish Landrace and Yorkshire pigs," Animal, vol. 4 no. 8, pp. 1319-1324, DOI: 10.1017/s1751731110000418, 2010.

[10] M. T. Nikkilä, K. J. Stalder, B. E. Mote, "Genetic parameters for growth, body composition, and structural soundness traits in commercial gilts," Journal of Animal Science, vol. 91 no. 5, pp. 2034-2046, DOI: 10.2527/jas.2012-5722, 2013.

[11] S. Y. Huang, H. L. Tsou, M. T. Kan, W. K. Lin, C. S. Chi, "Genetic study on leg weakness and its relationship with economic traits in central tested boars in subtropical area," Livestock Production Science, vol. 44 no. 1, pp. 53-59, DOI: 10.1016/0301-6226(95)00063-q, 1995.

[12] W. Laenoi, M. Uddin, M. Cinar, "Quantitative trait loci analysis for leg weakness-related traits in a Duroc x Pietrain crossbred population," Genetics Selection Evolution, vol. 43 no. 1,DOI: 10.1186/1297-9686-43-13, 2011.

[13] N. Rangkasenee, E. Murani, R. Brunner, "KRT8, FAF1 and PTH1R gene polymorphisms are associated with leg weakness traits in pigs," Molecular Biology Reports, vol. 40 no. 4, pp. 2859-2866, DOI: 10.1007/s11033-012-2301-9, 2013.

[14] L. M. Boyden, J. Mao, J. Belsky, "High bone density due to a mutation in LDL-receptor-related protein 5," New England Journal of Medicine, vol. 347 no. 12, pp. 943-944, DOI: 10.1056/nejm200209193471216, 2002.

[15] M. Kato, M. S. Patel, R. Levasseur, "Cbfa1-independent decrease in osteoblast proliferation, osteopenia, and persistent embryonic eye vascularization in mice deficient in Lrp5, a Wnt coreceptor," The Journal of Cell Biology, vol. 157 no. 2, pp. 303-314, DOI: 10.1083/jcb.200201089, 2002.

[16] C. N. Bennett, K. A. Longo, W. S. Wright, "Regulation of osteoblastogenesis and bone mass by Wnt10b," Proceedings of the National Academy of Sciences, vol. 102 no. 9, pp. 3324-3329, DOI: 10.1073/pnas.0408742102, 2005.

[17] A. Jackson, B. Vayssière, T. Garcia, "Gene array analysis of Wnt-regulated genes in C3H10T1/2 cells," Bone, vol. 36 no. 4, pp. 585-598, DOI: 10.1016/j.bone.2005.01.007, 2005.

[18] D. Huangfu, K. V. Anderson, "Signaling from Smo to Ci/Gli: conservation and divergence of hedgehog pathways from Drosophila to vertebrates," Development, vol. 133 no. 1,DOI: 10.1242/dev.02169, 2006.

[19] I. Zofkova, P. Nemcikova, P. Matucha, "Trace elements and bone health," Clinical Chemistry and Laboratory Medicine, vol. 51 no. 8, pp. 1555-1561, DOI: 10.1515/cclm-2012-0868, 2013.

[20] L. Tandon, G. V. Iyengar, R. M. Parr, "A review of radiologically important trace elements in human bones," Applied Radiation and Isotopes, vol. 49 no. 8, pp. 903-910, DOI: 10.1016/s0969-8043(97)10102-6, 1998.

[21] J. Aaseth, G. Boivin, O. Andersen, "Osteoporosis and trace elements--an overview," Journal of Trace Elements in Medicine and Biology, vol. 26 no. 2-3, pp. 149-152, DOI: 10.1016/j.jtemb.2012.03.017, 2012.

[22] A. J. Hall, B. M. Margetts, D. J. P. Barker, "Low blood manganese levels in Liverpool children with Perthes’ disease," Paediatric and Perinatal Epidemiology, vol. 3 no. 2, pp. 131-136, DOI: 10.1111/j.1365-3016.1989.tb00505.x, 1989.

[23] G. P. Staley, "Congenital skeletal malformations in Holstein calves associated with putative manganese deficiency," Journal of the South African Veterinary Association, vol. 65 no. 2, pp. 73-78, 1994.

[24] D. M. Medeiros, B. Stoecker, A. Plattner, D. Jennings, M. Haub, "Iron deficiency negatively affects vertebrae and femurs of rats independently of energy intake and body weight," The Journal of Nutrition, vol. 134 no. 11, pp. 3061-3067, DOI: 10.1093/jn/134.11.3061, 2004.

[25] M. Parelman, B. Stoecker, A. Baker, D. Medeiros, "Iron restriction negatively affects bone in female rats and mineralization of hFOB osteoblast cells," Experimental Biology and Medicine (Maywood), vol. 231 no. 4, pp. 378-386, DOI: 10.1177/153537020623100403, 2006.

[26] J. Diaz-Castro, M. R. López-Frías, M. S. Campos, "Severe nutritional iron-deficiency anaemia has a negative effect on some bone turnover biomarkers in rats," European Journal of Nutrition, vol. 51 no. 2, pp. 241-247, DOI: 10.1007/s00394-011-0212-5, 2012.

[27] P. Ammann, V. Shen, B. Robin, Y. Mauras, J.-P. Bonjour, R. Rizzoli, "Strontium ranelate improves bone resistance by increasing bone mass and improving architecture in intact female rats," Journal of Bone and Mineral Research, vol. 19 no. 12, pp. 2012-2020, DOI: 10.1359/jbmr.040906, 2004.

[28] J. A. Staessen, H. A. Roels, D. Emelianov, "Environmental exposure to cadmium, forearm bone density, and risk of fractures: prospective population study. Public health and environmental exposure to cadmium (PheeCad) study group," The Lancet, vol. 353 no. 9159, pp. 1140-1144, DOI: 10.1016/s0140-6736(98)09356-8, 1999.

[29] W. J. Boyle, W. S. Simonet, D. L. Lacey, "Osteoclast differentiation and activation," Nature, vol. 423 no. 6937, pp. 337-342, DOI: 10.1038/nature01658, 2003.

[30] T. J. Chambers, "Regulation of the differentiation and function of osteoclasts," The Journal of Pathology, vol. 192 no. 1,DOI: 10.1002/1096-9896(2000)9999:9999<::aid-path645>3.3.co;2-h, 2000.

[31] S. L. Teitelbaum, "Bone resorption by osteoclasts," Science, vol. 289 no. 5484, pp. 1504-1508, DOI: 10.1126/science.289.5484.1504, 2000.

[32] G. A. Rodan, T. J. Martin, "Therapeutic approaches to bone diseases," Science, vol. 289 no. 5484, pp. 1508-1514, DOI: 10.1126/science.289.5484.1508, 2000.

[33] W. Bi, W. Huang, D. J. Whitworth, "Haploinsufficiency of Sox9 results in defective cartilage primordia and premature skeletal mineralization," Proceedings of the National Academy of Sciences, vol. 98 no. 12, pp. 6698-6703, DOI: 10.1073/pnas.111092198, 2001.

[34] D. M. Bell, K. K. H. Leung, S. C. Wheatley, "SOX9 directly regulates the type-II collagen gene," Nature Genetics, vol. 16 no. 2, pp. 174-178, DOI: 10.1038/ng0697-174, 1997.

[35] P. Smits, P. Li, J. Mandel, "The transcription factors L-Sox5 and Sox6 are essential for cartilage formation," Developmental Cell, vol. 1 no. 2, pp. 277-290, DOI: 10.1016/s1534-5807(01)00003-x, 2001.

[36] U. M. Ayturk, C. M. Jacobsen, D. C. Christodoulou, "An RNA-seq protocol to identify mRNA expression changes in mouse diaphyseal bone: applications in mice with bone property altering Lrp5 mutations," Journal of Bone and Mineral Research, vol. 28 no. 10, pp. 2081-2093, DOI: 10.1002/jbmr.1946, 2013.

[37] B. S. Yoon, D. A. Ovchinnikov, I. Yoshii, Y. Mishina, R. R. Behringer, K. M. Lyons, "Bmpr1a and Bmpr1b have overlapping functions and are essential for chondrogenesis in vivo," Proceedings of the National Academy of Sciences, vol. 102 no. 14, pp. 5062-5067, DOI: 10.1073/pnas.0500031102, 2005.

[38] H.-S. Seo, R. Serra, "Deletion of Tgfbr2 in Prx1-cre expressing mesenchyme results in defects in development of the long bones and joints," Developmental Biology, vol. 310 no. 2, pp. 304-316, DOI: 10.1016/j.ydbio.2007.07.040, 2007.

[39] J. L. Carrington, P. Chen, M. Yanagishita, A. H. Reddi, "Osteogenin (bone morphogenetic protein-3) stimulates cartilage formation by chick limb bud cells in vitro," Developmental Biology, vol. 146 no. 2, pp. 406-415, DOI: 10.1016/0012-1606(91)90242-u, 1991.

[40] Y. Sasano, J.-X. Zhu, M. Tsubota, "Gene expression of MMP8 and MMP13 during embryonic development of bone and cartilage in the rat mandible and hind limb," Journal of Histochemistry & Cytochemistry, vol. 50 no. 3, pp. 325-332, DOI: 10.1177/002215540205000304, 2002.

[41] S. R. Jafarzadeh, D. T. Felson, "Updated estimates suggest a much higher prevalence of arthritis in United States adults than previous ones," Arthritis & Rheumatology, vol. 70 no. 2, pp. 185-192, DOI: 10.1002/art.40355, 2018.