Differences in the expression profiles of lncRNAs and mRNAs in partially injured anterior cruciate ligament and medial collateral ligament of rabbits

Animals and partially injured models of ACL and MCL

Three-month-old healthy male New Zealand white rabbits (n = 6) were purchased from Qingdao Kangda Biological Technology Co, Ltd. (Qingdao, China). Rabbits were bred at the Department of Experimental Animals of China Medical University (Shenyang, China) at 23 ± 2 °C with 12-hour light-dark cycles. The Medical Ethics Committee China Medical University approved the animal experiments on October 22, 2020 (approval number CMU2020310).

Animals were randomly divided into two groups with three rabbits in each group. Rabbits in the control group were bred under normal conditions and labeled N1, N2, and N3 respectively. The following treatment was done for the injured group: (1) Surgery to establish bilateral partially injured models of ACL and MCL: A longitudinal incision was made on the medial side of the knee joint of the rabbit, with an incision length of approximately 1.5 cm. The skin was incised layer by layer, and the subcutaneous tissue and joint capsule were bluntly separated to fully expose the ACL and MCL. ACL and MCL were laterally punctured by a 20ml needle at the proximal 1/3 of the femur and subsequently lacerated by a needle, resulting in blunt laceration of most (approximately2/3) of the ligament. The wounds were cleaned with physiological saline and the joint capsule and skin were sutured. (2) Postoperatively, ceftriaxone sodium 75 mg/kg was administered intramuscularly once daily for 3 days to combat the infection. Rabbits in the injured group were labeled as I1, I2, and I3. Rabbits in both groups were sacrificed to collect ligaments after deep anesthesia at 1 week after modeling. For all ligaments, snap-frozen in liquid nitrogen within 15min of surgery and stored at −80 °C. Two normal ACLs and MCLs and two injured ACLs and MCLs were randomly selected for RNA-Seq, and the rest of ligaments were used for qRT-PCR.

RNA extraction, library construction, and RNA-Seq

Total RNA of each sample was extracted from the ligament tissue as per the instruction manual of the TRlzol Reagent (Life technologies, Carlsbad, CA, USA). RNA concentration was checked using Nanodrop2000 (Thermo Fisher Scientific, Waltham, MA, USA), and RNA integrity was checked using Agient2100 (Agilent Technologies Inc., Santa Clara, CA, MA, USA) and LabChip GX (PerkinElmer Inc., Waltham, MA, USA). RNA contamination was monitored on agarose gels. The cDNA libraries were constructed using the Ribo-off rRNA Depletion Kit after the samples passed the quality control, and then sequenced using Illumina NovaSeq 6000 platform.

RNA-seq analysis

The high-quality reads were aligned to the Oryctolagus cuniculus reference genome (v2.Oryctolagus_cuniculus.v2.genome.fa) using HISAT2 software. StringTie (1.3.1) was used to calculate FPKMs of both lncRNAs and coding genes in each sample after completion of the comparison analysis. Differential expression analysis of two groups was performed using normalized count data which were derived using the DESeq2 package. Fold Change (FC) indicates the ratio of expression between two groups, and False Discovery Rate (FDR) is obtained by correcting the p-value for the significance of the difference (p-value) using Benjamini and Hochberg’s approach. Genes with FDR < 0.01 and the absolute value of FC ≥ 2 were assigned as differentially expressed.

Gene functional annotation

Gene function was annotated based on the following databases:Nr (NCBI non-redundant protein sequences); Pfam (Protein family); KOG/COG (Clusters of Orthologous Groups of proteins); Swiss-Prot (a manually annotated and reviewed protein sequence database); KEGG (Kyoto Encyclopedia of Genes and Genomes); GO (Gene Ontology). GO enrichment analysis was performed using the topGO R package, and Kobas software was used to test the statistical enrichment in the KEGG pathway.

Construction of PPI network

