Whole Genome Expression Profiling of Semitendinosus Tendons from Children with Diplegic and Tetraplegic Cerebral Palsy

1. Introduction

Tendons are fibrous connective tissues that connect skeletal muscle to bone and transmit the mechanical force generated by muscle contraction to the bones, driving the movement of the skeletal system. Tendons are primarily composed of cells and extracellular matrix (ECM), principally comprised of collagen type I, which accounts for about 65-80% of the dry weight and is organized into parallel fiber bundles arranged into a three-dimensional network surrounded by a thin connective tissue sheath known as the epitenon [1]. Additionally, the ECM contains proteoglycans (e.g. aggrecan, versican, decorin, biglycan, fibromodulin, lumican, keratocan, osteoglycin, syndecans), glycoproteins (e.g. elastin, cartilage oligomeric matrix protein (COMP/TSP5), tenascin-C, lubricin, tenomodulin, fibronectin), glycosaminoglycans, and smaller quantities of other collagens with different functions. Aggrecan and versican are large aggregating proteoglycans, whose main function is to provide compression resistance to the tendon by increasing water content [2]. Decorin, biglycan, fibromodulin, lumican, and osteoglycin are small leucine-rich proteoglycans, which play a role in fibril assembly and tendon integrity [3,4,5], while syndecans are transmembrane heparan sulfate proteoglycans whose roles in tendons have remained elusive [6]. Thus, while collagen 1 fibers are responsible for the tensile strength of tendons, proteoglycans are regulating their viscoelastic properties. Tenoblasts and tenocytes constitute 90-95% of tendon cells and are fibroblast-like cells responsible for the production of collagen and other ECM components [1]. Tenocytes respond to mechanical load of tendons by altering the expression of ECM proteins, thereby modulating tendon structure, composition, and mechanical properties [7].

Cerebral palsy (CP) is a group of disorders of the development of movement and posture, causing activity limitations, which are attributed to non-progressive brain lesions that occurred in the developing fetus or infant [8]. CP is the most common chronic movement disorder in children with a prevalence of 1.5 to 4 per 1000 live births [9,10,11,12,13,14,15], and one of the most disabling and costly chronic disorders of children and adults. Although CP is the result of a non-progressive lesion of the developing brain, it leads to progressive alterations of the musculoskeletal system, which can manifest with different clinical presentations, such as spasticity, dyskinesia, dystonia, hypertonia, often leading to the development of fixed contractures, i.e. a permanent shortening and stiffening of the muscle-tendon unit, resulting in the loss of joint motility (reviewed in [16,17]). Patients are diagnosed in respect to the anatomical distribution of their deformity as monoplegic, hemiplegic, diplegic, or tetraplegic (reviewed in [16]). Involvement of only one limb is referred to as monoplegia, unilateral involvement of an upper and lower limb is referred to as hemiplegia, predominant lower limb involvement is referred to as diplegia, involvement of the lower limbs and one upper limb is referred to as triplegia, and involvement of all four limbs and the trunk is referred to as quadriplegia. Various treatment strategies, including physiotherapy, pharmacological interventions, neurectomy, and orthosis, are used to reduce spasticity, limit pain, and prevent contractures [16,17]. Nevertheless, most CP patients eventually develop fixes contractures, necessitating tendon lengthening surgery to restore the joint range of motion and relieve symptoms. In previous studies, our coworkers found that tendons of CP patients show ECM remodeling [18] and altered expression of selected genes related to collagen turnover [19,20]. However, the effect of CP on genome-wide gene expression has not been studied and it remains poorly understood how the exposure of tendons to the spasticity-induced increased mechanical loading and functional demands in CP patients affect tendon homeostasis at the molecular level. In the present study we performed total RNA-sequencing (RNA-Seq) analysis on semitendinosus tendons from diplegic and tetraplegic CP patients undergoing tendon-lengthening surgery compared to age-matched controls undergoing cruciate ligament (ACL) reconstruction.

2. Materials and Methods

2.1. Sample collection and preparation

