Integrated Analysis of Transcriptome Expression Profiles Reveals miRNA-326–NKX3.2-Regulated Porcine Chondrocyte Differentiation

The porcine body length trait is an essential factor affecting meat production and reproductive performance. It is evident that the development/lengthening of individual vertebrae is one of the main reasons for increases in body length; however, the underlying molecular mechanism remains unclear. In this study, RNA-seq analysis was used to profile the transcriptome (lncRNA, mRNA, and miRNA) of the thoracic intervertebral cartilage (TIC) at two time points (1 and 4 months) during vertebral column development in Yorkshire (Y) and Wuzhishan pigs (W). There were four groups: 1- (Y1) and 4-month-old (Y4) Yorkshire pigs and 1- (W1) and 4-month-old (W4) Wuzhishan pigs. In total, 161, 275, 86, and 126 differentially expressed (DE) lncRNAs, 1478, 2643, 404, and 750 DE genes (DEGs), and 74,51, 34, and 23 DE miRNAs (DE miRNAs) were identified in the Y4 vs. Y1, W4 vs. W1, Y4 vs. W4, and Y1 vs. W1 comparisons, respectively. Functional analysis of these DE transcripts (DETs) demonstrated that they had participated in various biological processes, such as cellular component organization or biogenesis, the developmental process, the metabolic process, bone development, and cartilage development. The crucial bone development-related candidate genes NK3 Homeobox 2 (NKX3.2), Wnt ligand secretion mediator (WLS), gremlin 1 (GREM1), fibroblast growth factor receptor 3 (FGFR3), hematopoietically expressed homeobox (HHEX), (collagen type XI alpha 1 chain (COL11A1), and Wnt Family Member 16 (WNT16)) were further identified by functional analysis. Moreover, lncRNA, miRNA, and gene interaction networks were constructed; a total of 55 lncRNAs, 6 miRNAs, and 7 genes formed lncRNA–gene, miRNA–gene, and lncRNA–miRNA–gene pairs, respectively. The aim was to demonstrate that coding and non-coding genes may co-regulate porcine spine development through interaction networks. NKX3.2 was identified as being specifically expressed in cartilage tissues, and it delayed chondrocyte differentiation. miRNA-326 regulated chondrocyte differentiation by targeting NKX3.2. The present study provides the first non-coding RNA and gene expression profiles in the porcine TIC, constructs the lncRNA–miRNA–gene interaction networks, and confirms the function of NKX3.2 in vertebral column development. These findings contribute to the understanding of the potential molecular mechanisms regulating pig vertebral column development. They expand our knowledge about the differences in body length between different pig species and provide a foundation for future studies.


Introduction
The evolution of the vertebrate body length has primarily been achieved by changes in the axial skeleton. In many instances, body proportions are reflected by the number of vertebrae per region and the morphology of the vertebrae [1][2][3][4]. For example, body elongation is achieved by increasing the number of vertebrae, the elongation of individual vertebrae, or a combination of both [5]. In mammals, the vertebral formula shows developmental constraint [6]. The number of specific vertebrae in most mammalian species is relatively conserved, but domestic pigs are among the few mammals where variation exists with regard to the number of vertebrae. Several studies have reported quantitative trait loci (QTL) for the number of vertebrae in pigs by genome scan [7][8][9][10][11][12]. Two major QTL have been discovered in these studies, and causative genetic variation has been identified in germ cell nuclear factor (NR6A1) [8] and Vertnin (VRTN) [9]. However, research on how body elongation is achieved by extending individual vertebrae during vertebral column development remains scarce.
Two distinct developmental processes form bone in mammals. (1) Intramembranous ossification is the developmental process of craniofacial skeletal elements. (2) Endochondral ossification is the developmental process of most of the skeletal elements in the body [13]. The vertebral column is formed through endochondral ossification. The development of the vertebral column is a complex process involving the mesenchymal condensation of undifferentiated cells, hyperplasic and hypertrophic growth, and the mineralization of chondrocytes. The proliferation and hypertrophy of chondrocytes are the main driving forces of vertebral elongation. This process is tightly regulated by various factors, including transcription factors (SRY-box transcription factor 9 (SOX9) [14], RUNX family transcription factor 2 (RUNX2) [15], and NKX3.2 [16]), fibroblast growth factors (FGFs) [17], bone morphogenetic proteins (BMPs) [18], and the Wnt [19] or Notch [20] signaling pathway. Defects in these factors often lead to skeletal dysplasia and short stature.
The transcribed genome makes up 70-90% of the whole genome in mammals; only around 2% of it produces proteins. Among the non-coding RNAs (ncRNAs) that lack protein-coding potential, miRNAs are a class of small endogenous ncRNAs (18-24 nt) that regulate gene expressions by the translational repression or degradation of their targets. There is growing evidence that miRNA-mediated mechanisms also play critical roles in bone development [21,22]. For example, miR-196 has been shown to target multiple Hox genes, making it potentially important for the development of vertebrae and the axial skeleton [23]. The loss of miR-140 in mice causes growth defects in the endochondral bones, resulting in dwarfism [24]. miR-140-5p regulates zebrafish embryonic bone development via targeting BMP-2 [25]. LncRNAs are another class of regulatory ncRNAs, with sizes ranging from 200 bp to 100 kb [26]. The evidence suggests that lncRNAs serve, by various mechanisms, as versatile regulators of diverse aspects of biology. The involvement of lncR-NAs in the differentiation of mesenchymal stem cells into osteoblasts has been revealed [27]. Although an increasing number of lncRNAs are being characterized by high-throughput sequencing studies, the regulatory mechanisms involving these lncRNAs during the bone development of pigs or other species have rarely been described.
This study analyzed the non-coding RNA and mRNA expression profiles obtained through next-generation RNA sequencing (RNA-seq) in the TIC of two pig breeds at 1 and 4 months. The functions and interactions between NKX3.2 and miRNA-326 were verified in chondrocyte models in vitro. This study aimed to obtain deep insights into coding and non-coding RNA expression regulation in the TIC associated with vertebral column development in pigs and to identify the specific mechanisms underpinning the regulation of vertebral lengthening. Furthermore, one advantage of this work is that we investigated for the first time the lncRNA, miRNA, and gene crosstalk that mediates the vertebral column development in pigs.