The sequences of the DEGs were blast (blastx) to the genome of a related species (the protein-protein interaction of which exists in the STRING database: http://string-db.org/) to get the predicted PPI of these DEGs. Cytoscape software (version 3.8.2; https://cytoscape.org) was then used to build a PPI network and identify hub genes.

Construction of lncRNA-miRNA-mRNA network

Competitive endogenous RNA (ceRNA) (Salmena et al., 2011) has gained attention in recent years as a new mode of transcriptional regulation. LncRNAs can be used as a miRNA sponge to compete with miRNAs and inhibit the regulation of miRNAs on the target gene, thus indirectly regulating gene expression (Xiao et al., 2019). Our study obtains candidate ceRNA relationship pairs by targeting relationships of miRNAs that satisfy the following conditions. (1) the number of identical miRNAs between ceRNAs should be greater than 5. (2) the p-value of the hypergeometric test is less than 0.01, and the corrected FDR value is less than 0.01. (3) The results of co-expression analysis are taken into account to obtain a ceRNA co-expression network where both ceRNAs are co-expressed with each other. Based on the ceRNA network, the relationship pairs that differed in all three for each differential combination were extracted and visualized using Cytoscape software. The top 50 nodes in each group were identified using the cytoHubba plug-in in Cytoscape and Degree method mapped with Cytoscape.

Quantitative real-time PCR

Total RNA of each sample was extracted from the ligament tissue according to the instruction manual of the TRlzol Reagent (Ambion, Massachusetts, USA). The PrimeScript II 1st Strand cDNA Synthesis Kit (Takara Bio inc., Kusatsu, Japan) and SYBR Premix ex Taq II (Takara) were used for reverse transcription and quantitative real-time PCR (qRT-PCR) assay. Polymerase chain reaction primer sequences are listed in Table S1 (Sangon Biotech Co., Ltd., Shanghai, China). β-Actin was applied as an internal control. For quantitative analysis of differential expression, data were processed by the 2^−ΔΔCt method and the expression of all lncRNAs and mRNAs were presented as relative fold changes to β-actin, which was repeated three times for each sample.

Statistical analysis

All statistical data are expressed as the mean ± standard error of mean (SEM). Differences were analyzed by t-test using Graphpad Prism 8 statistical software. P-value < 0.05 was considered statistically significant.

Data availability statement

All sequences are accessible via SRA using the accession PRJNA751479. Data is available at NCBI SRA, accession numbers: SRR15858594 to SRR15858601.

Statistical results of RNA sequencing data

The statistical results of RNA-seq data are shown in Table S2. A total of 31,217 genes were detected, and 10,903 new genes were discovered, 2,466 of which were functionally annotated. On the basis of the comparison results, gene expression analysis was carried out. Based on the expression of genes in different samples, 1,233 DEGs were identified, and functional annotation and enrichment analysis were performed. A total of 26,242 lncRNAs and 1,263 DELs were identified, and the obtained lncRNAs were divided into four groups with the following classification and percentage (Fig. S1): lincRNA (40%), antisense_lncRNA (8.8%), intronic_lncRNA (48.8%), and sense_ lncRNA (2.4%).

DEGs screening, gene ontology, and Kyoto encyclopedia of genes and genomes pathway analysis

Differences between normal and injured groups

In total, 267 mRNAs were differentially expressed in the injured ACL compared with normal controls, of which 165 were up-regulated, and 102 were down-regulated (Figs. 1A and 1B). DEGs were enriched for several GO pathways such as cellular adhesion, positive regulation of angiogenesis, cellular response to tumor necrosis factor, elastic fiber assembly, extracellular exosomes, proteinaceous extracellular matrix, extracellular matrix, cellular adhesion molecule binding, and fibroblast growth factor binding (Figs. 1C to 1E). The KEGG pathway enriched analysis results showed that DEGs were mainly enriched in the PI3K-Akt signaling pathway, cell cycle, PPAR signaling pathway, ECM-receptor interaction, Cellular senescence, AMPK signaling pathway, p53 signaling pathway, Ferroptosis, TNF signaling pathway, and the IL-17 signaling pathway, which were related to cell proliferation and apoptosis or inflammation (Fig. 1F).

A total of 726 mRNAs were differentially expressed in the injured MCL compared with normal controls, of which 372 were up-regulated, and 354 were down-regulated (Figs. 2A and 2B). DEGs were enriched for several GO pathways such as positive regulation of ERK1 and ERK2 cascade, inflammatory response, acute-phase response, blood vessel remodeling, phagocytosis, engulfment, positive regulation of phagocytosis, leucocyte migration involved in the inflammatory response, extracellular exosomes, cell surface, proteinaceous extracellular matrix, lysosome, extracellular matrix, cellular adhesion molecule binding, collagen binding, cytokine receptor activity, and collagen receptor activity (Figs. 2C to 2E). The KEGG pathway enriched analysis results showed that DEGs were mainly enriched in Cytokine-cytokine receptor interaction, Cell adhesion molecules, and Phagosome, which were related to cell proliferation and apoptosis (Fig. 2F).

Differences between the ACL and MCL groups

In the normal group, 420 mRNAs were differentially expressed in ACL compared to MCL, of which 178 were up-regulated and 242 were down-regulated (Figs. 3A and 3B). DEGs were enriched for several GO pathways such as extracellular matrix organization, cell adhesion, positive regulation of ERK1 and ERK2 cascade, blood vessel development, regulation of fibroblast growth factor receptor signaling pathway, extracellular exosome, proteinaceous extracellular matrix, cell surface, extracellular matrix, collagen type I trimer, collagen type IX trimer, collagen type II trimer, extracellular matrix structural constituent, growth factor activity, cell adhesion molecule binding, platelet-derived growth factor binding, and extracellular matrix binding (Figs. 3C to 3E). The KEGG pathway enriched analysis results showed that DEGs were mainly enriched in PI3K-Akt signaling pathway, Focal adhesion, ECM-receptor interaction, cell adhesion molecules, ErbB signaling pathway, FcγR-mediated phagocytosis, Leukocyte transendothelial migration, and TNF signaling pathway, which were related to cell proliferation, differentiation, migration apoptosis, and inflammation (Fig. 3F).

In the injured group, 162 mRNAs were differentially expressed in ACL compared with MCL, of which 82 were up-regulated and 80 were down-regulated (Figs. 4A and 4B). DEGs were enriched for several GO pathways such as elastic fiber assembly, positive regulation of mesenchymal stem cell differentiation, proteinaceous extracellular matrix, cell surface, elastic fiber, collagen type II trimer, and extracellular matrix structural constituent (Figs. 4C to 4E). The KEGG pathway enriched analysis results showed that DEGs were mainly enriched in ECM-receptor interaction, Focal adhesion, PI3K-Akt signaling pathway, and Cell adhesion molecules, which were related to cell proliferation and apoptosis (Fig. 4F).

DELs screening, gene ontology, and Kyoto encyclopedia of genes and genomes pathway analysis

Differences between normal and injured groups

In total, 329 lncRNAs were differentially expressed in the injured ACL compared with normal controls, of which 209 were up-regulated and 120 were down-regulated (Figs. S2A and S2B). Target genes of DELs were enriched for several GO pathways such as positive regulation of angiogenesis, positive regulation of inflammatory response, proteinaceous extracellular matrix, extracellular matrix, growth factor activity, chemokine activity, and cell adhesion molecule binding (Figs. S2C to S2E). The KEGG pathway enriched analysis results showed that target genes of DELs were mainly enriched in cytokine-cytokine receptor interaction, the cAMP signaling pathway, cGMP-PKG signaling pathway, chemokine signaling pathway, JAK-STAT signaling pathway, and the ErbB signaling pathway, which were associated with cell proliferation, differentiation, and migration apoptosis (Fig. S2F).

A total of 609 lncRNAs were differentially expressed in the injured MCL compared with normal controls, of which 263 were up-regulated, and 346 were down-regulated (Figs. S3A and S3B). Target genes of DELs were enriched for several GO pathways such as positive regulation of cell proliferation, cell division, positive regulation of angiogenesis, growth, platelet-derived growth factor receptor signaling pathway, extracellular exosomes, focal adhesion, proteinaceous extracellular matrix, growth factor activity, and cell adhesion molecule binding (Figs. S3C to S3E). The KEGG pathway enriched analysis results showed that target genes of DELs were mainly enriched in PI3K-Akt signaling pathway, cAMP signaling pathway, Focal adhesion, cell adhesion molecules, ECM-receptor interaction, and ErbB signaling pathway, which were associated with cell proliferation, differentiation, and migration apoptosis (Fig. S3F).

Differences between the ACL and MCL groups

In the normal group, 470 lncRNAs were differentially expressed in ACL compared to MCL, of which 166 were up-regulated and 304 were down-regulated (Figs. S4A and S4B). Target genes of DELs were enriched for several GO pathways such as positive regulation of ERK1 and ERK2 cascade, positive regulation of JNK cascade, regulation of cellular metabolic process, extracellular exosome, focal adhesion, proteinaceous extracellular matrix, stress fiber, growth factor activity, and cell adhesion molecule binding (Figs. S4C to S4E). The KEGG pathway enriched analysis results showed that target genes of DELs were mainly enriched in the PI3K-Akt signaling pathway, MAPK signaling pathway, Focal adhesion, cGMP-PKG signaling pathway, ECM-receptor interaction, and Ferroptosis, which were associated with cell proliferation, differentiation, and migration apoptosis (Fig. S4F).

In the injured group, 205 lncRNAs were differentially expressed in ACL compared to MCL, of which 116 were up-regulated and 89 were down-regulated (Figs. S5A and S5B). Target genes of DELs were enriched for several GO pathways such as acute-phase response, regulation of ERK1 and ERK2 cascade, and proteinaceous extracellular matrix (Figs. S5C to S5E. The KEGG pathway enriched analysis results showed that target genes of DELs were mainly enriched in the cAMP signaling pathway, focal adhesion, cGMP-PKG signaling pathway, ECM-receptor interaction, and the ErbB signaling pathway, which were associated with cell proliferation, differentiation, and migration apoptosis (Fig. S5F).

Identifying candidate genes for PCR-validation

DEGs related to the altered pathway were further searched in the PubMed literature database. “ Gene symbol ”(gene name) was used as a query keyword and searched in the“ title/abstract ”. Three randomly selected up- and down-regulated genes associated with tissue damage repair were validated in each group. Finally, since very little information was available on DELs, three randomly selected DELs in each of the top five up-regulated and top five down-regulated DELs were chosen for further validation (Fig. 5). All the results of qRT-PCR are corresponded to RNA-seq data.

Construction of lncRNA-miRNA-mRNA network

Based on the ceRNA regulatory network, the top 50 nodes in each group (Figs. 6A to 6D) were identified, including lncRNAs, miRNAs, and mRNAs(only 45 nodes in the normal group of ACL compared to the injured group). Each DEGs could be associated with one or more miRNAs and lncRNAs, which contained some genes related to the tissue damage repair process, such as COL7A1, LIF, FGD5, SLC20A1, and their associated lncRNAs and miRNAs were searched for, thus their networks were displayed with Cytoscape (Figs. 6E to 6G). The ceRNA network suggested that lncRNAs could not only participate in the regulation of biological processes independently but also act as miRNA sponges together to affect downstream gene expression and thus participate in the processes after ACL and MCL injury.

Construction of PPI network

We used PPI analysis and online STRING platform to detect the interaction effects between proteins in each group of DEGs, and then using cytoHubba plug-in and MCC method, we identified the top 30 ranked hub genes in each group (Figs. 7A to 7D), which contained several genes related to the tissue damage repair process, such as FGFR2, HBEGF, EPHA2, CSF1, MMP2, MMP9, SOX5, LOX.