A cohort of 30 children were included in the study: 12 children with diplegic cerebral palsy (10 males and 2 females; mean age: 13.9 years), 13 children with tetraplegic cerebral palsy (8 males and 5 females; mean age: 12.6 years), and 5 typically developing pediatric patients as control patients (4 males and 1 female; mean age: 17.9). The characteristics of the subjects are listed in Table 1. Semitendinosus tendon biopsies were obtained during surgery of cerebral palsy patients undergoing tendon lengthening or control patients undergoing ACL reconstruction. Biopsies were immediately frozen on dry ice and stored at –80°C for further processing.

2.2. RNA extraction

For extraction of tendon RNA, tendons were first powderized with a hammer on dry ice after which TRIzol Reagent (Thermo Fisher Scientific) was added and tendons were homogenized using a TissueLyser II (Qiagen). RNA was subsequently extracted using the RNA-ZOL Direct Clean-up Kit (Fisher Molecular Biology) following the instructions of the manufacturer. After DNase I treatment, RNA purity and concentration were evaluated using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific) and a Qubit Fluorometer (Thermo Fisher Scientific). RNA quality was analyzed on an Agilent Tape Station 4200 system and the RNA integration number (RIN) varied from 2.7 – 9.1.

2.3. RNA-Seq analysis