Vertebral Development in Pigs
The vertebral length and scleromere size of the thoracic and lumbar vertebrae were measured in Yorkshire and Wuzhishan pigs after slaughter. As the pigs grew, the vertebral length of the 4-month-old pigs was significantly longer than that of the 1-month-old pigs, in both breeds (Y4 vs. Y1, W4 vs. W1, p < 0.01, Figure 1A). However, the Yorkshire pigs showed much faster growth, and thus, ended up with much longer spines than the Wuzhishan pigs by the age of 4 months (Y4 vs. W4, Figure 1A). The scleromere size of the thoracic and lumbar vertebrae also showed the same pattern, with faster growth in the Yorkshire than in the Wuzhishan pigs (Y4 vs. Y1, W4 vs. W1, p < 0.01, Figure 1B). Interestingly, it was also observed that the thoracic vertebrae were always significantly shorter than the lumbar vertebrae in both pig breeds (lumbar vertebrae vs. thoracic vertebrae, p < 0.01, Figure 1B). Moreover, frozen sections were made in order to observe the morphology of the chondrocytes in the TIC ( Figure S1). Larger hypertrophic chondrocytes were observed in Y4, which is consistent with the results described above. of vertebral lengthening. Furthermore, one advantage of this work is that we investigated for the first time the lncRNA, miRNA, and gene crosstalk that mediates the vertebral column development in pigs.

Vertebral Development in Pigs
The vertebral length and scleromere size of the thoracic and lumbar vertebrae were measured in Yorkshire and Wuzhishan pigs after slaughter. As the pigs grew, the vertebral length of the 4-month-old pigs was significantly longer than that of the 1-month-old pigs, in both breeds (Y4 vs. Y1, W4 vs. W1, p < 0.01, Figure 1A). However, the Yorkshire pigs showed much faster growth, and thus, ended up with much longer spines than the Wuzhishan pigs by the age of 4 months (Y4 vs. W4, Figure 1A). The scleromere size of the thoracic and lumbar vertebrae also showed the same pattern, with faster growth in the Yorkshire than in the Wuzhishan pigs (Y4 vs. Y1, W4 vs. W1, p < 0.01, Figure 1B). Interestingly, it was also observed that the thoracic vertebrae were always significantly shorter than the lumbar vertebrae in both pig breeds (lumbar vertebrae vs. thoracic vertebrae, p < 0.01, Figure 1B). Moreover, frozen sections were made in order to observe the morphology of the chondrocytes in the TIC ( Figure S1). Larger hypertrophic chondrocytes were observed in Y4, which is consistent with the results described above.

RNA-Seq Analysis Identifies DETs in the TIC
Among the 89-114 million raw reads from the RNA sequencing, 97% were kept after trimming and filtering, and approximately 80% of the clean reads were mapped to the reference genome (Table 1). With the read counts mapped to each gene, the correlations of the gene expressions among the different individuals were analyzed. The gene expressions among the three biological replicates in each group showed relatively high correlation coefficients (0.845-0.949), suggesting a high similarity of gene expression between the pigs in the same breed at the same age ( Figure S2). In addition, the correlation of the gene expression between the two species at the same age was also higher than that between the two ages within the same breed, indicating that age has a more significant effect than breed on gene expression regulation. After strict filtering ( Figure S3A) and removal of the potential coding transcripts identified using the coding-non-coding index (CNCI), coding potential calculator (CPC), Pfam scan (PFAM), and phylogenetic codon substitution frequency (PhyloCSF) ( Figure S3B), 1735 lncRNAs were obtained. After differential expression analysis, 161 DE lncRNAs were detected in the TIC of the Yorkshire pigs, among which 119 were up-regulated and 42 were down-regulated in Y4 compared to Y1. A total of 275 DE lncRNAs were identified in the TIC of the Wuzhishan pigs, with 148 up-regulated genes and 127 down-regulated genes in W4 compared to W1. Eighty-six DE lncRNAs were detected in the TIC of the 4-month-old pigs in both breeds, including fifty-eight that were up-regulated and twenty-eight that were down-regulated in Y4 compared to W4. One hundred and twenty-six DE lncRNAs were identified in the TIC of the 1-month-old pigs in both breeds, including fifty that were up-regulated and seventy-six that were down-regulated in Y1 compared to W1 (Figure 2A). two ages within the same breed, indicating that age has a more significant breed on gene expression regulation. After strict filtering ( Figure S3A) and rem potential coding transcripts identified using the coding-non-coding index (CN potential calculator (CPC), Pfam scan (PFAM), and phylogenetic codon subst quency (PhyloCSF) ( Figure S3B), 1735 lncRNAs were obtained. After different sion analysis, 161 DE lncRNAs were detected in the TIC of the Yorkshire pi which 119 were up-regulated and 42 were down-regulated in Y4 compared to of 275 DE lncRNAs were identified in the TIC of the Wuzhishan pigs, with 14 lated genes and 127 down-regulated genes in W4 compared to W1. Eigh lncRNAs were detected in the TIC of the 4-month-old pigs in both breeds, inclu eight that were up-regulated and twenty-eight that were down-regulated in Y4 to W4. One hundred and twenty-six DE lncRNAs were identified in the TIC month-old pigs in both breeds, including fifty that were up-regulated and seven were down-regulated in Y1 compared to W1 (Figure 2A).  In the Yorkshire pigs, 1478 DEGs were identified, among which 1051 gen higher expression and 427 genes showed lower expression in Y4 than in Y1. A to DEGs, including 1365 with relatively higher and 1278 with relatively lower e were identified in the TIC of the Wuzhishan pigs in W4 compared to W1. A t DEGs were detected in the TIC of the 4-month-old pigs in both breeds, with 24 higher and 159 relatively lower expression genes in Y4 than in W4. A total of were seen in the TIC of the 1-month-old pigs in both breeds, including 253 higher and 497 relatively lower expressions in Y1 than in W1 ( Figure 2B). Y4 vs. Y1 revealed 74 DE miRNAs in the TIC, among which 57 were up and 17 were down-regulated. Fifty-one DE miRNAs, including twenty-three up-regulated and twenty-eight that were down-regulated, were detected in the on W4 vs. W1. Thirty-four DE miRNAs were identified in the TIC of the 4-mon in both breeds, with twenty-two up-regulated and twelve down-regulated g on Y4 vs. W4. Y1 vs. W1 also revealed 23 DE miRNAs in the TIC of the 1-mon in both breeds, including 15 that were up-regulated and 8 that were down-regu In the Yorkshire pigs, 1478 DEGs were identified, among which 1051 genes showed higher expression and 427 genes showed lower expression in Y4 than in Y1. A total of 2643 DEGs, including 1365 with relatively higher and 1278 with relatively lower expression, were identified in the TIC of the Wuzhishan pigs in W4 compared to W1. A total of 404 DEGs were detected in the TIC of the 4-month-old pigs in both breeds, with 245 relatively higher and 159 relatively lower expression genes in Y4 than in W4. A total of 750 DEGs were seen in the TIC of the 1-month-old pigs in both breeds, including 253 somewhat higher and 497 relatively lower expressions in Y1 than in W1 ( Figure 2B). Y4 vs. Y1 revealed 74 DE miRNAs in the TIC, among which 57 were up-regulated and 17 were down-regulated. Fifty-one DE miRNAs, including twenty-three that were upregulated and twenty-eight that were down-regulated, were detected in the TIC based on W4 vs. W1. Thirty-four DE miRNAs were identified in the TIC of the 4-month-old pigs in both breeds, with twenty-two up-regulated and twelve down-regulated genes based on Y4 vs. W4. Y1 vs. W1 also revealed 23 DE miRNAs in the TIC of the 1-month-old pigs in both breeds, including 15 that were up-regulated and 8 that were down-regulated ( Figure 2C).  Figure S4B). Interestingly, among the shared DE lncRNAs, DEGs, and DE miRNAs, there were 3 DE lncRNAs, 17 DEGs, and 1 DE miRNA that were up-regulated in Y4 vs. W4 and down-regulated in Y1 vs. W1, while 4 DEGs and 4 DE miRNAs were up-regulated in Y1 vs. W1 and down-regulated in Y4 vs. W4. Subsequently, these DE lncRNAs, DEGs, and DE miRNAs in the TIC based on Y4 vs. Y1, W4 vs. W1, Y4 vs. and W4, and Y1 vs. W1 were integrated into the analysis to reveal the potential regulatory pathways in vertebral column development.

Functional Analysis of DETs
To better understand the biological function of the non-coding RNAs, their target genes were predicted. Previous studies have shown that lncRNAs regulate the expression of neighboring protein-coding genes through cis-acting mechanisms [28,29]. In addition, lncRNAs can also regulate the expression of genes on other chromosomes by trans-acting mechanisms [30,31]. In the comparison groups of Y4 vs. Y1 and W4 vs. W1, 958, 2736, and 1785 genes were predicted by bioinformatic analysis to be the targets of 73 unique DE lncRNAs in Y4 vs. Y1; 187 unique DE lncRNAs in W4 vs. W1; and 88 shared DE lncRNAs in Y4 vs. Y1 and W4 vs. W1, respectively, including 157, 433, and 182 cis-target genes and 820, 2303, and 1658 trans-target genes, of which 19, 86, and 55 genes were regulated by both cis-acting and trans-acting mechanisms. For the 50 unique DE miRNAs in Y4 vs. Y1, 27 unique DE miRNAs in W4 vs. W1, 24 shared DE miRNAs in Y4 vs. Y1 and W4 vs. W1, and 5908, 5381, and 4914 target genes were predicted by miRanda for the 43, 10, and 12 miRNAs that were more abundant in the 4-month-old pigs than in the 1-month-old pigs, and 3148, 4058, and 3253 target genes were predicted for the 7, 17, and 9 miRNAs that were less abundant in the 4-month-old pigs than in the 1-month-old pigs. In addition, 711 and 470 target genes were predicted for two DE miRNAs (up-regulated in Y4 vs. Y1 and down-regulated in W4 vs. W1) and one DE miRNA (up-regulated in W4 vs. W1 and down-regulated in Y4 vs. Y1), which came from 24 shared DE miRNAs in Y4 vs. Y1 and W4 vs. W1.
In the comparison groups of Y4 vs. W4 and Y1 vs. W1, 828, 1214, and 219 genes were predicted by bioinformatics analysis to be the targets of 67 unique DE lncRNAs in Y4 vs. W4; 107 unique DE lncRNAs in Y1 vs. W1; and 19 shared DE lncRNAs in Y4 vs. W4 and Y1 vs. W1, respectively, including 147, 189, and 58 cis-target genes and 693, 1051, and 164 trans-target genes, of which 12, 26, and 3 genes were regulated by both cis-acting and trans-acting mechanisms. For the 23 unique DE miRNAs in Y4 vs. W4, 12 unique DE miRNAs in Y1 vs. W1, and 11 shared DE miRNAs in Y4 vs. W4 and Y1 vs. W1, 4185, 4920, and 2758 target genes were predicted by miRanda for the 17, 7, and 4 miRNAs that were more abundant in the 1-and 4-month-old Yorkshire pigs than in the Wuzhishan pigs, and 2477, 2629, and 1180 target genes were predicted for the 6, 5, and 2 miRNAs that were less abundant in the 1-and 4-month-old Yorkshire pigs than in the Wuzhishan pigs. In addition, 235 and 2301 target genes were predicted for one DE miRNA (up-regulated in Y4 vs. W4 and down-regulated in Y1 vs. W1) and four DE miRNAs (up-regulated in Y1 vs. W1 and down-regulated in Y4 vs. W4), which came from 11 shared DE miRNAs in Y4 vs. W4 and Y1 vs. W1, respectively.
GO and KEGG analyses were performed on the identified DEGs and target genes of the DE lncRNAs and DE miRNAs to evaluate the biological functions and potential regulatory pathways in vertebral column development. In the comparison group of Y4 vs. Y1 and W4 vs. W1, 495 (specific to Y4 vs. Y1), 1660 (specific to W4 vs. W1), and 983 (shared in Y4 vs. Y1 and W4 vs. W1) DEGs were assigned to 422, 511, and 784 significant (p < 0.05) GO terms. In the comparison groups of Y4 vs. W4 and Y1 vs. W1, 301 (specific to Y4 vs. W4), 647 (specific to Y1 vs. W1), and 103 (shared in Y4 vs. W4 and Y1 vs. W1) DEGs were assigned to 197, 563, and 391 significant (p < 0.05) GO terms. These GO terms were mainly concentrated in the biological process and molecular function categories, including system development, spindle elongation, cell differentiation, cartilage development, and bone maturation ( Figure S5A). The top GO terms related to bone development are shown in Table 2. The results of the KEGG analysis revealed that the commonly enriched pathways were involved in the cAMP, TGF-beta, and Wnt signaling pathways ( Figure S5B). The GO analysis (p < 0.01) suggested that the cis-and trans-target genes of the DE lncRNAs were primarily involved in the regulation of single organism processes, bone development, and bone maturation ( Figure S5C). The KEGG pathway analysis indicated that the target genes were involved in Wnt, TGF-beta, and Notch signaling pathways ( Figure S5D). The significant GO terms (p < 0.01) of the DE miRNA target genes predominantly included protein binding, developmental processes, skeletal system development, ossification, and axis elongation ( Figure S5E). The KEGG analysis significantly enriched the TGF-beta, Wnt, and Notch signaling pathways ( Figure S5F).