Indexed sequencing libraries were generated from 140 ng of tendon RNA using the KAPA RNA HyperPrep with RiboErase (Roche), which allows for library construction from partially degraded samples. Paired-end multiplexed sequencing of libraries to generate reads of 150 bp (PE 2x250) was performed on a NovaSeq 6000 instrument (Illumina Inc.). 22-40 million reads per sample were obtained. RNA-Seq samples were first demultiplexed and FASTQ files were created from BCL files using bcl2fastq v2.17.1.14 software (Illumina Inc.). Quality control and assessment were performed on FASTQ files using FastQC v0.11.9; http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Paired-end reads of 150 bp were aligned to the Gencode Human reference genome (build GRCh38) using STAR v2.7.2b [21]. Annotation of genes was performed using the biomaRt package v2.54 [22] based on annotations from Ensembl release 110. The raw read counts were normalized with TMM implemented in the edgeR v3.40 package [23] in R/Bioconductor [24] and lesser expressed genes were filtered out with the filterByExpr (min.count=6) function. Differential expression analysis of read counts was performed using the voom, lmFit, and eBayes (robust=T) functions of the limma v3.46 package [25] in R. Significantly differentially expressed genes were chosen based on a false discovery rate (FDR) < 0.1 and |log2 fold change (FC)| ≥ 0.4. The Sva v3.38 package [26] was used to estimate artifacts (n.sv=1) and correct the CPM values. Hierarchical clustering of significantly modulated genes was performed using hclust and dist functions on sva corrected log2CPM in R. Clustering was performed with the ward.D2 method and Euclidean distance to generate a heatmap using pheatmap (https://cran.r-project.org/web/packages/pheatmap/index.html) and scaling the rows. Principal components analysis (PCA) was performed using the PCA function of the FactoMineR v2.7 [27] package in R. The Venn Diagrams were obtained using the R package VennDiagram (https://cran.r-project.org/web/packages/VennDiagram/index.html).

Raw data files for the RNA-Seq analysis have been deposited in the NCBI BioProject database (https://www.ncbi.nlm.nih.gov/bioproject/) under accession code PRJNA1004310.

2.4. Gene ontology and pathway analysis

For gene enrichment analysis, the gene ontology (GO) Biological Process 2021 and Reactome 2022 databases within EnrichR [28] were used to find significant enriched GO terms and pathways with an adjusted p value ≤ 0.05.

3. Results

Total RNA-Seq analysis was performed on semitendinosus tendons from 12 diplegic, 13 tetraplegic and 5 control patients. Unfortunately, due to muscle contamination 5 control samples had to be excluded. Principal component analysis (PCA) (Figure 1A) and hierarchical clustering showed clustering of the control group distant from the CP groups (Figure 1A,B). 1170 genes (676 upregulated and 494 downregulated) were differentially expressed in tetraplegic vs. control patients, 112 in diplegic vs. control patients (61 upregulated and 51 downregulated), and 101 (70 upregulated and 31 downregulated) in tetraplegic vs. diplegic patients. As illustrated in the Venn diagram in Figure 1C, 101 genes (58 upregulated and 43 downregulated) were commonly altered in diplegic and tetraplegic patients, while 1068 and 11 genes were differentially expressed specifically in tetraplegic vs. control patients and diplegic vs. control patients, respectively. Thus, most genes that were differentially expressed in diplegic patients were altered also in tetraplegic patients. The full lists of differentially expressed (DE) genes are shown in Table S1-S3. For genes that were differentially expressed in tetraplegic vs. control patients, gene ontology analysis for biological processes as well as Reactome pathway analysis were performed (Figure 2 and Tables S4–S7). For upregulated genes, there was an enrichment for genes involved in collagen formation, ECM organization, and categories related to vesicle-mediated transport (Figure 2A,B, left), while for downregulated genes, genes related to cytokine and apoptotic signaling, including TNF-related apoptosis-inducing agent (TRAIL) signaling, were most significantly affected (Figure 2A,B, right). The differentially expressed genes in diplegic patients did not show significant enrichment for genes involved to any specific process. Differentially expressed genes in tetraplegic CP patients vs. control patients involved in ECM organization, vesicle-mediated transport, cytokine signaling, and apoptotic signaling are listed in Figure 3. Many genes encoding collagens or involved in collagen formation and organization were upregulated, including proteoglycans, glycoproteins, integrins, and growth factors (Figure 3A). Increased expression of several genes encoding metallopeptidases (MMPs), involved in ECM turnover, was also found, including MMP2, which is the most expressed MMP in tendons. Furthermore, many genes involved in membrane trafficking were upregulated (Figure 3B). Among the downregulated genes were many genes involved in cytokine signaling, including interferon signaling, interleukin signaling, and chemokines (Figure 3C) as well as a number of genes related to TRAIL signaling and apoptosis (Figure 3D).

4. Discussion

The present study is the first to report genome-wide expression profiling on tendons from CP patients. Consistent with our previous findings [19], tetraplegic CP patients showed increased expression of genes encoding collagens and ECM proteins. The most upregulated genes were the collagen 1-encoding genes COL1A1 and COL1A2, which were 15- and 8-fold upregulated, respectively. Furthermore, genes encoding both fibrillar (COL5A1/2, COL11A1/2) and non-fibrillar collagens (COL6A1/2/3 (beaded filament forming collagen), COL8A1/2 (hexagonal network forming collagen), COL10A1 (hexagonal network forming collagen), COL14A1 (fibril-associated collagen with interrupted triple-helix (FACIT)), COL15A1 (FACIT)) were highly upregulated as shown in Figure 3A [29].

Collagen is composed of a triple helix, consisting of three α chains with the repeating (Gly-X-Y)n amino acid sequence, where X and Y are frequently proline or hydroxyproline but can be any amino acid [29]. While some collagens contain three identical α chains, others contain two or three different chains referred to as α1, α2, and α3. In tendons, fibrillar collagens are bundled into fibrils, which form larger fibers. During collagen biosynthesis, collagen undergoes post-translational modification by hydroxylation of proline residues into hydroxyproline by prolyl 3-hydroxylases (P3H1, P3H2, P3H3, P3H4) and prolyl 4-hydroxylases (a tetramer composed of two α (P4H1, P4H2, P4H3) and one β (P4HB) subunit) as well as hydroxylation of lysine into hydroxylysine by lysyl hydroxylases (LH1, LH2, LH3), encoded by procollagen-lysine,2-oxoglutarate 5-dioxygenases (PLOD1, PLOD2, PLOD3) [30,31]. Hydroxylysine residues are often subject to glycosylation and P4HB is responsible for the formation of disulfide bands between the α chains. Peptidyl-prolyl cis-trans isomerases (FKBP10, FKBP14, PPIB) catalyze the cis-trans isomerization of peptide bonds, which is a rate limiting step for triple helix formation. After secretion, lysyl oxidases (LOX, LOXLs) oxidize hydroxylysine residues to induce cross-link formation [32], which provide the collagen fibrils with mechanical stability and tensile strength as well as contribute to the stiffness of the collagen fibril. In tendons of tetraplegic CP patients, genes involved in all steps of collagen formation and cross-linking were upregulated, including P3H3, P3H4, P4HA1, P4HA3, P4HB, PLOD1, FKBP10, FKBP14, PPIB, LOXL1, LOXL2, and LOXL3, while only PPIB was increased in diplegic CP patients. This suggests that collagen synthesis is strongly induced in tetraplegic CP patients.

In tetraplegic CP patients, genes encoding the small leucine-rich proteoglycans biglycan (BGN), keratocan (KERA), and osteoglycin (OGN) were highly upregulated by 2.4-, 7.4- and 4.7-fold, respectively. Biglycan was recently shown to be involved in the maintenance of tendon structure and mechanics in mature tendons [33] and osteoglycin has been implicated in collagen fibrillar organization and tendon mechanical function [5,34], while the role of keratocan has remained elusive [35]. Decorin, which constitutes about 80% of the total amount of small leucine-rich proteoglycans in tendons, was not affected. This may be explained by a recent study, which demonstrated that while decorin plays a more important role than biglycan (the second most abundant small leucine-rich proteoglycan in tendons) in the modulation of collagen fibril structure and viscoelastic mechanics during tendon development, biglycan is more important for the maintenance of tendon structure and mechanical properties during homeostasis in mature tendons [33]. While the expression of none of the genes encoding small leucine-rich proteoglycans was altered in diplegic CP patients, SDC2, encoding the transmembrane heparan sulfate proteoglycan syndecan-2, was upregulated in both diplegic and tetraplegic CP patients by 1.7- and 1.6-fold, respectively.

Several genes encoding glycoproteins, including tenascin-C (TNC), fibronectin 1 (FN1), tenomodulin (TNMD), secreted protein acidic and rich in cysteine (SPARC), and thrombospondin 4 (THBS4) were upregulated in tetraplegic CP patients. In particular, SPARC and THBS4 were strongly upregulated by 6.0- and 5.6-fold, respectively. SPARC has been shown to be essential for load-induced tendon tissue maturation and homeostasis, affecting ECM composition and tendon biomechanical properties [36], while tenomodulin plays a role in tenocyte proliferation and collagen fibril maturation [37]. Thrombospondin 4, encoded by THBS4, which was 2.9-fold upregulated, is involved in the organization of collagen fibrils [37], while tenascin C and fibronectin, whose genes were 3.3- and 2.2-fold upregulated, respectively, both contribute to the mechanical stability of the ECM [38,39]. Among the glycoproteins, only TNC was upregulated (2.2-fold) in diplegic CP patients. In tetraplegic CP patients, alterations in the expression of genes encoding different subunits of the large multidomain heterotrimeric glycoprotein laminin was also found. Laminin is located in the basement membrane, where it interacts with collagen type IV, integrins, and dystroglycans, and play important roles in cell adhesion, differentiation, and migration [40]. LAMA1 (encoding laminin α1) and LAMB2 (encoding laminin β2) were 3.3- and 1.5-fold upregulated, respectively, while LAMC2 (encoding laminin γ2) was 3.3-fold downregulated.

Altered mRNA expression of several integrins was also found in the tendons of tetraplegic CP patients, including ITGA2, ITGB1, and ITGBL1, which were 3.5-, 1.4,- and 2.2-fold upregulated, respectively as well as ITGA3, which was 1.8-fold downregulated. Integrins are heterodimeric transmembrane receptors, which play a major role in linking the ECM to the cytoskeleton [41]. It is believed that integrins can sense and transmit mechanical stimuli from the ECM to tenocytes, thereby triggering intracellular signaling pathways leading to adaptive regulation of gene expression. In particular, mechanical stretch was found to activate the AKT/mTOR pathway via β1 integrin, thereby regulating collagen expression [42]. Also, GJA1, encoding connexin 43, a gap junction protein with an important role in the communications between tenocytes allowing for the passage of passage of free metabolites and ions [43], was 2.7-fold upregulated in tetraplegic CP patients.

The expression of number of genes encoding growth factors were increased in tendons of tetraplegic CP patients, including integrin like growth factor 1 (IGF1), integrin like growth factor binding protein 3 (IGFBP3), fibroblast growth factor 12 (FGF12), transforming growth factors (TGFB1, TGFB2, TGFB3), and vascular endothelial growth factor B (VEGFB), while FGF11 was downregulated (see Figure 3A). IGF1 has been shown to be required for adult tendon growth in response to increased mechanical overload through stimulation of tenocyte proliferation and protein synthesis [44]. Similarly, TGFβ was shown to promote collagen synthesis and matrix remodeling during tendon healing, thereby enhancing mechanical strength [45]. On the other hand, excessive TGFβ activation as a result of mechanical overload or repetitive mechanical loading can lead to tendinopathy [46]. VEGF promotes angiogenesis and increases vascular permeability and is important for neovascularization during tendon healing [47]. Also, the MKX gene encoding Mohawk, an atypical homeobox transcription factor involved in postnatal tendon maturation [48,49] and critical for the tendon response to mechanical stimuli [50], was 2.6-fold upregulated in tetraplegic CP patients.

Among the differentially expressed genes in the tendons of tetraplegic CP patients were many genes encoding metalloproteases, including matrix metalloproteinases (MMP2, MMP10, MMP14, MMP16) and members of the ADAMTS (a disintegrin and metalloproteinase with thrombospondin motifs) protein family (ADAMTS2, ADAMTS3, ADAMTS17, ADAMTSL2). Metalloproteases are important for the turnover of the ECM, which is a finely balanced, dynamic process of protein synthesis and degradation taking place at a low rate during homeostasis and at higher rates during conditions of inflammation, tissue damage, or increased mechanical load [51]. Metalloproteinases are mediating the proteolytic degradation of components of the ECM and are regulated by tissue inhibitors of metalloproteinases (TIMS). While ECM turnover allows for tendon adaptation in response to altered mechanical loading conditions, an imbalance in the activities of MMPs and TIMS can lead to pathological conditions [52]. MMP2, MMP14, and MMP16, whose mRNA levels were 3.5-, 2.0-, and 3.6-fold upregulated, respectively, are able to cleave fibrillar collagen I, while MMP10, whose gene was 7.4-fold downregulated, does not cleave collagen I [52], but degrades various components of the ECM, including proteoglycans, fibronectin, and collagen III-V [53] as well as activates other MMPs [54]. ADAMTS2 and ADAMTS3, which were 3.4- and 3.3-fold upregulated at the mRNA level, respectively, are involved in fibrillar collagen maturation though cleavage of the N-terminal propeptide of procollagens [55,56]. Also, PCOLCE, encoding procollagen C-endopeptidase enhancer, which enhances the activity of procollagen C proteinases cleaving the C-terminal propeptide of procollagens, was 2.6-fold upregulated [57]. In addition, mRNA encoding ADAMTS17, whose function in tendons has remained elusive, was 2.5-fold upregulated. On the others hand, transcript levels of ADAMTSL2, which belongs to the subfamily of ADAMTS-like proteins without catalytic activity and is involved in the modulation of microfibril formation, were 2.1-fold reduced [58]. Among the metalloproteases, only MMP2 was 2.1-fold upregulated in diplegic CP patients.

Consistent with increased protein turnover and remodeling in tetraplegic CP patients, Reactome analysis and gene ontology analysis for biological processes revealed an enrichment for upregulated genes involved in vesicle-mediated transport (Figure 2A,B, left). Among the downregulated genes, there was an enrichment for genes involved in interferon and cytokine signaling, including a number of chemokines (Figure 2A,B, right, Figure 3C). In particular, CCL5, CCL16, CXCL10, and CXCL11 were strongly downregulated by respectively 3.5-, 2.1-, 5.6-, and 4.4-fold, while CXCL16 was 1.7-fold upregulated. Thus, the higher mechanical load experienced by tendons of CP patients does not appear to be associated with inflammation. Furthermore, apoptosis does not seem to be activated as there was an enrichment for downregulated genes associated with TRAIL signaling and apoptotic processes (Figure 3D).

The results of our RNA-Seq analysis are consistent with increased collagen synthesis and turnover in tendons from tetraplegic CP patients as a consequence of the spasticity-induced chronic mechanical stimulation. Many of the changes in gene expression observed in CP patients are similar to those observed in the tendons of patients with tendinopathy, but there are important differences suggesting that the altered gene expression in tetraplegic CP patients may be adaptive rather than pathological as in tendinopathy. Like in tendinopathy patients the expression of many collagens was increased. However, while the mRNA expression of a number of collagens was 1.6-2.7-fold upregulated in tendinopathy patients, COL1A1, COL1A2, COL10A1, and COL11A1/2 were 5.3-15.2-fold upregulated in tetraplegic CP patients, suggesting strongly increased collagen synthesis. In regard to the expression of metalloprotease-encoding genes, increased expression of MMP2, MMP14, MMP16, ADAMTS2, ADAMTS3 as well as downregulation of MMP10 has also been found in tendinopathy [59,60]. Also, upregulation of genes encoding proteoglycans (BGN) glycoproteins (TNC, FN1, SPARC), integrins (ITGB1), and growth factors has been reported in tendinopathy patients [60,61,62]. On the other hand, we found no alterations in genes associated with Notch3, Wnt, and JAK/STAT signaling, which were reported to be increased in tendinopathy patients ([60]). In addition, although the role of inflammation in tendinopathy remains controversial [63,64,65,66], a number of genes encoding pro-inflammatory cytokines have been reported to be upregulated in tendinopathy, whereas many chemokines and interleukins were downregulated in tetraplegic CP patients [60,62] (see Figure 3C). Moreover, whereas increased expression of markers of apoptosis was found in tendinopathy [62], many apoptosis-related genes were downregulated in tetraplegic CP patients. In addition, an enrichment for upregulated genes associated with vesicle-mediated transport has not been reported in tendinopathy patients. Overall, our results indicate that tetraplegic CP patients show increased ECM turnover with an increased net synthesis of collagen to allow for reorganization of the ECM as an adaptive response to increased mechanical loads and functional demands as a result of spasticity. This is similar to observations after exercise training, where ECM remodeling has been found to result in reduced tendon stress and increased tendon stiffness and tensile strength [67,68,69,70,71,72,73], likely generating more resistant tendons to support higher loads. The absence of pro-inflammatory cytokines and markers of apoptosis in CP patients compared to tendinopathy patients could be due to the difference in stress stimuli, being sub maximal and continuous in the former and over maximal and discontinuous in the ladder. Although tendon remodeling in CP patients does not appear to be directly pathological, increased tendon stiffness and shortening in CP patients as a result of continuous mechanical stimulation is an important issue in clinical practice, where tendon extension surgery is often required to restore joint mobility. In clinal practice, various kinds of muscle relaxants and neuromuscular blocking agents are also often used to reduce spasticity and consequently decrease pathological stimulation of tendons and prevent the development of contractures [16,17]. In future studies it would be interesting to perform mechanical tension studies and histological analyses of patients with different severity of CP to determine whether they show increased tendon stiffness and cross-sectional area as observed in tendons from long distance runners [74,75] and after high resistance training [68,69,71,72,76]. The much fewer gene expression changes in diplegic compared to tetraplegic CP patients implies that tendons of patients with a less severe form of CP are not affected by the disease to an extent that induces major tendon remodeling, consistent with the fact that only intense exercise is associated with substantial tendon remodeling [68,74,76]. This is also consistent with our personal observations in the clinic, where increased tendon stiffness is generally more evident in tetraplegic CP patients compared to diplegic CP patients.

In conclusion, the present study provides new insights into the alterations in gene expression in CP patients as a consequence of the chronically increased mechanical loads due to muscle hyperactivity. A better understanding of the molecular alterations in the tendons of CP patients during the course of the disease may be beneficial for disease monitoring and lead to better treatment strategies to improve quality of life in CP patients.