LncRNA, miRNA, and Gene Interaction Network Construction
To generate an lncRNA, miRNA, and gene interaction network map, the essential genes with functions associated with bone development were identified by integrative GO enrichment analysis with the DEG and non-coding RNA target genes ( Table 2). Seven essential candidate genes (NKX3.2, WLS, GREM1, FGFR3, HHEX, COL11A1, and WNT16), which were mainly involved in the Wnt and Notch signaling pathway and skeletal development regulation, were further identified by combining functional analysis with literature mining. The interaction network among the lncRNA, miRNA, and essential candidate genes was constructed using Cytoscape v3.2.1 (Figure 3), showing the miRNA-mRNA, lncRNA-mRNA, and lncRNA-miRNA-mRNA interactions. In total, the interaction networks contained 55 lncRNAs (21 annotated and 34 novel lncRNAs), 6 miRNAs, and 7 genes. The results prove that the genes, lncRNAs, and miRNAs were either more or less abundant in the TIC of the Yorkshire and Wuzhishan pigs, based on a comparison of the 4-month-old pigs with the 1-month-old pigs, and that they may contribute together to regulate vertebral column development by the interaction networks.
GO enrichment analysis with the DEG and non-coding RNA target genes (T essential candidate genes (NKX3.2, WLS, GREM1, FGFR3, HHEX, COL11A1 which were mainly involved in the Wnt and Notch signaling pathway and opment regulation, were further identified by combining functional analy ture mining. The interaction network among the lncRNA, miRNA, and essen genes was constructed using Cytoscape v3.2.1 (Figure 3), showing the m lncRNA-mRNA, and lncRNA-miRNA-mRNA interactions. In total, the in works contained 55 lncRNAs (21 annotated and 34 novel lncRNAs), 6 m genes. The results prove that the genes, lncRNAs, and miRNAs were eithe abundant in the TIC of the Yorkshire and Wuzhishan pigs, based on a com 4-month-old pigs with the 1-month-old pigs, and that they may contribu regulate vertebral column development by the interaction networks. Figure 3. LncRNA, microRNA, and key candidate gene interaction network constr nal, prismatic, and circular symbols represent the lncRNA, miRNA, and mRNA, r represents an up-regulation, and green represents a relative down-regulation com month-old pigs.

Quantitative PCR Validation of Differentially Expressed Genes
To validate the results of the RNA-seq analysis, six mRNAs, six lncR miRNAs from the TIC were randomly selected and quantified using qPCR results confirm that the expression patterns of the DE mRNAs, lncRNAs were consistent with their expression levels calculated from the RNA-seq d Figure 3. LncRNA, microRNA, and key candidate gene interaction network construction. Hexagonal, prismatic, and circular symbols represent the lncRNA, miRNA, and mRNA, respectively; red represents an up-regulation, and green represents a relative down-regulation compared to the 1-month-old pigs.

Quantitative PCR Validation of Differentially Expressed Genes
To validate the results of the RNA-seq analysis, six mRNAs, six lncRNAs, and six miRNAs from the TIC were randomly selected and quantified using qPCR (Figure 4). The results confirm that the expression patterns of the DE mRNAs, lncRNAs, and miRNAs were consistent with their expression levels calculated from the RNA-seq data.

NKX3.2 Inhibits Chondrocyte Differentiation
Tissue expression profiles of the key candidate genes were examined in the heart, liver, spleen, lung, kidney, longissimus dorsi, thoracic intervertebral cartilage, lumbar intervertebral cartilage, and bone tissue of the Yorkshire pigs ( Figure 5). NKX3.2 was expressed explicitly in the cartilage tissues, including the thoracic and lumbar intervertebral cartilage, and the band was brighter in the thoracic intervertebral cartilage.
To explore whether NKX3.2 is involved in chondrocyte development, ATDC5 chondrogenic progenitor cells and pig primary chondrocytes were treated with Lipofectamine 2000 to overexpress NKX3.2 during culturing in a differentiation medium. Proteoglycan production was reduced by the overexpression of NKX3.2, as indicated by a reduction in Alcian blue ( Figure 6A,B). Moreover, NKX3.2 overexpression significantly promoted the gene expression of the chondrocyte marker collagen type II alpha 1 chain (COL2A1) ( Figure 6C,D), while inhibiting the hypertrophic chondrocyte marker collagen type X alpha 1 chain (COL10A1) ( Figure 6C,D). Interestingly, the results are consistent in both the ATDC5 cells ( Figure 6C) and the pig primary chondrocytes ( Figure 6D) Tissue expression profiles of the key candidate genes were examined in the heart, liver, spleen, lung, kidney, longissimus dorsi, thoracic intervertebral cartilage, lumbar intervertebral cartilage, and bone tissue of the Yorkshire pigs ( Figure 5). NKX3.2 was expressed explicitly in the cartilage tissues, including the thoracic and lumbar intervertebral cartilage, and the band was brighter in the thoracic intervertebral cartilage. To explore whether NKX3.2 is involved in chondrocyte development, ATDC5 chondrogenic progenitor cells and pig primary chondrocytes were treated with Lipofectamine 2000 to overexpress NKX3.2 during culturing in a differentiation medium. Proteoglycan production was reduced by the overexpression of NKX3.2, as indicated by a reduction in Alcian blue ( Figure 6A, B). Moreover, NKX3.2 overexpression significantly promoted the gene expression of the chondrocyte marker collagen type II alpha 1 chain (COL2A1) (Figure 6C,D), while inhibiting the hypertrophic chondrocyte marker collagen type X alpha 1 chain (COL10A1) (Figure 6C,D). Interestingly, the results are consistent in both the ATDC5 cells ( Figure 6C) and the pig primary chondrocytes ( Figure 6D), suggesting that the overexpression of NKX3.2 delayed the differentiation of chondrocytes into hypertrophic chondrocytes. Collectively, these data indicate that NKX3.2 can inhibit chondrocyte hypertrophy.  cartilage, and the band was brighter in the thoracic intervertebral cartilage. To explore whether NKX3.2 is involved in chondrocyte development, ATDC5 chondrogenic progenitor cells and pig primary chondrocytes were treated with Lipofectamine 2000 to overexpress NKX3.2 during culturing in a differentiation medium. Proteoglycan production was reduced by the overexpression of NKX3.2, as indicated by a reduction in Alcian blue ( Figure 6A, B). Moreover, NKX3.2 overexpression significantly promoted the gene expression of the chondrocyte marker collagen type II alpha 1 chain (COL2A1) (Figure 6C,D), while inhibiting the hypertrophic chondrocyte marker collagen type X alpha 1 chain (COL10A1) (Figure 6C,D). Interestingly, the results are consistent in both the ATDC5 cells ( Figure 6C) and the pig primary chondrocytes ( Figure 6D), suggesting that the overexpression of NKX3.2 delayed the differentiation of chondrocytes into hypertrophic chondrocytes. Collectively, these data indicate that NKX3.2 can inhibit chondrocyte hypertrophy. The mapping of the lncRNA, miRNA, and gene interaction networks showed that miRNA-326 was predicted to target NKX3.2 ( Figure 3) and had two possible binding sites (binding site I and binding site Ⅱ) with NKX3.2 3'UTR ( Figure 7A). To confirm the miRNA-

NKX3.2 Was the Direct Target of miRNA-326
The mapping of the lncRNA, miRNA, and gene interaction networks showed that miRNA-326 was predicted to target NKX3.2 ( Figure 3) and had two possible binding sites (binding site I and binding site II) with NKX3.2 3 UTR ( Figure 7A). To confirm the miRNA-326 target on NKX3.2, the binding site I mutation plasmid (MT1), binding site II mutation plasmid (MT2), and double mutations vector in binding site I and binding site II (MT12) were constructed to disrupt putative base pairing between miRNA-326 and NKX3.2 3 UTR. The wild-type plasmids of binding site I and binding site II (WT) and MT1 or MT2 or MT12 were co-transfected into HeLa cells with miRNA-326 mimics or NC mimics. The mimics of miRNA-326 significantly decreased the luciferase activity of WT, MT1, and MT2 but did not affect the activity of MT12 ( Figure 7B). These results indicate that miRNA-326 could randomly bind binding site I and binding site II of NKX3.2 3 UTR.

NKX3.2 Was the Direct Target of miRNA-326
The mapping of the lncRNA, miRNA, and gene interaction networks showed that miRNA-326 was predicted to target NKX3.2 ( Figure 3) and had two possible binding sites (binding site I and binding site Ⅱ) with NKX3.2 3'UTR ( Figure 7A). To confirm the miRNA-326 target on NKX3.2, the binding site I mutation plasmid (MT1), binding site Ⅱ mutation plasmid (MT2), and double mutations vector in binding site I and binding site Ⅱ (MT12) were constructed to disrupt putative base pairing between miRNA-326 and NKX3.2 3'UTR. The wild-type plasmids of binding site I and binding site Ⅱ (WT) and MT1 or MT2 or MT12 were co-transfected into HeLa cells with miRNA-326 mimics or NC mimics. The mimics of miRNA-326 significantly decreased the luciferase activity of WT, MT1, and MT2 but did not affect the activity of MT12 ( Figure 7B). These results indicate that miRNA-326 could randomly bind binding site I and binding site Ⅱ of NKX3.

miRNA-326 Promotes Chondrocyte Differentiation
To determine whether miRNA-326 played a role in the chondrocyte differentiation, ATDC5 chondrogenic progenitor cells and pig primary chondrocytes were transfected with ssc-miRNA-326 mimics during culturing in a differentiation medium. In contrast to the treatment with NC, the overexpression of miRNA-326 resulted in an increase in proteoglycan production ( Figure 8A,B) and COL10A1 expressions ( Figure 8C,D) in both the ATDC5 cells and the pig primary chondrocytes. COL2A1 expression decreased upon miRNA-326 overexpression ( Figure 8C,D). Collectively, these data suggest that miRNA-326 can promote chondrocyte hypertrophy.

miRNA-326 Promotes Chondrocyte Differentiation
To determine whether miRNA-326 played a role in the chondrocyte differentiation, ATDC5 chondrogenic progenitor cells and pig primary chondrocytes were transfected with ssc-miRNA-326 mimics during culturing in a differentiation medium. In contrast to the treatment with NC, the overexpression of miRNA-326 resulted in an increase in proteoglycan production ( Figure 8A,B) and COL10A1 expressions ( Figure 8C,D) in both the ATDC5 cells and the pig primary chondrocytes. COL2A1 expression decreased upon miRNA-326 overexpression ( Figure 8C,D). Collectively, these data suggest that miRNA-326 can promote chondrocyte hypertrophy.

Discussion
The vertebral column is a metameric array of tightly interconnected skeletal elements, the vertebrae. It constitutes the supportive yet flexible backbone of all vertebrates. During embryogenesis, the vertebral column is derived from somites, which are primary segments sculpted from the presomitic paraxial mesoderm during somitogenesis [32]. Vertebral column elongation is achieved by somite segmentation and individual vertebra development. In pigs, the differences in body length among different breeds are mainly attributed to the number of thoracolumbar vertebrae, whereas in the same species, the number of vertebrae is fixed, and the changes in body length are primarily caused by the elongation of individual vertebrae. Our results show that the vertebral column of the Yorkshire pigs was significantly longer (p < 0.01) than that of the Wuzhishan pigs at the ages of both 1 and 4 months ( Figure 1A); this is because Yorkshire pigs have 21 to 23 thoracolumbar vertebrae, while Wuzhishan pigs have 19 thoracolumbar vertebrae, which is mainly due to the lack of thoracic vertebrae. Both pig breeds showed increasing vertebra scleroderma size as they grew ( Figure 1B), but the vertebra scleromere in the Yorkshire pigs became much more significant than that in the Wuzhishan pigs by the age of 4 months. Interestingly, the size of the lumbar vertebra scleromere was much longer (p < 0.01) than the thoracic vertebra scleromere ( Figure 1B) in both pig breeds. Moreover, the frozen sections of TIC also showed larger hypertrophic chondrocytes in Y4. This result is consistent with that described above. To better understand the potential molecular regulation activity underlying this variation of vertebral development, we explored the regulation of gene expression by RNA-seq analysis.
RNA-seq analysis allows for the rapid exploration of the critical genes associated with particular phenotypes or essential biological processes. Although previous RNA-seq studies have shown transcriptome changes in the cartilage of piglet tibia [33], studies on cartilage tissues are still lacking. Therefore, we conducted this RNA-seq study to analyze the non-coding RNAs and mRNA expression profiles in the TIC of the pig. A large number of DEGs that play essential roles in the regulation of bone development were identified. Key candidate genes were further identified through a combination of functional analysis with literature mining, including NKX3.2 [34][35][36][37][38][39][40], WLS [41][42][43], GREM1 [44][45][46], FGFR3 [47][48][49][50][51], HHEX [52][53][54], COL11A1 [55][56][57], and WNT16 [58,59]. These genes were mainly enriched in the BMP, FGF, Wnt, and Notch signaling pathways, suggesting that they are involved in the development of the pig vertebral column. In addition to the coding genes, the target genes of the lncRNAs or miRNAs were also enriched in these pathways.
WLS (also known as GPR177) is a chaperone protein that is required for the secretion of both canonical and non-canonical Wnt ligands [41]. The loss of WLS specificity in cartilage resulted in delayed chondrocyte differentiation, disorganized chondrocyte arrangement, and impaired endochondral bone formation [42]. Recent studies have shown that skeletal homeostasis was altered with decreased bone formation within 1-2 weeks of conditional WLS knockout in adult mice (5 months old) [43]. GREM1 is an antagonist of bone morphogenetic proteins and is expressed in osteoblasts. GREM1 null mice exhibited developmental skeletal abnormalities, decreased weight and body fat, and shortened femoral length [44]. The latest research suggests that the hyperexpression of GREM1 specifically blocks chick phalanx development by inhibiting PFR activity [45] and causing severe defects of skeletogenesis in Trim28 MKO mice [46]. HHEX, a homeodomain protein, is a transcription factor that participates in cell proliferation, differentiation, and migration. A study in mouse chondrogenic cell line ATDC5 showed that HHEX protein expression and subcellular localization were associated with chondrocyte maturation [52]. Recent studies have shown that HHEX negatively regulates osteoclast differentiation [53]. Our recent study also found that a TA haplotype on the HHEX promoter was significantly associated with body length in pigs [54]. COL11A1 is normally found in cartilage, which is responsible for the formation of cartilage structures in the early stages of skeleton development. Li et al. showed that chondrodysplasia in mice is associated with mutations in COL11A1 [55]. A novel deleterious heterozygous mutation in COL11A1 was associated with severe skeletal dysplasia [56]. A novel missense variant (G to A in exon 45) of COL11A1 resulted in short-limb skeletal dysplasia [57]. WTN16, a ligand of the Wnt signaling pathway, is associated with bone density, cortical thickness, bone strength, and fracture risk [58,59]. Wnt16 −/− zebrafish exhibited significant deformities in the head, spine, and tail and decreased bone mineral density and trabecular bone [59].
Non-coding RNAs (ncRNAs) have evolved in eukaryotes as epigenetic regulators of gene expression. The most abundant regulatory ncRNAs are miRNAs and lncRNAs. Each class of ncRNAs regulates gene expression through distinct mechanisms [60]. Therefore, to better understand the complex mechanism of vertebral column development, we constructed the lncRNA-mRNA, miRNA-mRNA, and lncRNA-miRNA-mRNA interaction networks for those key DEGs we identified (Figure 3). The lncRNA-target regulatory interaction networks, comprising 55 DE lncRNAs (of which 21 were annotated and 34 were novel) and the essential candidate genes, were constructed ( Figure 3). FGFR3 was predicted to be a putative target of ALDBSSCT0000008143, LNC000074, and LNC000799. FGFR3 is mainly expressed in cartilage and negatively regulates bone growth [47,48]. A recent study suggested that mice with a conditional knockout of FGFR3 in the chondrocytes displayed overgrowth of bone with a significantly increased bone mass at both 1 and 4 months of age [49]. A study on postmenopausal osteoporosis mice showed that FGFR3 activation inhibited the ability of bone regeneration and bone mineralization [50]. Human FGFR3-related skeletal dysplasia was also summarized in a recent review [51]. In the present study, FGFR3 showed a significantly lower expression in the 4-month-old pigs than in the 1-month-old ones ( Figure 4A), which is consistent with that reported above. ALDBSSCT0000008143 was predicted to regulate the expression of FGFR3 via cis-acting mechanisms and was expressed at a lower level in the 4-month-old than in the 1-month-old pigs ( Figure 4B). A previous study suggested that lncRNAs could activate the transcription of nearby genes via cis-acting mechanisms by recruiting remodeling factors to local chromatin [61]. Therefore, it is reasonable to speculate that ALDBSSCT0000008143 may target FGFR3 via cis-acting mechanisms to regulate vertebral column development based on the positive relationship between ALDBSSCT0000008143 and FGFR3 observed in this study.
To understand how the interactions between miRNAs and their target genes regulate vertebral column development, interaction networks between DE miRNAs and key candidate genes were constructed. A total of six miRNAs and seven genes formed eight miRNA-gene pairs. Of the above pairs, the integrated analyses identified anti-correlated ssc-miR-326-NKX3.2 pairs (Figure 3). NKX3.2, also known as Bapx1, is a transcription factor that belongs to the NK family of homeobox genes. It has been suggested that it plays a role in inhibiting chondrocyte hypertrophy and maintaining the chondrocyte phenotype during chondrogenic differentiation [34,35]. In recent years, the mutation and function of NKX3.2 have been extensively studied in mouse [36] and zebrafish models [37,38], as well as in human [39] skeletal abnormalities. In this study, NKX3.2 was explicitly expressed in cartilage tissues ( Figure 5). This suggests that NKX3.2 may play an essential role in the development of the porcine vertebral column. The function of NKX3.2 in inhibiting chondrocyte hypertrophy, and thus, delaying chondrocyte differentiation was further identified in vitro ( Figure 6). These results are similar to those of previous studies. Interestingly, in vivo cartilage-specific NKX3.2 overexpression resulted in overall size reduction and shorter vertebral columns [40].
Furthermore, the targeted relationship between miRNA-326 and NKX3.2 was demonstrated by a dual-luciferin assay. miR-326 was also further identified to promote chondrocyte hypertrophy by regulating the expression of NKX3.2 (Figures 7 and 8). Recently, studies have shown that miR-326 acts as a negative modulator for sonic hedgehog (Shh) signaling [62], which is a crucial regulatory pathway for bone morphogenesis [63]. Collectively, these results indicate that the interaction of miR-326 and NKX3.2 may jointly influence the development of vertebral columns by regulating chondrocyte differentiation. Thus, further investigations into the signaling pathways in which NKX3.2 is involved in regulating vertebral column development are needed, by ChIP-seq or RNA-seq analysis or a combination of both. The study of this pathway can be applied not only to the breeding of body length traits in pig production, the improvement of pork yield, and an increase in the efficiency of pig production, but it can also provide an ideal small-sized pig model for human physiological and pathological research.
Recent reports have suggested that lncRNAs can potentially interact with miRNAs and then regulate the expression of genes. This study systematically analyzed the complex interactions among miRNAs, lncRNAs, and their target genes and provided lncRNA-miRNA-gene interaction networks. For example, LNC000173 (Gene ID: XLOC_007104) was observed to have a significantly lower expression in the 4-month-old pigs than in the 1-month-old pigs ( Figure 4B) and was predicted to compete with miR-326 to bind with NKX3.2. Therefore, it was speculated that LNC000173 might regulate the porcine vertebral column development via the pig's LNC000173-miR-326-NKX3.2 interaction. The lncRNA-miRNA-mRNA interaction networks revealed a setting for the competition of endogenous RNAs that was possibly involved in the vertebral column development process. However, this hypothesis requires further study for verification.

Animal Sample Collection and RNA Isolated
Yorkshire pigs at the ages of 1 month (Y1, n = 3) and 4 months (Y4, n = 3) were obtained from the Beijing pig breeding farm (Beijing, China). Wuzhishan pigs at the ages of 1 month (W1, n = 3) and 4 months (W4, n = 3) were obtained from the Animal Sciences and Veterinary Institute, Hainan Academy of Agricultural Sciences (Haikou, China). The TIC was sampled after the animals were slaughtered. The tissue samples were frozen in liquid nitrogen and stored at −80 • C before RNA isolation.
According to the manufacturer's instructions, the total RNA was isolated from each TIC tissue using TRIZOL ® reagent (Invitrogen, Carlsbad, CA, USA). The RNA quality was examined on 1% agarose gel. The RNA purity was checked using the NanoPhotometer ® spectrophotometer (Implen, Westlake Village, CA, USA). The RNA concentration was measured using the Qubit ® RNA Assay Kit in the Qubit ® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA). The RNA integrity was assessed using the RNA Nano 6000 Assay Kit in the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).

Library Preparation for lncRNA Sequencing and Data Analysis
Sequencing libraries were prepared using 3µg of RNA per sample. The samples were indexed with the NEBNext ® Ultra™ Directional RNA Library Prep Kit for Illumina (NEB, Ipswich, MA, USA), following a procedure described previously [64]. The libraries were sequenced on a HiSeq 2500 platform (Illumina, San Diego, CA, USA), and 125 bp paired-end reads were generated. Raw reads were processed through in-house Perl scripts to remove reads containing adapter, poly-N, and low-quality reads. The remaining reads were mapped to the porcine reference genome using TopHat (v2.0.9). Mapped reads from each sample were assembled by both Scripture (beta2) [65] and Cufflinks (v2.1.1) [66] in a reference-based approach. The assembled transcripts were evaluated using five criteria to identify lncRNAs: (1) transcripts with exon number ≥2 were retained; (2) transcripts with length > 200 bp were retained; (3) known no-lncRNA annotations were removed; (4) transcripts with fragments per kilobase of exon per million fragments mapped (FPKM) < 0.5 were removed; (5) CNCI (v2) [67], CPC (0.9-r2) [68], PFAM (v1.3) [69], and PhyloCSF (v20121028) [70] were used to distinguish the mRNAs from the lncRNAs ( Figure S3A). The transcripts predicted to have coding potential by all four tools described above were excluded, and those without coding potential were classified as lncRNA candidates ( Figure S3B). The transcripts excluded above were used as candidate mRNAs. PHAST v1.3 was used for conservation analysis for the coding genes and lncRNAs [71].
Cuffdiff (v2.11) was used to calculate the FPKM of the lncRNAs and coding genes in each sample. Transcripts with a p-adjust of <0.05 were assigned as DETs. We searched the coding genes 100 kb upstream and downstream of the lncRNA as the cis-target genes. The potential trans-target genes of the lncRNA were identified by the Pearson correlation coefficients (r > 0.95 or r < −0.95) of their expression levels.

Library Preparation for microRNA Sequencing and Data Analysis
Three µg of RNA per sample was used as the input material for the small RNA library. The sequencing libraries were generated using NEBNext ® Multiplex Small RNA Library Prep Set for Illumina ® (NEB, Ipswich, MA, USA), following a procedure described previously [64]. The levels of miRNA expression were estimated using transcript per million (TPM) through the following normalization formula: normalized expression = mapped reads/total reads × 1,000,000 [72].
Differential expression analysis of the miRNAs among the groups was performed using the DESeq R package (1.8.3). The p-values were adjusted using the Benjamini and Hochberg methods. miRNAs with a corrected p-value of < 0.05 were classified as significantly differential expressions.

GO and KEGG Enrichment Analysis
Gene Ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analyses of the DEGs or target genes of the DE lncRNAs and DE miRNA were conducted as described previously [64]. The GO terms and KEGG pathways with p-values less than 0.05 were considered significantly enriched by the DEGs.

Cell Isolation, Culture, and Histological Analysis
The HeLa and ATDC5 cell lines were purchased from Xiehe Medical University (Beijing, China) and Otwo Biotech (Shenzhen, China), respectively, and the primary chondrocytes from the TIC of newborn piglets were isolated as previously described [74]. All the cells were cultured in a medium that contained DMEM, 10% FBS, and 1% penicillin/streptomycin.
For chondrocyte differentiation, the ATDC5 cells and pig primary chondrocytes were grown in a medium containing 1% insulin-transferrin-selenium (ITS, Gibco). To perform Alcian blue staining, the cells were first fixed in 10% NBF for 30 min and then incubated with a 1% Alcian blue/3% acetic acid solution for 30 min. The cultures were then rinsed with 70% ethanol, followed by ddH 2 O, and air-dried before imaging.
For transfection, the ATDC5 cells and pig primary chondrocytes were transfected with 2 µg of the overexpression vectors or empty plasmid using Lipofectamine 2000 (Invitrogen) in each well of a 6-well plate. After being transfected for two days, the cells were harvested for RNA extraction.
For the tissue slices, the cartilage tissue samples from the TIC were embedded in OCT tissue freezing medium (Leica, Germany) and serially sectioned (6 µm thickness). The obtained cartilage tissue sections were observed under a 10 × 20 camera microscope.

Plasmid Construction and Dual-Luciferase Reporter Assay
For the NKX3.2 overexpression plasmids, the CDS sequence was cloned into the EcoRI and BamHI (NEB, Beverly, MA, USA) restriction sites of the expression vector pEGFP-N1. The 3 UTR of the NKX3.2 fragment containing the predicted target site was cloned into the psiCHECK2 vector with AsiSI and PmeI (NEB, Beverly, MA, USA). The primer was synthesized by Sangon (Shanghai, China) to introduce the mutation of some bases at the putative binding site of miRNA-326 in NKX3.2. All the primers used are listed in Supplementary Table S1.
HeLa cells were seeded into 24-cell plates and transfected 24 h later. The miRNA-326 was co-transfected with 100 ng of psicheck2-NKX3.2-fragment and 40 pmol of miRNA-326 mimics or negative control (NC) mimics (RiboBio, Guangzhou, China). Thirty-six hours after transfection, the Firefly and Renilla luciferase activities were measured using the Dual-Glo luciferase kit (Promega). The assays were repeated three times.

Measurement of Gene Expression
Semi-quantitative reverse transcription PCR (SqRT-PCR) was performed for the tissue expression profiles as previously described [75], and β-actin was an endogenous control gene.
For validating the differentially expressed mRNAs, lncRNAs, and miRNAs, six lncR-NAs, six mRNAs, and six miRNAs (Table 3) were randomly selected. Quantitative PCR (qPCR) was performed using the SYBR Green Master Mix (TIANGEN, Beijing, China) according to a previously described procedure [64]. The relative expression was normalized using the 2 −∆∆Ct method with β-actin (genes and lncRNAs) and U6 small nuclear RNA (miRNAs) as an internal reference. First, the cycle threshold (CT) values of the target gene were normalized with the CT values of the internal reference gene for all the samples and calibration samples. The calibration sample was a mixture of all the individual sample cDNAs. Second, the ∆CT values of the calibration samples were used to normalize the ∆CT values of the test samples. Finally, the expression-level ratio was calculated. Table 3. Differentially expressed lncRNAs, mRNAs, and microRNAs between different ages in two pig breeds.

Genes
Transcript ID Gene ID/Name Y4_vs_Y1_ log2FC Each biological duplicate consisted of three technical replicates. The miRNA primers were designed and synthesized by RiboBio (Guangzhou, China), who treat primer sequences as confidential information. All the primers used for this section are listed in Supplementary Table S2.

Statistical Analysis
The statistical analysis results are shown as the means ± standard error (SE). The significant differences among the groups were identified through one-way analysis of variance (ANOVA). A comparison of the different groups was made using Duncan's new multiple range test in the R package agricolae (version 1.2-8), R version 3.3.3. The differences were considered statistically significant at p-values < 0.05.

Conclusions
In this study, we profiled the expressions of lncRNAs, miRNAs, and mRNAs at twotime points of TIC development (1 and 4 months old) in Yorkshire and Wuzhishan pigs. Our analysis suggests that several lncRNAs, miRNAs, and genes are involved in critical biological processes associated with bone development. We also constructed the lncRNA-miRNA-gene interaction networks based on the transcription profiles derived from pig TIC. Moreover, the interaction between miRNA-326 and NKX3.2 was identified as being involved in the regulation of chondrocyte hypertrophy differentiation. Therefore, our results contribute to the understanding of the underlying molecular mechanisms of pig body length variation and provide new insights and valuable information for the study of body length variations in other animals from the perspective of vertebral column development.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/ijms24087257/s1. Author Contributions: M.F. designed the study; Q.X., Y.L., Z.C., X.L., Q.T., K.W. and S.T. collected the samples; Q.X. generated the data. M.F., Q.X. and Y.L. analyzed the data. Q.X., M.F. and J.Z. wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement:
The experimental procedures were approved by the Animal Welfare Committee in the State Key Laboratory for Agrobiotechnology at China Agricultural University (approval number XK257), and this study was carried out in strict accordance with the guidelines and regulations established by this committee.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data are available under project number PRJNA949733.