**Preface to "Intervertebral Disc Regeneration II"**

This book is devoted to the thriving research on the regeneration of intervertebral discs. It is a colorful mixture of different applied science fields with the shared aim of reducing the number of lower back pain surgeries worldwide. I thank all contributors for submitting their manuscripts to this Special Issue who made this project a success.

> **Benjamin Gantenbein** *Editor*

## *Editorial* **New Developments on Growth Factors, Exosomes, and Single Cell RNA-Sequencing for Regeneration of the Intervertebral Disc**

**Benjamin Gantenbein 1,2**


#### **1. Introduction**

Low back pain (LBP) is the number one cause of disability worldwide, with incidences increasing exponentially [1–3]. A recent study estimates that by the year 2050, an increase of 200 million people is expected, with a current peak at 619 million people [1]. This Special Issue targets the specific niche of finding innovative ways to address the clinical problem of LBP, which is often induced by the prolapse of the spinal column caused by genetic or epigenetic factors. Intervertebral disc (IVD) degeneration is often believed to be the root cause of chronic pain. Future research aims to understand the contribution of metabolic factors such as nutrition, besides other risk factors such as smoking and *Diabetes mellitus* [1]. This Special Issue provides a heterogeneous snapshot of recent applied research on IVD and LBP, ranging from cell biology studies to artificial intelligence in diagnostics. In the following subsections, I provide a short overview and summarize the seven articles' main findings.

#### **2. Wet Laboratory Studies in the Second Special Issue**

Two original studies focused on the single-cell transcriptomics to characterize phenotypes in rats in an IVD degeneration model [4] or performed a pathway study of ~100 genes using qPCR gene arrays using total RNA extractions from human donors suffering from idiopathic skeletal hyperostosis (DISH). The first study established an in vivo IVD degeneration model in 8-week-old Sprague Dawley rats that underwent surgery for retroperitoneal exposure using a 27 Gauge needle of the L4-L6 lumbar spine. Rohanifar et al. (2022) [4] then allowed the rats to recover for two and eight weeks postoperatively. Then, they digested the single cells from the tissue using mild digestion protocols, extracted and sequenced the total RNA and compared the next-generation sequencing (NGS) data relative to untreated controls. Rohanifar et al. (2022) confirmed that the nucleus pulposus (NP) mainly expressed key markers such as CD24 [5] as well as aggrecan and collagen type 2 [6–8]. In the outer part of the IVD, the cells mainly expressed collagen type 1, as previously identified in rats [6–8] but also in other species such as bovine and human [9]. Furthermore, their data allowed us to distinguish IVD cells from lymphoid, endothelial and myeloid cells [4]. The needle-punctured groups subsequently had a significantly higher amount of myeloid cells and lymphocytes than controls. These data are logical since needle puncture causes inflammation [4]. The second transcriptome study is on DISH patients. DISH is also known as Forestier's disease [10–13]. To the best of my knowledge, the study by Gantenbein et al. (2021) [14] is the first to elucidate on the possible phenotypic changes and deregulations of DISH cells compared to the phenotype of IVD cells isolated from trauma patients. This comparison has its limitations, of course, as trauma cells are not necessarily "healthy" cells. However, it most likely that these cells are still in a better state than cells

**Citation:** Gantenbein, B. New Developments on Growth Factors, Exosomes, and Single Cell RNA-Sequencing for Regeneration of the Intervertebral Disc. *Appl. Sci.* **2023**, *13*, 7346. https://doi.org/ 10.3390/app13137346

Received: 12 June 2023 Accepted: 19 June 2023 Published: 21 June 2023

**Copyright:** © 2023 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

isolated from degenerated IVDs. Nevertheless, by comparing the expression levels of bone morphogenic protein pathway cytokines and their inhibitors, i.e., Gremlin, Noggin and Chordin, my group found that these were dys-regulated [15], although these changes were not significant. However, DISH-IVD, in contrast to IVD obtained from trauma, showed a significant up-regulation of early growth response 2 (EGR2) interleukin 6 (IL6), and insulin-like growth factor 1 (IGF1) tended to be up-regulated [14]. IGF1 has been proposed by other authors to be a possible marker in serum samples of DISH patients [16,17].

The third study in the Special Issue is that of Bischof et al. (2021) [18], who focused on cell culture and tissue-specific progenitor cells. They presented a study on the optimization of culture conditions on the so-called nucleus pulposus progenitor cells (NPPC). In this study, primary IVD cells were isolated from human disc tissue with a mild digestion protocol [19,20]. After reaching confluence in monolayer recovery, these mixed IVD cells, mainly from the nucleus pulposus (NP), were then trypsinized and sorted with a surface marker named Tie2 (or CD202b), which stands for angiopoeitin receptor-1. The Tie2+ cells and the Tie2- cells were then further cultured in normoxic and hypoxic (i.e., 2%) conditions in presence of Angiopoetin-1 (Ang-1) or Ang-2 at increasing doses [21]. However, the results of the study [18] did not produce the expected results, namely, that Tie2+ were stimulated and would proliferate faster compared to Tie2− cells. Despite this, it seemed very clear that hypoxia, i.e., at 2% oxygen level, was the most important factor for higher cellular metabolic activity. Their conclusion that hypoxia is beneficial for these NPPC is in agreement with other studies performed by the team of Daisuke Sakai from Tokai University School of Medicine [22,23].

#### **3. Clinical/Radiological Studies in the Second Special Issue**

There are two radiological studies published in this Special Issue. Firstly, Landauer and Trieb (2022) [24] provide radiological evidence that the lumbosacral transitional vertebrae (LSTV) are valid as a model for IVD regeneration. They scanned 1482 patients radiologically, and their LSTV were then classified according to Castellvi classification type II–IV [24]. Additionally, magnetic resonance scans (MRIs) were also obtained from selected patients. The authors concluded that the reduced or absent mobility in the LSTV segments led to an overload of the adjacent segments in these patients.

The second study, by Kim et al. (2022) [25], is on the usage of natural language processing (NLP), which is defined as understanding, analyzing, and extracting meaningful information from text (natural language) by computer science [26]. This research targets a highly significant area of research, which is "big data". It is obvious that artificial intelligence (AI) will be necessary to make full use of all the available clinical data and to help surgeons to take decisions with the assistance of fast data processing. In this approach, the authors tested their NLP pipeline on a balanced sample of 300 X-ray, 300 CT, and 300 MRI reports. When evaluating their NLP model performances, four parameters—recall, precision, accuracy, and a so-called "F1" score (the harmonic mean of precision and recall [27])—were greater than 0.9 for all 23 radiologic findings.

#### **4. A Review on Secretomes of the IVD**

Extracellular vesicles (EVs) have long attracted the attention of the regenerative community. The importance of this topic is underlined by a current "wave" of comparable articles that also focus on the usage of EVs to regenerate the IVD [28–31]. This is not altogether surprising as regulatory hurdles toward proving non-toxicity and patient safety have recently been introduced by authorities such as the Federal Drug and food Agency (FDA) and the label of the Conformité Européene (CE) [32]. This applies in particular to cellular applications. Thus, secretomes, or so-called conditioned media, have been the focus of recent IVD-related research [29,33]. The review by Tilotta et al. (2023) provides a valuable insight into the field of EVs and a summary of their characterization [29].

#### **5. Conclusions**

Overall, this Special Issue offers a good insight into the heterogeneity of IVD research and the recent findings not only from clinics, but also from biologists and engineers. I hope that this Special Issue will give scientists an overview of this highly translational and applied fast growing research field. There is still yet further research to come to help to find possible "cures" for affected patients. With the recent prognosis by Ferreira et al. (2023) [1] warning of an increase of one-third more LBP patients in the next 50 years, the scientific community is urged to find better treatment options and also especially early diagnostic tools to foresee critical cases to come.

**Acknowledgments:** B. Gantenbein wishes to acknowledge his funding sources, namely, the Marie Skłodowska Curie International Training Network (ITN) "disc4all" (https://disc4all.upf.edu, accessed on 6 June 2023) grant agreement #955735 (https://cordis.europa.eu/project/id/955735, accessed on 6 June 2023), the iPspine project, and the Swiss National Science Foundation projects #310030E\_192674/1 (https://data.snf.ch/grants/grant/192674, accessed on 6 June 2023) and project #40B2-0\_211510/1 (https://data.snf.ch/grants/grant/211510, accessed on 6 June 2023).

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

**Veronica Tilotta, Gianluca Vadalà \*, Luca Ambrosio, Fabrizio Russo, Claudia Cicione, Giuseppina Di Giacomo, Rocco Papalia and Vincenzo Denaro**

> Laboratory for Regenerative Orthopaedics, Department of Orthopaedic and Trauma Surgery, Campus BioMedico University of Rome, Via Alvaro del Portillo 200, 00128 Rome, Italy; v.tilotta@unicampus.it (V.T.); l.ambrosio@unicampus.it (L.A.); fabrizio.russo@unicampus.it (F.R.); c.cicione@unicampus.it (C.C.); g.digiacomo@unicampus.it (G.D.G.); r.papalia@unicampus.it (R.P.); denaro@unicampus.it (V.D.) **\*** Correspondence: g.vadala@unicampus.it

> **Abstract:** Low back pain (LBP) is one of the most frequent symptoms associated with intervertebral disc degeneration (IDD) and affects more than 80% of the population, with strong psychosocial and economic impacts. The main cause of IDD is a reduction in the proteoglycan content within the nucleus pulposus (NP), eventually leading to the loss of disc hydration, microarchitecture, biochemical and mechanical properties. The use of mesenchymal stem cells (MSCs) has recently arisen as a promising therapy for IDD. According to numerous reports, MSCs mediate their regenerative and immunomodulatory effects mainly through paracrine mechanisms. Recent studies have suggested that extracellular vesicles (EVs) extracted from MSCs may be a promising alternative to cell therapy in regenerative medicine. EVs, including exosomes and microvesicles, are secreted by almost all cell types and have a fundamental role in intercellular communication. Early results have demonstrated the therapeutic potential of MSCs-derived EVs for the treatment of IDD through the promotion of tissue regeneration, cell proliferation, reduction in apoptosis and modulation of the inflammatory response. The aim of this review is to focus on the biological properties, function, and regulatory properties of different signaling pathways of MSCs-derived exosomes, highlighting their potential applicability as an alternative cell-free therapy for IDD.

> **Keywords:** extracellular vesicles; exosomes; mesenchymal stem cells; miRNA; low back pain; intervertebral disc degeneration; regenerative medicine

#### **1. Introduction**

Low back pain (LBP) is one of the most common musculoskeletal symptoms affecting up to 80% of adults aged between 40–80 years. Chronic LBP is one of the main causes of disability and loss of work ability with a strong impact on patients' quality of life; additionally, it poses a significant socioeconomic burden on public health worldwide [1]. According to a recent global report on the prevalence of LBP, it has been estimated that two thirds of adults are afflicted with LBP at least once in their lives. In addition, the financial burden of LBP in the United States is estimated to exceed \$100 billion per year, depending on the inevitable effect of changes in work status, prolonged hospitalization, and increased outpatient visits [2–4].

The majority of LBP episodes are associated with intervertebral disc degeneration (IDD). The intervertebral disc (IVD) acts as a shock absorber between the vertebrae while allowing for flexion–extension, lateral flexion, and rotation movements. The IVD is a complex structure connecting two adjacent vertebrae and is composed of the annulus fibrosus (AF), nucleus pulposus (NP) and cartilaginous endplates (CEP) [5]. The AF is a fibro-cartilaginous ring encircling the NP and it consists of concentric lamellae of type I collagen fibers associated with an interlamellar matrix, in which scattered cells with fibroblast-like morphology and phenotype are present [6]. The NP presents a gelatinous

**Citation:** Tilotta, V.; Vadalà, G.; Ambrosio, L.; Russo, F.; Cicione, C.; Di Giacomo, G.; Papalia, R.; Denaro, V. Mesenchymal Stem Cell-Derived Exosomes: The New Frontier for the Treatment of Intervertebral Disc Degeneration. *Appl. Sci.* **2021**, *11*, 11222. https://doi.org/10.3390/ app112311222

Academic Editor: Benjamin Gantenbein

Received: 18 October 2021 Accepted: 25 November 2021 Published: 26 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

5

extracellular matrix (ECM) rich in proteoglycans (mainly aggrecan), type II collagen and other proteins such as elastin and fibronectin. The high level of water content helps maintain disc height and supports the functions of the IVD [7]. The NP cell population has a notochordal origin: notochordal cells (NCs) act as progenitor cells and play an important role in stimulating ECM synthesis by NP cells (NPCs), although in adult human NP tissue the number of NCs diminishes drastically [8]. To preserve IVD architectural features, it is essential to maintain disc cell metabolism and physiological biomechanics, but the avascular and relatively acellular environment of the NP leads to a limited regenerative capacity against stressful and degenerative stimuli. IDD is an age-related disease and is mainly linked to loss of NP hydration and reduction in the proteoglycan content [9]. Current approaches to treat IDD are based on conservative or surgical procedures that aim to relieve the symptoms; however, they are not able to change the natural history of the disease [10].

Recent advances in the understanding of IVD biology have led to an increased interest in the development of novel biological treatments. Growth factors, gene therapy, stem cell transplantation and tissue engineering might support intervertebral disc regeneration [11–16]. The transplantation of mesenchymal stem cells (MSCs), including bone marrow-derived (BM-MSCs) and adipose tissue (Ad-MSCs), has been evaluated in different clinical trials. Although encouraging results have been achieved, there are still several impediments related to the survival and differentiation of transplanted MSCs due to the particularly hostile microenvironment of degenerated IVDs [17]. Moreover, costs related to cell expansion under Good Manufacturing Practice (GMP) conditions and to the quality control guidelines for Advanced Therapy Medical Products (ATMP) limits the large-scale application of a cell therapy strategy.

In the last decade, the interest in extracellular vesicles (EVs) and their therapeutic potential has grown regarding different applications in several fields, including oncology, neurology, cardiology, and orthopedics [18,19]. EVs are endogenous nanoparticles released by different cell types and are involved in various physiological mechanisms including apoptosis, tissue regeneration, inflammation, ECM synthesis, cell proliferation and immunomodulation [20]. Based on morphology, biogenesis, size, and content, EVs are generally distinguished in apoptotic bodies, microvesicles and exosomes [21]. One of the contemporary challenges in regenerative medicine is the development of new biology-based therapeutic strategies employing EVs. Indeed, MSCs-derived exosomes have recently demonstrated their potential benefits in preclinical models of cardiovascular diseases, graft-versus-host disease and osteoarthritis [19], while their role in IDD still remains unclear [22].

The aim of this review is to focus on the biological properties, function, and regulatory properties of different signaling pathways of MSCs-derived exosomes, highlighting their potential applicability as an alternative cell-free therapy for IDD.

#### **2. Intervertebral Disc Degeneration**

IDD is an age-related chronic and progressive process, characterized by a gradual loss of proteoglycans and water content in the NP with a reduction in the disc's capacity to withstand compressive loads [17].

The IVD is an avascular organ, and the supply of nutrition occurs by diffusion from its surrounding vasculature through the CEPs [23]. However, IVD aging and mechanical overload progressively foster endplate calcification and capillary obliteration, thus reducing nutrient delivery to the IVD itself [24]. Due to the lack of blood flow, a hypoxic microenvironment is consequently established with higher O2 concentrations at the surface of the AF (19.5%) and a reduction towards the center of the NP (0.65%). However, low oxygen concentration has not been shown to be detrimental for NPC viability [25]. Another factor involved in making the IVD microenvironment hostile is the low glucose concentration, which ranges from about 5 mM at the periphery to about 0.8 mM at the center of the IVD [26].

pH is another crucial parameter in the IVD microenvironment. Local production of lactic acid by IVD maintains a pH between 2 and 6, with significant effects on cell survival and matrix synthesis, as well as on the expression of proinflammatory cytokines and pain-related factors [27]. The main change in NP during IDD is in the imbalance between anabolic factors, such as insulin-like growth factor 1 (IGF-1), growth differentiation factor 6 (GDF-6), transforming growth factor β (TGF-β), bone morphogenetic proteins (BMPs) and catabolic enzymes, such as matrix metalloproteinases (MMPs) and a disintegrin and metalloproteinase with thrombospondin motifs (ADAMTS) [28].

Aggrecan, the most abundant proteoglycan in the NP matrix, binds to sulphated glycosaminoglycans (GAGs, e.g., keratan and chondroitin sulfate), and consequently retains hydration and osmotic pressure within the IVD. Negatively charged GAGs regulate the ion balance of the ECM [29]. In a healthy state, osmolarity can range from 430 mOsm/L (isosmotic pressure) to about 496 mOsm/L (hyperosmotic pressure) [30]. In degenerative conditions, due to the progressive loss of proteoglycans, the osmolarity of the IVD decreases until it reaches values of about 300 mOsm/L [31]. Tao et al. reported that high osmolarity can not only decrease the viability and proliferation of disc cells but can also influence the expression levels of collagen type II, SOX-9 and aggrecan in progenitor NPCs [32].

Excessive mechanical loading may be another crucial factor that compromises the biomechanical characteristics of IVD. Static compression stress induces apoptosis, senescence, mitochondrial damage, and matrix catabolism due to a reduction in collagen and aggrecan expression with an increase in MMPs [33]. On the contrary, dynamic compressive loading often triggers an anabolic response and improves the transport of nutrients and growth factors within the NP [34].

The inflammatory response may also contribute to the onset and progression of the degenerative process [35]. The nutritional decrease and accumulation of degraded matrix products promote macrophage-mediated production of pro-inflammatory cytokines such as interleukin (IL)-1β and tumor necrosis factor-α (TNF-α) through activation of the nuclear factor kappa-light-chain-enhancer of activated B cells (NF-kB) pathway [36]. IL-1β and TNF-α are involved in the progression and initiation of IDD by regulating inflammatory response, IVD cell proliferation, senescence, apoptosis, pyroptosis, autophagy, ECM destruction, and oxidative stress [37]. TNF-α is able to amplify the inflammatory response by stimulating IVD cells to synthesize additional cytokines, including IL-6, IL-8, IL-17, IL-1β and numerous chemokines [38]. IL-1β is expressed at higher levels in degenerated IVDs with an increased production of genes and proteins for the IL-1 receptor. Therefore, although both cytokines are involved in the pathogenesis of IDD, IL-1β may play a more significant role, and could thus be considered an ideal therapeutic target [39]. Tissue aging increases oxidative stress, mitochondrial dysfunction, DNA damage, cell senescence and apoptosis, thus compromising ECM homeostasis [40]. Besides, the viability of NPCs has an important impact on ECM composition. The declining number of cells, their phenotype modification and the inability to withhold compression stresses lead to the deterioration of the AF lamellar architecture and the reduction in ECM synthesis [41,42]. The low metabolic activity of cells and the avascular nature of IVD may explain the inability to self-repair following injury and degeneration [43]. However, recent research has reported the presence of progenitor stem cells in the NP, which may be physiologically responsible, to a certain extent, for endogenous repair of intrinsic insults to IVD integrity.

The natural history of IDD may progress to disc herniation, degenerative spondylolisthesis, spinal instability, and stenosis associated with neurological symptoms such as radiculopathy and/or myelopathy, and surgical intervention is often required [44]. However, the available treatment approaches are not able to arrest IDD; they are only able to relieve accompanying symptoms. Therefore, strenuous research efforts are being made to develop effective strategies that can achieve IVD regeneration. Therefore, the potential use of EVs appears to be a promising therapeutic tool that may support endogenous repair and re-establish the proteoglycan content, thus improving disc hydration and biomechanics.

#### **3. MSC Interaction with NP Cells**

Human MSCs are multipotent cells present in almost all tissues such as bone marrow, adipose tissue, skin, lung, synovial membrane, dental pulp, nasal olfactory mucosa, muscle, periosteum, corneal limbus, peripheral blood, endometrial and menstrual blood, uterine cervix, and in fetal/neonatal tissues [9]. Stem cell transplantation represents a promising approach to restore the homeostasis of the degenerated environment of IVD due to the ability of the stem cells to differentiate towards NP-like cells and their immunomodulatory and anti-catabolic effects [14]. In recent decades, several preclinical studies have evaluated the interaction between MSCs and NPCs in vitro [44–46].

Cunha et al. tested the systemic delivery of BM-MSCs in a rat IVD and reported a diminished hypoxic response, risk of herniation and synthesis of cytokines [45]. In an experimental study, Wangler et al. described the homing of MSCs as a potential strategy to prevent the onset of the degenerative cascade in IVDs. They have also evaluated the effect of human MSCs homing in bovine and human disc cells expressing Tie2 (also known as angiopoietin-1 receptor TEK tyrosine kinase). Tie2+ cells represent a progenitor cell population with discogenic differentiation potential. The results showed a proliferative response and reduction in dead cells in both tissues [46].

Our group was one of the first to report the crosstalk mechanism between MSCs and NPCs. In a 3D culture system, where a short distance may improve communication among cells, a paracrine effect of MSCs on NPCs collected from a degenerative disc was demonstrated. Gene expression analysis on both cell types showed that MSCs exerted a trophic effect on NPCs while acquiring a chondrogenic phenotype due to the communication with NPCs [47]. In another similar in vitro study, we highlighted the significant increase in GAG content and proliferation of NPCs when co-cultured with MSCs [48].

In the last decade, intradiscal stem cell therapy for IDD has provoked increasing interest in the field. Since cell isolation and implantation procedures can be safely performed in a sterile environment and with low immunogenicity, both autologous and allogeneic transplants can be used [49]. In a pilot study by Orozco et al., 10 patients with chronic LBP were injected with autologous BM-MSCs for the treatment of lumbar disc degeneration. Three months after surgery, LBP and disability were significantly reduced. Despite the small number of patients, this study demonstrated that intradiscal BM-MSC transplantation was a safe procedure with a promising role in the treatment of IDD [50]. Similarly, a randomized controlled trial including 24 patients with lumbar IDD reported a significant reduction in pain in the lumbar region and reduced disability one year after BM-MSC allogeneic transplant, whereas the control group, who received a sham infiltration within the paravertebral musculature did not show any improvement [51].

Earlier research has proved that MSCs secrete bioactive molecules and trophic factors such as soluble proteins, nucleic acids, lipids and EVs, and highlights their anti-apoptotic, immunomodulatory, angiogenetic, anti-fibrotic functions and differentiation capacity, probably mediated by paracrine mechanisms [52–54]. The products released by MSCs in vitro and in vivo constitute their secretome [55]. In a bovine ex vivo model, Teixeira et al. investigated the immunomodulatory paracrine effect of BM-MSCs in the proinflammatory IDD microenvironment. While MSCs viability was not influenced, cell migration increased, collagen type II and aggrecan levels remained stable while inflammatory cytokines IL-6, IL-8 and TNF-α were downregulated by MSCs in IL-1β-stimulated IVDs [56]. Furthermore, MSCs have been shown to migrate from their niches to damaged sites in response to injury, immune or physicochemical signals. This physiological mechanism of MSCs recruitment or homing may play a role in promoting the repair of different musculoskeletal tissues [57].

Interestingly, recent studies investigated the cross-talk between MSCs and NPCs through trophic factors [58]. Shim et al. investigated paracrine interactions in a co-culture system by isolating MSCs, AF and NPCs from the same donor. The characterization revealed an increase in mRNA levels expression of ECM genes (SOX9, COL2A1 and ACAN), the downregulation of MMPs and ADAMTS, as well as the inhibition of the production of pro-inflammatory factors compared to monocultures [59]. Li et al. also employed a co-culture system to understand the influence of BM-MSCs on TNF-α-induced NPC degeneration and reported an increase in cell proliferation and type II collagen production, a decrease in MMP-9 expression, and reduced TGFβ/NF-κB signaling in pro-senescent NPCs. Premature senescence and aging of NPCs were restored in the coculture system [60]. Even though the regenerative, anabolic, and inflammatory-modulating properties of MSCs have been demonstrated in several investigations, many aspects are still unclear.

The harsh microenvironment with the absence of vasculature, hypoxia, low glucose concentration, acidity, hyperosmolarity and continuous mechanical stimulation plays a crucial role for cell survival and their regenerative potential [61]. Furthermore, MSC implantation may further impair the balance between the nutrient supply and demand of resident disc cells, thus causing cell death. Miguélez-Rivera et al. investigated the effectiveness of a conditioned medium (CM) derived from the culture of MSCs in a rat IDD model. MSC-derived CM has shown immunomodulatory and anti-inflammatory properties in a detrimental environment by down-regulating the expression levels of pro-inflammatory cytokines such as IL-1β, IL-6, IL-17 and TNF [62].

Although MSC-based therapy may be considered an attractive and a safe option, there are still unsolved issues. In addition to the abovementioned shortcomings, MSCs require a long and expensive process for in vitro expansion according to the minimum standard to guarantee controlled conditions in production processes and clinical procedures [63].

#### **4. Extracellular Vesicles**

EVs were discovered in 1987 by Johnstone et al. during the in vitro maturation of reticulocytes in culture [64]. Secreted EVs are a group of cell-derived membranous structures with biological potential in numerous human diseases (e.g., cancer, neurodegenerative and musculoskeletal disorders). Collection, pre-processing, separation, concentration, and characterization of EVs have been widely studied. In 2018, the International Society for Extracellular Vesicles (ISEV) proposed the Minimal Information for Studies of Extracellular Vesicles (MISEV) as an update of general knowledges presented in MISEV 2014 guidelines to allow the reproducibility of results among investigators [65].

#### *4.1. Size and Morphology*

The determining factors that discriminate the different classes of EVs are size, morphology, origin and content (Table 1) [66].


**Table 1.** General features of exosomes, microvesicles and apoptotic bodies.

ILVs = intraluminal vesicles; MVEs = multivesicular endosomes.

It is possible to distinguish three main EV categories: apoptotic bodies, microvesicles and exosomes (Figure 1) [67].

**Figure 1.** Extracellular vesicles distinguished by size: apoptotic bodies (800–5000 nm), microvesicles (50–1000 nm) and exosomes (30–100 nm). Created with BioRender.com (accessed on 10 October 2021).

Apoptotic bodies have a size range of 800–5000 nm with a heterogeneous morphology. They also include smaller vesicles known as apoptotic cell-derived microparticles or apoptotic blebs released by apoptotic cells [66].

Microvesicles were originally described as "platelet dust", as they were first isolated from platelets [68]. They have a diameter of 50–1000 nm, a spherical structure, and are expelled through exocytosis of the plasma membrane [69].

The last class of vesicles are smaller than those mentioned above, have a size range of 30–100 nm, a cup-shaped morphology and are called exosomes [70]. However, some studies have shown that the cup-shaped morphology could be an experimental artefact and the overlapping range of size may lead to unclear terminology [71]. Further isolation and characterization protocols are expected to establish an appropriate classification and nomenclature [72].

#### *4.2. Content and Biogenesis*

Cell type and the physiological or pathological state of the parent cell may influence its EV content [73]. To date, the interest in EVs has grown exponentially due to their abundant presence in body fluids, their wide range of regulatory functions and their capacity as vehicles for DNA, mRNAs, miRNAs, proteins, and lipids involved in intercellular communication (Figure 2) [74].

**Figure 2.** Exosomal composition in terms of bioactive molecules: lipids, proteins, and nucleic acids. Created with BioRender.com (accessed on 10 October 2021).

The protein content of EVs is often used as a marker to recognize EV sub-populations. Exosomes are rich in endosome-associated proteins (e.g., Rab GTPases, SNAREs, annexins, and flotillin), proteins involved in EV biogenesis (e.g., Alix and Tsg101), transmembrane proteins composed of four domains such as tetraspanins (CD63, CD81, CD9 and CD82), major histocompatibility complex (MHC) molecules and cytosolic proteins such as heat shock proteins [20]. EVs also present a lipidic component. Exosomes are enriched in cholesterol, sphingomyelin, phosphatidylcholine, and phosphatidylethanolamine, as well as proteins associated to lipid rafts such as glycosylphosphatidylinositol-anchored proteins and flotillin. Nevertheless, the lipid content of microvesicles is less known than exosomes and is primarily involved in their budding from the plasma membrane [75]. Regarding nucleic acids, it was reported that microvesicles and exosomes carry single-stranded DNA (ssDNA), double-stranded DNA (dsDNA), genomic DNA, mitochondrial DNA, and reverse-transcribed complementary DNAs [76,77]. On the other hand, EVs can deliver mRNA, mRNA fragments, long non-coding RNA (lnRNA), miRNA, ribosomal RNA (rRNA) and fragments of tRNA, vault- and Y-RNA [78]. More importantly, miRNAs may act on the behavior of recipient cells by regulating gene expression, immune response, and additional physiological processes [79,80].

Exosomes and microvesicles have two different pathways of biogenesis. They are constantly produced by viable cells, while apoptotic bodies are released exclusively during the apoptotic phase of the cell cycle [81]. Whereas the releasing mechanisms of apoptotic bodies from the cell membrane have been widely elucidated, those related to microvesicles have been investigated only recently. Microvesicles are formed by fusion with the internal side of the plasma membrane and requires several rearrangements within the membrane itself such as reorganization of lipid and protein components and Ca2+ levels [82,83]. It has been reported that a reduction in cholesterol in microvesicles may compromise their generation in activated neutrophils [84]. In addition, the small GTPase family Rho and the protein kinase associated with Rho (ROCK) induce the biogenesis of microvesicles in various populations of cancer cells by regulating the activity of cytoskeletal components such as actin [85].

Exosomes are defined as intraluminal vesicles (ILVs) formed by the invagination of the endosomal membrane during the maturation of multivesicular endosomes (MVEs). This process leads to the acquisition of cargo and emission of vesicles in the lumen through particular sorting machineries (Figure 3) [86].

**Figure 3.** Visualization of MVEs by electron microscopy (scale bar = 1 μm). Biogenesis process of ILVs by budding into early endosomes. MVEs release exosomes by fusion with the plasma membrane or fuse with lysosomes. ILVs = intraluminal vesicles; MVEs = multivesicular endosomes. Created with BioRender.com (accessed on 10 October 2021).

The discovery of the endosomal sorting complex required for transport (ESCRT) machinery was important as it seems to be a fundamental component in several cell types for the generation of MVEs and ILVs [87,88]. The first step is the enrolling of ESCRT-0 subunits, which select cargo in a ubiquitin-dependent manner; ESCRT-I and ESCRT–II are involved in budding endosomal membrane formation, and ESCRT-III is deputed to ILV scission and recycling of the ESCRT machinery [89]. Tumor susceptibility gene 101 (TSG101) and ALIX are two important members of the ESCRT complex discovered by Théry et al. in a proteomic analysis of mouse dendritic cell-derived exosomes [90].

Exosomes can also be formed through an ESCRT-independent mechanism. The sphingolipid ceramide plays a role by creating a spontaneous negative curvature in the membrane leaflets of MVEs to shape ILVs [91]. The formation and the sorting of ILV cargo can also be regulated by the syndecan–syntenin–ALIX pathway. Indeed, syndecan heparan sulphate proteoglycans are connected to their cytoplasmic adaptor syntenin, which interacts with ALIX hence promoting endosomal membrane budding [92].

Finally, exosomal surface proteins belonging to the tetraspanine family, such as CD63 CD81, CD82 and CD9 or heat shock protein (Hsp70), may control ESCRT-independent exosomes biogenesis [93,94]. Intracellular vesicular trafficking, fusion with the plasma membrane, and consequently, exosome secretion are driven by Rab GTPases proteins [95]. Mechanisms regulating the sorting of nucleic acids in exosomes are still not known.

In summary, microvesicles bud directly from the plasma membrane, while exosomes are firstly incorporated into the endosomal vesicles and then released into the extracellular milieu by fusion of MVEs with the plasma membrane. Once released, EVs are absorbed by neighboring cells by endocytosis, phagocytosis, pinocytosis, or fusion of membranes. Ligands on the extracellular vesicle membrane can also bind to a receptor on target cells to thereby induce an intracellular signal [96].

#### **5. Exosome Therapy for IVD Regeneration**

Emerging cell therapies have demonstrated promising results by supplementing damaged IVDs with stem cells harvested from different sources [97]. Despite this, there is a growing interest in cell-free therapy. Recent research has proposed that exosomes derived from different tissues such as bone marrow and adipose tissue could represent an emerging tool in regenerative medicine [98]. The use of MSC-derived exosomes may offer several advantages. In fact, they do not require tissue engraftment but present high biocompatibility, low immunogenicity and toxicity, less oncogenic potential and substantial economic benefits compared to cell therapy. Such aspects may promote the use of MSCs-derived exosomes for IVD regeneration as an attractive alternative to intradiscal cell therapy [99].

However, before the exploitation and application of exosomes, it is necessary to isolate them efficiently from biological fluids or supernatants from in vitro cell cultures [100,101]. Ultracentrifugation is the most used technique and is often coupled with ultrafiltration to improve the purity of the sample, although this provides the risk of nanoparticle deformation [102,103]. Other sophisticated approaches include tangential flow filtration, dimensional exclusion chromatography and polymeric precipitation [104–106]. Several methods are available for exosome characterization. Transmission electron microscopy (TEM) and scanning electron microscopy (SEM) are often employed to verify their morphology. Nanoparticle tracking analyzers (NTA) measure the diameter size and specific EV markers (CD9, CD63, CD81, TSG101 and ALIX) can be identified through flow cytometry [107].

The role of exosomes has been evaluated with regards to ECM synthesis, metabolic dysfunction, anti-apoptotic processes, inflammation, and wound healing (Table 2) [108,109].


**Table 2.** Summary of available evidence on the role of EVs in IDD.

AF = annulus fibrosus; BM-MSCs = bone marrow-derived mesenchymal stem cells; CEP = cartilaginous endplate; ECM = extracellular matrix; ER = endoplasmic reticulum; EVs = extracellular vesicles; GAG = glycosaminoglycan; hMSCs = human mesenchymal stem cells; HUVECs = human umbilical vein endothelial cells; IDD = intervertebral disc degeneration; IL = interleukin; IVD = intervertebral disc; MMP = matrix metalloproteinase; MSC = mesenchymal stem cell; NGF = nerve growth factor; NP = nucleus pulposus; NPC = nucleopulpocyte; ROS = reactive oxygen species.

The beneficial effect of MSC-derived exosomes on NPCs in acidic conditions has been recently demonstrated. BM-MSC-exosomes increased proliferation and attenuated the apoptosis of NPCs and enhanced the expression of ECM components [115]. Moreover, exosomes are able to inhibit IL-1β-induced inflammation and apoptosis in AF cells and they promote cell proliferation [116]. On the other hand, Liu et al. identified various similarities between CEP stem cells and MSCs, including their capacity to communicate with IVD cells by secreting extracellular vesicles [121]. Indeed, IVD repair may be promoted by CEP stem cell exosomes by inhibiting NPC apoptosis and encouraging the migration of CEP stem cells into the IVD and differentiate into NPCs via autocrine mechanisms [117,118]. AF-derived exosomes have also showed a protective effect in IDD [120]. Aging and degenerative tissues usually accumulate in endoplasmic reticulum (ER) advanced glycation end-products (AGEs) related to inflammation, cell death, metabolic dysfunction, and stress [122,123]. Liao et al. found that the injection of MSCs-derived exosomes in an in vivo rat IDD model modulated ER stress by protecting NPCs against cell death and ameliorating IDD through an attenuation of AGE-induced ER stress mediated by AKT and ERK pathways [110]. Exosomes may also induce differentiation of stem cells towards different phenotypes, although the underlying mechanisms are still unknown. It was reported that the inhibition of Notch1 may promote exosome-induced differentiation of MSCs into NP-like cells in vitro [111]. Lu et al. analyzed exosomes secreted by both NP cells and BM-MSCs in a co-culture system. They highlighted that NPC-derived exosomes induced BM-MSC differentiation towards a NP-like phenotype, whereas BM-MSC-derived exosomes promoted NPC proliferation and ECM production [112].

Recently, a study by Xia et al. described the protective role of exosomes in a rabbit IDD model. Exosomes inhibited H2O2-induced NLRP3 inflammasome activation, drastically reduced the catabolic action of MMP-3, MMP-13, IL-6 and inducible nitric oxide synthase (iNOS), repressed reactive oxygen species (ROS) generation and mitochondrial dysfunction was restored [113]. Hingert et al. have reported the potential abilities of MSC-derived EVs in a 3D in vitro model to improve cell viability, proliferation, ECM production, apoptosis, chondrogenesis and cytokine secretion in NPCs harvested from patients with IDD and chronic LBP [114]. Furthermore, an actual challenge is the potential of engineered EVs to lead NPC towards a healthy NP-like phenotype, decrease catabolism and inflammation, and improve GAG accumulation both in vitro and in vivo.

However, the clinical application of exosomes is still premature, and appropriate strategies will be necessary to exploit exosomes as drug delivery vehicles with high specificity, non-cytotoxic effects and low immunogenicity [67].

#### **6. The Role of miRNAs in IDD**

Exosomes may be promising diagnostic biomarkers due to their ability to incorporate several molecules including lipids, RNAs and proteins targeting specific recipient cells or tissues [124]. Interestingly, many studies have shown that several miRNAs can be predictive markers for IDDs. miRNAs are single-stranded and small non-coding RNAs that can bind to the 3- -untranslated region (3- -UTR) of target mRNA molecules, leading to their translational repression or degradation [125]. They represent key biomolecular players that regulate the expression of several target genes and are involved through various pathways in cell proliferation, apoptosis, immunomodulation, ECM anabolism and catabolism (Figure 4) [126].

**Figure 4.** Potential therapeutic role of exosomes in IDD through the modulation of ER stress, stimulation of production of ECM components, anti-apoptotic and anti-inflammatory effects on NP cells. BM-MSCs = bone marrow-derived stem cells; ECM = extracellular matrix; ER = endoplasmic reticulum; NP = nucleus pulposus; NPC = nucleus pulposus cell. Created with BioRender.com (accessed on 10 October 2021).

The anti-apoptotic effects of exosomes from MSCs-derived CM were described in a study conducted by Cheng et al. These results seem to be partially mediated by the transfer of exosomal miR-21 associated with the downregulation of homologous phosphatase and tensin (PTEN), which leads to the activation of the phosphoinositide 3-kinase (PI3K)/Akt pathway and the inhibition of TNF-α-induced NPC apoptosis [127]. The regulatory role of miR-210 was also investigated. Indeed, miR-210 upregulation inhibited apoptosis in human NPCs, suggesting its involvement in the etiology of IDD [128]. It has hypothesized that increased miR-223 expression may be a predictive marker of chronic lumbar radicular pain, as it would be released from the NP following disc herniation [129]. Similarly, the transfection of miR-146a in bovine NPCs decreased IL-1β-driven inflammatory and catabolic responses both in vitro and ex vivo [130]. Chen et al. have reported that miR-494-5p was upregulated in NP tissues of a mouse IDD model, while tissue inhibitors of metalloproteinases-3 (TIMP-3) were downregulated. The reduced expression of this miRNA promoted viability and attenuated the apoptosis and senescence of NPCs [131]. The upregulation of miR-142-3p or miR-199 delivered by MSC-derived exosomes may protect NPC from injury through inhibition of mitogen-activated protein kinase (MAPK) signaling [132,133]. In addition, miR-199a decreased MMP-2, MMP-6, TIMP-1 and apoptotic cells in a mouse model of IDD, while the expression of collagen type II, aggrecan, increased following the exposition to BMSC EVs in vitro [134]. Among other miRNAs previously identified, miR-27a was involved in the regulation of apoptotic pathways and in the prevention of IDD through the ability to suppress ECM degradation targeting MMP-13 [135,136]. Similarly, under TNF-α stimulation, exosomes extracted from BM-MSCs were able to secrete miR-532-5p. This exosomal miRNA has been linked to a reduction in apoptosis in NP cells and matrix degradation by targeting Ras association domain family 5 (RASSF5), attenuating disc degeneration [137]. Zhang et al. suggested the potential effect of MSC-derived exosomes on pyroptosis activation and the role of miR-410 as a specific mediator of NPC pyroptosis [138].

Activating transcription factor 6 (ATF6) is a transcription factor identified as the target gene of miR-31-5p. When miR-31-5p was upregulated, the expression of ATF6 and other proteins involved in oxidative stress-induced ER-stress was suppressed, resulting in a protective effect in IDD [139].

Another family of molecules delivered by exosomes is represented by circular noncoding RNAs (circRNAs). CircRNAs may act as post-transcriptional regulators and interact with miRNAs as miRNA sponges. In a recent study, the role of exosomal circRNA\_0000253 was explored, suggesting that it may competitively absorb miRNA-141-5p and downregulate SIRT1 promoting IDD progression [140].

Further research is required to optimize the effects of MSCs-derived exosomes and to confirm the association of their miRNAs with the pathogenesis of IDD. It will be crucial to know the active components responsible for their actions, underlying mechanisms, pharmacokinetic aspects, and safety of this approach. Table 3 summarizes the different miRNA studied in IVD research.


**Table 3.** Summary of available evidence on the role of miRNA in IDD.

ATF = activating transcription factor; BM-MSCs = bone marrow-derived mesenchymal stem cells; CEP = cartilaginous endplate; ECM = extracellular matrix; EVs = extracellular vesicles; GREM = gremlin; HOX = homeobox; IL = interleukin; IRAK = interleukin 1 receptor associated kinase; IDD = intervertebral disc degeneration; IVD = intervertebral disc; MAPK = mitogen-activated protein kinase; MLK = mixed-lineage protein kinase; MMP = matrix metalloproteinase; NLRP = nucleotide-binding oligomerization domain, leucine rich repeat and pyrin domain containing; NPC = nucleus pulposus cell; PI3K = phosphoinositide 3-kinase; PTEN = phosphatase and tensin homolog; RASSF = Ras association domain family; TGF = transforming growth factor; TIMP = tissue inhibitor of metalloproteinases; TNF = tumor necrosis factor; TRAF = tumor necrosis factor receptor-associated factor.

#### **7. Conclusions**

In recent years, the transplantation of adult MSCs has been considered a promising approach for IVD tissue regeneration due to their capacity to secrete growth factors, cytokines, and other molecules with paracrine activity. However, the success of MSC-based therapies depends on the survival rate and functionality of transplanted cells within the harsh IVD microenvironment. A growing number of clinical trials have opted for cell therapy, even though there are still many doubts about their routine clinical application.

The use of cell-free therapies could circumvent many obstacles. Numerous studies have recently shown that EVs are released by many cell types including MSCs. These small vesicles can transfer functional proteins, RNA, miRNAs, and lipids among cells, hence promoting intercellular communication, proliferation, inhibition of apoptosis and tissue damage repair. Exosomes represent a promising alternative to traditional stem cell-based therapies. Interestingly, exosomal miRNAs may be predictive biomarkers that regulate the balance between regeneration and degeneration mechanisms in IDD. However, due to many unknown pathophysiological processes in which exosomes are involved, further studies are necessary to confirm the ability of these nanoparticles to promote IVD health and homeostasis.

**Author Contributions:** Conceptualization, V.T., G.V.; methodology, V.T.; data curation, V.T., C.C., G.D.G.; writing—original draft preparation, V.T.; writing—review and editing, G.V., L.A., F.R., C.C., G.D.G.; supervision, R.P., V.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

**Acknowledgments:** The support by the ACTIVE (ClinicalTrials.gov Identifier: NCT04759105), DREAM (ClinicalTrials.gov Identifier: NCT05066334), RESPINE (ClinicalTrials.gov Identifier: NCT03737461) and iPSspine (European Union's Horizon 2020 Research and Innovation Programme under grant agreement no. 825925) projects are gratefully acknowledged.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Effects of Growth Factor Combinations TGF**β**3, GDF5 and GDF6 on the Matrix Synthesis of Nucleus Pulposus and Nasoseptal Chondrocyte Self-Assembled Microtissues**

**Shani Samuel 1,2, Emily E. McDonnell 1,2 and Conor T. Buckley 1,2,3,4,\***


**Abstract:** There has been significant interest in identifying alternative cell sources and growth factor stimulation to improve matrix synthesis for disc repair. Recent work has identified nasoseptal chondrocytes (NC) as a possible alternative cell source with significant matrix-forming abilities. While various growth factors such as members of the TGFβ superfamily have been explored to enhance matrix formation, no consensus exists as to the optimum growth factor needed to induce cells towards a discogenic phenotype. This study assessed both nucleus pulposus (NP) and NC microtissues of different densities (1000, 2500 or 5000 cells/microtissue) stimulated by individual or combinations of the growth factors TGFβ3, GDF5, and GDF6. Lower cell densities result in increased sGAG/DNA and collagen/DNA levels due to higher nutrient availability levels. Our findings suggest that growth factors exert differential effects on matrix synthesis depending on the cell type. NP cells were found to be relatively insensitive to the different growth factor types examined in isolation or in combination. Overall, NCs exhibited a higher propensity to form extracellular matrix compared to NP cells. In addition, stimulating NC-microtissues with GDF5 or TGFβ3 alone induced enhanced matrix formation and may be an appropriate growth factor to stimulate this cell type for disc regeneration.

**Keywords:** microwell; GFD5; GDF6; TGFβ3; microtissues; intervertebral disc

#### **1. Introduction**

Cell-based therapies may hold significant potential as a treatment strategy for the repair and regeneration of the intervertebral disc (IVD). However, it remains challenging to identify an appropriate cell source that can be obtained with minimal donor site morbidity and produce a matrix with high levels of proteoglycan containing predominately collagen type II and low levels of collagen type I. The most commonly explored cell types for IVD regeneration include mesenchymal stem cells (MSCs) [1,2], articular chondrocytes [3] and disc-derived cells [4]. Studies have demonstrated that reimplantation of extracted NP cells can retard degenerative changes, and the injection of autologous nucleus pulposus (NP) cells has been clinically tested in humans with positive outcomes [5]. DiscGenics Inc. (Salt Lake City, Utah) is actively investigating their propriety technology IDCT (injectable discogenic cell therapy) which contains a mixture of progenitor cells derived from allogeneic "discogenic" cells and a viscous biomaterial [6]. However, the use of NP cells is limited due to the matrix-forming capacity of expanded NP cells derived from degenerated tissue [7] and the number of nucleus pulposus (NP) cells that can be isolated from a degenerated

**Citation:** Samuel, S.; McDonnell, E.E.; Buckley, C.T. Effects of Growth Factor Combinations TGFβ3, GDF5 and GDF6 on the Matrix Synthesis of Nucleus Pulposus and Nasoseptal Chondrocyte Self-Assembled Microtissues. *Appl. Sci.* **2022**, *12*, 1453. https://doi.org/10.3390/ app12031453

Academic Editor: Benjamin Gantenbein

Received: 13 December 2021 Accepted: 27 January 2022 Published: 29 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

disc is insufficient [8] to meet the requirements for successful treatment without significant culture expansion.

In contrast, stem cells can be obtained from various sources (bone marrow or adipose tissue) and have been shown to survive and proliferate after implantation into the disc [9]. However, cell leakage at the injection site has been shown to result in osteophyte formation [10]. Mature chondrocytes produce a matrix similar to NP cells with nasoseptal chondrocytes (NCs) being recently explored as an alternative non-disc cell source for IVD regeneration. NCs have a higher cell yield and sulphated glycosaminoglycan (sGAG) synthesis [11] when compared to articular chondrocytes (ACs) and MSCs [12]. Recently, Borrelli et al. also demonstrated that NCs resulted in higher matrix synthesis levels in comparison to NP cells [13]. The use of NCs for IVD regeneration may present many advantages compared to NP cells, including a high cell yield [11], the ability to produce tissue in an age-independent manner [14,15], enhanced matrix synthesis [13] and can be obtained with minimal donor site morbidity facilitating their autologous use. These intrinsic properties make NCs an attractive non-disc cell source for IVD regeneration.

Growth factors regulate IVD homeostasis, including extracellular matrix (ECM) synthesis and degradation, cell differentiation, or apoptosis [16]. Researchers have widely investigated the potential of various growth factor-mediated induction of cells to mimic the expression profile resembling native disc tissue. Studies have shown that the activation of transforming growth factor beta (TGFβ) signalling pathways delays IVD degeneration by increasing ECM synthesis [17,18] and hence members of the TGFβ superfamily may be potential candidates for driving discogenic differentiation and phenotype. TGFβ3 has been shown to maintain the phenotype of disc cells in organ culture [19] and when encapsulated with MSCs, they induced IVD regeneration in vivo [20]. Growth differentiation factor 5 (GDF5) and GDF6 are members of the bone morphogenetic protein (BMP) family. GDF5 and GDF6 play important roles in the development of bones, limb joints, skull and axial skeleton [21] and are also expressed in developing cartilage, tendons and ligaments [22,23]. In GDF5/6-knockout mouse models, the vertebral column showed severe lateral curvature and reduced staining of the IVDs indicating lower proteoglycan content [22]. These results suggest that GDF5 and GDF6 are required for normal development and maintenance of the IVD. In pellet culture, GDF5 results in reduced expression of the catabolic enzyme MMP13 in human chondrocytes [24]. There is also evidence to suggest that GDF-5 confers protection against NP cell apoptosis, promotes the synthesis of the main components of the extracellular matrix and can inhibit the activation of the NF-κB signalling pathway, thereby down-regulating the expression of inflammatory cytokines [25]. Inflammatory cytokines play an important role in the pathogenesis of disc degeneration by promoting matrix breakdown and recruitment of immune cells [26]. GDF6 supplementation has been shown to increase both proteoglycan and collagen production in 3D alginate bead cultures of human NP and AF cells [27]. Clarke et al. previously tested the individual effects of TGFβ3, GDF5, and GDF6 on discogenesis of bone marrow-derived MSCs and adipose-derived stem cells (ADSCs) reporting that GDF6 induced differentiation of these cell types towards an NP-like phenotype to a greater extent compared to other growth factors [28]. Furthermore, the authors found optimal expression of genes COL2 and ACAN when stimulated with 10 ng/mL TGFβ3 and 100 ng/mL GDF5 and GDF6. However, they did not investigate the effect of combining these growth factors. The overall objective of this study was to assess and compare the extent of matrix deposited by NP and NC microtissues stimulated by individual or combinations of the growth factors TGFβ3, GDF5, and GDF6.

#### **2. Materials and Methods**

#### *2.1. Nucleus Pulposus and Nasoseptal Chondrocyte Isolation and Expansion*

All tissue was sourced from a local abattoir and dissected within 24 h. NP cells were isolated from the IVDs of porcine spines (3 female donors; 4 months old). Briefly, NP tissues were harvested aseptically from the central nucleus pulposus of the IVD avoiding the annulus fibrosus region, washed with phosphate-buffered saline (PBS) and minced. Tissue fragments were placed in T-25 flasks containing Low Glucose-Dulbecco's Modified Eagle Medium (LG-DMEM) with 10% foetal bovine serum (FBS) and penicillin (100 U/mL) streptomycin (100 μg/mL) (PenStrep) and cultured in a humidified atmosphere at 37 ◦C and 5 %CO2. After ~7 days, cells had migrated from the tissue fragments, the flasks were washed with PBS and the cells were expanded to 80% confluence. NCs were isolated from porcine nasal septa, washed with PBS and minced. To isolate NC cells, minced tissue was digested in serum free LG-DMEM containing PenStrep and 3000 U/mL collagenase type II (Gibco, Invitrogen, Ireland) at a ratio of 10 mL/g of minced tissue. Digestion was performed under constant rotation for 3 h at 37 ◦C and subjected to physical agitation using a tissue dissociator (GentleMACSTM, Miltenyi Biotech, Surrey, UK) as previously described [29]. Digested tissue/cell suspensions were passed through a 40 μm cell strainer to remove tissue debris. Cell yield and viability in both cases were determined with a hemocytometer and trypan blue exclusion. Cells were seeded at an initial density of <sup>5</sup> × <sup>10</sup><sup>3</sup> cells/cm<sup>2</sup> in T-175 flasks in LG-DMEM supplemented with 10% FBS and PenStrep. Cultures were expanded to passage two in a humidified atmosphere at 37 ◦C and 5 %CO2.

#### *2.2. Fabrication of PDMS Microwell Moulds and Microtissue Formation*

Concave PDMS (polydimethylsiloxane) microwells were employed to generate selfassembled microtissues using lower numbers of cells than traditional pellet cultures or alginate bead models using a similar approach as described previously [30]. The microwell layout was designed using SolidWorks software (Solid Solutions Management Ltd, Dublin, Ireland). Moulds were designed with 100 microwells (10 × 10 array) containing 500 μm hemispherical concave microwells, with 1390 μm centre–centre distance and a depth of 1750 μm (Figure 1). A Formlabs Form 2 printer (Formlabs GmbH, Berlin, Germany) was used to 3D print the clear resin master stamps. After printing, the master stamps were washed in isopropyl alcohol for 20 min to remove excess residue and subsequently cured under an ultraviolet lamp (4 W, Uvitec, Cambridge, UK) for 1 h. 3D printed stamps were inserted into 12-well plates containing approximately 1 mL of PDMS (Sylgard® 184, Sigma-Aldrich, Arklow, Ireland), degassed in a vacuum oven and cured at 80 ◦C for 4 h. After cooling, the master stamps were removed carefully with the use of a tweezers. The PDMS microwell moulds formed were washed with ethanol and phosphate-buffered saline (PBS, pH 7.4), sterilised by UV-irradiation for 3 h and allowed to dry overnight in a laminar hood. The morphology and the diameter of the PDMS microwells formed were characterised by scanning electron microscopy (SEM) (SEM, SUPRA 35 V P, Carl Zeiss, Oberkochen, Germany). For microtissue formation, the moulds were placed in 12-well plates and seeded with different densities of NP and NC cells (1000, 2500, and 5000 cells/microwell or 100,000, 250,000, 500,000 cells/mould) and centrifuged at 300× *g* for 3 min to initiate microtissue formation. The plates were then cultured in a humified atmosphere at 37 ◦C and 5 %O2 in 1 mL of media supplemented with or without growth factors as described below. Preliminary experiments investigated the formation of microtissues with less than 1000 cells/microtissue. However, these were found to be inconsistent and lacked cohesion and integrity, making them unsuitable for further experiments.

#### *2.3. Growth Factor Stimulation*

Microtissues were cultured in 1 mL of LG-DMEM with PenStrep, 1.5 mg/mL bovine serum albumin (BSA; Sigma-Aldrich), 1× insulin-transferrin-sodium selenite (ITS; Sigma-Aldrich), 40 μg/mL L-proline, 100 nM Dexamethasone, 50 μg/mL L-ascorbic acid 2 phosphate, 4.7 μg/mL linoleic acid and were either supplemented with no growth factor (control) or with 10 ng/mL TGFβ3 [31]. Microtissues were maintained at 37 ◦C for 7 days under low oxygen (5 %O2) conditions. After 7 days of culture, microtissue size was evaluated through image analysis using Image J software (NIH, Bethesda, MA, USA). Initial results showed that lower cell density microtissues (1000 cells/microtissue) resulted in higher matrix production. Hence, this cell density was chosen for further experiments

and they were individually supplemented with either 100 ng/mL of GDF5, 100 ng/mL of GDF6 [28], 10 ng/mL TGFβ3 + 100 ng/mL GDF5, 10 ng/mL TGFβ3 + 100 ng/mL GDF6 or 10 ng/mL TGFβ3 + 100 ng/mL GDF5 + 100 ng/mL GDF6.

**Figure 1.** Fabrication of PDMS microwells for microtissue formation. (**A**) Design of the master stamp with negative pattern containing 10 × 10 microwells and the dimensions of the individual microwell. (**B**) 3D printed master stamp. (**C**) PDMS moulds with a 10 × 10 array for microtissue culture. (**D**) Microwells exhibited uniformity in size and shape as confirmed by SEM imaging.

#### *2.4. In-Silico Modelling of Microtissue Nutrient Microenvironment*

The in silico nutrient model of microtissues was created using COMSOL Multiphysics 5.6 (COMSOL Ltd., Cambridge, UK). Oxygen concentration at the microtissue boundary was dependent on multiple local environmental factors such as external level at the air/media interface, diffusion rate within the media, and the volume of media used. Glucose concentration and pH level at the microtissue boundary was dependent on the initial concentration of the media and the volume of media. Therefore, the in silico model was based on the in vitro geometry of the media-filled microwell containing an idealised radially symmetric microtissue with an initial cell seeding density of 1000, 2500 or 5000 cells. Since oxygen kinetics occur on a faster timescale than cell cycle kinetics the radius of the microtissue was assumed to be constant. The steady-state nutrient microenvironment was governed by coupled reaction-diffusion equations as described previously for disc cells [32]. Briefly, the metabolic rates were modelled as being dependent on local oxygen and pH levels by employing Michaelis–Menten equations derived and published previously [32–35]. Results for oxygen, glucose and pH levels were predicted and displayed as concentration contour plots through the midplane of the microtissues and graphically as a function of normalised radial distance through the microtissues.

#### *2.5. Live/Dead Analysis*

Cell viability was assessed using the LIVE/DEAD® viability/cytotoxicity assay kit (Invitrogen, Biosciences, Dublin, Ireland). Briefly, microtissues were washed with PBS and incubated in live/dead solution containing 2 μM calcein AM (live cell membrane, abs/em = 494/517 nm) and 4 μM ethidium homodimer-1 (dead cell DNA, ex/em = 528/617 nm; both from Cambridge Bioscience, Cambridge, UK) in PBS for 1 h. Samples were then washed in PBS, imaged with a Leica SP8 scanning confocal microscope at 515 and 615 nm channels and assessed using Leica Application Suite X (LAS X) Software (Version 3.5.5.19976).

#### *2.6. Quantitative Biochemical Analysis*

After 7 days in culture, microtissues were flushed from the moulds with PBS and frozen at −80 ◦C for further analysis. Microtissues were digested with 3.88 U/mL papain in 0.1 M sodium acetate, 5 mM L-cysteine-HCl and 0.05 M ethylenediaminetetraacetic acid (EDTA) (pH 6.0) (all from Sigma-Aldrich) at 60 ◦C under constant rotation for 18 h. DNA content was quantified using the Quant-iT PicoGreen® dsDNA (Invitrogen) assay. sGAG content was quantified using the dimethylmethylene blue dye-binding assay at pH 1.35 (Blyscan, Biocolor Ltd., Carrickfergus, Northern Ireland), with a chondroitin sulphate standard. Total collagen content was determined by measuring the hydroxyproline content. Samples were hydrolysed at 110 ◦C for 18 h in 12 M HCl, assayed using a chloramine-T assay [36] and

the collagen content determined using a hydroxyproline:collagen ratio of 1:7.69. Samples of media supernatants were also analysed for both sGAG and collagen content. DNA and sGAG contents were normalised based on the mass of all the microtissues (~100) from the entire microwell. The sGAG/collagen ratio was determined by dividing sGAG (μg) by collagen (μg).

#### *2.7. Histology and Immunohistochemistry*

Microtissues were treated with 4% paraformaldehyde (4 ◦C, 12 h) and washed in PBS. Microtissues were encapsulated in 2% agarose (Sigma-Aldrich) to facilitate handling, transferred to a cassette, dehydrated in a series of graded alcohols and finally wax embedded. Sections of 5 μm were stained with 1% alcian blue 8GX in 0.1 M HCl to assess sGAG deposition and with picrosirius red to evaluate collagen distribution. Collagen types I and II were assessed using immunohistochemistry techniques. Sections were treated with chondroitinase ABC (37 ◦C, 1 h) (Sigma-Aldrich), and non-specific sites were blocked using 5% BSA. Sections were incubated at 4 ◦C overnight with collagen type I (Abcam, Cambridge, UK) or collagen type II (Santa Cruz, Heidelberg, Germany) primary antibodies. The secondary antibody (Anti-Mouse igG biotin conjugate, Sigma-Aldrich) was applied for 1 h followed by incubation (45 min) with ABC reagent (Vectastain PK-400, Vector Labs, 2BScientific Ltd., Oxfordshire, UK). DAB peroxidase (Vector Labs, UK) was used as a developer.

#### *2.8. Statistical Analysis*

Statistical analysis was performed using GraphPad Prism (version 9.3.1, GraphPad Software, San Diego, CA, USA) software with 3–4 samples analysed for each experimental group. Two-way ANOVA was used for analysis of variance with Tukey's multiple comparisons test to compare between groups. Three technical replicates from three different porcine donors (biological replicates) were analysed. Results are displayed as mean ± standard deviation, with significance accepted at a level of *p* < 0.05.

#### **3. Results**

#### *3.1. Design and Fabrication of PDMS Microwells*

The master stamp prototype designed using SolidWorks software (Figure 1A) was fabricated using a Form 2 printer (Figure 1B). The printed stamps had a smooth and clear surface and were used to form the negative PDMS moulds. The PDMS moulds formed had 100 microwells (10 × 10 array) with 500 μm hemispherical concave microwells, with 1390 μm centre-centre distance and depth of 1750 μm (Figure 1C). The microwells had a uniform size with an average diameter of 504 ± 24 μm based on SEM images (Figure 1D).

#### *3.2. Morphology and Viability of Microtissues*

In terms of cell viability, results revealed high viability in low cell density microtissues with more dead cells being detected in the microtissues formed using higher cell numbers (Figure 2A). Moreover, the microtissues in the control group were loosely packed, exhibiting an irregular surface with loose arrangement of cells at the peripheral layers while the microtissues formed in the presence of TGFβ3 exhibited a smooth surface with compact organisation of the interior cells. The low cell density microtissues in the control group were also easily dissociated during the flushing process. The size of the microtissues increased with increasing cell numbers in both control and TGFβ3 treated groups for both NP and NC microtissues (Figure 2B). As expected, DNA content increased with higher content observed in higher cell density microtissues. NC microtissues exhibited increased DNA content (which indicates higher cell content) compared to NP microtissues, with a significant difference observed for 2500 and 5000 cells/microtissue cultured in the presence of TGFβ3 (Figure 2C).

**Figure 2.** Morphology and viability of microtissues formed using nucleus pulposus (NP) and nasoseptal chondrocyte (NC) cells. (**A**) Representative brightfield and fluorescent Live/Dead® images demonstrating the morphology and viability of microtissues in microwells seeded with 1000, 2500 and 5000 cells/microwell, after 7 days (**B**) Diameter of microtissues were analysed using ImageJ (**C**) DNA content (μg) was determined in digested microtissues from moulds at each timepoint (day 0 or day 7). Data represent mean ± SD of at least three independent experiments performed in triplicate. ! (*p* < 0.05) indicates significance compared to 1000 cells/microtissue. # (*p* < 0.05) indicates significance compared to the control group for the same cell number. \$ (*p* < 0.05) indicates significance compared to day 0 for the same cell number. \* (*p* < 0.05) indicates significance compared to the NP microtissues for the same experimental group and cell number.

#### *3.3. ECM Synthesis and Optimising Cell Density of TGFβ3 Stimulated Microtissues*

Total sGAG content increased with increasing cell densities, with NC microtissues achieving significantly higher sGAG levels for all cell densities compared to NP microtissues cultured in the presence of TGFβ3 (Figure 3A). sGAG normalised to DNA showed increased sGAG/DNA for lower cell density microtissues and the NC 1000 cells/microtissue in TGFβ3 achieved a significantly higher level compared to NP 1000 cells/microtissue (Figure 3B). A similar trend was observed for total collagen content with higher cell density microtissues resulting in higher collagen accumulation (Figure 3C). Interestingly, when normalised to DNA content, NP 1000 cells/microtissue in TGFβ3 was significantly higher compared to NC 1000 cells/microtissue in TGFβ3 (Figure 3D). Strong alcian blue staining in lower cell density NC microtissues in TGFβ3 confirmed the higher synthesis of sGAG (Figure 3E). Picrosirius red staining for collagen also corroborated the collagen biochemical data with dense staining observed in NP 1000 cells/microtissue supplemented with TGFβ3. Comparatively, low cell density microtissues resulted in better ECM synthesis on a per cell basis and in silico modelling was performed to determine the nutrient microenvironment within the microwells to elucidate how this may be influencing matrix synthesis.

#### *3.4. In-Silico Modelling of Microtissue Nutrient Microenvironment*

An in silico modelling analysis was performed based on the size of the microtissues for each of the cell densities investigated. As expected, an inverse relationship between metabolite concentration and cell density was observed and correlated well with the biochemical observations, whereby higher ECM on a per cell basis was observed with higher oxygen, glucose and pH levels. In terms of oxygen, minimum levels of 2.7 %O2, 1.8 %O2 and 0.9 %O2 for 1000, 2500 and 5000 cells/microwell were observed, respectively (Figure 4A,B). For glucose, minimum levels of 4.6 mM, 4.1 mM and 3.4 mM (Figure 4C,D) and for pH, 7.3, 7.2 and 7.0 were observed for 1000, 2500 and 5000 cells/microwell respectively (Figure 4E,F). Minimal gradient effects from the centre to the periphery of microtissues were observed.

#### *3.5. Viability and ECM Synthesis of Microtissues in the Presence of Different Growth Factors and Combinations*

Based on the observations of lower cell density microtissues resulting in higher cell viability and matrix production, a cell density of 1000 cells/microtissue was chosen for further experiments. Microtissues formed by NCs were more compact compared to NP microtissues. No significant difference in viability was observed among the microtissues cultured in the presence of various growth factors, with high cell viability observed on day 7 (Figure 5A). Weak alcian blue staining for sGAG was observed in NP microtissues supplemented with GDF5 and GDF6 growth factors in isolation. Increased positive staining was observed by co-stimulation with TGFβ3. For NC microtissues, GDF5 had a positive effect with minimal staining observed for GDF6. TGFβ3 co-stimulation resulted in intense staining and was comparatively stronger than NP microtissues (Figure 5B). Both NP and NC microtissues stained positive for collagen in the presence of GDF5, with higher deposition of collagen observed in the presence of GDF6 and when co-stimulated with TGFβ3 (Figure 5C).

#### *3.6. Quantitative Biochemical Analysis of Microtissues in the Presence of Different Growth Factors*

Both cell types responded positively to growth factor treatment relative to the no growth factor control (red solid line). For NP microtissues, no significant differences in DNA content were observed after 7 days irrespective of the growth factor type or combination. In contrast, a higher DNA content was observed in NC microtissues when cultured in the presence of GDF5 and all other combinations of growth factors relative to GDF6 (Figure 6A). Total sGAG was highest in NP microtissues stimulated with TGFβ3 + GDF5, although there were no significant differences between the groups investigated. Total sGAG levels in the NC microtissues were similar for all growth factor groups and combinations and were significantly higher compared to GDF6. Overall, sGAG levels for

NC microtissues were significantly higher compared to NP microtissues (Figure 6B). When normalised by DNA content (sGAG/DNA), no differences were found between growth factor groups for a specific cell type, although NC microtissues still exhibited higher content compared to NP cells (Figure 6C).

**Figure 3.** Biochemical and histological staining for sGAG and collagen of nucleus pulposus (NP) and nasoseptal chondrocyte (NC) microtissues stimulated with TGFβ3 on day 7. (**A**) Total sGAG (μg) content (**B**) Total sGAG/DNA (**C**) Total Collagen (μg) (**D**) Total Collagen/DNA. Data represent mean ± SD of at least three independent experiments performed in triplicate. ! (*p* < 0.05) indicates significance compared to 1000 cells/microtissue. # (*p* < 0.05) indicates significance compared to the control for the same cell number on day 7. \$ (*p* < 0.05) indicates significance compared to day 0 for the same cell number. \* (*p* < 0.05) indicates significance compared to the NP microtissues for the same experimental group and cell number. (**E**) Histological evaluation with alcian blue staining and picrosirius red to identify sGAG collagen content, respectively.

**Figure 4.** In silico modelling of the nutrient microenvironment within microtissues as a function of cell density. Predicted contour plots and radial profile through the midplane containing 1000, 2500 and 5000 cells for (**A**,**B**) oxygen (%O2) (**C**,**D**) glucose (mM) and (**E**,**F**) pH. The radial distance was normalised by the radius of the microtissue for each cell density.

In terms of total collagen content, for NP microtissues, a synergistic effect was observed through co-stimulation of growth factors (TGFβ3 + GDF5, TGFβ3 + GDF6 and TGFβ3 + GDF5 + GDF6) when compared to GDF5 and GDF6 alone. A similar synergistic effect was observed for NC microtissues when cultured in the presence of TGFβ3 + GDF5. However, the collagen content in NC microtissues treated with TGFβ3 + GDF6 and TGFβ3 + GDF5 + GDF6 were significantly higher than those treated with GDF5 alone. NC microtissues cultured in the presence of GDF6, TGFβ3 + GDF5 and TGFβ3 + GDF5 + GDF6 exhibited significantly higher collagen content compared to NP microtissues exposed to similar conditions (Figure 6D). Total collagen/DNA was significantly higher in NP microtissues cultured in the presence of TGFβ3 + GDF6, and no significant difference was observed when compared to TGFβ3 (red dashed line). However, for NC microtissues, GDF6-treated groups exhibited significantly higher concentrations compared to TGFβ3

(red dashed line), GDF5 and TGFβ3 + GDF5 + GDF6 treated groups (Figure 6E). In terms of sGAG/collagen level, which is a surrogate measure of NP-like matrix, with a higher ratio being desirable, no significant difference—irrespective of the growth factor treatment—was observed for NP microtissues. NC microtissues treated with GDF5 showed significantly higher concentrations compared to GDF6, TGFβ3 + GDF5 and TGFβ3 + GDF6. They also exhibited higher ratios compared to NP microtissues when cultured in the presence of TGFβ3 (red dashed line), GDF5, TGFβ3 + GDF5, TGFβ3 + GDF6 and TGFβ3 + GDF5 + GDF6 growth factor combinations (Figure 6F).

**Figure 5.** Live/Dead® and histological staining of nucleus pulposus (NP) and nasoseptal chondrocyte (NC) microtissues (1000 cells/microtissue) treated with no growth factor (control) or stimulated with TGFβ3, GDF5, GDF6, TGFβ3 + GDF5, TGFβ3 + GDF6, TGFβ3 + GDF5 + GDF6 on day 7. (**A**) Live/Dead® fluorescent images demonstrating cell viability. Histological evaluation with (**B**) alcian blue staining to identify sGAG and (**C**) picrosirius red to detect collagen.

**Figure 6.** Biochemical analysis of nucleus pulposus (NP) and nasoseptal chondrocyte (NC) microtissues stimulated with GDF5, GDF6, TGFβ3 + GDF5, TGFβ3 + GDF6, TGFβ3 + GDF5 + GDF6 on day 7. (**A**) DNA (μg), (**B**) Total sGAG (μg), (**C**) Total sGAG/DNA, (**D**) Total Collagen (μg), (**E**) Total Collagen/DNA and (**F**) sGAG/Collagen ratio. & (*p* < 0.05) indicates significance compared to other growth factors in the same experimental group, \* (*p* < 0.05) indicates significance compared to the NP microtissues for the same growth factor group. Solid red line indicates no growth factor treatment group (control); red dashed line indicates TGFβ3 stimulation group.

#### *3.7. Immunostaining for Collagen Type of Microtissues in the Presence of Different Growth Factors*

Irrespective of the growth factor used, immunostaining showed limited presence of collagen type I in both NP and NC microtissues (Figure 7A). NP microtissues also showed weak positive staining for collagen type II while NC microtissues stimulated with TGFβ3, GDF5 and TGFβ3 + GDF5 growth factors had higher and more dispersed collagen type II content (Figure 7B).

**Figure 7.** Immunostaining for collagen types I and II of nucleus pulposus (NP) and nasoseptal chondrocyte (NC) microtissues stimulated with TGFβ3, GDF5, GDF6, TGFβ3 + GDF5, TGFβ3 + GDF6, TGFβ3 + GDF5 + GDF6 on day 7. (**A**) Collagen type I and (**B**) Collagen type II.

#### **4. Discussion**

In this study, a microwell array system was used to promote the formation of selfassembled homogenous microtissues using two different cell types, nucleus pulposus cells and nasoseptal chondrocytes, and were subsequently assessed in response to growth factor stimulation. Growth factors TGFβ3, GDF5, GDF6 and a combination of these growth factors were assessed for their ability to induce matrix synthesis. The choice of growth factors selected was based on previous research in the field. For example, TGFβ3 is known to support the maintenance of the phenotype of disc cells in a rat lumbar organ culture [19] and MSCs supplemented with TGFβ3 support their differentiation towards an NP-like phenotype [37]. In addition, murine IVD explants treated exogenously with TGFβ3 have been shown to result in an up-regulation of aggrecan [38]. Studies have also shown that cells primed with TGFβ3 promoted higher sGAG and collagen synthesis [31,39] and also mitigated the detrimental effect of the harsh acidic microenvironment [40]. GDF5 (also known as cartilage-derived morphogenetic protein, CDMP1) and GDF6 (CDMP2) are members of the bone morphogenetic protein (BMP) family and are key regulators of cellular condensation and chondrogenesis [41]. GDF5 binds with higher affinity to BMP receptor type-1B (BMPR-IB), which controls the primary stages of cellular condensation and promotes cartilaginous tissue formation [42,43]. GDF6 has been shown to promote chondrogenesis and it positively regulates growth and maintenance of articular cartilage [44]. Both GDF5 and GDF6 have been shown to induce the expression of NP associated genes in MSCs and adipose-derived stem cells (ASCs) [45] with a healthy NP phenotype exhibiting stabilised expression of hypoxia inducible factor HIF-1α, the glucose transporter glut-1, the PG aggrecan (ACAN), type II collagen (COL2A), the signalling factor sonic hedgehog (SHH), the transcription factor Brachyury [T], keratins KRT18, KRT19, carbonic anhydrase CA12, and CD24 [46].

We selected the growth factor concentrations of TGFβ3, GDF5, and GDF6, based on previous work stimulating bone marrow-derived MSCs and adipose-derived stem cells (ADSCs) which reported optimal expression of genes COL2 and ACAN when stimulated with 10 ng/mL TGFβ3 and 100 ng/mL GDF5 and GDF6 [28].

This study confirms that the initial cell seeding density of microtissues plays an important role in regulating ECM deposition. Although the amount of sGAG and collagen produced were highest in high cell density microtissues, sGAG and collagen synthesised on a per cell basis (per DNA) were highest in low density microtissues (1000 cells/microtissues). Based on the in silico models, higher cell density microtissues experienced lower nutrient levels and this in turn could alter the metabolism of the cells [32,47], and inhibit or impede matrix synthesis [39,48]. Overall, higher levels of matrix synthesis were observed for NC microtissues compared to NP microtissues, and it was significantly pronounced for NC microtissues stimulated with GDF5 or in combination with TGFβ3. In agreement with previous studies, microtissues stimulated with TGFβ3 resulted in increased proteoglycan synthesis compared to the controls [49,50]. TGFβ3 is widely known to modulate cell proliferation, differentiation and ECM synthesis and in particular supports chondrogenic differentiation of various cell types including mesenchymal stem cells [37,51–54], through the activation of the Wnt signalling pathway [55,56]. It has also been reported that Smad2/3 cooperatively is one of the important signalling pathways stimulated by TGFβ3 that helps in developing and maintaining a chondrocytic phenotype [57].

Given that increased ECM synthesis was observed in low cell density microtissues, we next explored the effect of growth factors—GDF5, GDF6 in isolation or in combination with TGFβ3. A synergistic effect on DNA content (indicative of increased cell number), total sGAG and total collagen was observed for both cell types with higher values observed in groups stimulated with the addition of TGFβ3 to GDF5 or GDF6, and this was evident for both cell types. NP microtissues stimulated with TGFβ3 + GDF5 or TGFβ3 + GDF6 or TGFβ3 + GDF5 + GDF6 were comparable to TGFβ3 alone while the NC microtissues cultured in the presence of GDF5 alone were equally effective to TGFβ3 and had a synergistic effect when combined (TGFβ3 + GDF5). When these effects were evaluated on a per cell basis, the effects were diminished, with only a marginal increase observed in microtissues that were co-stimulated with growth factors. This was evident in the NC microtissues, but was not explicit in NP microtissues. Interestingly, Coleman et al. showed that MSCs co-stimulated with TGFβ3 (10 ng/mL) and GDF5 (150 ng/mL) resulted in a significant increase in sGAG/DNA by day 7 compared to TGFβ3 alone [58]. In this study, microtissues cultured in the combination of those growth factors did not result in any significant increase in proteoglycan synthesis on a per cell basis compared to TGFβ3 alone, irrespective of the cell type. This could be due to the lower concentration of GDF5 used (100 ng/mL) or could be attributed to the difference in the cell type or culturing conditions with lower oxygen and glucose concentrations. Another study by Jenner et al. showed that GDF5 stimulation caused an increase in cell numbers but collagen/DNA was significantly high in MSC scaffolds treated with TGFβ1 compared to GDF5 [59]. This trend was also observed in this study for NC microtissues, whereby TGFβ3 stimulation resulted in a significant increase in collagen/DNA compared to GDF5 only.

This study showed that GDF6 alone had little effect on DNA content (an indicator of cell number) for both NP and NC microtissues. This is in agreement with Bobacz et al., who reported no increase in cell numbers when articular chondrocytes were cultured in GDF6-supplemented media [60]. Gulati et al. showed that GDF6 at a concentration of 400 ng/mL resulted in a significant increase in proteoglycan accumulation in NP cells compared to the non-treated controls [27]. Similar observations were made in this study with both NP and NC microtissues exhibiting significantly higher sGAG/DNA in the presence of GDF6 compared to the non-treated controls.

sGAG/collagen was quantified to estimate the potential of the microtissues stimulated with various growth factors to produce the appropriate matrix type normally found in the disc. A higher ratio of sGAG composition of the ECM in comparison to the collagenous matrix is a predominant characteristic of NP tissue [45]. While no significant differences were observed for NP microtissues irrespective of the growth factor used, supplementation of NC microtissues with GDF5 or TGFβ3 alone or with the combination of growth factors (TGFβ3 + GDF5 + GDF6) yielded a significantly higher concentration. This implies that GDF5 or TGFβ3 alone could produce similar effects to the combination of growth factors indicating their potential to induce NC microtissues to form an NP-like matrix.

Recent research has shown that a combination of TGFβ1 and GDF5 significantly enhanced glycosaminoglycan content of human-derived MSCs, suggesting that this combination is optimal for MSC to NP cell induction [61]. Clarke et al. previously investigated the effect of TGFβ3, GDF5 and GDF6 growth factors on the discogenic differentiation of bone marrow- and adipose-derived MSCs (AD-MSCs). After 14 days, they observed that GDF6 resulted in increased expression of NP phenotypic genes KRT8, 18, and 19, FOX1 and CAXII compared to either TGFβ3 or GDF5 alone in both cell types with greater effects observed for AD-MSCs. However, they did not evaluate the effects of combining the growth factors [28].

In this study, GDF6 resulted in a higher sGAG/collagen ratio in NP microtissues, although this was not found to be statistically significant. On the other hand, GDF5 stimulation resulted in the highest sGAG/collagen ratio in NC microtissues, which was significantly higher than GDF6, TGFβ3 + GDF5 and TGFβ3 + GDF6. Interestingly, TGFβ3 only stimulation had a comparable effect to combinations with GDF5/GDF6. For NP microtissues, TGFβ3 also had a similar effect to other growth factors used in isolation and in combination, thereby demonstrating the potency of TGFβ3 in promoting matrix synthesis. The differential effects induced by the different growth factors, despite all belonging to the TGFβ superfamily, could be due to the difference in the signalling pathways activated. TGFβ3 is recognised by receptor type II, which activates the SMAD 2/3 signalling pathway, while GDF5 and GDF6 utilise the BMP receptor II, which activates the SMAD 1/5/8 pathway [28,62]. Activation of these distinct pathways ultimately results in different downstream signals, which may explain the differential effects observed in this study.

There are several limitations associated with the present study worth highlighting. First, longer term evaluation beyond 7 days may result in larger differences in microtissue maturation. However, this study is valuable as it demonstrates that microtissues stimulated with growth factors can produce a significant amount of ECM in a short time frame of 7 days, which would be compatible with a priming or conditioning strategy prior to implantation in the challenging disc microenvironment. We have previously shown that bone marrow stem cells (BMSCs) respond to growth factor stimulation (TGFβ3) under acidic conditions typically found in the degenerated disc (pH 6.8) [63]. A growth factor supplementation strategy with any of the growth factors assessed in this work may also prove to be beneficial. In addition, we have previously demonstrated that priming of BMSCs, NP and chondrocyte tissues with growth factors prior to a challenge or insult, confers better protection against acidic conditions, possibly due to the ECM being produced during the priming phase, thereby providing a protective niche [31,40,64]. In addition, further investigations on the gene expression profile of NCs would be valuable to identify whether typical NP markers are up-regulated.

#### **5. Conclusions**

In this study, we have developed and evaluated self-assembled microtissues stimulated with various growth factors and combinations. In general, while NP cells did respond to growth factor stimulation, NP cells were found to be relatively insensitive to the different growth factor types examined in isolation or in combination in this experimental setting. In contrast, NCs, which are more easily isolated in a less invasive manner and are more easily expandable, may prove to be an alternative cell type for NP repair and these cells could be primed with either GDF5 or TGFβ3 to enhance matrix synthesis.

**Author Contributions:** S.S. provided substantial contribution to study design, data acquisition data analysis and presentation, interpretation of data, drafting of the article, revising it critically and final approval. E.E.M. developed and performed the computational modelling, analysis, interpretation and presentation. C.T.B. is the overall project funding holder, takes responsibility for the integrity of the work as a whole from inception to finalised article, provided substantial contributions to study design, data presentation, interpretation of data, drafting of the article, revising it critically, and final approval. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by a Science Foundation Ireland Career Development Award (15/CDA/3476).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Single Cell RNA-Sequence Analyses Reveal Uniquely Expressed Genes and Heterogeneous Immune Cell Involvement in the Rat Model of Intervertebral Disc Degeneration**

**Milad Rohanifar 1, Sade W. Clayton 2, Garrett W.D. Easson 2, Deepanjali S. Patil 1, Frank Lee 1, Liufang Jing 1, Marcos N. Barcellona 1, Julie E. Speer 1, Jordan J. Stivers 2, Simon Y. Tang <sup>2</sup> and Lori A. Setton 1,2,\***


**Abstract:** Intervertebral disc (IVD) degeneration is characterized by a loss of cellularity, and changes in cell-mediated activity that drives anatomic changes to IVD structure. In this study, we used singlecell RNA-sequencing analysis of degenerating tissues of the rat IVD following lumbar disc puncture. Two control, uninjured IVDs (L2-3, L3-4) and two degenerated, injured IVDs (L4-5, L5-6) from each animal were examined either at the two- or eight-week post-operative time points. The cells from these IVDs were extracted and transcriptionally profiled at the single-cell resolution. Unsupervised cluster analysis revealed the presence of four known cell types in both non-degenerative and degenerated IVDs based on previously established gene markers: IVD cells, endothelial cells, myeloid cells, and lymphoid cells. As a majority of cells were associated with the IVD cell cluster, sub-clustering was used to further identify the cell populations of the nucleus pulposus, inner and outer annulus fibrosus. The most notable difference between control and degenerated IVDs was the increase of myeloid and lymphoid cells in degenerated samples at two- and eight-weeks post-surgery. Differential gene expression analysis revealed multiple distinct cell types from the myeloid and lymphoid lineages, most notably macrophages and B lymphocytes, and demonstrated a high degree of immune specificity during degeneration. In addition to the heterogenous infiltrating immune cell populations in the degenerating IVD, the increased number of cells in the AF sub-cluster expressing *Ngf* and *Ngfr*, encoding for p75NTR, suggest that NGF signaling may be one of the key mediators of the IVD crosstalk between immune and neuronal cell populations. These findings provide the basis for future work to understand the involvement of select subsets of non-resident cells in IVD degeneration.

**Keywords:** intervertebral disc degeneration; single-cell RNA sequencing; cell type

#### **1. Introduction**

The intervertebral disc (IVD) provides for motion and flexibility in the spine [1]. Degeneration of the IVD is characterized by a loss of cellularity, changes in composition, and loss of hydration that is manifested as changes in disc height and MRI signal intensity. These features of IVD degeneration in the lumbar and cervical spine can contribute to mal-alignment of the adjoining vertebral bodies and affect the integrity and health of the IVD and adjacent neural structures, leaving the IVD vulnerable to progressive damage during the loading conditions of daily living [2–4]. The cell-mediated changes affecting IVD growth, extracellular matrix synthesis, and metabolism can affect IVD proteolytic degradation, and are of great interest for identifying the relevant cellular and molecular targets to support IVD repair and regeneration.

The cell density of the nucleus pulposus (NP) region is at the lowest in the IVD, and are the principal contributors to the loss of cellularity with age. NP cells have long been considered to share features with chondrocytes in the adult IVD, due to their roundedness and high expression of multiple chondrocyte markers including SOX9, type II collagen,

**Citation:** Rohanifar, M.; Clayton, S.W.; Easson, G.W.; Patil, D.S.; Lee, F.; Jing, L.; Barcellona, M.N.; Speer, J.E.; Stivers, J.J.; Tang, S.Y.; et al. Single Cell RNA-Sequence Analyses Reveal Uniquely Expressed Genes and Heterogeneous Immune Cell Involvement in the Rat Model of Intervertebral Disc Degeneration. *Appl. Sci.* **2022**, *12*, 8244. https:// doi.org/10.3390/app12168244

Academic Editor: Benjamin Gantenbein

Received: 3 July 2022 Accepted: 18 July 2022 Published: 18 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and aggrecan upon bulk mRNA transcriptional profiling [5–9]. NP cells are derived from the embryonic notochord and express transcriptionally unique markers that reflect their notochordal origins, including CD24, cytokeratins, and the brachyury and FOX transcription factors [10–13]. Cells of the annulus fibrosus (AF) are derived from mesenchyme and express some overlapping and some differing molecular markers from adjacent NP cells. Indeed, there have been numerous studies designed to identify the unique cellular phenotypes of NP and AF cells, as well as non-resident cells involved in degeneration and repair, using tools of bulk RNA transcription or RNA-sequencing, proteomic profiling, flow cytometry, and protein assays [5,14–18]. The improved understanding of cell-specific molecular changes is crucial for identifying important protein targets that drive progress degeneration and molecular targets for stem cell-mediated IVD repair.

In recent years, single-cell RNA sequencing (scRNA-seq) of IVD cells, both from anatomically distinct regions and for comparisons with pathology, has been used to characterize the phenotype and abundance of distinct cell populations in the IVD in human and animal tissue sources [5,15,17,19–21]. When cells are isolated rapidly and immediately prepared for RNA-sequencing, this technique has the potential to reveal the endogenous RNA expression of resident cells and thus identify differences in cell sub-populations in the native IVD [5,15,17,19–22]. Results of sc-RNA-seq from rat and bovine IVD tissues, and from non-degenerate human IVD tissues, rely on mapping cells with similar transcriptomic profiles into "clusters" in order to name and number discrete cell populations in the non-degenerate IVD. While the number of clusters is highly variable in the absence of a consensus-based approach, it is most commonly observed that the majority of IVD cells map to a few clusters based on their expression of common genes (50–99% of isolated cells). Studies have used this approach to identify small but meaningful stem cell or progenitor cell populations in the IVD [15,20], or to identify relationships between identified clusters and primary cell functions of matrix regulation, stress responses, cell cycle, and more [22].

Only a few studies have used sc-RNA-seq to better understand the progression of IVD degeneration at the cellular level. In donor human tissues, sc-RNA-seq has been used to identify genes associated with cells of the degenerative IVD with some findings that corroborate prior work including the unique associative expression for CTGF, S100A1/A2, and TNF receptors in cells from degenerated IVDs [5,22]; both human studies and a study of a rodent bipedal IVD degeneration model also present surprising findings for novel gene involvement that have not yet been confirmed in other studies or species. In some cases, sc-RNA-seq identifies cells that express macrophage markers at high levels [5] consistent with prior literature showing increased involvement of macrophages in animal models and human tissues with IVD degeneration [23]. The IVD is an aneural, alymphatic, and avascular structure upon the termination of development and growth, however, it can remain immunologically separate from the host over the lifetime of an individual. The presence of macrophage markers in degenerated IVD tissues and cells suggests integrity of the IVD may be disrupted such that cells of the systemic circulation may infiltrate and reside in the IVD. These cells can be expected to be minor and difficult to localize via immunostaining or RT-PCR; for this reason, sc-RNA-seq provides the potential to identify the unique phenotypes of these relatively small cell populations in the IVD.

For these reasons, we performed scRNA-seq analysis on cells obtained from rat lumbar IVD tissues at two time points following induction of IVD degeneration via stab injury [24–28]. In brief, we identified four major cell types in the non-degenerate and degenerative IVD tissues for 2 weeks following surgery, based on the expression of "classical" cell-specific markers for NP and AF cells as identified previously [5]. The cell clustering scheme was held fixed between control and IVD degeneration conditions to estimate the differences in cell numbers and major gene expression levels for respective cell populations within each cluster. Further, we used sub-clustering to identify discrete native cell types within the IVD, and gene set analysis to identify the numerically minor cells of the degenerated IVD and their changes from 2 to 8 weeks after induced IVD degeneration. This study is innovative for utilizing sc-RNA-seq to reveal both temporal and spatial changes to cell

presence in the IVD following onset of disc degeneration, and for presenting data on a specimen-specific basis.

#### **2. Materials and Methods**

#### *2.1. IVD Tissue Preparation and Single-Cell Isolation*

Rats (n = 8, male Sprague-Dawley, 8-weeks-old) underwent surgery for retroperitoneal exposure of the L4-L6 lumbar spine in protocols approved by the Washington University IACUC. The L4-5 and L5-6 IVDs each received a single unilateral stab injury via a 27G needle. The rats were allowed to recover postoperatively for either 2 (n = 4) or 8 weeks (n = 4) (Figure 1). Following the respective recovery periods, the animals were euthanized and the IVD tissues were isolated and pooled from L4-5 and L5-6 (LDP = Lumbar disc puncture), or from L2-3 and L3-4 (CON = control), with care to remove peripheral muscle, ligaments, and attachments by microdissection. The isolated IVD tissue was digested in medium containing 0.2% type 2 collagenase (Worthington Biochemical, Lakewood, NJ, USA) and 0.3% pronase (Roche, Basel, Switzerland) for <4 h total at 37 ◦C and 5% CO2. After digestion, a cell pellet was obtained by centrifuging the medium for 10 min at 400 rcf. The media was aspirated, and cells were resuspended in PBS (Sigma Aldrich, St. Louis, MO, USA), followed by filtering through a 70 μm filter to get rid of undesired cell debris. The flowthrough was again centrifuged for 10 min at 400 rcf, and the resulting cell pellet was resuspended in PBS. The cells isolated from pooled IVDs for each rat were kept together and on ice and considered as a separate sample for a total of n = 16 samples in total (n=4 for CON and LDP at 2 weeks post-surgery; n = 4 for CON and LDP at 8 weeks post-surgery). Samples of IVD cells were immediately transported on ice to the sequencing facility. More than 80% of cells were determined to be viable. All analyses described below were performed on individual rat cell samples unless otherwise indicated.

**Figure 1.** Flowchart depicting an overview of the single-cell RNA-seq experimental design.

#### *2.2. cDNA Library Generation and Single-Cell RNA-Sequencing (scRNA-Seq) with Data Standardization*

The input samples were submitted to the Washington University Genome Technology Access Center to obtain and sequence the cDNA libraries (10xGenomics, 3'v3.1; Illumina NovaSeq S4) according to established protocols. For CON and LDP samples, only cells with <10% mitochondrial features, gene counts of 500–3000 and UMI (unique molecular identifier) counts between 500–30,000 were determined to be of sufficient quality to include in the analyses (Partek FlowTM, Partek Inc., St. Louis, MO, USA). Under these selection criteria, the total number of cells from rat samples was 42,560 and 51,365 for CON and LDP at 2 weeks post-surgery, 27,064 and 28,293 for CON and LDP at 8 weeks post-surgery, respectively.

#### *2.3. Dimensionality Reduction and Unsupervised Clustering*

In order to identify resident cell populations, we first performed Principal Components Analysis (PCA) on pooled data from n = 4 CON samples at 2 weeks, and then again for the 8 weeks samples. PCA was performed to reduce dimensionality of the data for CON at each timepoint with selection of optimum number of principal components after normalization; this process was repeated for the n = 4 LDP samples at 2 and at 8 weeks. Cluster analysis was then performed with the iterative procedure of K-means clustering on CON cells at 2 weeks to identify optimal separation of cell clusters. A total of five major cell clusters were identified as providing for optimal separation of the data for CON samples. The percentage of cells mapped to each cluster was evaluated for each tissue sample at both 2 and 8 weeks post-surgery as a test of inter-subject variability and appropriateness of the clustering scheme. Annotated cluster data were visualized with uniform manifold approximation and projection (UMAP) analyses.

#### *2.4. Cell Identification*

The nomenclature classifying cell types within each identified cluster for CON samples at 2 weeks was based on a set of "marker genes" based on reports within the literature (see Table 1). For CON samples, we first validated the presence of known cell subsets based on prior scRNA-seq studies of intervertebral cells, endothelial cells, myeloid and lymphoid cells of the rat IVD [5,15]. Violin plots of gene expression levels for all cells within each cluster were used to define the majority cell type in each cluster.

**Table 1.** Genes examined for definition of cell subtypes in the IVD of rat samples (i.e., marker genes). Percentages are numbers of total cells assigned to each cluster averaged across four rats. CON—cells from control (unoperated) IVDs; LDP—cells from lumbar disc punctured IVDs at 2 or 8 weeks following surgery.


As a majority of cells in the CON populations were identified with a cluster named intervertebral disc cells, we sought to further sub-cluster this population in a post-hoc process separate from the data for cells mapped to the other clusters. IVD cells for CON samples at 2 weeks were subjected to PCA to select an optimal number of principal components, followed by K-means clustering on pooled data to identify further separation of cells in the intervertebral disc clusters based on the lowest Davies–Bouldin index.

The above-described process for clustering and sub-clustering was repeated for CON samples at 8 weeks post-surgery following the same cluster nomenclature and marker genes. We sought to retain five major cell clusters for samples from the LDP population at both 2 and 8 weeks post-surgery, in order to directly compare cell numbers and gene expression levels against those of the CON population. In large part, the nomenclature for these five cell clusters did not differ from that of CON populations; however, the numbers of cells mapped to each cluster and gene expression for marker and novel genes varied substantially.

#### *2.5. Gene Specific Analysis*

Differential gene expression analyses were performed with the GSA toolbox in Partek FlowTM to identify the differentially expressed genes within each rat sample between CON and LDP at each of 2 and 8 weeks post-surgery. The genes exhibiting the greatest fold-change were identified (log2 (fold change)) for each sample at 2 weeks or 8 weeks with a statistical significance of false discovery rate (FDR) *p*-value < 0.001. This set of genes identified as most highly differentially expressed differed among samples with some overlap as described in the results. Genes with the highest fold-difference between CON and LDP that are common across samples at 2 or 8 weeks are presented here.

#### *2.6. Gene Set Analyses*

Gene set analyses were performed using the Gene Ontology resource (https://www. uniprot.org/help/gene\_ontology (accessed on and after 1 September 2021) (UniProt Database) to predict the physiological roles and molecular processes of the upregulated gene networks. The top 50 differentially regulated genes for the respective cluster or subclusters were used to define the relevant cellular processes.

#### *2.7. Statistical Analyses*

Optimal separation of cells into a minimum number of distinct clusters was identified via the minimal Davies-Bouldin index, as indicated above for CON cells at 2 weeks and for IVD cells when separated into sub-clusters. Data for the magnitude of gene expression for cell-specific markers, e.g., genes defining B and T lymphocytes, were tested for differences between CON and LDP populations at each of 2 and 8 weeks post-surgery by the Student's *t*-test or the Mann–Whitney test as appropriate (GraphPad PrismTM). Likewise, similar comparisons were made for the expression of *Ngf* and *Ngfr* (p75NTR) in matched subclusters between CON and LDP populations. Chi-squared tests were used to compare the proportions between matched cell clusters/subclusters across CON and LDP populations. Gene set enrichment analyses were performed by matching genes to defined pathways based on known cell functions and reporting overrepresented cell numbers in groups via a Fisher's exact test Partek FlowTM.

#### **3. Results**

#### *3.1. Cell Populations Identified in IVD Tissues from 2 Weeks Post-Surgery*

Samples in control IVD tissues 2 weeks post-surgery. Unsupervised K-means clustering was used for the CON samples to group the cells into clusters with equal variance. Results show that five clusters separated all cells into discrete populations. There were identified four named clusters: intervertebral disc cells, endothelial cells, myeloid, and lymphoid cells (Figure 2a); cluster 5 was annotated as "other cells" due to the lack of common identification markers. Identification of cell populations within each cluster was based on known gene markers that were pre-selected as described in methods (Figure 2b). The identified clusters had a distribution of cell counts that was very similar across samples for the four CON samples derived from four separate rats (Figure 2c). Each sample showed a similar cell distribution profile where intervertebral disc cells contributed to most of the cluster (~87.7%), followed by endothelial cells (~7.3%), myeloid (~1.0%) and lymphoid cells (~0.2%).

**Figure 2.** (**a**) UMAP plot representation of five cell types within CON group at two weeks after surgery. (**b**) Expression levels of marker genes within intervertebral disc cells, endothelial cells, myeloid, and lymphoid cells. (**c**) Pie chart showing the distribution of each cell analyzed for the four rat samples in the CON group. (**d**) The expression level of marker genes for three subclusters of intervertebral disc cells: nucleus pulposus, inner annulus fibrosus, and outer annulus fibrosus.

The selected marker genes were chosen based on *a priori* knowledge for identifying the discrete cell clusters (Figure 2b). There was some commonality in extracellular matrix genes across clusters, primarily for *Col2a1* and *Comp*; this observation is consistent with the prevalent expression of these extracellular matrix proteins in different tissues and suggests that relying on non-extracellular matrix genes as differentiating markers may be preferable. Markers selected as lymphoid and myeloid markers proved to be effective in separating cells into clusters with more separation.

As more than 85% of all cells mapped to the "intervertebral disc cells" cluster, we sought to sub-cluster this population to identify the known cell populations. The IVD cell cluster was partitioned into three distinctively adjacent subclusters that were classified as nucleus pulposus (NP), inner annulus fibrosus (IAF), and outer annulus fibrosus (OAF) based on relative expression levels of classical IVD cell markers (Figure 2d). Subcluster 1 was identified to be NP cells based on a relatively higher expression for *Acan* and *Col2a1* while also expressing lower *Col1a1* compared to other IVD cells [8,29–32]. Additionally, these cells highly express *Cd24*, a classical marker for NP cells [13,33–35], previously identified markers of the chondrocyte lineage including S100B, and a novel marker for bovine NP tissues, PDIA4 [19,36,37]. Cells mapping to the IAF cluster were similar in the pattern of expression to that of NP cells, differing with modestly increased *Col1a1* expression and decreased presence of *Cd24*. Finally, OAF cells were identified based on the expression of a number of previously identified markers for this cell type including *Igfbp5*, *Igfbp6*, *Lum*, *Myoc*, and *Fibin* (Figure 2d) [8,19,29,31,32,38].

Samples in degenerated IVD tissues. Samples harvested from the intervertebral discs of the LDP rats at 2 weeks post-surgery underwent unsupervised K-means clustering as applied to the CON group. Five clusters provided for optimal separation of the LDP cells based on known gene markers (Figure 3a,b). The IVD cell cluster contributed to the greatest proportion of the cells (~82.4%), followed by endothelial (~6.4%), myeloid (1.4%), and lymphoid clusters (~8.2%). The myeloid and lymphoid clusters exhibited the most notable difference 2 weeks after LDP (Figure 3c).

Differential gene expression between LDP and CON samples on a subject-specific basis. To directly assess differences in mRNA levels for cells from the CON (L2-3 and L3-4 levels) and LDP (L4-5 and L5-6 levels) groups within each rat sample, we performed specimen-specific differential gene expression (Figure 4a,b). Differential gene analyses demonstrated that between 700–900 genes were more highly expressed in LDP cells at 2 weeks after surgery, as compared to their counterparts amongst the CON cell population (FDR *p*-value < 0.001, Figure 4c). Analyses also showed between 250 and 1200 genes to be more highly expressed in CON cells as compared to the LDP cell population at 2 weeks after surgery (FDR *p*-value < 0.001, Figure 4c). While there was variability in differentially expressed genes across samples, we identified a subset of 25 commonly upregulated genes in LDP populations of cells from two samples with the highest total numbers of sequenced cells (i.e., samples #3 and #4); expression levels for these genes against their cell identification are shown in the heat map in Figure 4d (representative results shown for sample #4). Of these 25 genes with upregulated expression levels in LDP compared to CON cells, the majority (22/25) were associated with cells of the lymphoid or myeloid cell clusters (Figure 4d).

**Figure 3.** (**a**) UMAP plot representation of five cell types within the LDP group at two weeks after surgery. (**b**) The expression level of marker genes within intervertebral disc cells, endothelial cells, myeloid, and lymphoid cells in LDP groups at two weeks after surgery. (**c**) Percentage of all analyzed cells mapped to each cell type for CON and LDP groups at two weeks after surgery. Inset provided to highlight apparent differences in CON and LDP groups for immune cell populations.

Immune cell involvement in differences between LDP and CON gene expression patterns. We determined that immune cells were contributing to the upregulation genes due to LDP. Then, we classified the immune cell sub-types as either the myeloid or lymphoid clusters. There were significant increases in the number of cell mapping to the myeloid and lymphoid clusters for LDP samples as compared to CON samples at 2 weeks after surgery (*p* < 0.00001; Chi-square test), indicating an increase of infiltrating immune cells post-injury. The lymphoid cluster showed the largest difference in cell number between CON and LDP conditions, where the LDP samples had a cell count larger than 4.6× that of CON samples (*p* < 0.00001; Figure 5a). For the myeloid cells, we performed statistical analysis on the expression level of the pan-macrophage markers (*Cd4*, *Cd14*, *Cd68*, and *Lyz2*), M1 polarization (*Tlr2*, *Tlr4*, *Cd80*, and *Cd86*), and M2 polarization (*Arg1*, *Mrc1* (*Cd206*), *Mrs1* (*Cd204*)) cell surface markers (Figure 5b–d). The number of cells expressing pan-macrophage markers (*Cd4*, *Cd14*, *Cd68*, and *Lyz2*), as well as their respective levels of expression, indicated increased macrophages and heightened activity (Figure 5b). Furthermore, the number of cells expressing the canonical M1 cell surface markers *Tlr2*, *Tlr4*, and *Cd80* were not

different between groups, but the number of *Cd86*+ cells (*p* < 0.00001) and their expression were higher in the LDP group (Figure 5c) [39]. Despite the similarity in the number of cells expressing M2 polarization markers between LDP and CON groups, including *Arg1*, *Mrc1* (*Cd206*), and *Mrs1* (*Cd204*), the expression levels were either indistinguishable (*Arg1*, *Cd206*) or lower (*Cd204*) compared with the CON group (Figure 5d) [39,40].

**Figure 4.** Volcano plot depicting specimen-specific differenti gene expression between CON and LDalP cells at two weeks after surgery. Two representative samples are shown—(**a**) sample 3 (left) and (**b**) sample 4 (right) (Red dots represent genes expressed at higher levels in the LDP group while blue dots represent genes with higher expression levels in CON group) (**c**) Venn diagram represents the total number of commonly upregulated and downregulated genes for LDP and CON cells in across multiple samples (i.e., samples 3 and 4 here). (**d**) Heatmap showing the pattern for a subset of commonly upregulated genes with LDP compared to CON groups as a function of cell grouping (*x*-axis).

**Figure 5.** The identification of immune cell subtypes in the myeloid and lymphoid clusters. (**a**) Diagram showing the differences in cell numbers mapped to the myeloid and lymphoid clusters between the control (CON) and degenerated (LDP) samples at two weeks after surgery. (**b**) Panmacrophage marker expression and cell count differences between CON and LDP groups were used to confirm the identification of cells in the myeloid cluster. Macrophage polarization was assessed by identifying the presence of either (**c**) M1 or (**d**) M2 polarization cell surface markers, and the number of cells that expressed them. The presence of (**e**) B and (**f**) T lymphocytes was identified in the lymphoid cluster by measuring pan and subtype-specific marker expression levels and cell count changes. Numbers underneath each plot denote the number of cells expressing each marker. A Student's *t*-test or a Mann–Whitney test was run to determine significance of marker expression levels. \* = < 0.05, \*\* = <0.01, \*\*\* = <0.001, \*\*\*\* = <0.0001.

For the lymphoid clusters, we likewise performed statistical analysis for the B lymphocytes (*Cd19*, *Cd79b*, *Blnk*, *Cd38*, *Cxcr4*) and T lymphocytes (*Cd3g*, *Cd4*, *Cd8b*) markers by assessing the number of cells that expressed them (Figure 5e,f). B cell receptor and co-receptor, *Cd79B* and *Cd19* in LDP samples, showed greater levels of gene expression compared to CON cells, and were expressed in the majority of the lymphoid cells (Figure 5e). Markers of B lymphocyte subtypes, *Blnk*, *Cd38*, and *Cxcr4*, all showed significantly greater gene expression in LDP samples concurrent with large increases in cell numbers between CON and LDP at 2 weeks after surgery (Figure 5e). As for T lymphocytes, the gene expression and numbers of cells expressing pan T cell markers, *Cd3g*, the T cell receptor, were both significantly increased (Figure 5f). A small number of cells expressed *Cd4* or *Cd8b* subtype markers with significant increases in gene expression (Figure 5f).

Differential *Ngf* expression in IVD cells in LDP populations. The number and marker expression of cells mapping to the IVD subclusters were evaluated between CON and LDP populations at 2 weeks after surgery. The number of cells expressing *Ngf* and *Ngfr* (p75NTR) and their relative expression were used to compare the IVD subclusters including the NP, IAF and OAF (Figure 5a,b). The Gene ontology functions ranked from top to bottom by enrichment score derived from Fisher's exact test. The five most highly enriched processes identified by the Gene Set analyses were related to immune responses and associated regulatory processes (Figure 6c), consistent with the known infiltration of monocytes into the injured IVD post-surgical insult [41–43]. Genes associated with the top three regulated processes were upregulated in LDP samples and that are large drivers of the observed results are *Cd19*, *Cd37*, *Rag1*, *Blnk*, and *Myb*, as T- and B-cell associated marker genes.

While few differences were observed in relative expression levels of *Ngf* across the subclusters for the LDP group, the number of cells expressing *Ngf* increased following injury (*p* < 0.00001; Figure 6a). In CON cells from the NP cell sub-cluster, just 7.9% of cells expressed *Ngf* which increased to 12.9% following injury (*p* < 0.00001). *Ngf* is believed to be important for promoting neurite infiltration into the otherwise avascular and aneural tissues of the NP [38]. In OAF cells, the incidence of *Ngf* expressing cells increased from 1.0% of cells in CON samples at 2 weeks after surgery to 3.8% of cells in LDP samples, a greater than three-fold increase (*p* < 0.00001). Additionally, the gene expression of *Ngf* was significantly higher for all cells in the OAF subcluster of LDP samples (*p* < 0.00001). Likewise, the expression of *Ngfr*, which encodes for p75NTR—a low-affinity receptor for *Ngf*, was elevated in the OAF of LDP samples. Further, *Ngfr* was expressed in 2.3% of CON cells in the OAF sub-cluster increasing to 7.6% of cells expressed in the LDP injury group (*p* < 0.00001; Figure 6b).

**Figure 6.** (**a**) *Ngf* and (**b**) *Ngfr* (p75NTR) expression and cell count counts for samples from CON and LDP populations at two weeks after surgery. Numbers underneath each plot denote the number of cells expressing each marker. A Student's *t*-test or a Mann–Whitney test was run to determine significance of marker expression levels. \*\*\*\* = <0.0001. (**c**) Representative plot of gene set enrichment analysis of upregulated genes corresponding to one LDP sample. Gene ontology functions are ranked from top to bottom by enrichment score derived from Fisher's exact test; The <sup>−</sup>*logp value* <sup>10</sup> is represented by color where green denotes the highest level of statistical significance. The size of each node represents the number of genes in each gene ontology category.

#### *3.2. Cell Populations Identified in IVD Tissues from 8 Weeks Post-Surgery*

Samples in control IVD tissues. As described for samples harvested at 2 weeks postsurgery, unsupervised K-means clustering was applied to the CON group at 8 weeks post-surgery. The four identified clusters had a distribution of cell counts that was very similar across samples for the four CON samples at 8 weeks post-surgery (Figure 7a); note that the absence of a 5th identified cluster was related to a cluster (previously "endothelial

cell cluster") with no mapped cells within. Violin plots of marker genes expression levels for all cells within each cluster were used to define the majority cell type in each cluster (Figure S1). The IVD cells contributed to most of the cluster (~83.7%) of all CON cells (8 weeks), followed by myeloid (~1.2%) and lymphoid cells (~8.4%). These differences from samples at 2 weeks post-surgery reflect a natural aging of the spine in the rat.

**Figure 7.** UMAP plot representation of all cell types within (**a**) CON and (**b**) LDP groups at 8 weeks after surgery. (**c**) Volcano plot depicting genes differentially expressed between sample-specific LDP relative to CON populations for one sample at 8 weeks after surgery (Red dots represent genes expressed at higher levels in LDP while blue dots represent genes with higher expression levels in CON group) (**d**) Heatmap showing the pattern for a subset of commonly upregulated genes with LDP compared to CON groups at 8 weeks after surgery as a function of cell grouping (*x*-axis).

Samples in degenerated IVD tissues. Samples harvested from the LDP segments of the rats at 8 weeks post-surgery underwent unsupervised K-means clustering as described above, with the observation that five clusters provided for optimal separation of the data (Figure 7b). The IVD cells contributed to most of the cluster (~87.4%), followed by endothelial cells (~2.3%), myeloid (0.8%) and lymphoid cells (~5.1%). Identification of cell populations within each cluster was based on known gene markers (Figure S2). These differences in sample numbers reflect the remodeling response to LDP at 8 weeks after surgery.

Differential gene expression between LDP and CON samples on a subject-specific basis. Differential gene expression between CON and LDP cells at 8 weeks after surgery demonstrated that 86 genes were significantly upregulated in LDP compared to CON, while 128 genes were more highly expressed in CON than in LDP samples (representative plot shown for sample #8 in Figure 6c). A subset of 25 genes was identified as commonly upregulated in LDP compared to CON samples and used to generate a heat map, showing that the majority of genes identified were associated with cells of the immune system, including *Blnk* and *Cd79b* used to identify immune cells, along with *Myb*, *Rag1*, and *Cecr2*. These findings point to similarities in upregulated genes at the 2 and 8 weeks timepoints after LDP surgery; however, the identification of a lower number of upregulated genes and lesser fold effect changes suggest that the dramatic differences observed at 2 weeks after surgery in LDP cells may reflect an acute response to the injury.

The heatmap shows that the upregulated genes mostly belong to the lymphocyte cluster (Figure 7d). Therefore, we explore the presence of B and T lymphocyte-specific markers within the lymphoid clusters. B cell receptor and co-receptor, *Cd79b* and *Cd19* in LDP samples, showed greater levels of gene expression compared to CON cells, and were expressed in the majority of the lymphoid cells (Figure 8b). *Blnk*, *Cd38*, and *Cxcr4* all showed significantly greater gene expression levels in LDP samples concurrent with large increases in cell numbers between CON and LDP at 8 weeks after surgery (Figure 8b). As for T lymphocytes, the only gene expression is related to *Cd3g* with zero expression level for *Cd4* and *Cd8b* in both CON and LDP (Figure 8c). These data suggest that the lymphocytes cluster is almost composed of B cells by less sign of T cells with respect to 2 weeks.

Gene set analyses revealed the five most highly significant processes identified point to B cell activation, receptor signaling, and differentiation (Figure 8a). Genes associated with the top three regulated processes that were found to be upregulated in LDP samples and that are large drivers of the observed results are *Cd19*, *Cd7bB*, *Rag1*, *Blnk*, *Dppr4*, *Il7r*, *Lef1*, and *Myb*.

**Figure 8.** (**a**) Representative plot of gene set enrichment analysis of upregulated genes corresponding to one LDP sample. Gene ontology functions are ranked from top to bottom by enrichment score derived from Fisher's exact test; the <sup>−</sup>*logp value* <sup>10</sup> is represented by color where green denotes the highest level of statistical significance. The size of each node represents the number of genes in each gene ontology category. (**b**) B and (**c**) T lymphocytes were identified in the lymphoid cluster of one rat sample taken at 8 weeks after surgery (sample #8) by measuring pan and subtype-specifi marker expression levels and cell count changes. Numbers below each plot denote the number of cells expressing each marker.

#### **4. Discussion**

Cells were isolated from the degenerated IVDs from a widely used rat model of lumbar disc puncture [24–28,44]. Single-cell RNA-seq was then performed on over 149,282 nondegenerate and degenerate IVD cells harvested from two timepoints following surgery to allow the progression of IVD degeneration. In brief, we identified four major cell types in the non-degenerate and degenerative IVD tissues ofc 2 weeks following surgery, based on the identification of genes previously identified as "classical" cell-specific markers [15]. For example, IVD cells are closely clustered around classical "chondrocyte-like" markers including *Acan*, *Col2a1*, *Sox9*, *Comp*, and *Fmod*. The cell classifications here in the master clustering and sub-clustering schemes for IVD cells from CON populations are consistent with a subset of those identified previously for NP cells; differences in cell identities were not noted between CON populations at 2 and 8 weeks following surgery, only differences in cell numbers mapping to each cluster.

Initial results that showed over 80% of all cells were identified as IVD cells at both 2 and 8 weeks, similar to a prior study showing that up to 99% of all sequenced cells mapped to a common transcriptomic profile [5]. The IVD cluster was then sub-clustered into NP, IAF, and OAF cells. For the markers chosen here for cell identification, the analysis revealed the absence of any singular genes that are specific to each subcluster. Rather, IVD cell type identification requires a set of genes that are either more or less highly expressed when compared to their spatial counterparts. Classical markers of IVD cells such as *Acan*, *Col1a1*, and *CoL2a1* and some NP-specific markers such as *Cd24* were useful to support the sub-clustering scheme; however, expression levels exhibited a continuous distribution from inner to outer IVD. Prior studies have sometimes attempted to anatomically isolate NP tissues from adjacent inner AF prior to performing sc-RNA-seq analyses, yet yield the same outcome [5,20,22,45], indicating the spatial overlapping nature of these tissues and the transition between cell populations in the IVD exists in a continuum. Instead, we found additional gene markers including *Igfbp5*, *Igfbp6*, *Lum*, *Myoc*, and *Fibin* useful as supplements of the "classical" markers used to separate outer AF from NP populations within the IVD [5].

In the current study, effort was taken to evaluate sample-specific variability in the clustering results for the CON populations at each time point (e.g., Figure 2c). Studies that pool multiple tissue levels across animals to procure sufficient cells for a single sample for RNA-sequencing may be susceptible to one sample disproportionately driving the clustering outcome. While we had differences in cell yields between rats, the clustering scheme was valued for its ability to produce repeatable results across animals (e.g., 82–92% of all cells mapping to IVD cluster. Approaches that define correlation analyses for clustering schemes as used by Calio and co-workers [16] would be a useful addition to all sc-RNA-seq analyses that employ multiple tissues from multiple subjects.

With LDP, differences were observed in the numbers of cells associated with each identified cluster. We observed that the LDP samples contained 1.5 times more myeloid cells than CON tissue and that macrophages were the most prominent immune subtype present in the myeloid cluster (Figure 3c; Figure 5b). Macrophages have been shown to infiltrate the herniated or injured disc in multiple species [23,42,46–50]. Accordingly, *Cd4* is a pan-macrophage marker in rats [51], *Cd68* is involved in antigen presentation [52,53], *Cd14* is essential for macrophage phagocytosis [54], and *Lyz2* is important for bacteriolysis [55]. In the LDP group, there were increased proportions of *Cd4*+ (*p* < 0.00001), *Cd14*+ (*p* < 0.05), *Cd68*+ (*p* < 0.00001) and *Lyz2*+ (*p* = 0.05) cells. Analysis of gene expression levels for each marker showed few differences in expression levels for these macrophage markers between CON and LDP groups, despite the higher numbers of cells expressing each marker in the LDP samples compared to CON (Figure 5b).

Macrophages can also undergo polarization into the classically activated pro-inflammatory phenotype, M1, or the alternative anti-inflammatory phenotype, M2 [39,56]. They have a remarkable range of cellular functions that can be elicited during tissue injury and repair dependent upon macrophage polarization within the affected tissue [57]. M1-type

macrophages are classically activated and differentiation is stimulated by IFNγ to secrete pro-inflammatory cytokines such as IL6, IL12, and TNFα. M2-type macrophages are alternatively activated macrophages, and their differentiation is stimulation by IL4 produced by Th2 T lymphocytes. Their roles are to secrete anti-inflammatory cytokines and factors such as Arginase 1, IL10, and aid in wound healing, tissue repair, and regeneration. The myeloid cluster was highly enriched for both M1 macrophages and M2 macrophages at 2 weeks post IVD injury (Figure 5c,d). While both M1 and M2 macrophage markers are present, there is a significantly higher expression for *Cd86* (M1) and significantly lower expression for *Arg1* (M2) genes, suggesting a shift towards M1, or pro-inflammatory state with LDP at 2 weeks post-injury. This shift away from the M2 phenotype is consistent with the observations in the human degenerate nucleus pulposus cells that exhibited declining M2 proportion with increasing degeneration severity. Single-cell transcriptomic analyses of the human IVD cells show also the increasing activation of immune recruitment pathways during degeneration, demonstrating the robust interactions between the immune system and the degenerating IVD.

Our study here also identified the increased activation of B- and T-lymphocytes two weeks into the degenerative process. The majority of the top 25 commonly upregulated genes were in the myeloid and lymphoid clusters, and gene analysis led to the intriguing discovery that over 80% percent of those genes were highly enriched in lymphoid cell lineages. Notably, B cell-associate gene markers were the most represented group (40%: *Pika3p1*, *Rdac2*, *Blnk*, *Syk*, *Ikzf3*, *Mzb1*, *Cd19*, *Mybl1*) while the T cell-associated gene markers were slightly less prominent (25%: *Cyfip2*, Lef1, *Mybbp1a*, *Prag1*, *Myb*). Several gene markers were expressed in both T- and B- cells (35%: *Ptprcap*, *Sash3*, *Rag1*, *Top2a*, *Mybl2*, *Cd53*, *Laptm5*) (Figure 4d). B-lymphocyte specific markers included *Blnk*, *Rac2*, and *Pik3ap1*, which are involved in receptor signaling; *Syk*, *Cd72*, *Ikzf*, and *Cd19*, which are involved in differentiation and proliferation; and *Mzb1*, *Mybl1*, which are involved in antibody production and secretion. Due to the large presence of lymphocyte genes, we checked for the presence of B and T lymphocyte specific markers within the lymphoid clusters to confirm these immune cell types were present. B cell co receptor and receptor, *Cd19* and *Cd79b*, were present along with *Blnk*, which plays a role in B cell maturation; *Cd38*, which is expressed on B cell progenitors and plasma cells [58,59]; and *Cxcr4*, which plays a role in B cell activation and trafficking (Figure 5e) [60]. T lymphocytes cell receptor, *Cd3g*, had increased gene expression levels and an increase in cell number. Interestingly, only a small number of cells expressed *Cd4*, helper T cells, or *Cd8b*, cytotoxic T cells, emerged as identifiable in the LDP samples (Figure 5f).

Macrophages are capable of expressing cytokines that activate T-cells so that we examined T-cell specific markers for the incidence and genotype of T-cells. Notable T cell-specific markers were *Cyfip2*, needed for T cell adhesion; *Lef1*, required for IL17A expressing T cell maturation; and *Myb*/*Mybbp1a*, which are involved in T cell differentiation. Five genes were highly expressed in both B and T lymphocytes and regulate a range of functions such as cell adhesion and migration, *Cd53*; antigen receptor signaling, *Ptprcap* and *Sash3*; and VDJ recombination, *Rag1*, and *Top2a*. Though T-lymphocytes have been found in degenerating human IVDs [5], the presence of B-lymphocytes, which accounts for the majority of the highly expressed genes in the LDP population (e.g., *Cd72*, *Blnk*, *Syk*, and *Rag1*) has yet to be described during degeneration (Figure 5e,f). The involvement of immune cells and hematopoietic cells in the LDP tissue suggests that cell infiltration and differentiation were promoted by the surgical stab injury [14]. The role of lymphocytes after tissue injury had not been studied extensively; recent studies in bone fracture healing show clear roles for lymphocytes in mediating the tissue response after injury [61]. In bone fracture healing, lymphocytes were observed to infiltrate the injury site in a temporal fashion where they initially infiltrate the tissue early during the inflammation process, by 3 days post-injury, then leave and return after the peak inflammation window has passed and tissue proliferation has begun, around 14 days post-injury [61]. Lymphocytes function in the inflammation process by recruiting antigen-presenting cells, secreting cytokines to help mediate the inflammatory

response, and secreting growth factors to aid in regeneration [61,62]. We also observed that lymphocytes persisted up to 8 weeks post injury (Figure 8a,b; Table 1). B-lymphocytes act as antigen-presenting cells and dampen the inflammatory response by promoting Treg survival and secreting IL10 [62]. The lymphoid cluster contained *Cd79b*/*Cd38* positive mature B-lymphocytes [59], and the presence of *Cd79b*/*Cxcr4* positive B-lymphocytes indicates activation and trafficking (Figure 5e,f) [60]. T-lymphocytes, in contrast, modulate antigen presenting macrophages, and regulate the inflammatory response with specific subtypes such as *Cd4* positive Tregs [63]. *Cd4* positive T-lymphocytes are associated with tissue repair after injury due to eliciting immunomodulatory effects, while cytotoxic *Cd8* positive T cells have been shown to have a negative effect on repair [61]. We observe *Cd3g* broadly expressed in T-lymphocytes, *Cd3g/Cd4* expression in T-lymphocytes helper cells, and *Cd3g/Cd8* expression in cytotoxic T-lymphocytes (Figure 5e,f). Our data here reveal a novel role for B-lymphocytes in concert with T-lymphocytes as part of the immune cascade during degeneration that persists in both our model here and in humans. Moreover, the immune markers, such as *Myb*, *Blnk*, and *Cd53*, remain similar over time for CON and LDP, suggesting that this is a sustained, progressive response due to degeneration rather than the acute trauma of the injury.

In addition to the described differences in cell number and gene expression between cells within the clusters of control and degenerative IVDs, we observed changes in *Ngf* levels and expression patterns for IVD cell sub-clusters with IVD degeneration. Both *Ngf* and *Ngfr* (encoding for p75NTR) were increased in expression levels in OAF cells, and the number of OAF cells expressing *Ngfr* was also increased. Additionally, NGF is proangiogenic that upregulates VEGF, and indirectly recruits endothelial cells [64]. Sustained NGF sensitizes innervating primary afferent neurons to produce hyperalgesia [65] via upregulated NGF receptor expression in ganglion neurons [66], and promotes the expression of substance P, calcitonin gene-related peptide, and neuronal ion channels that may lead to long-term adaptive effects on nociceptors [67,68]. Moreover, NGF binding to its receptor can regulate collateral sprouting of sympathetic fibers to dorsal root ganglion neurons that may relate to the maintenance of sensitization or chronic pain. Indeed, a prior study of this LDP model in the rat revealed elevated protein expression for the NGF receptors TrkA and p75NTR in the innervating dorsal root ganglion neurons [44]. Secreted NGF can also drive the growth and differentiation of both B- and T-lymphocytes [69,70] in a pattern consistent with the elevated expression of these cell populations in the LDP samples.

A limitation of the current study is the absence of PCR confirmation of identified targets as performed in some prior studies [5,15,22]. As our interest was in numerically minor cell types involved in IVD degeneration, rather than phenotyping cells derived from discrete regions of the IVD, both immunolabeling and bulk-PCR proved difficult with remaining sample-specific tissues. Future work would be needed to design a surgical approach that made use of more lumbar levels, or other strategies such as flow cytometry, to support the collection of tissues and cells so that bulk-PCR or immunolabeling could be used to confirm the current findings. As is, the clear corroboration of clustering schemes that supported the use of "classical" cell-specific markers in the rat as used by Wang and co-workers [15], and the consistency observed in clustering results across four distinct rat subjects provides some measure of confidence and consistency to support the results provided here.

In summary, this study of single-cell RNA-seq in a rat model of IVD degeneration demonstrated changes in the numbers of cells associated with immune cell populations, with increased infiltration of B-lymphocyte and macrophage subsets into the IVD, particularly in early degeneration. While the surgical stab of the IVD undoubtedly provoked the release of soluble factors [71] that support endothelial cell infiltration, neurite invasion, and the activation of infiltrating monocytes that lead to neuronal sensitization and the identification of B- and T-lymphocytes as shown here, the results also point to a regionspecific increase in the number of cells expressing *Ngf* and *Ngfr* (p75NTR) and associated transcriptional expression that may be one of the key findings that support the development

of progressive IVD degeneration following stab injury to the IVD in this model and to annulus injury in the human. Most notably, there was a large presence of B-lymphocytes, which accounts for the majority of the highly expressed genes in the LDP population. These findings provide the basis for future work to investigate immune cell cross-talk and sub-type specificity during IVD degeneration.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/10.3390/app12168244/s1, Figure S1: Expression level of marker genes within intervertebral disc cells, myeloid, and lymphoid cells from CON samples at 8 weeks post-surgery. Figure S2: Expression level of marker genes within intervertebral disc cells, endothelial cells, myeloid, and lymphoid cells from LDP samples at 8 weeks post-surgery.

**Author Contributions:** M.R., J.E.S., S.W.C., G.W.D.E., D.S.P. and F.L. contributed to study design, and analyzed data. D.S.P., M.N.B., J.J.S. and L.J. performed the experiment and gathered the data. L.A.S., J.J.S. and S.Y.T. designed the research study and contributed to data interpretation. M.R., S.W.C., D.S.P., J.J.S., S.Y.T. and L.A.S. wrote the manuscript and contributed to data interpretation. All authors have read and agreed to the published version of the manuscript.

**Funding:** Funding was received from the National Institutes of Health (R01AR070975, R01AR078776, R01074441, T32EB028092), The Orthopedic Research and Education Foundation (resident research grant P20-03413), the National Science Foundation (Graduate Research Fellowship No. DGE-1745038), and the Rita Levi-Montalcini Postdoctoral Fellowship in Regenerative Medicine. We also acknowledge the Genome Technology Access Center (GTAC) at Washington University School of Medicine for help with genomic analysis (NIH P30 CA91842, NIH UL1TR002345). Figure 1 was created with Biorender.com (accessed on 1 June 2022).

**Institutional Review Board Statement:** The animal study protocol was approved by the Washington University Institutional Animal Care and Use Committee under Animal Welfare Assurance # A-3381- 01 (Protocol # 20-0059).

**Data Availability Statement:** The RNA sequencing dataset generated during this study was deposited in the NCBI Gene Expression Omnibus repository (accession number: GSE211407).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Influence of Angiopoietin Treatment with Hypoxia and Normoxia on Human Intervertebral Disc Progenitor Cell's Proliferation, Metabolic Activity, and Phenotype**

**Muriel C. Bischof 1, Sonja Häckel 2, Andrea Oberli 1, Andreas S. Croft 1, Katharina A. C. Oswald 2, Christoph E. Albers 2, Benjamin Gantenbein 1,2,\* and Julien Guerrero <sup>1</sup>**


**Abstract:** Increasing evidence implicates intervertebral disc (IVD) degeneration as a major contributor to low back pain. In addition to a series of pathogenic processes, degenerated IVDs become vascularized in contrast to healthy IVDs. In this context, angiopoietin (Ang) plays a crucial role and is involved in cytokine recruitment, and anabolic and catabolic reactions within the extracellular matrix (ECM). Over the last decade, a progenitor cell population has been described in the nucleus pulposus (NP) of the IVD to be positive for the Tie2 marker (also known as Ang-1 receptor). In this study, we investigated the influence of Ang-1 and Ang-2 on human NP cell (Tie2+, Tie2<sup>−</sup> or mixed) populations isolated from trauma patients during 7 days in normoxia (21% O2) or hypoxia (≤5% O2). At the end of the process, the proliferation and metabolic activity of the NP cells were analyzed. Additionally, the relative gene expression of NP-related markers was evaluated. NP cells showed a higher proliferation depending on the Ang treatment. Moreover, the study revealed higher NP cell metabolism when cultured in hypoxia. Additionally, the relative gene expression followed, with an increase linked to the oxygen level and Ang concentration. Our study comparing different NP cell populations may be the start of new approaches for the treatment of IVD degeneration.

**Keywords:** intervertebral disc (IVD); nucleus pulposus progenitor cells (NPPCs); angiopoietin-1 receptor (also known as Tie2); fluorescence-activated cell sorting; hypoxia; normoxia

#### **1. Introduction**

#### *1.1. Low Back Pain and Intervertebral Disc Degeneration (IVDD)*

Intervertebral disc degeneration is a significant and multifaceted cause of low back pain [1]. The so-called "discogenic pain" is described as a major cause of low back pain in 26–42% of cases [2]. Yet, it is known that the IVD undergoes degenerative changes earlier in life than do other tissues of major organ systems [1]. Whilst for the process of aging and degeneration, the IVD cells enter a state of cell senescence, showing an altered and often more catabolic metabolism [3], oftentimes the nucleus pulposus (NP) is primarily affected [3]. Moreover, the IVD possesses a very limited self-healing potential [4,5]. Conventional treatment strategies, such as analgetics, physical therapy, or surgical interventions, mainly address the symptoms of LBP, while no regeneration mechanisms in the IVD are activated [6]. The combination of tissue engineering strategies with stem cell-based therapies and approaches with native cells have gained significant momentum, being a desirable and promising novel treatment opportunity for degenerated IVDs since they aim to restore the properties of the IVD and, thus, their biological function [7].

**Citation:** Bischof, M.C.; Häckel, S.; Oberli, A.; Croft, A.S.; Oswald, K.A.C.; Albers, C.E.; Gantenbein, B.; Guerrero, J. Influence of Angiopoietin Treatment with Hypoxia and Normoxia on Human Intervertebral Disc Progenitor Cell's Proliferation, Metabolic Activity, and Phenotype. *Appl. Sci.* **2021**, *11*, 7144. https:// doi.org/10.3390/app11157144

Academic Editor: Ming Pei

Received: 12 July 2021 Accepted: 30 July 2021 Published: 2 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### *1.2. Hypoxia vs. Normoxia in the Intervertebral Disc (IVD)*

Another critical aspect in the IVD cell biology is that the NP is avascular, so its cell population does not receive a blood supply (nutrients diffuse through the vertebral endplates); therefore, oxygen levels in this central area are very low [8,9], causing a hypoxic environment [10]. As consequence, these NP cells may have adapted to these conditions. An increase in oxygen concentration during IVDD could be considered a pathological condition, unlike that observed in other tissues or wound healing [11].

#### *1.3. Nucleus Pulposus Progenitor Cells (Also Known as Tie2+ Cells)*

A cell source of current major interest is the tissue-specific Tie2 (also known as TEK receptor tyrosine kinase or angiopoietin-1 receptor or CD202b) positive NP progenitor cell (NPPC) population. These cells are perfectly adapted to the hypoxic environment and can therefore uphold their regenerative potency [12]. Notably, NPPCs are attributed with having potential in regenerative medicine because of their remarkable in vitro multipotency capability through their osteogenic, chondrogenic, and adipogenic lineages in addition to their self-renewal potential [13]. The harsh IVD microenvironment forms an obstacle to exogenous cells (such as bone marrow-derived mesenchymal stem cells) to survive and exhibit a tissue-regenerative function. Thus, NPPCs might pose a promising target to stimulate and induce tissue repair [14]. However, a major challenge in the NPPC therapy approach is caused by the decrease in the percentage of the Tie2<sup>+</sup> NPPC population with aging and the extent of disc degeneration [1]. The methodology for optimized NP cell isolation (by cell sorting) and cell expansion are already described [6,15,16].

#### *1.4. Tie2 Receptor and Its Ligands; Angiopoietin-1 and Angiopoietin-2*

The receptor tyrosine kinase Tie2 is mainly expressed in endothelial cells and is essential for embryonal and adult angiogenesis and vascular maturation and maintenance [17,18] and its also expressed by early hematopoietic cells and subsets of monocytes express Tie2 [19]. Ang-1–4 are the growth factor ligands of the Tie2 receptor. In this project, we focused on Ang-1 and -2, which are mainly produced by vascular smooth muscle cells, throughout the entire vascular system, and endothelial cells respectively [20–22]. The binding of Ang-1 to its extracellular receptor Tie2 results in autophosphorylation and subsequent activation of downstream signaling pathways including the PI3-kinase/Akt [18,19,23]. The Tie2 receptor activation mediates stabilization, remodeling, and repair of the vascular system [18,24]. The inactivation of this pathway protects endothelial cells from apoptosis and inhibits inflammatory responses [18]. Cell-based Ang-1 gene therapy may represent a new strategy for treating various endothelium disorders, including several severe human pulmonary diseases [25,26]. Ang-2 fulfills a dual function by acting as an agonist or antagonist in a context-dependent manner [17]. while showing similar biological outcomes to Ang-1. However, it has to be emphasized that they are weaker in potency and kinetics. As a consequence, Ang-2, thus, may act as a partial Ang-1 agonist [17].

#### *1.5. The Role of Tie2/Ang-1 and Ang-2 Signaling Pathway in Human NPCs and NPPCs*

Sakai et al. [1] suggested that Tie2 is a sensitive marker of aging and the degeneration of IVDs, and therefore is strongly relevant in IVDD. In this study, the elementary role of the Tie2/Ang-1 signaling was already demonstrated. Endogenous Ang-1 provides an innate anti-apoptotic effect on human NP cells and is crucial for their survival. Additionally, enhancement of the proliferation of Tie2+ populations was achieved when NP cells were cultured with Ang-1. More recently, it was discovered that with aging and increasing degeneration, the imbalance between the Ang-2/Ang-1 ratio significantly rises [23]. Downregulation by Ang can be associated with suppressing cell adhesion and viability and promoting the apoptosis of NPCs; therefore, the blockage of Ang-2 may represent another novel therapeutic approach to prevent NPC loss in IVD.

#### *1.6. Hypothesis and Aims*

This study aimed to investigate the effect of angiopoietin-1 and angiopoietin-2 in combination with culture under hypoxia and normoxia on three NP cell populations' (mixed, Tie2<sup>−</sup> and Tie2+) proliferation, metabolism, and phenotype in the context of IVD regenerative medicine. The present work was undertaken to explore the hypothesis that, as NP tissue represents a heterogenous cell population, cells will react differently, in a dose-dependent effect to Ang-1 and Ang-2 stimulation. Secondly, we hypothesized that the variation in the measured parameters within those three cell populations is modified under normoxic and hypoxic conditions with or without stimulation with Ang-1 and Ang-2.

#### **2. Materials and Methods**

#### *2.1. NPC Isolation*

NP cells were isolated from IVD. The IVDs were obtained from human donors undergoing spinal surgery after trauma (Table 1). Written and informed consent was obtained from all participants. The tissue was dissected into nucleus pulposus tissue (NP), annulus fibrosus (AF) tissue, cartilaginous endplate (CEP) tissue, and cut into smaller fragments (0.3 cm3). To avoid drying out the sample fragments, a sufficient amount of phosphatebuffered salt solution (PBS) was given over the tissue fragments.


**Table 1.** Overview of the sex and age of the patient when the sample was collected, IVD tissue location, cell's passage, and IVD status from all patients used in this project. Abbreviations: Th, thoracic; L, lumbar; S, sacrum; P, passage.

> Subsequently, NP tissue samples were centrifuged for 5 min at 500× *g* (Eppendorf Centrifuge 5810, Vaudaux-Eppendorf, Schönenbuch, Switzerland). After centrifugation, PBS was removed and the samples were incubated with sterile 0.19% pronase solution (#10165921001; Roche Diagnostics, Basel, Switzerland) and digested on a shaker for 60 min at 12 rpm.

> After one hour of digestion, tissue digested samples were centrifuged again for 5 min at 500× *g*. After removing the supernatant, the tissues were washed with PBS and digested with sterile collagenase type II (255 U/mg; #LS004174; Worthington Biochemical, Inc., London, UK). A solution of 0.025% collagenase type II was used, diluted in LG-DMEM (low glucose Dulbecco's Modified Eagle Medium (1 g/L; LG-DMEM; #31600-083; Thermo Fisher Scientific, Basel, Switzerland) +10% FBS (fetal bovine serum; #10500-064; Gibco Life Technologies, Basel, Switzerland). Then, the tissue samples were transferred into T75 flasks and left on a shaker at 12 rpm for overnight digestion.

> The following day, cells and the digested tissues were filtered through a 100 μm strainer (#431752; Falcon; Becton, Dickinson and Company, Allschwil, Switzerland) to remove debris and washed once again with PBS. The cells were then resuspended in a LG-DMEM complete medium (CM) containing 0.22% sodium hydrogen carbonate (#31437

500G-R; Sigma-Aldrich, Buchs, Switzerland), 10% FBS, 1 mM sodium pyruvate (#11360-039; Thermo Fisher Scientific), 1% penicillin/streptomycin/glutamine (#10378-016; P/S/G, 100 units/mL, 100 μg/mL and 292 μg/mL, respectively; Thermo Fisher Scientific), and 10 mM HEPES buffer solution (#15630-056; Thermo Fisher Scientific). The NP cell suspension was then transferred into T150 culture flasks (#90151; TPP, Trasadingen, Switzerland). Furthermore, cells from the AF and CEP tissues were frozen.

#### *2.2. NPC Expansion*

The cells were expanded with LG-DMEM (CM) supplemented with 2.5 ng/mL of FGF-2 (#100-18B, PeproTech, London, UK) [13]. The medium was changed twice per week. To obtain a sufficient number of cells for future experiments, NP cells were passaged until 100% confluency was reached, which corresponded to about 5–8 days of culture.

#### *2.3. NPC Sorting*

As soon as the NP cells reached 90% confluency, most of them were used to be sorted into a Tie2<sup>+</sup> cell population and a Tie2<sup>−</sup> cell population with fluorescence-activated cell sorting (FACS) as previously described [6,15,16]. In short, freshly isolated NPCs were incubated with an anti-human Tie2 PE-conjugated monoclonal mouse antibody (#FAB3131P, clone 83715, R&D systems) with 2 μL/106 cells for 30 min on ice and protected from light in 100 μL of FACS buffer ([PBS] with 0.1% bovine serum albumin (BSA; #A7030-100G; Sigma-Aldrich), P/S (100 units/mL and 100 μg/mL; Thermo Fisher Scientific), and 0.5 M ethylenediaminetetraacetic acid (EDTA; #03610; Fluka, Buchs, Switzerland)). Propidium iodide (Sigma-Aldrich) was used to exclude dead cells. Cell sorting was performed by FACS Diva III (BD Biosciences) based on the procedure previously reported [6]. The mouse IgG1 PE-conjugated antibody (#IC002P, clone 11711, R&D systems) was used as an isotype control to set the gate for sorting. The NP cell population prior to the sorting process was considered a mixed NP cell population, and the NP cell population after the sorting process was considered Tie2+ or Tie2<sup>−</sup> NP cell populations.

#### *2.4. NPC Seeding and Incubation with Ang-1 or Ang-2*

In our study, 104 cells/cm2 NPCs were seeded on 12-well plates (#92012, TPP) and incubated for 7 days, either within a control medium (CM; LG-MEM) or supplemented medium (LG-MEM with either angiopoietin-1 or angiopoietin-2 (#130-06 and #130-07, respectively, Peprotech) at a concentration of 10 ng/mL, 50 ng/mL, and 100 ng/mL). For each donor, the plates were prepared as duplicates, such that one plate per condition could be tested under a normoxic condition (21% of O2) and the second plate was tested under 2% O2 hypoxia (Invivo2 400 #0612NCF05, Ruskinn Technology Ltd., London, UK).

#### *2.5. Resazurin Sodium Salt Cell Activity Assay*

The resazurin sodium salt (#R7017-1G, Sigma-Aldrich) assay was performed at a concentration of 10 μg/mL added to the respective medium for 4 h at 37 ◦C under normoxic or hypoxic conditions, respectively. All test conditions were run in triplicate (#655101\_100, Greiner Bio-One, St. Gallen, Switzerland). The fluorescence signal was quantified with a SpectraMax M5 (Molecular Devices, distributed by Bucher Biotec, Switzerland) microplate reader with excitation and emission wavelengths of 560 nm and 600 nm, respectively. The resulting values were normalized to the respective day 7 control sample CM without any added growth factors and incubated in normoxia or hypoxia, respectively.

#### *2.6. Papain Digestion and DNA Quantification*

The quantity of DNA present in each condition was measured for each sample after overnight papain digestion. The absolute amount of DNA was determined, using the CyQUANTTM Cell Proliferation Assay Kit (#7026 Molecular Probes, Fisher Scientific) according to the manufacturer's instruction. All test groups were analyzed in triplicate.

#### *2.7. RNA Extraction and cDNA Synthesis*

The total RNA was isolated from monolayer cultures using the GenEluteTM Mammalian Total RNA Miniprep Kit (#RTN70-1KT; Sigma-Aldrich), according to the manufacturer's instruction. The following four steps were performed: releasing of the RNA from cells, binding of the RNA to the affinity column, washing to remove contaminants, and elution of the purified RNA. The purity of the RNA was assessed by measuring the A260/A280 ratio.

Reverse transcription of total RNA into double-stranded cDNA was performed with the High Capacity Reverse Transcription cDNA Kit (Thermo Fisher Scientific, Inc., cat #4368814), including 10X RT Buffer, 10X Random Primers, 25X dNTP Mix, and MultiScribe Reverse Transcriptase (50 U/mL), which were combined into a mastermix with additional RNase-free water and the addition of about 500 ng of total RNA.

#### *2.8. Quantitative Polymerase Chain Reaction (qPCR)*

To study the relative gene expression, cDNA (2 μL) was mixed with the PCR reaction solution (2X SYBR green master mix, #B21202; Bimake) containing 0.25 μM specific primers (Table 2). A real-time quantitative polymerase chain reaction (qPCR) was performed using the CFX96TM Real-Time System (#185-5096 C1000 Touch Thermal Cycler; Bio-Rad Laboratories) under standard thermal conditions. The qPCR program consisted of a twostep cycling protocol: an initial denaturation at 95 ◦C for 5 min, followed by 39 PCR cycles at 95 ◦C for 10 sec, and an annealing temperature of 61 ◦C for 30 sec. To control the specificity of the amplification products, a melting curve analysis was carried out for each reaction. The results were expressed relative to gene expression levels on day 0. Gene expression was calculated by the 2−ΔΔCt method [27], and the results were normalized to the reference genes glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and 18S ribosomal RNA (18S) expression. All reactions were performed in triplicate.


**Table 2.** Overview of human genes used for qPCR using glyceraldehyde-3-phosphate dehydrogenase (*GAPDH*) and *18S* ribosomal RNA (*18S*) as a reference gene. Genes analyzed for the NPCs: receptor tyrosine kinase (*TEK*), collagen type II (*COL2*), hypoxia-inducible factor-1α (*HIF1A*), aggrecan (*ACAN*). (F = forward,R=reverse).

> We first assessed the effect of Ang-1/2 at different concentrations in hypoxia and normoxia on the relative expression of two genes that are associated with ECM production in the IVD, namely aggrecan (ACAN). We also examined the HIF1A gene, which is the crucial mediator of the adaptive response of cells to hypoxia. Furthermore, to estimate the expression of the Tie2 receptor, regarding the effect of the Ang-1/2 proteins and the microenvironments (normoxia vs. hypoxia) in the NP cells mixed population, and the NP Tie2<sup>−</sup> and NP Tie2+ cell populations, we analyzed the TEK gene.

#### *2.9. Statistical Analysis*

Statistics were performed using GraphPad Prism (version 8.0.1 for Windows, Graph-Pad Software; San Diego, CA, USA). All data were initially tested for normal distribution using the Shapiro–Wilk test. Due to lack of normality a nonparametric distribution was then assumed. Statistical significance was determined using the Kruskal–Wallis's test followed by the Sidak's multiple comparison test. Results were presented as mean ± standard deviation (SD). A *p*-value < 0.05 was considered significant. Asterisks within figures denote the degree of statistical relevance observed: \* < 0.05; \*\* < 0.01; \*\*\* < 0.001; \*\*\*\* < 0.0001.

#### **3. Results**

#### *3.1. Effect of Oxygen Tension and Angiopoietin-1/2 on NP Cell's Proliferation*

To explore the role of Ang-1 and Ang-2 in cell proliferation on the three NP cells populations (namely, Tie2+, Tie2<sup>−</sup>and mixed populations) we assessed the amount of DNA under normoxic and hypoxic conditions for the NP cells mixed, NP Tie2− cells, and NP Tie2+ cells populations (Figure 1).

**Figure 1.** Amount of DNA measured of NP cells populations under normoxia and hypoxia assessed by Hoechst DNA assay on Day 7. Ang-1 and Ang-2 were administrated in three concentrations: 10, 50 and 100 ng/mL. The control Day 0 samples were used as reference. (**a**) Amount of DNA measured of NP cells mixed population (N = 7); (**b**) amount of DNA measured of NP cells Tie2<sup>−</sup> population (N = 3); (**c**) amount of DNA measured of NP cells Tie2+ population (N = 2). Bars represent mean ± SD. A *p*-value < 0.05 was considered significant. Asterisks within figures denote the degree of statistical relevance observed: \* < 0.05; \*\* < 0.01; \*\*\*\* < 0.0001.

Concerning the NP cells mixed population (Figure 1a), under normoxic conditions, the control Day 7 and the stimulated samples showed a significant decrease in their amount of DNA in a range of 35 to 47%, compared to the control Day 0 sample. However, all the samples did not show any significant differences in their amount of DNA, compared to the Day 0 control sample in hypoxia. Moreover, over the whole NP cells mixed population, no concentration-dependent or protein-dependent trends were apparent.

In the NPCs Tie2− population (Figure 1b), a slight upregulation in the amount of DNA was observed in every sample, in normoxia and hypoxia, compared to the Day 0 sample. In normoxia, the amount of DNA was increased up to 1.5 folds and in hypoxia, it was increased 1.25- to 1.6-fold. However, no concentration-dependent or O2-dependent effect or trend was detectable for any sample condition, treated or either untreated with Ang-1 or Ang-2.

The entire collection of samples treated in the NP cells Tie2+ population (Figure 1c) showed an increased amount of DNA, compared to the Day 0 sample. However, no concentration-dependent stimulation or protein-dependent effect was observable either in normoxia. In the hypoxic condition, only the sample treated with Ang-1 at (10 ng/mL) showed a significant increase in their amount of DNA compared to the control.

#### *3.2. Effect of Oxygen Tension and Angiopoietin-1/2 on NP Single Cell's Metabolism*

Further, we examined the metabolic activity/DNA ratio at Day 7, which is the most relevant parameter, as it refers to single-cell metabolic activity (Figure 2).

**Figure 2.** Metabolic activity on single-cell level of NP cells populations in normoxia and hypoxia on Day 7. Ang-1 and Ang-2 were administrated in three concentrations: 10, 50 and 100 ng/mL. The control Day 7 samples were used as reference. (**a**) Metabolic activity on single-cell level of NP cells mixed population; (**b**) metabolic activity on single-cell level of NP cells Tie2<sup>−</sup> population; (**c**) metabolic activity on single-cell level of NP cells Tie2+ population. Bars represent mean <sup>±</sup> SD. A *p*-value < 0.05 was considered significant. Asterisks within figures denote the degree of statistical relevance observed: \* < 0.05; \*\* < 0.01; \*\*\* < 0.001; \*\*\*\* < 0.0001.

In the NP mixed population under normoxic conditions, the stimulated samples showed a 2- to 3-fold increase in the metabolism/DNA ratio, compared with the Day 7 control sample. In summary, for normoxia, a stimulation effect was detectable, but not specifically in concentration- or protein-dependency (Figure 2a). Furthermore, in hypoxic conditions, for the NP mixed population, the stimulated samples all presented an increase in their ratio metabolism/DNA, as they were 4 to 8 times higher than in normoxic conditions. Hence, no clear trend in concentration- or protein-dependency was observable. Concerning the Ang-1 stimulated samples, the metabolism/DNA ratio increased from 10 to 50 ng/mL concentrations and then decreased for the 100 ng/mL sample, with a significant increase for Ang-1 at 50 and 100 ng/mL. Regarding the Ang-2 treated samples, the metabolism/DNA ratio decreased from 10 to 50 ng/mL and then increased again for the Ang-2 for the 100 ng/mL treated sample, with a significant increase for Ang-2 at 10 and 100 ng/mL.

For the NP Tie2− cells population under normoxic condition, the ratios of the Ang-2 at 10 ng/mL and the Ang-1 at 50 ng/mL samples were 3- to 3.5-fold increased, compared to the Day 7 control sample (Figure 3b). For the other Ang-1 or Ang-2 treated samples, the ratios stayed at an equal level to the Day 7 control sample, or they even decreased. In hypoxia, all the stimulated samples showed a decrease in their metabolism/DNA ratio, except for the Ang-2 at the 100 ng/mL sample, whose ratio stayed at a similar level to the Day 7 control sample. as Additionally, for the mixed population, no trend in the effect of the concentration, nor the protein, was assessable. The Ang-1 treated samples presented an increase in their metabolism from the 10 to 50 ng/mL samples and then a decrease to the 100 ng/mL sample. The Ang-2 treated samples showed a decrease from the 10 ng/mL sample to the 50 ng/mL sample and then an increase to the 100 ng/mL sample. However, no significant increase or decrease was detected.

In the NP Tie2<sup>+</sup> cells population, we observed a similar ratio for all the samples in normoxia, either stimulated with Ang-1/2 or not (Figure 2c). The same picture was present in hypoxia. All samples showed similar or slightly elevated ratios, compared to the day 7 control sample, except for the Ang-2 at the 10 ng/mL sample, which significantly decreased in its ratio, compared to the control sample. Therefore, a significant increase in stimulation with Ang-1 at 50 ng/mL and Ang-2 at 10, 50 and 100 ng/mL was detected.

#### *3.3. Effect of Oxygen Tension and Angiopoietin-1/2 on ECM Related-Genes* 3.3.1. Analysis of Aggrecan Relative Gene Expression

Concerning the relative gene expression of *ACAN* for the NP cells mixed population, we did not observe any significant difference for normoxia, compared to the control condition (Figure 3a). No effects were detected in a dose-dependent manner for angiopoietin on the relative gene expression nor for specific variation in the use of Ang-1 or Ang-2. Nonetheless, all treated samples, either with Ang-1 or Ang-2, showed an equal amount of *ACAN* as the control samples, or at least a 2- to 3-fold increase in their ACAN relative expression. However, Ang-2 samples at 100 ng/mL in normoxia presented a downregulation in their *ACAN* gene expression.

**Figure 3.** Relative gene expression of the three NP cells populations, comparing the effect of Ang1/2 at different concentrations in normoxia and hypoxia on the ACAN gene expression. The control (Day 0) samples were used as reference (set at 1). Ang-1 and Ang-2 were administrated in three concentrations: 10, 50 and 100 ng/mL. (**a**) ACAN relative gene expression of NP cells mixed population; (**b**) ACAN relative gene expression of NP cells Tie2− population; (**c**) ACAN relative gene expression of NP cells Tie2+ population. Bars represent mean <sup>±</sup> SD.

Like the NP cells mixed population, the control Day 0 and control Day 7 samples showed non-significant differences in their relative expression of *ACAN* in normoxia and hypoxia for the NP cells Tie2− population (Figure 3b). However, a remarkable effect in a dose-dependent manner was detectable for the Ang-1/2 sample at 10 ng/mL in normoxia, the Ang-1 sample at 100 ng/mL in normoxia and hypoxia, and the Ang-2 sample at 100 ng/mL in hypoxia.

The NP cells mixed, NP Tie2<sup>−</sup> cells, and NP Tie2<sup>+</sup> cells populations showed no differences in the gene expression of *ACAN,* compared to the control Day 0 (Figure 3c). We detected no specific variation in the use of Ang-1 or Ang-2 in normoxia and hypoxia but detected a dose-dependent increase in relative gene expression for hypoxia. In normoxia, we detected no effect in a dose-dependent manner of angiopoietin on the relative gene expression or specific variation in the use of Ang-1 or Ang-2. All the stimulated samples in normoxia were downregulated in their *ACAN* gene expression by 40 to 75%, compared to the control Day 0 sample.

#### 3.3.2. Analysis of Collagen Type II Relative Gene Expression

The relative gene expression of *COL2* for every sample in normoxia and hypoxia increased compared to the control Day 0 sample (Figure 4a). The peak values in *COL2* relative gene expression compared to the control Day 0 sample were reached on a significant level in the Ang-1 and Ang-2 at 10 ng/mL treated samples in normoxic condition. All stimulated hypoxia samples and the control Day 7 sample showed a 2- to 2.5-fold increase in their *COL2* relative gene expression. Thus, we detected no effect in a dose-dependent manner of angiopoietin on the relative gene expression or specific variation in the use of Ang-1 or Ang-2 in hypoxia.

**Figure 4.** Relative gene expression of the three NP cells populations, comparing the effect of Ang1/2 at different concentrations in normoxia and hypoxia on the COL2 gene expression. The control (Day 0) samples were used as reference (set at 1). Ang-1 and Ang-2 were administrated in three concentrations: 10, 50 and 100 ng/mL. (**a**) COL2 relative gene expression of NP cells mixed population; (**b**) COL2 relative gene expression of NP cells Tie2− population; (**c**) COL2 relative gene expression of NP cells Tie2+ population. Bars represent mean <sup>±</sup> SD. A *<sup>p</sup>*-value < 0.05 was considered significant. Asterisks within figures denote the degree of statistical relevance observed: \*\* < 0.01; \*\*\*\* < 0.0001.

The *COL2* relative gene expression for the Ang-1/2 at 10 ng/mL and Ang-1 at 100 ng/mL sample in normoxia followed the same patterns concerning the peak values as described above (Figure 4b). The NP Tie2− population showed a similar relative gene expression, compared to the NP cells mixed population for the control day 7 sample and the Ang-2 at 100 ng/mL sample in normoxia and hypoxia. As showed previously in the NP cells mixed population, the effect of the stimulation of the 10 ng/mL concentration for Ang-1 and Ang-2 in normoxia is noteworthy. Additionally, the Ang-1 at 100 ng/mL sample reacted remarkably to its stimulation in normoxia and hypoxia. For all the other samples, no trend in the Ang-1 and Ang-2 administration was detectable.

Concerning the NP Tie2+ population in normoxia, all the stimulated samples, as well as the control Day 7 sample, showed a 1.5- to 2.5-fold increase in their *COL2* relative gene expression, compared to the normoxia control Day 0 sample (Figure 4c). In hypoxia only, the control Day 7 sample with a 3.5-fold increase and the Ang-1/2 at 10 ng/mL sample with a 2- to 2.25-fold increase in their *COL2* relative gene expression showed a difference, compared to their control Day 0 sample. Interestingly, the 10 ng/mL concentration for Ang-1 and Ang-2 did not provoke a remarkable difference in the gene expression for normoxia as it did in the NP cells mixed population and the NP cells Tie2− population, compared to the normoxia control Day 0 sample. Summing up, the unstimulated control Day 7 sample in hypoxia reached the highest *COL2* gene expression in the NP cells Tie2+ population.

#### *3.4. Effect of Oxygen Tension and Angiopoietin-1/2 on Oxygen Level Related-Genes* Analysis of Hypoxia-Inducible Factor 1-Alpha Relative Gene Expression

The *HIF1A* relative gene expression in NP cells mixed population presented a similar trend, compared to the *COL2* relative gene expression in the NP cells mixed population. As in the *COL2* relative gene expression of NP cells mixed and NP Tie2− cells populations, the 10 ng/mL stimulation in normoxia of the Ang-1 and Ang-2 samples reached, also in this population, the highest amount of relative gene expression (Figure 5a). All stimulated hypoxia samples and the control Day 7 sample range in a 10% decrease and increase around the baseline in their *COL2* relative gene expression. Thus, we detected no effect in a dose-dependent manner of angiopoietin on the relative gene expression or specific variation in the use of Ang-1 or Ang-2 in hypoxia. However, the stimulation in normoxia clearly induced a higher potency than in hypoxia for 4 out of the 6 samples, whereas, again, no clear trend in concentration- or protein-dependency could be evaluated.

**Figure 5.** Relative gene expression of the three NP cells populations, comparing the effect of Ang1/2 at different concentrations in normoxia and hypoxia on the HIF1α gene expression. The control (Day 0) samples were used as reference (set at 1). Ang-1 and Ang-2 were administrated in three concentrations: 10, 50 and 100 ng/mL. (**a**) HIF1α relative gene expression of NP cells mixed population; (**b**) HIF1α relative gene expression of NP cells Tie2− population; (**c**) HIF1α relative gene expression of NP cells Tie2+ population. Bars represent mean <sup>±</sup> SD. A *<sup>p</sup>*-value < 0.05 was considered significant. Asterisks within figures denote the degree of statistical relevance observed: \* < 0.05.

The normoxia associated samples of the Ang-1 and Ang-2 at 10 ng/mL concentrations reached a 12-fold and 15-fold increase in their *HIF1A* relative gene expression, compared to the control Day 0 sample (Figure 5b). Remarkably, similar to the NP cells mixed population, the NP Tie2− cells population showed a similar pattern in relative gene expression, compared to the *COL2* relative gene expression within the NP Tie2− cells population. It is also worth mentioning that the Ang-2 at 50 ng/mL and the Ang-1 at 100 ng/mL samples achieved the highest effect in stimulation in hypoxia for both *COL2* and *HIF1A* relative gene expression within the NP cells Tie2− population.

In the NP cells Tie2<sup>+</sup> population (Figure 5c), the incubation with Ang-1 and Ang-2 in their different concentrations does not seem to have any effect on the relative *HIF1A* gene expression. Notably, all of the stimulated samples showed downregulation in their gene expression, compared to the control conditions at Day 0.

#### *3.5. Effect of Oxygen Tension and Angiopoietin-1/2 on NP Progenitor Cells Related-Genes* Analysis of Angiopoietin-1 Receptor Relative Gene Expression

In the NP cells mixed population, all stimulated samples in normoxia, as well as in hypoxia, showed an increase in their gene expression (Figure 6a). With a stimulation dose of 10 ng/mL for Ang-1/2, the gene expression in normoxia reached two peak values: the Ang-1 at 10 ng/mL sample presents a statistically relevant 400-fold increase and the Ang-2 at 10 ng/mL sample, a 200-fold upregulation of their *TEK* gene expression, compared to their control Day 0 sample. Thus, in this population, again, the 10 ng/mL dose in normoxia reached the highest values in gene expression. However, in hypoxia, no trend in gene expression upon stimulation by Ang-1 or Ang-2 and their different concentrations were detectable.

**Figure 6.** Relative gene expression of the three NP cells populations, comparing the effect of Ang-1/2 at different concentrations in normoxia and hypoxia on the TEK gene expression. The control (Day 0) samples were used as reference (set at 1). Ang-1 and Ang-2 were administrated in three concentrations: 10, 50 and 100 ng/mL. (**a**) TEK relative gene expression of NP cells mixed population; (**b**) TEK relative gene expression of NP cells Tie2− population; (**c**) TEK relative gene expression of NP cells Tie2+ population. Bars represent mean <sup>±</sup> SD. A *<sup>p</sup>*-value < 0.05 was considered significant. Asterisks within figures denote the degree of statistical relevance observed: \* < 0.05.

Concerning the NP cells Tie2− population (Figure 6b), no clear trends in dosedependency nor protein-dependency are detectable. In normoxia, for all the NP cells Tie2+ samples treated with Ang-1 or Ang-2 and including the Day 7 control sample, the *TEK* relative gene expression was increased up to 2-fold, compared to the control Day 0 equivalent (Figure 6c). Therefore, normoxia upregulated the gene expression slightly. Yet, the stimulation with Ang-1 or Ang-2 in its different doses did not have any effect in normoxia, as the gene expression stays in a similar range for all normoxia samples. In addition, in the hypoxic condition, the stimulation did not have any upregulating effect. This missing reaction to stimulation in hypoxia was further confirmed by the control day 7 sample, which reached the maximum value of *TEK* relative gene expression for this population.

#### **4. Discussion**

Our study showed that the three different cells population studied here and present in the NP tissue act/behave differently in hypoxic and normoxic conditions. The increase in the O2 tension significantly decreases the NP cells mixed population proliferation, while, conversely, we observed an increase in proliferation of the NP cells Tie2<sup>+</sup> population. Concerning the cell's activity, the decrease in the O2 tension induced an increase in the cell's metabolism for the NP cells mixed and Tie2+ populations.

Additionally, the three NP cells population showed variation in their phenotype/genotype, according to the presence of Ang-1 or Ang-2 in combination with normoxic or hypoxic conditions. To our knowledge, no other studies are comparing the influence of the oxygen tension combined on three distinct cells population, well-characterized in the human NP tissue, including NP progenitor cells. Moreover, to the best of our knowledge, as of today, there is no literature published about the different components of the heterogeneous NP

tissue cell's population, including NPPCs isolated from relatively young and healthy IVDs in relation to the activation of the receptor-ligand system Tie2/Ang.

#### *4.1. Ang-1/2 and Its Effect on DNA Measurement in Normoxic and Hypoxic Conditions*

Evidence has arisen in recent years suggesting that a tissue renin–angiotensin system (tRAS) is involved in the progression of various human diseases [28], including IVD degeneration [29]. Numerous reports and studies showed the positive effects of pathologic tRAS pathway inhibition and protective tRAS pathway stimulation on the treatment of cardiovascular, inflammatory, and autoimmune disease and the progression of neuropathic pain. Cell proliferation, senescence, and aging are known to be related to RAS pathways. In this context, the DNA concentration measured in our study in all the NP cells mixed and Tie2+ population samples cultured within normoxic conditions were statistically different from the controls at Day 0. None of the samples in the NP cells Tie2− population showed statistically significant differences. The amount of DNA in the NP cells mixed population decreased, whereas the amount was increased in the NP cells Tie2+ population on Day 7 for the control and the stimulated samples in normoxia, compared to their control Day 0 equivalents. In hypoxia, no effect on the proliferation was detected, except for the NP cells Tie2+ population treated with Ang-1 at 10 ng/mL.

In the NP cells Tie2− population, a slight upregulation in the amount of DNA was seen in every sample, in normoxia and hypoxia, compared to the Day 0 sample. The reason for this discrepancy between the NP cells mixed and NP Tie2− cells populations is not clear. As the NP cells mixed population consists of 95 to 99% of the exact cell types as the Tie2− cells population, the 1 to 5% difference in cell type should not provoke the observed difference [1]. Furthermore, both samples underwent the same schedule of incubation and medium change, which should result in a more similar outcome. However, a possible cell-to-cell interaction or co-culture between the Tie2+ and Tie2<sup>−</sup> cells within the mixed population could influence the proliferation as shown in other tissue and cell types [30,31].

All samples treated in the NP cells Tie2+ population showed an increased amount of DNA, compared to the Day 0 sample and a trend of higher amounts of DNA in normoxia, compared to hypoxia. A possible explanation for the lower amount of DNA detected within the cells population incubated under hypoxic conditions, compared to the normoxic condition, could be the potential differentiation within the chondrogenic lineage of the cells in hypoxia, reducing, at the same time, cell proliferation. As described by Lyu et al. [12], an explanation could be that hypoxia can exert either positive or negative effects on IVD progenitor cells. In particular, hypoxia can inhibit proliferation and induce apoptosis but might also promote the chondrogenic differentiation of IVD progenitor cells. Huang et al. [32] identified tissue-specific intervertebral disc progenitor cells (DPCs) from healthy Rhesus monkeys; in particular, they found that they are sensitive to low oxygen tension and undergo apoptosis under hypoxic conditions due to their inability to induce/stabilize hypoxia-inducible factors (HIF). These results confirm our findings. Our findings and the findings of the previously mentioned studies contradict the publication of Tekari et al. [13], who found a better cell survival of Tie2<sup>+</sup> cells cultured in hypoxia (2%), which were based on young-aged bovine cells.

We identified in our results a trend for an increased amount of DNA in normoxia for the NP Tie2+ population. These are contradictory findings to several studies, in which it was shown that an oxemic shift would be expected to lead to a failure in progenitor cell activation [13,33]. The variation in oxygen tension is possibly mediated by alterations in the vascular supply to the endplate cartilage or even the annulus fibrosus. In turn, the cells no longer reside in their native associated microenvironment and are associated with a reduced cell proliferation capacity [34–36].

Concerning the amount of DNA, we detected either no effect in a dose-dependent manner of angiopoietin or specific variation in the use of Ang-1 or Ang-2 in normoxia and hypoxia for the three tested populations. Sakai et al. [1] demonstrated selective enhancement of the proliferation of Tie2<sup>+</sup> GD2<sup>−</sup> CD24<sup>−</sup> (i.e., Tie2 single-cell sorted positive cells (T/sp)) and of Tie2<sup>+</sup> GD2+ CD24<sup>−</sup> (i.e., Tie2 single-positive cells (TG/dp)) by using a co-culture with a cells feeder layer (murine stromal cells overexpressing human Ang-1 (AHESS-5). A marked Ang-1 dependent growth resulted in increased proliferation of the T/sp population by 3.1 ± 0.4 times and for the TG/dp population by 3.2 ± 0.4 times. Unfortunately, no information about the administrated Ang-1 dose or incubation time was published. Nevertheless, these results could not be confirmed with the NP cells Tie2+ population in this study. A possible explanation could be that the Tie2 receptor was not expressed anymore over in vitro culture time, as cells needed expansion before starting the experiment. Summing up, at the time of incubation with Ang-1 and Ang-2, the cells were supposedly not Tie2+ progenitor cells anymore, and therefore could not react to any stimulation, as the receptor for our ligands was not present anymore.

Ligand-deficient Ang1−/<sup>−</sup> or receptor-deficient Tie2−/<sup>−</sup> is embryonically lethal for transgenic mice. Both deletions showed common features in their phenotype, and both failed to develop a fully functional cardiovascular system [4]. Vasculogenesis proceeds normally, but remodeling and maturation of the vessels are defective [24]. Ang-2-deficient mice seem phenotypically normal at birth but die within 14 days or develop normally throughout adulthood. On the other hand, systemic embryonic Ang-2 overexpression results in an embryonic lethal phenotype of Ang-2 transgenic mice [18].

Furthermore, Wang et al. [37] determined the role of Ang-2 in the pathology of IVDD. They verified the expression of Ang-2 in human degenerative NPC and showed that Ang-2 expression was considerably higher in severe IVDD than in mild IVDD. Additionally, they demonstrated that NP cells treated with exogeneous Ang-2 resulted in marked downregulation of type II collagen and aggrecan and a significant increase in the expression of MMP13, ADAMTS4, and the pro-inflammatory cytokine IL-1β, both on protein and mRNA levels. The resulting dysregulation in ECM degradation in NPCs upon Ang-2 exposure contributed to the pathogenesis of IVDD. Conclusively, the study reveals Ang-2 inhibition as a potential novel therapeutic target for IVDD treatment.

#### *4.2. Ang-1/2 and Its Effect on Metabolic Activity Per DNA in Normoxic & Hypoxic Conditions*

Out of all analyses, in cell metabolism relative to DNA amount, the difference between oxygen tension (≤ 5% vs. 21%) was in favor of hypoxia for every single sample tested in our three populations of interest, i.e., NP cells mixed, NP cells Tie2−, and NP cells Tie2<sup>+</sup> populations. In our study, the samples in the hypoxic environment all achieved a higher metabolism, compared to the normoxic environment. This is in accordance with previous studies upon oxygen levels in NPC cultures, such as those of Feng et al. [38], Mwale et al. [39], Stoyanov et al. [40], and Gantenbein et al. [31]. The authors demonstrated that hypoxia enhances NPCs phenotype, which results in greater ECM components production and, indeed, a higher metabolism.

Additionally, as the NP cells mixed population and NP Tie2− cells population only differ in 1 to 5% of their cell population (only depletion of Tie2<sup>+</sup> cells), this could explain that stimulated samples showed comparable values in their metabolic activity [6]. In contrast, the NP cells Tie2<sup>+</sup> population responded in higher potency to the normoxic and hypoxic environments, compared to the two others population (i.e., mixed and Tie2−). On one hand, the upregulation of the cell's metabolism would therefore seem logical, as this population expresses the Tie2 receptor and thus can react to the stimulation (i.e., presence/absence of Ang) by its ligands Ang-1 and Ang-2. However, on the other hand, the stimulation effect remains in question, as the treated samples do not show any trend either in a dose- or protein-dependent manner.

According to the findings of Wang et al. [23], Ang-2 plays a role in suppressing cell adhesion, decreasing cell viability, and promoting NP cell apoptosis. Concerning the cell's metabolism, in our study, the findings for Ang-2 could not be confirmed, which supports the aforementioned lack of effect upon Ang-1 and Ang-2 stimulation. Thus, the more reliable and convincing explanation for the upregulation of fluorescence intensity (i.e., cell's metabolism) in the NP cells Tie2+ population in normoxia and hypoxia, compared to

the NP cells mixed and NP cells Tie2+ populations, could potentially be that progenitor cells show a higher metabolism, compared to mature cells [41].

Another aspect that has to be taken into consideration by comparing the findings of Wang et al. [23] with our study is that the cells from our study were isolated from trauma discs and not from degenerated IVDs. Thus, the cells in Wang's study were already stimulated by Ang-1, due to potential blood vessel infiltration and inflammation in the degenerated disc as well as by Ang-2, as the Ang-2/Ang-1 ratio increases in degenerated discs for the same reason [37]. Whether such contradictory findings could be related to differences in IVD degeneration pathways between injuries (trauma), natural aging, or the degree of IVD degeneration is still unclear.

#### *4.3. Ang-1/2 and Its Effect on the ECM Gene Expression in Normoxia & Hypoxia*

Looking at the gene expression levels, we did not observe significant results between normoxia and hypoxia for the majority of the samples. The effect of hypoxia on gene expression associated with chondrogenic lineage and ECM production in the IVD is already well documented in the literature [38,39]. At least, in the NP cell Tie2<sup>+</sup> population for ACAN, we detected a trend of upregulation in hypoxia, whereas for *COL2,* it was not the case. If Ang-1 and/or Ang-2 had a stimulatory effect, this would be the case for the NP cells Tie2<sup>+</sup> population, as this population has the Tie2 receptor for the binding of the Angligands. It is, therefore, astonishing and unexpected to find the highest upregulations for ACAN in normoxia in the NP cells mixed and Tie2− populations for Ang-1/2 at 10 ng/mL and Ang-1 at 100 ng/mL. A possible explanation for these particular increases could be related to any cross-reactions of the ligands Ang-1 and Ang-2 with another receptor. This hypothesis is difficult to defend in the sense that there are only two or three out of the six treated samples in normoxia that show an effect upon stimulation.

However, the main explanation for the stimulation of Ang-1 and Ang-2 not having any effect is the outcome for the NP cells Tie2+ population. The stimulation did not show any differences upon Ang-1 or Ang-2 protein stimulation, including their different concentrations. In the case of the NP cells Tie2<sup>+</sup> population, several explanations can be considered. First, there could have been a cross-reaction of the Tie2 receptor with other ligands, which subsequently blocked the binding of Ang-1/2 to the receptor and thus the downstream pathways. Second, it would be possible that the Tie2 receptor was not expressed after the expansion culture phase. Briefly, at the time of incubation with Ang-1 and Ang-2, the cells were not Tie2+ progenitor cells anymore. To confirm the progenitor state of the cell right before the start of the experiment, another FACS analysis could have been performed. Third, the first FACS analysis itself could already have stressed the cells and favored the loss of their progenitor attributes.

The hypothesis for the lack of effect of Ang stimulation is encouraged by the examination of the TEK gene. It is eminently important to see that, similar to *ACAN* and *COL2*, we did miss out on any specific reaction in the relative gene expression of *TEK* to either Ang-1, Ang-2, including their three different doses in the Tie2+ cells population. This is unexpected for the same reasons as already mentioned for *ACAN* and *COL2*. In the NP mixed and NP Tie2− cells populations, the relative gene expression in *TEK* and *HIF1A* do show the same patterns as *COL2* and *ACAN*. The peak values in normoxia deserve special mention, such as Ang-1/2 at 10 ng/mL and Ang 1 at 100 ng/mL. However, again, it was unexpected to find such upregulation in normoxia of the genes within the NP mixed and NP cells Tie2− populations. As we find the same pattern several times, there must be attractive or favorable attributes for these doses in normoxia. Up to that point, explanations remain elusive. As already mentioned before, the data must be interpreted with caution, due to the high donor variation. The downregulation of *HIF1A* in hypoxia in the NP cells mixed, NP cells Tie2<sup>−</sup> and Tie2<sup>+</sup> populations can be explained, as hypoxia progressively decreases *HIF1A* expression at the mRNA level, due to the mRNA-destabilizing protein tristetraprolin [42].

However, Huang et al. [32] detected that the level of relative gene expression of *HIF1A* for DPCs of the NP remained constant under hypoxia, compared to normoxia. They suggest that incompetency of the DPCs of the NP to adapt to oxygen tension may be in part due to the inability to upregulate HIF expression. Their findings are consistent with our results of the *HIF1A* relative gene expression of the NP cells Tie2+ population.

#### *4.4. Study Weaknesses and Limitations*

One of the main limitations of this study was the relatively small sample size for the NP cells Tie2<sup>+</sup> population. Sample expansion of the NP progenitor cells after a sorting process was a challenge, due to the relatively low percentage who represent this population in the NP tissue, and future studies should recruit IVD material from bigger cohorts and multiple populations and etiologies. Moreover, previous studies demonstrated the difficulties to expand, in vitro, this population [16,43]. Additionally, the "controls" used in the study were not truly healthy IVDs, as they were traumatic discs. In this respect, a comparison of fresh cadaveric IVD material of undegenerated nature would be essential for any future investigations. Moreover, the 2D culture was shown to not be the most accurate system for culture and expansion of NP progenitor cells [16]. Previous studies, such as that of Sakai et al., dealt with only Tie2<sup>+</sup> cells from degenerated discs; Wang et al. investigated the effect of Ang-2 in degenerated discs. It remains unclear whether the Tie2<sup>+</sup> cells from the degenerated or the trauma discs react in a comparable way upon Ang stimulation. Thus, it remains unclear whether we can compare our results (cells isolated from trauma patients) to the results from publications based on degenerated discs. We suspect that the results from degenerated discs are potentially biased, as discs subjected to degeneration present a higher Ang-2/Ang-1 ratio than healthy discs [23], and the Ang-1 release should be increased as the degenerated discs present blood vessels ingrowth.

Another strong limitation of the study is that we only used the Tie2− population as a control population. Previous studies [1,37] investigated Tie2 blocking antibodies as a negative control. For our study, a second negative control with a Tie2 blocking antibody would have made sense, as we observed, especially in the qPCR analysis, the effects in normoxia for Ang-1/2 at 10 ng/mL and Ang-1 at 100 ng/mL. Alternatively, other small inhibitory molecules could have been tested, which are known to interfere and block the Ang-1 and/or Ang-2 receptors. These inhibitors could have been tested by adding these in equimolar concentrations or by titration approaches. These are important negative controls that were not considered here.

#### *4.5. Outlook and Clinical Relevance*

In this study, we used human NP cells isolated from trauma patients' IVDs. We demonstrated the significant impact of culturing NP cells under hypoxic conditions to improve the proliferation, metabolism and gene expression of NP progenitor cells. However, no impact or effects were observed concerning the implementation of Ang-1 or Ang-2 in our culture medium.

The ensuing findings could directly impact clinical and therapeutic strategies in contrast to preclinical studies in other species that often need to be verified in humans first. Additionally, we provided novel evidence that culture within a hypoxic environment can drastically influence the proliferation, metabolism, and gene expression of the different cell populations present in human NP tissue. Moreover, we demonstrated that the hypoxic conditions affect more the progenitor cell populations from the human NP tissue (also known as Tie2+) than the two other populations (mixed and Tie2−). All those findings could lead to a new way to culture NP progenitors cells prior to clinical applications. Several studies have demonstrated Ang-1 and Ang-2 act in a concentration-dependent manner [17,44–47]. Beyond the crucial role in blood vessel development, recent studies suggest a broader role for Tie2, and the importance of Ang-1 and Ang-2 seem to be evident, as they are ligands (agonists) of the Tie2+ receptor and, therefore, expressed in NPPC. Nevertheless, previous studies, such as those of Sakai et al. and Wang et al., demonstrated

the effects of Ang-1 on cell proliferation and of Ang-2 on cell viability and cell adhesion. Their findings prompt us to clarify the favorable methods, adjust and improve our study design, and to further investigate the effects of Ang-1 and Ang-2 on NP cells. Despite our results, Ang-1 and Ang-2 may constitute a novel therapeutic target for the treatment of IVDD. A better understanding of the influence of Ang-1 and Ang-2 on the NPPCs would be of eminent importance in the context of IVD tissue engineering and regenerative medicine.

Other parameters that could have been taken into consideration as well are the GAG content, and the measurement of inflammatory cytokines. Moreover, HIF1A should be investigated at the protein level with Western blot analysis or immunohistology, as the HIF1A mRNA expression could have a fast turnover, especially in prolonged hypoxia [42].

#### **5. Conclusions**


**Author Contributions:** Conceptualization, M.C.B., B.G. and J.G.; methodology, M.C.B. and A.O.; formal analysis, M.C.B. and J.G.; investigation, M.C.B. and J.G.; resources, C.E.A. and B.G.; data curation, M.C.B. and J.G.; writing—original draft preparation, M.C.B. and J.G.; writing—review and editing, M.C.B., B.G., A.S.C., S.H., K.A.C.O. and J.G.; visualization, J.G. and B.G.; supervision, J.G. and B.G.; project administration, J.G. and B.G.; funding acquisition, B.G. and C.E.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** Financial support was received from the iPSpine H2020 project #825925 entrusted to B.G, www.iPSpine.eu (accessed on 1 August 2021), https://cordis.europa.eu/project/id/825925 (accessed on 1 August 2021). Further funds were received from the Swiss Society of Orthopedics (SGOT), the clinical trials unit of Bern University Hospital, and by a Eurospine Task Force Research grant #2019\_22 entrusted to C.E.A.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Ethics Committee of University of Bern (# 2019-00097).

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The cytometer equipment was from FACSlab core facility of the University of Bern (Available online: www.facslab.unibe.ch (accessed on 1 August 2021).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **EGR2, IGF1 and IL6 Expression Are Elevated in the Intervertebral Disc of Patients Suffering from Diffuse Idiopathic Skeletal Hyperostosis (DISH) Compared to Degenerative or Trauma Discs**

**Benjamin Gantenbein 1,2,\*,†, Rahel D. May 1,†, Paola Bermudez-Lekerika 1,2, Katharina A. C. Oswald 2, Lorin M. Benneker <sup>2</sup> and Christoph E. Albers <sup>2</sup>**


**Abstract:** Diffuse idiopathic skeletal hyperostosis (DISH) is characterised by ectopic ossification along the anterior spine and the outer intervertebral discs (IVD). However, the centre of the IVD, i.e., the nucleus pulposus, always remains unaffected, which could be due to the inhibition of the bone morphogenetic protein (BMP) pathway. In this study, we investigated the transcriptome for the BMP pathway of DISH-IVD cells versus disc cells of traumatic or degenerative discs. The disc cells originated from nucleus pulposus (NP), annulus fibrosus (AF) and from cartilaginous endplate (CEP) tissue. Here, ninety genes of the transforming growth factor β-BMP signalling pathway were screened by qPCR. Furthermore, the protein expression of genes of interest was further investigated by immune-staining and semi-quantitative microscopy. IVDs of three DISH patients were tested against three control patients (same disc level and similar age). Early Growth Response 2 (*EGR2*) and Interleukin 6 (*IL6*) were both significantly up-regulated in DISH-IVD cells compared to controls (12.8 ± 7.6-fold and 54.0 ± 46.5-fold, respectively, means ± SEM). Furthermore, Insulin-like Growth Factor 1 (IGF1) tended to be up-regulated in DISH-IVD donors, i.e., 174.13 ± 120.6-fold. IGF1 was already known as a serum marker for DISH and other rheumatoid diseases and is confirmed here to play a possible key role in DISH-IVD.

**Keywords:** low back pain; diffuse idiopathic skeletal hyperostosis; RNA extraction; qPCR; TGFβpathway; immune-histochemistry; Interleukin 6; Insulin-like growth factor 1; early growth response 2

#### **1. Introduction**

Diffuse idiopathic skeletal hyperostosis (DISH), also known as Forestier's disease, was first described by Forestier and Rotes-Querol in 1950 [1]. The authors characterised this disorder by spinal stiffness, osteophytosis, and aberrant new-bone formation along the spine's anterior lateral aspect. The underlying causes are yet unknown and thus are called "idiopathic" [2]. DISH is generally characterised by ossification of ligamentous insertions (fibro-osteosis) and ossification or calcifications of tendons, ligaments, fasciae, joint capsule and annulus fibrosus (AF) fibres [3–5]. Furthermore, DISH increases the amount of normal cancellous and cortical bone. Moreover, a tendency of post-operative heterotopic bone formation has been observed [6]. The incidence of DISH, however, remains unknown, but it seems that factors of metabolic, endocrine, genetic, and epigenetics seem to be important

**Citation:** Gantenbein, B.; May, R.D.; Bermudez-Lekerika, P.; Oswald, K.A.C.; Benneker, L.M.; Albers, C.E. EGR2, IGF1 and IL6 Expression Are Elevated in the Intervertebral Disc of Patients Suffering from Diffuse Idiopathic Skeletal Hyperostosis (DISH) Compared to Degenerative or Trauma Discs. *Appl. Sci.* **2021**, *11*, 4072. https://doi.org/10.3390/ app11094072

Academic Editor: Roger Narayan

Received: 26 February 2021 Accepted: 25 April 2021 Published: 29 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

in the development of this condition [4,7]. This disease mainly affects the elderly and is more frequent in men (65%) than in women, and the prevalence also increases with age (48–85) [8,9]. DISH patients often suffer from severe pain and need special treatment. Hence, the disease is often co-diagnosed with osteoarthritis [10,11].

Recently, Kuperus et al. (2016) [12] investigated the relationship between DISH and IVD degeneration. These authors found that partial or complete bone bridges were often present in DISH cadaveric IVDs that were only rarely found in healthy cadaveric controls. However, the disc height and degeneration level of DISH and control specimens were comparable. This specific feature also distinguishes the DISH condition from spondylosis [11].

To date, only very few molecular pathway studies have been conducted on the effect of DISH on IVDs. ElMiedany et al. (2000) proposed that the new formation of bone is mainly driven by abnormal cell growth/activity in the bony-ligamentous region [13]. Previous studies on the metabolomics of this disease found that patients suffering from symptomatic DISH disease have elevated levels of growth hormone (GH) and insulin-like growth factor 1 (IGF1) compared to patients of the same age and gender [9,14–19]. Thus, the molecular mechanism of DISH and the specific characteristics of underlying cell signalling remain obscure. Bone morphogenic proteins (BMPs), such as BMP2, BMP4, BMP6, BMP7, and BMP9, have shown their potential to induce the osteogenic differentiation of stromal cells and osteoblastic lineage cells in vitro and in vivo [20–23]. Furthermore, the recent literature has summarised and questioned the supra-physiological usage of BMP2 for the application of improved spinal fusion [24]. Thus, recently it has been hypothesised that BMP inhibitors could be the main reason for the lacking effect of BMP2 for improved spinal fusion [25–28]. In this context, it is of uttermost interest to analyse the transcriptome of DISH-IVDs in more detail, as all the surrounding tissue starts to ossify but not the IVD itself, which seems to be resistant to osteogenic signalling [29].

Hence, the current study aimed to compare the gene expression at the BMP pathway of cells of IVDs from patients suffering from DISH disease and compare it with cells originating from traumatic or degenerative IVDs. We assumed that trauma discs, especially from younger patients, contain presumably "healthier" cells that could be defined as a "control". Here, we investigated whether DISH disease is linked to alterations in BMP/TGF*β* signalling pathways to gain further insights into how the DISH-IVD was affected by this idiopathic disease at the molecular level. We hypothesised that IVDs of patients diagnosed with DISH demonstrate specific gene expression changes affecting the TGF*β*/BMP pathway compared to degenerative IVDs or trauma patients.

#### **2. Materials and Methods**

#### *2.1. IVD Donor Materials and Cell Isolation*

Human IVDs were obtained from patients suffering from DISH disease. As controls, traumatic or degenerative IVDs of patients undergoing spinal surgery were used (Table 1). Patients gave their written consent, and procedures were approved by the ethics committee of the Canton of Bern. Tissues were processed within 24 h after surgery.

Human IVD tissue was separated into nucleus pulposus (NP), annulus fibrosus (AF) and cartilaginous endplate (CEP) tissue by an experienced surgeon or processed as mixed IVD. Cell isolation from the native extracellular matrix was performed as described earlier [28]. Briefly, tissue was enzymatically digested by 1 h in pronase (Roche, Basel, Switzerland) followed by overnight incubation in collagenase type 2 (Worthington, London, United Kingdom). Subsequent expansion of NP cells (NPC), AF cells (AFC), and CEP cells (CEPC) followed in proliferation medium (low-glucose (1g/L) Dulbecco's Modified Eagle Medium (LG-DMEM) supplemented with 10% FBS and 1% P/S).

IVD cells were obtained from six DISH donors aged from 65 to 84 years (72.3 ± 3.1 years) and three "control" donors, which were lacking symptoms of DISH, suffering from degenerated or traumatic discs aged from 17 to 66 years (45 ± 14.5 years).



#### *2.2. Analysis of Specific Gene Expression in Human Primary IVD Cells with Quantitative Polymerase Chain Reaction (qPCR)*

Cells were cultured at passage 1–2 (P1–2) in monolayer, and then trypsinised before 2 Mio of cells were lysed directly with TRI. The extraction of total RNA was performed as described in the methods before [28,30]. Prior qPCR RNA integrity and purity were checked on selected samples by the Experion™ Automated Electrophoresis System (Bio-Rad, Reinach, Switzerland). Any possible remaining DNA was digested by using DNase I (AMP-D1, Sigma-Aldrich, Buchs, Switzerland) for 15 min. Reverse transcription was performed by using all-in-one cDNA Synthesis SuperMix (Bimake distributed by LuBio-Science GmbH, Lucerne, Switzerland). Real-time PCR was performed using SYBR Green PCR master mix on a CFX96touch qPCR system (all from Bio-Rad). Ninety genes (Table 2) were tested with PrimePCR™ Pathway plate 96 well (PrimePCR™ TGFβ BMP Signalling Pathway Plus H96, 90 genes, Bio-Rad, US). The qPCR was run using a two-step protocol with an annealing temperature of 60 ◦C for 30 s and 95 ◦C for 5 s for 40 cycles (according to the recommendations of the manufacturer). Amplicon specificity was ensured by performing a melting curve analysis after cycling. CFX manager software version 3.1 (Bio-Rad) was used for the quantification analysis of normalised relative gene expression. The number of PCR cycles needed for each sample to reach the threshold level was recorded as the C*q* value [31]. Values were normalised relative to two reference genes (18S and ACTB). DISH and control cells were compared in scatter plots and RNA transcripts that differed more than four times were indicated in red (up-regulated in DISH patients) or in green (down-regulated in DISH patients).

**Table 2.** List of genes tested in PrimePCR™ TGFβ BMP Signalling Pathway Plus H96, 90 genes, Bio-Rad, US.



#### **Table 2.** *Cont.*


#### **Table 2.** *Cont.*

#### *2.3. Statistical Analysis and Visualisation of Heat Maps and Clustering*

Relative gene expression was quantified between DISH and "control" patients. Statistical significance was tested using a non-parametric Mann–Whitney U test due to relatively low sample sizes using Prism 7.0d for Mac OS X (GraphPad, La Jolla, CA, USA). Nonparametric distribution was assumed. Consequently, statistical significance was determined using the Mann–Whitney U test. A *p*-value < 0.05 was considered to be significant. Relative gene expression data were normalised using CFX manager software 3.1. Heat mapping and hierarchical clustering was performed in CFX manager based on the degree of similarity of expression for different samples, and a heat map was generated alongside as implemented in CFX manager software. Pairwise scatter plot analysis was performed between DISH tissues versus "control" tissues and relative gene expression differences > 4.0 were highlighted (Supplementary Figure S1).

#### *2.4. Immunocytochemistry*

Cells at P1/P2 of frozen DISH donor 1 and control donor 1 (Table 1) were thawed and seeded at a density of 8000 cells/mL in glass chamber slides with removable 12 well silicone chambers for immunofluorescence staining (growth area per well 1.9 cm2) (Vitaris, Baar, Switzerland). Cells were then cultured in LG-DMEM supplemented with 10% FBS and 1% P/S until they reached 70% confluency. Cells were then fixed in 4% buffered paraformaldehyde (Sigma-Aldrich) for 15 min and stored in 70% ethanol at 4 ◦C prior to staining. Fixed cells were permeabilized with 0.2% Triton (Sigma-Aldrich) for 15 min. Again, the cells were washed three times with Tris-buffered saline (TBS) for 5 min. Subsequently, cells were blocked with 10% FCS in 0.025% Tris-buffered saline

containing Tween 20 (TBS-T) for 30 min at room temperature. After blocking, cells were incubated overnight at 4 ◦C using a primary antibody (AB) in blocking solution and under continuous agitation (see Table 3). Secondary ABs were applied in 1:200 dilution in 0.025% TBS-T for 3 h at room temperature (Alexa Fluor® 555 Rabbit Anti-Goat IgG (H+L), Alexa Fluor® 555 Goat anti-mouse SFX-kit, Alexa Fluor® 488 Goat anti-Rabbit SFX Kit, Molecular Probes, Life Technologies, Basel, Switzerland). Cells were viewed under a microscope (Leica DMI4000 (Biosystems inc., Muttenz, Switzerland) and LAS AF Software).

**Table 3.** Primary antibodies used for immunohistochemistry.


#### **3. Results**

#### *3.1. Prime PCR of DISH and Traumatic/Degenerative Discs*

PrimePCR™ screening for 90 genes of the BMP pathway was initially run. Most strikingly, these data revealed a significant mean up-regulation of *EGR2* (12.8 ± 7.6-fold, *p* = 0.0068, overall mean ± SEM of all comparisons) and of interleukin 6 (*IL6*) (54.0 ± 46.5-fold, *p* = 0.0005) in DISH-IVDs relative to the "control" discs (Figure 1, Supplementary Figure S1). Furthermore, *IGF1* tended to be up-regulated in DISH-IVD donors (174.1 ± 120.6-fold, *p* = 0.1704), although this was non-significant. Growth and Differentiation Factor 5 (*GDF5*) and *GDF6* both tended to be down-regulated in the DISH-IVDs (i.e., −11.5 ± 10.0, *p* = 0.26 and −3.7 ± 3.1 *p* = 0.30, respectively, Figure 1, Supplementary Figure S1). However, these differences were not significant.

**Figure 1.** Normalised relative gene expression in DISH-IVDs relative to IVDs of "control" patients. The relative gene expression of IL6 and EGR2 was significantly up-regulated: \*\*\* *p* < 0.01. The scatterplots of the pairwise comparisons of the 90 PrimePCR™ gene array is given in Supplementary Figure S1. Means ± SEM.

Clustering the gene expression profiles of the 19 analysed samples of the three control patients and the six DISH patients revealed that not all samples clustered together as expected (Figure 2). However, it seems also evident that the hierarchical phenogram generated two main clusters of phenotypes with a similar TGFβ pathway pattern (Figure 2). Clade 1 consisted of mainly DISH samples, and clade 2 consisted mainly of control samples. However, samples of DISH patients 3 and 5 (i.e., D P3 IVDC and D P5 AFC and EPC) showed a phenotype more similar to the cells of the included trauma patients. On the other hand, the dendrogram confirmed that cells originating from different tissue types of the IVD from the same patient were grouped closely together, which was expected.

**Figure 2.** Clustergram displaying hierarchy based on the degree of similarity of expression for different samples and normalised gene expression heat map of relative gene expression at 88 genes of PrimePCR™ of the TGFβ-BMP Signalling Pathway Plus H96. Red indicates up-regulations, green indicates down-regulations and black indicates no regulations. The lighter the shade of colour, the greater the relative normalised expression difference. If no normalised expression value could be calculated, this is shown in black with a white X. Relative gene expression was estimated using ACTB and GAPDH as the two reference genes. Samples encoded in orange are samples categorised as trauma patients = control (C), and samples in light blue are from patients that were diagnosed with DISH (D). Tissue types were colour-coded in dark blue for CEP, in yellow for NP, and in red for AF.

#### *3.2. Immunocytochemistry of IGF1, IL6, EGR2, BMP2, GDF5, and GDF6*

Gene expression differences, which showed a more than four times up- or downregulation in PrimePCR™ between DISH and "control" patients (Figure 1), were further evaluated using immune-staining of isolated cells grown in monolayer. BMP2 showed a stronger signal in the AFC in DISH than in "control" patients, as found in the relative gene expression (Figure 3). Furthermore, the EGR staining was more intense in the AFC and NPC in DISH compared to control patients' staining, as was also found in PrimePCR™. In particular, the AFC showed a strong immunocytochemistry staining for EGR2 and BMP2 compared to the AFC of the control patient. Additionally, IGF1 staining of DISH was more intense compared to control patients, as the NPC and the AFC showed almost no IGF1 protein expression by immunocytochemistry staining. GDF5 and GDF6 were overall more strongly stained in control cells compared to DISH cells.

**Figure 3.** Fluorescence microscopy of immune-histochemistry of nucleus pulposus (NPC), annulus fibrosus (AFC) and cartilaginous endplate cells (CEPC) of DISH versus "control" patient 1 using human GDF5, human IL6, human IGF1, human BMP2, and human EGR2, and human GDF6 antibody and the secondary antibody labeled with the fluorescent dye Alexa 555 (red), and Alexa 488 (green), respectively.

#### **4. Discussion**

In this study, we report, for the first time, on the phenotypic changes in IVD cells from DISH patients compared to IVD cells from supposedly "healthier" cells from younger and mid-aged trauma patients (Table 1). We observed the consistent up-regulation of IGF-1, and EGR2, which seems to be a hallmark of the DISH-IVD cells [4,29], and cells from degenerative discs showed a different BMP gene-related transcript and protein expression (Figures 1–3). One of the main limitations of this study was the relatively small sample size. Thus, sample collection from a single hospital was a challenge, and future studies should recruit IVD material from bigger cohorts and multiple populations and etiologies. Additionally, the "controls" used in the study were not truly healthy IVDs, as they were either degenerative or traumatic discs. In this respect, a comparison of fresh cadaveric IVD material of undegenerated nature would be essential for any future investigations. Acquiring more DISH and matching healthy IVDs are required to support our first findings regarding the phenotypic status of the IVD during DISH disease. Moreover, it should be investigated whether DISH affects the overall phenotype of IVD cells and their tissues. To answer this particular question, complete transcriptome analyses using bulk RNA nextgeneration sequencing would definitively shed new light on the better understanding of DISH pathology and into yet-unidentified modified pathways of this disease.

Since the discovery of DISH, various studies, especially concerning the metabolic conditions of these patients, were conducted. Nevertheless, this disease is still poorly understood. In the current study, we investigated the changes at the transcript and the protein expression level of IVDs from patients diagnosed with DISH. We observed significant up-regulation of EGR2 and IL6 and a trend of higher IGF1 expression at the transcript level. These changes could be confirmed qualitatively at the protein level using immunohistochemistry. For GDF5 and GDF6, we noticed, for both genes, a down-regulation in gene expression in the IVD of those patients, though this was not statistically significant. However, these differences were not apparent at the protein level (Figure 3). This gene expression profile points to an involvement in inflammation in DISH, recently identified in a review by Mader et al. (2021) [29].

Previous research by Levi et al. (1996) [32] has shown that EGR2-/- mice develop skeletal abnormalities, i.e., reducing length and thickness in newly formed bone and calcified trabeculae. These phenomena were explained by the activation of EGR2 in a subpopulation of hypertrophic chondrocytes of the growth plate and differentiating osteoblasts (OBs) [32]. Furthermore, Leclerc et al. (2008) [33] showed up-regulation of EGR2 expression during the late stage of OB differentiation in-vitro. Furthermore, Zaman et al. (2012) [34] found that EGR2 expression was up-regulated in the osteosarcoma-derived UMR106 cell line through IGF1 stimulation. The cause of the observed up-regulation of EGR2 and its role in DISH-IVD cells would need to be further explored.

DISH has been associated with the change in inflammatory mediators [35]. Interleukin 1 (IL1) and IL6 are increased in DISH patients through the activation of nuclear factor kappa-B (NF-κB) ligand, and these cytokines stimulate the proliferation of OBs and bone deposition [29,35,36]. Two other studies [37,38] have also suggested the stimulating effect of IL6 on bone formation and turnover. Kang et al. (1996) [39] further demonstrated an increased concentration of IL6 in herniated lumbar discs, which may indicate an association with degradation of the IVD. We could most interestingly confirm the involvement of IL6 in our data, as it was up-regulated in all analysed DISH-IVDs (Figure 1, and Supplementary Figure S1).

Several studies have shown the possible overlap between DISH, obesity and Diabetes mellitus, with a pathogenic effect of insulin, IGF1 and GH [40–43]. Yakar et al. (2002) [44] demonstrated an association of circulating blood serum IGF1 and bone mass in mice with reduced IGF1 serum levels. These mice showed abnormalities in bone formation including decreased bone mineral density and cortical thickness [44]. However, IGF1 is not only present in the serum but also readily available in the skeletal environment. The skeletal homeostasis is driven by the release of IGF1 during bone modelling and remodelling, as an essential mitogen and chemotactic player in the skeleton [45]. It also could have been shown that IGF1 stimulates IVD cell proliferation and matrix synthesis in-vitro [46].

GDF5, GDF6, and BMP2 are members of the TGFβ superfamily, with GDF5 and GDF6 being involved in the growth and development of cartilage and bone tissues [47]. Moreover, in the IVD, GDF5 seems to play a role in altering homeostasis, as several studies showed its beneficial effect of inducing cell proliferation, proteoglycan stimulation and COL2 expression [48,49]. Furthermore, GDF6 plays a role in spinal column development, especially in IVD development and homeostasis [50,51]. Both GDF5 and GDF6 have recently been shown to promote mesenchymal stromal cells towards a discogenic phenotype [52–56]. Interestingly, both GDF5 and GDF6 were down-regulated in several cell types of DISH-IVDs compared to controls. This might indicate a DISH-related change of the discogenic phenotype in these patients.

With increasing rates of obesity and type 2 Diabetes, the prevalence of DISH is also expected to rise in the future, possibly also affecting younger patients below 50 years of age [29,57]. This development is leading to an urgency to further understand and investigate this disease and its cell signalling.

#### **5. Conclusions**


**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/10.3390/app11094072/s1: Figure S1: TGFβ pathway transcriptome analysis of three donors with age- and disc level matched controls.

**Author Contributions:** B.G. and R.D.M. were responsible for the conception and design of the study. R.D.M. performed the experiments, data analysis and prepared the figures. R.D.M. drafted the manuscript; B.G. did extensive editing and prepared all of the figures; P.B.-L. and K.A.C.O. reviewed and edited the manuscript; C.E.A. and L.M.B. provided donor materials and clinical input and contributed to editing. All authors participated in the interpretation of the findings and approved the final version of the manuscript. B.G. and C.E.A. provided funding. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported from a start-up grant from the Center for Applied Biotechnology and Molecular Medicine (CABMM) to C.E.A. and B.G, and by funds from the Marie Skłodowska Curie International Training Network (ITN) "disc4all" (https://cordis.europa.eu/project/id/955735, accessed on 28 April 2021). Further funds were received from the clinical trials unit of Bern University Hospital and by a Eurospine Task Force Research grant #2019\_22.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** All primary human cells presented in this article were obtained with written consent from the donors. The study was conducted in accordance with the Declaration of Helsinki, and the Ethics Committee approved the protocol of the Canton of Bern (Ref. 2019-00097).

**Data Availability Statement:** All data presented in this manuscript can be obtained upon request from the corresponding author.

**Acknowledgments:** We thank Andrea Oberli, Selina Steiner, and Eva Roth for laboratory assistance. The microscopes were provided by the microscope core facility of the University of Bern (www.mic. unibe.ch, accessed 28 April 2021).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**


#### **References**


### *Article* **Diagnostic Limitations and Aspects of the Lumbosacral Transitional Vertebrae (LSTV)**

**Franz Landauer <sup>1</sup> and Klemens Trieb 1,2,\***


**Abstract:** The regeneration of an intervertebral disc can only be successful if the cause of the degeneration is known and eliminated. The lumbosacral transitional vertebrae (LSTV) offer itself as a model for IVD (intervertebral disc) regeneration. The aim of this work is to support this statement. In our scoliosis outpatient clinic, 1482 patients were radiologically examined, and ambiguous lumbosacral junction underwent MRI examination. Patients with Castellvi classification type II–IV were included and the results are compared with the current literature in PubMed (12 October 2022). The LSTV are discussed as a possible IVD model. A total of 115 patients were diagnosed with LSTV Castellvi type II–IV. A Castellvi distribution type IIA (n-55), IIB (n-24), IIIA (n-20), IIIB (n-10) and IV (n-6) can be found. In all, 64 patients (55.7%) reported recurrent low-back pain (LBP). Scoliosis (Cobb angle >10◦) was also confirmed in 72 patients (58 female and 14 male) and 56 (75.7%) had unilateral pathology. The wide variation in the literature regarding the prevalence of the LSTV (4.6–35.6%) is reasoned by the doubtful diagnosis of Castellvi type I. The LSTV present segments with reduced to absent mobility and at the same time leads to overload of the adjacent segments. This possibility of differentiation is seen as the potential for a spinal model.

**Keywords:** lumbosacral; transitional; vertebrae; intervertebral disc; scoliosis; Castellvi classification; lumbalgia

#### **1. Introduction**

In the introduction, the state of knowledge of the lumbosacral transitional vertebrae (LSTV) is shown on the basis of the current literature. The degeneration of the intervertebral discs is the focus of the work. The literature is reviewed for links between degeneration and regeneration of the intervertebral discs. Attention is paid to the clinical symptoms of lumbar complaints and radiological criteria.

The lumbosacral transitional vertebrae (LSTV) are a congenital malformation of the spine. It is characterized by enlargement of the L5 transverse process(es), with potential pseudoarticulation or fusion with the sacrum. It leads to altered rotational movement of the lower spine, which gives rise to low-back pain (LBP). A diagnosis of Bertolotti syndrome is given to patients experiencing pain caused by the presence of the LSTV. Estimates of the prevalence of the LSTV in the general population vary widely in the literature, ranging from 4.0% to 35.9%. Apazidis et al. have found that type IA is the most common, with a prevalence of 14.7%. However, type I is generally considered to be clinically insignificant [1].

In the Castellvi classification, only type I uses a dimension of >19 mm. Castellvi types II–IV, on the other hand, are defined by a structural change as pseudarthrosis or ossification. This is important because a differentiation between normal structure without altered biomechanics and pathological structure with expected symptoms of discomfort cannot be justified by a measurement of 19 mm. The scattering of the prevalence, which is repeatedly found in the literature, is only one indication of this (Figures 1 and 2).

**Citation:** Landauer, F.; Trieb, K. Diagnostic Limitations and Aspects of the Lumbosacral Transitional Vertebrae (LSTV). *Appl. Sci.* **2022**, *12*, 10830. https://doi.org/10.3390/ app122110830

Academic Editor: Benjamin Gantenbein

Received: 2 October 2022 Accepted: 24 October 2022 Published: 25 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

97

**Figure 1.** Castellvi classification.

**Figure 2.** Castellvi IIIA (dotted area shows pathology).

The topic of the inaccuracy of the representation of the L5 transverse process will therefore be discussed in detail at this point. In Figure 1, Castellvi type I is represented as a sketch in such a way that a slight to more severe restriction of the mobility of the affected segment can be assumed (Figure 1). As reflected in the review of the following international literature, this problem has not yet been conclusively addressed. The problem of the imaging quality of the LSTV is presented in a value-neutral manner, even though this topic is repeatedly addressed in the further course. Depending on the pelvic shape and pelvic tilt, the radiological ap image of the lumbosacral junction leads to superimpositions in the case of an enlarged L5 transverse process, which do not always permit measurement and assignment with certainty. The measurement of the L5 transverse process (>19 mm) is thus not guaranteed. As outlined in Figure 3, the imaging result is significantly influenced by the radiological beam path. The imaging technique according to Fergusion with a tube tilt of 30◦–40◦ cranially and caudally suggests improved imaging (Figure 3). However, the dimension of 19 mm remains a variable dependent on many factors (patient body size, tube–object distance, object–film distance, etc.). The pelvic setting in relation to the spine is not only crucial for visualizing pathology, but also an important criterion for possible causes of complaints. Especially the expected disc load in the context of the LSTV is crucial. If now again the question of the cause of complaints in Castellvi type I arises, it becomes clear that a cause differentiation between the LSTV and a pathological pelvic setting is currently hardly possible. No conclusive answer can be found in the literature either. An interaction between the LSTV and a hip complaint can be assumed in individual cases. We therefore always recommend a simultaneous hip examination (Figures 3 and 4).

**Figure 3.** Course of the X-ray beam.

**Figure 4.** Pelvic setting.

The availability of an MRI examination significantly improves the diagnosis. However, even this examination technique must be specific to the question of the LSTV. Initial referrals of our own to various institutes yielded astonishing results that rarely allowed a diagnosis to be made. This ranged from imaging of only the sagittal and axial slices, which is often sufficient for disc assessment, to imaging that did not portray the L5 transverse process. The combination with a 1.5 Tesla device and a delivered slice thickness of 3 mm also does not allow a more accurate Castellvi type I assessment. These experiences led to arranging for standardized images, as described in the method below, after consultation with one MRI institute. These findings of our own and the still sparse literature at the beginning of the evaluation period led to the decision not to pursue Castellvi type I further for the time being (Figure 5).

**Figure 5.** Representation of the L5 transverse process(es) in MRI or CT.

In a special issue on the subject of "Intervertebral Disc Regeneration", the question of the ligamentous situation at the lumbosacral junction must also be considered. The findings available in the literature to date are presented, and the question of the quality of the imaging is part of the problems described above. In particular, the question of mobility restriction of the affected segment and the associated intervertebral disc stress becomes the focus of consideration. This makes it obvious that the Castellvi classification type I is not suitable for differentiating between a normal segment L5/S1 with a pronounced L5 transverse process and a pathological change. This illustration should clarify why the assessment of Castellvi type I was excluded in this study, although these changes are crucial not only for disc degeneration but even more so for the outcome of disc regenerative treatment.

In the work of Min Tang, a prevalence of 15.8% for the LSTV is diagnosed in 5860 radiographs. In all, 44.8% of them are classified as type I [2]. Regarding the question of diagnosis of the LSTV comparing X-ray and MRI, a good agreement was established very early [3].

In the work of the last 10 years, the problem of classification of the LSTV is increasingly worked on. Thus, the imaging quality and reliable assignment of the LSTV, especially of type I, is becoming the focus of discussion.

Farshad et al. advocate measuring segmental differences. Measurement methods in the sagittal slices of MRI images are described for identification of LSTV. However, the variability of the vertebral bodies makes it difficult to accurately define L5 and S1 [4].

Of importance is the total number of vertebral bodies diagnosed with 23 (sacralization) or 25 (lumbarization) with a normal value of 24 vertebral bodies. Thus, a total imaging of the spine is necessary [5].

Quantitative measurements are used to try to find new parameters for the diagnosis and differentiation of the LSTV. Especially LSTV type I with a height of the transverse process of >19 mm forms an uncertainty in the presentation and therefore in the diagnosis. The previously mentioned large differences in the prevalence of the LSTV are partly due to this [6].

In the article by Konin et al., attention to the identification and correct numbering of the LSTV and recognition of imaging findings related to the development of LBP are comprehensively addressed [7].

In previous work, the iliolumbar ligament is recommended for identification of L5/S1 [8]. However, attempts to identify L5/S1 using the iliolumbar ligament have not been successful to date [9].

The iliolumbar ligament in its anatomical structure and as a possible cause of complaints is currently coming into focus. In an MRI study by McGrath et al., unilateral pathology (types IIA and IIIA) shows a thinner ligament than on the unaffected side. Bilateral pathology (types IIB, IIIB, and IV) shows no side difference and no difference from the healthy control group. However, the differentiation of LSTV type I is not answered [10].

For a special issue, "Intervertebral Disc Regeneration II", the cause of LBP and disc degeneration is the focus of the LSTV.

The relationships between IVD degeneration (intervertebral disc degeneration) and the LSTV are differentiated in the literature. In LSTV type I, the highest numbers are reported in some papers; however, no pathological effect on the disc segment is described. In contrast, in LSTV type II–IV, the adjacent cranial disc segment is reported as the segment most exposed to degeneration. The distinction between unilateral and bilateral pathology forms another focus, which is also associated with the question of scoliotic change. A clear prediction cannot be derived so far.

Pathologies of adjacent structures as a cause of complaints are also increasingly subject to clarification. The segmental nervous assignment was examined thereby, too.

The disc degeneration is concentrated on the cranially adjacent segment L4/5. In contrast, the disc of the affected segment of the LSTV type II–IV is mechanically protected by the restriction of movement. However, this general statement still lacks a differentiated assessment dependent on the type of the LSTV. The biomechanical workup of motion restriction and the associated disc loading or unloading are at the beginning of scientific interest [11].

A paper by Ashour et al. states that the incidence of disc herniation is decreased in the segment of the LSTV. In contrast, the incidence of disc herniation or spondylolisthesis in the cranial segment L4-5 is increased in the presence of the LSTV. Therefore, the LSTV are considered a risk factor for degenerative disc changes in the level above the transitional vertebra [12]. In the LSTV type II group, compression of the L5 nerve due to a bony or synovial cause of nearthrosis has also been reported [13]. The question of nerve innervation in the LSTV has also been raised by intraoperative EMG studies, and no pathological distribution was found compared with the normal spine [14]. A significant association between lumbar stenosis and the LSTV is described in a retrospective cross-sectional study [15]. The bony structures show altered trabecular distribution in the LSTV compared to healthy spine in a cadaveric study [16]. In a prospective cadaver study by Golubovsky et al. the segmental mobility of the affected segment and adjacent segments was measured and described as follows:

The results of this study suggest that patients with Bertolotti syndrome have significantly altered spinal biomechanics and may develop pain due to increased loading forces at the LSTV joint during ipsilateral lateral flexion and axial rotation. In addition, increased motion in the upper levels in the presence of the LSTV can lead to degeneration over time [17].

After outlining the pathological changes and movement limitations, the focus turns to the possible causes of LBP:

There are different answers to the question of whether LSTV type I causes symptoms, although diagnostic uncertainty remains as described above. There is also no clear answer to the distinction between unilateral and bilateral changes as the cause of symptoms. However, the more stable the constructs of the LSTV are (type II–IV), the more likely degeneration of the cranially adjacent segment is responsible for the complaints.

Lumbar complaints (LBP) in the presence of the LSTV were originally noted by Mario Bertolotti in 1917 and are therefore referred to as "Bertolotti syndrome" [18].

Quinlan et al. found Bertolotti syndrome in 769 MRI scans of LBP patients in 4.6% of the patients studied. In patients with LBP younger than 30 years, 11.4% showed LSTV [19].

A high frequency of the LSTV was found in young male patients with LBP, and it was mostly subtype I. LSTV type I and associated disc and facet degeneration were notable in this group [20].

A review of the literature on the problem of discomfort in the LSTV (Bertolotti syndrome) shows that comparable changes are found with and without discomfort.

LSTVs pose a diagnostic dilemma for the treating physician, as they may remain undetected on plain radiographs and even on advanced imaging; moreover, even if the malformation is recognized, patients with the LSTV may be asymptomatic or have nonspecific symptoms, such as low-back pain with or without radicular symptoms.

Special attention should be given to younger patients under 30 years of age with chronic lumbar symptoms. Here, diagnostic exclusion or confirmation of the LSTV is required [21].

In lumbalgia, differences in muscle mass distribution are also associated with lordosis of the lumbar spine [22]. The LSTV are associated with a reduction in muscle volume and an increase in muscle degeneration of both the lumbar and trunk muscles compared to healthy individuals [23]. A differential relationship between lumbalgia and the type of the LSTV can be established. However, the segment affected by the LSTV must be distinguished from the adjacent segment as the cause of the symptoms [24]. Unilateral LSTV result in asymmetric biomechanical stresses. The side with the extra L5/S1 joint bears a greater proportion of the load, resulting in lateral tilt with scoliotic curvature. This local loading leads to greater wear and unilateral muscle activity. In addition, asymmetric motion may also affect disc degeneration [25]. The LSTV may also be associated with spina bifida occulta as another pathological change [26,27].

The topic of scoliosis is addressed in 16 papers. In a study in the Turkish Army, the LSTV were diagnosed in 12.9% of 3132 men [28]. Another regional assignment of prevalence to the southern European population is reported by Vinha A. et al. with 24.9% for the group LSTV Castellvi type II, III and IV [29]. In the work of Garg B. et al., a prevalence of LSTV in 198 patients with scoliosis is found in 18.2% [30].

The change in segmental mobility is measured in a paper by Becker L. et al. and summarized as follows: Patients with the LSTV show reduced range of motion in the transitional segment and a significantly increased distribution of motion to the cranially adjacent segment on flexion–extension radiographs. The increased proportion of mobility in the cranial adjacent segment may explain the higher rates of degeneration within the segment [31]. The influence of the LSTV on pelvic alignment with acetabular orientation is also described by the same group [32].

The impact of the LSTV on surgical interventions is also scientifically discussed but is not the focus of this article. An improvement of the patients' quality of life is described by surgical interventions. This concerns both the lysis of LSTV and the fusion of the segment depending on the degenerative changes [33]. However, surgical problems due to the altered anatomy, especially of the ventral vascular plant, must also be considered in surgical planning [34].

The LSTV have a wide range in diagnosis, especially in type I, which is reflected in the widely scattered prevalence data. The biomechanical limitations and the possible causes of the symptoms are not yet sufficiently differentiated. In particular, there is little information about the intervertebral disc in the LSTV. The aim of this paper is to investigate the diagnosis of the LSTV in relation to LBP. The prevalence of the LSTV in the literature and in our own work should be documented with respect to the problem of diagnosis of LSTV, especially of Castellvi type I.

#### **2. Materials and Methods**

The present article follows the current PRISMA criteria from 2020. The structure of the contribution is shown in a flow diagram [35] (Table 1).

**Table 1.** Structure of the contribution.

```
(1) Introduction:
```
Presentation of current literature of LSTV: Prevalence, diagnosis, Castellvi classification, LBP, Bertolotti syndrome, scoliosis, etc.

Flow diagram to Disc: PubMed (12 October 2022):

```
LSTV n = 150
 Prevalence n = 50
              Castellvi classification n = 51
                           LBP n = 59
                                        Bertolotti syndrome n = 19
                                                     Disc n = 45
                                                                  Scoliosis n = 16
```
Surgery *n* = 80

(2) Materials and Methods:


#### (3) Results:

Results of own LSTV-data (2014–2021) (gender, age, Castellvi classification, total number of vertebral bodies, low-back pain, scoliosis, concomitant pathologies, family history). Castellvi classification

total number of vertebral bodies

low-back pain (Bertolotti syndrome)

#### scoliosis

gender

age

concomitant pathologies family history

(4) Discussion of the literature and own data of LSTV:


(5) Summary:

LSTV as an experimental model for intervertebral disc degeneration/regeneration

Diagnosis of LSTV: In our scoliosis outpatient clinic 1482 patients were examined by one consultant orthopaedic surgeon (2014–2021) with more than 25 years of experience in that field.

Each radiological, structurally ambiguous assignment of the lumbosacral junction was submitted to further MRI examination with a 3.5 Tesla device (sagittal, coronal, axial and adjustment to the segment L5/S1 with a slice thickness of 3 mm, Siemens Medical Group).

Patients with Castellvi I classification (height of L5 transverse process >19 mm) are not included in the study. This decision was made in the initial phase of the study, because in our own unpublished studies the imaging inaccuracy allows an exact measurement of the L5 transverse process only in the Ferguson image. However, this additional X-ray is not permitted for reasons of radiation hygiene. On the other hand, 19 mm should be evaluated differently for a very tall person than for a short person.

Only patients with Castellvi classification type II–IV are included in the study.

#### **3. Results**

In 1482 spine patients examined, 115 patients were diagnosed with LSTV Castellvi II–IV (male n-32, female n-83).

In the Castellvi IIA (n-55) and IIB (n-24) groups, spondylolysis of the adjacent segment was seen in two patients and lumbarization in one patient, respectively. In patients with Castellvi IIIA (n-20), spondylolysis of the adjacent segment was found in three patients and there was lumbarization in one patient. In classification IIIB (n-10), spondylolisthesis (Meyerding 2) was seen. Six patients showed Castellvi type IV.

A total of 64 (55.7%) of all patients reported recurrent LBP. Of these, 23 patients showed orthograde lumbosacral transition and 41 patients showed scoliosis.

In unilateral pathology (IIA, IIIA, n-75), 39 complained of symptoms, and in bilateral pathology (IIB, IIIB, IV, n-40), 26 patients did. Only five patients belong to group IIIA and thus have neither unilateral nor pseudarthrosis as cause of pain.

In lumbalgia, a clear age-dependent differentiation is evident. In children and adolescents <16 years, 26.7% complained of chronic complaints. In the age group 16–30 years, it was 35.5%, and in >30 years it was 88.6% who reported chronic complaints.

Five patients were therefore already in psychological treatment. In the anamnesis, up to 30 visits to the doctor for lumbalgia can be recorded up to the time of the diagnosis of LSTV. Three patients with recurrent lumbalgia have already undergone spinal surgery without the LSTV being known to the patients.

In 72 patients, LSTV Castellvi II–IV and scoliosis was confirmed (58 female and 14 male). The curvature extent of scolioses showed a mean Cobb angel of 24.3◦ with a range from 11◦–55◦ (n-42) in children and in adults a mean of 32.4◦ with a range from 12◦–66◦ (n-30) in the adults.

Of the 72 patients with scoliosis, 56 showed unilateral pathology (IIA and IIIA) and 16 bilateral pathologies (IIB, IIIB, IV, Table 2).

A positive family history existed in 20 patients from the scoliosis group (*n* = 72). In one family, the mother and daughter were affected type IIA, and in another case it was twins with IIB with scoliotic component. Other malformations were raised: longitudinal malformation of the right upper extremity, spina bifida, other vertebral body malformations (n-3), a sacral cyst, syringomyelia (n-2) and a fork rib.

When comparing our own data with the results in the literature, it should be noted that in our own data only LSTV Castellvi type II, III and IV are collected. The data from the literature refer to the total sum of LSTV Castellvi type I–IV. If Castellvi type I is shown in the data, the value is given in parentheses (Table 3).


**Table 2.** Data of the own investigation.

**Table 3.** Own data in comparison to the literature.


#### **4. Discussion**

Castellvi type I:

In the initial phase of the examination, the Castellvi type I problem with its enlargement of the L5 transverse process of >19 mm becomes obvious. The standard radiograph of the lumbar spine, as well as the spinal X-ray for scoliosis clarification, does not result in a reliable orthograde representation of the L5 transverse process. Lordosis and pelvic tilt, as well as scoliosis, prevent standardized imaging. The magnification factor that inevitably occurs in a radiograph as a function of object–film distance and tube distance make the height indication of the L5 transverse process an unsuitable measurement for Castellvi I classification.

Only the Ferguson image (30◦ angle AP radiograph) serves as a standard method for detecting LSTV. Additional imaging is prohibited for radiation hygiene reasons. However, the problem of the magnification factor remains.

The alternative is an MRI examination, which provides additional information in all three planes in space. The standard images for disc assessment (coronal and sagittal) are not sufficient for the assessment of LSTV. The frequently used slice thickness of 3 mm and the quality of the device (<3.5 Tesla) form another uncertainty factor for diagnostics. If in the literature there are reported (in many cases not specified) investigations with a CT and MRI slice thickness of 0.6 mm–3 mm, this should not be considered an insignificant representation and therefore also not an inaccurate measurement [3,4,6,10].

However, the height of the L5 transverse process must also depend on the overall size of the patient, so that a measurement >19 mm is only of limited use for classification.

In conclusion, the shape of the L5 transverse process could be more indicative of pathology than the measurement alone.

In our own work, because of the intractability of the exact classification of Castellvi type I, we did not work on this type further. This does not mean that Castellvi type I cannot be responsible for complaints [8–10]. It should be emphasized that the scoliosis outpatient clinic, with its patients and the frequent oblique position of L5 compared to S1, presents an additional problem of presentation [1,2].

This raises the following questions:

Is pathology of the ligamentous apparatus (iliolumbar ligament) responsible for the shape of the L5 transverse process? If so, the shape of the L5 transvers process would be the expected identifying feature.

A standardized clarification of the ligamentous apparatus and appropriate classification must be required.

If it is a pathology with an effect on the intervertebral disc, then any interventional procedure on the intervertebral disc must also take into account a possible pathology due to LSTV. This assessment, at first sight very critical, becomes crucial when a new treatment modality is aimed at repairing the disc [37–39].

If the incidence of the LSTV indeed reaches values as reported in some literature between 4% and 35.9% and the LSTV leads to treatment failure, this would be crucial for success and failure for this new treatment modality [21].

#### LSTV type II:

In unilateral or bilateral pseudarthrosis, limitation of motion of the L5/S1 segment is assured. Therefore, the pseudarthrosis and the intervertebral disc must be differentiated as the cause of the complaint. It is to be expected that the pseudarthrosis cannot follow the age-related process of height reduction in the intervertebral disc without complaints. Corresponding complaints are thus to be expected. For this differentiation, further MRI criteria have to be defined.

We would like to point out another aspect to this.

As patients age, degenerative de novo curvatures or increases in curvature occur. This can lead to a contact zone of the L5 transverse process with the os sacrum (os ileum). If there is no previous X-ray from youth, the assignment as Castellvi II is also fraught with uncertainty. For the question of the IVD, which is the focus of this article, this may have no clinical significance. The literature provides the following evidence in this regard:

Type II LSTV had significant effects in promoting transitional and adjacent disc degeneration.

Thus, Castellvi type III/IV LSTV predisposed the adjacent spinal components to degeneration and protected the transitional discs [31,40].

#### LSTV type III:

In its unilateral form, stress loading is inevitable. In its bilateral form, stability of the L5/S1 segment is assumed. In both cases, L4/L5 overload is expected due to the lack of L5/S1 mobility [17,31].

#### LSTV type IV:

This is a bilateral form with a bony bridging and pseudarthrosis. In this case, an increased load on the adjacent segment L4/L5 may be assumed. Moreover, the pseudarthrosis may be responsible for discomfort since there is no symmetrical situation.

Disc degeneration:

In a paper by Ahmed Ashour et al., the incidence for disc prolapse is found to be decreased in the segment of LSTV. In contrast, the incidence for disc herniation or spondylolisthesis in the cranial segment L4-5 is increased in the presence of LSTV [12]. Hence, the LSTV are considered a risk factor for disc degenerative changes at the level above the transitional vertebra level. The IVD only makes sense (can be successful) if it does not affect the adjacent segment of the LSTV. Thus, a differential diagnosis of the LSTV is important before any IVD intervention.

Low-Back Pain (Bertolotti syndrome):

Explanation of own data:

A total of 64 patients of the total group (55.7%) reported recurrent LBP (orthograde lumbar spine *n* = 23, scoliosis *n* = 41).

The low number of LBP can be explained by the fact that the age limit of the patients is very low (scoliosis outpatient clinic). Additionally, the differentiation according to age shows a rapid increase in complaints in the patient distribution, and the dominance of scoliosis with 64.1% is striking.

Basis for an intervertebral disc model with different calculable limitation of motion:

A prospective cadaveric study measured the segmental range of motion of the affected segment and adjacent segments.

A biomechanical analysis of the tissue, joints and associated motion restriction (range of motion, resulting torques and the LSTV joint forces) and histological assessment of the segments was performed. The aim of this work is surgical therapy.

This study's results indicate that patients with Bertolotti syndrome have significantly altered spinal biomechanics and may develop pain due to increased loading forces at the LSTV joint with ipsilateral lateral bending and axial rotation. In addition, increased motion at superior levels when the LSTV are present may lead to degeneration over time.

However, the model can be adapted to address disc treatment issues.

The LSTV are a model given in nature for a differential study of motion restriction on the affected and adjacent disc segments [17,31].

The preceding illustration suggests the LSTV as a potential model for IVD research. The LSTV provide segments with reduced to absent mobility and at the same time an overloaded adjacent segment.

MRI imaging capabilities should differentiate between complaint origin in the LSTV and disc degeneration.

The model representation should also be suitable for further FEM calculation [36–40].

#### Disc model:

The work shown here offers the LSTV as a disc model for further investigation. As a model for the degeneration of the lumbosacral junction and thus clarification of the cause of the complaint, treatment options become delineable. It can thus lead to a further gain in knowledge for regenerative disc treatment. We are well aware of the limitations of the model, especially in the area of Castellvi type I. Knowing the limitations of imaging will be an advantage for every practitioner and will lead to the avoidance of frustrated therapy attempts in cases of mechanically persisting pathology.

Limitation of own data:

The data were collected in the spine outpatient clinic with a focus on scoliosis. This explains the high number of scolioses. This is clearly shown in Table 3 in comparison with the literature.

The lack of inclusion of LSTV Castellvi I cases prevents direct comparison with the prevalence data in the literature. A wide dispersion is shown, which cannot be narrowed down by our own work. This will be discussed in detail in the discussion.

Through cooperation with an MRI institute, an attempt was made to provide standardized imaging as compensation and thus to ensure a reliable assignment of the patients LSTV Castellvi II, III and IV.

#### **5. Conclusions**

This paper clarifies the diagnosis of the LSTV as a possible cause of LBP. From the prevalence of LSTV, in the literature and in our own work, it becomes clear that this is a clinically non-negligible number of patients. In the wide dispersion of prevalence, the problem of diagnosis of LSTV Castellvi type I becomes apparent.

Despite these aforementioned diagnostic problems, the LSTV offer the full spectrum of possible segmental movement limitations. This ranges from a slight limitation of the range of motion, without clinical relevance, to complete bony fixation. Thus, the cranially adjacent disc segment is also included in a possible functional model of the IVD.

This in vivo IVD model is thus suitable to provide further insights into the degeneration mechanism. Knowing the cause of degeneration is the basic prerequisite for successful regeneration treatment using the IVD. For the therapy decision of an IVD regeneration or a necessary alternative treatment procedure, the presented model can provide information in the future.

**Author Contributions:** Data curation, F.L.; Investigation, F.L.; Methodology, K.T.; Writing—original draft, F.L.; Writing—review & editing, K.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Ethical review and approval were waived for this study because only radiographs were analyzed.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**


#### **References**


**Yeji Kim 1,†, Chanyoung Song 2,†, Gyuseon Song 2, Sol Bi Kim 3, Hyun-Wook Han 2,\* and Inbo Han 3,\***


**Abstract:** A natural language processing (NLP) pipeline was developed to identify lumbar spine imaging findings associated with low back pain (LBP) in X-radiation (X-ray), computed tomography (CT), and magnetic resonance imaging (MRI) reports. A total of 18,640 report datasets were randomly sampled (stratified by imaging modality) to obtain a balanced sample of 300 X-ray, 300 CT, and 300 MRI reports. A total of 23 radiologic findings potentially related to LBP were defined, and their presence was extracted from radiologic reports. In developing NLP pipelines, section and sentence segmentation from the radiology reports was performed using a rule-based method, including regular expression with negation detection. Datasets were randomly split into 80% for development and 20% for testing to evaluate the model's extraction performance. The performance of the NLP pipeline was evaluated by using recall, precision, accuracy, and the F1 score. In evaluating NLP model performances, four parameters—recall, precision, accuracy, and F1 score—were greater than 0.9 for all 23 radiologic findings. These four scores were 1.0 for 10 radiologic findings (listhesis, annular fissure, disc bulge, disc extrusion, disc protrusion, endplate edema or Type 1 Modic change, lateral recess stenosis, Schmorl's node, osteophyte, and any stenosis). In the seven potentially clinically important radiologic findings, the F1 score ranged from 0.9882 to 1.0. In this study, a rule-based NLP system identifying 23 findings related to LBP from X-ray, CT, and MRI reports was developed, and it presented good performance in regards to the four scoring parameters.

**Keywords:** natural language processing; low back pain; lumbar spine imaging

#### **1. Introduction**

Low back pain (LBP) is defined as pain and discomfort localized under the ribs and above the inferior gluteal folds, with or without leg pain [1]. A total of 70–80% of adults experience LBP in some form during their lives [2]. More than 85% of people under 45 years of age experience at least one LBP symptom that requires medicine or interventional treatment [3]. LBP can be classified as acute or chronic, depending on its onset. Acute LBP usually occurs suddenly, and if the pain persists for more than 3 months, it can be called chronic LBP [4]. LBP is one of the most common causes of hospital visits and the second-leading cause of sick leave [5]. Because of its high direct and indirect costs, the health, social, and economic impacts on individuals, families, and society are significant [6]. Particularly, between 5% and 10% of chronic low back pain (CLBP) cases require high costs and long-term care [7]. It is important to develop tools to help patients with recently developed back pain accurately predict whether persistent pain will occur [8].

The most common methods to evaluate chronic LBP are X-radiation (X-ray) examinations, computed tomography (CT), and magnetic resonance imaging (MRI) [9]. The

**Citation:** Kim, Y.; Song, C.; Song, G.; Kim, S.B.; Han, H.-W.; Han, I. Using Natural Language Processing to Identify Low Back Pain in Imaging Reports. *Appl. Sci.* **2022**, *12*, 12521. https://doi.org/10.3390/app 122412521

Academic Editors: Valentino Santucci and Habib Hamam

Received: 17 October 2022 Accepted: 6 December 2022 Published: 7 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

radiological findings identified from reviewing those images are translated by clinicians and transferred to a radiology report. A radiology report is a formal interpretation of a radiological examination and contains radiological findings related to the clinical diagnosis and treatment decisions [10,11]. However, many radiologic reports use a free-text form language and remain unstructured, requiring a manual review to harvest clinical information. Manual extraction of information is labor-intensive and impractical for large-scale research. It is important to create a pipeline for extracting clinical information from the radiology reports.

Natural language processing (NLP) is defined as understanding, analyzing, and extracting meaningful information from text (natural language) by computer science [12]. NLP, a computer technology specializing in text processing, can be used to extract critical information from unstructured text in electronic health records [13]. The NLP system is classified into three approaches: the rule-based approach, the machine learning-based approach, and the hybrid types approach [12]. For example, a Bidirectional Encoder Representations from Transformers (BERT)-based NLP pipeline (a state-of-the-art deep learning model for language processing) was used to identify important findings in intensive care chest radiograph reports [14]. Emilien et al. used the NLP system developed by the BERT model. The free-text form records in the University Hospital of Amiens-Picardy in France was used to develop the BERT model. This BERT model learned contextual embeddings. In this research, the utilization of the BERT method in the medical care field was discussed [15]. Tan et al. developed a rule-based NLP system and a machine learningbased NLP system to identify lumbar spine imaging findings related to LBP and compared their performance [16].

In previous studies, there were some points that need to be improved. First, while the machine learning approach has been actively studied, the rule-based approach has not been studied much in the NLP system for extracting words associated with LBP. The set-up cost of the rule-based approach is lower than that of the machine learning approach. The training dataset size of the rule-based approach is also smaller than that of the machine learning approach. In previous studies, the rule-based approach presented high specificity, but moderate sensitivity [16]. This study developed an improved rule-based NLP system to that used in the previous research. Second, there is a need to develop the NLP system for extracting words associated with LBP in CT radiology reports. In previous studies, the NLP system was not trained and tested in the CT radiology reports [16–21]. In this study, the utilization of the NLP system was expanded to the CT radiology reports for extracting words associated with LBP.

This study aimed to develop an NLP system to recognize radiologic findings associated with LBP in X-ray, CT, and MRI radiologic reports. A rule-based approach of the NLP system was trained and evaluated with radiologic reports from patients who sought medical care for LBP. By developing an NLP system, it is possible to analyze the free-text from radiology reports unique to the medical institution. Unstructured clinical information from the radiology reports was processed into a structured data form. The structured form clinical data will be used for the data-drive research of LBP. This study laid the foundation for extracting the clinical data from free-text from radiology reports and using it in various clinical studies.

#### **2. Materials and Methods**

#### *2.1. Dataset*

This retrospective study of lumbar spine imaging reports enrolled patients who visited CHA University Bundang CHA Medical Center between January 2011 and December 2020. A dataset was assembled from 18,640 reports, and a sample of 300 X-ray, 300 CT, and 300 MRI reports was randomly extracted. With reference to the method of the NLP system for extracting words associated with LBP in previous research, 5% of these words were extracted from the original dataset and used for this study [17]. The dataset was selected randomly. The dataset was split into training and test datasets in an 8:2 ratio, maintaining the ratios of the three different modalities. The training dataset was used for developing the NLP system, and the test dataset was used to validate the information extraction performance. A consensus was reached through a discussion for any cases of discrepancies. This study was approved by the Institutional Review Board of the Bundang CHA Medical Center (IRB No. 2021-04-212).

#### *2.2. Inclusion Criteria*

The extracting of clinical information by the NLP system was carried out according to the inclusion criteria. The criteria were composed of radiologic findings in the form of words about spine pathology, structural problems of the spine, and spinal disease. A total of 23 radiologic findings known to be related to LBP from previous studies were established, and each report was annotated by two physicians specializing in LBP diagnosis [18,22] (Table 1). In the 23 radiologic findings, 7 items—disc extrusion, endplate edema or Type 1 Modic, any stenosis, central stenosis, foraminal stenosis, nerve root displaced/compressed and lateral recess stenosis represent the potentially clinically important findings [16].

**Table 1.** 23 Radiologic findings of the study.


#### *2.3. NLP System*

This NLP pipeline was developed using Python (3.8.10). Section and sentence segmentation from the radiology reports was performed using a rule-based approach, using regular expression with negation detection. A regular expression is a set of characters specialized for analyzing patterns in text.

The "re" module is one of the Python libraries used to implement regular expressions in Python. By using "re", syntactic analysis was performed to develop clinical information extraction rules. The training dataset was analyzed to construct a customized dictionary of terms associate with the 23 predefined LBP features. From the X-ray, CT, and MRI radiology reports, an NLP pipeline was developed to extract only terms relevant to LBP diagnosis. The pipeline extracted LBP information in three steps: text preprocessing, concept mapping, and summarizing.

First, in the text preprocessing step, tokenization was performed. The NLP refined the unstructured raw text to remove unexpected white spaces and punctuation and filtered out sentences unrelated to LBP. In this step, section segmentation, sentence segmentation, and normalization is performed. Second, LBP concepts were mapped by referring to the constructed dictionary and linguistic rules. Finally, in the summarizing step, the mapped concepts were summarized in a predefined format (Figure 1).

**Figure 1.** Summary of the natural language processing pipeline.

The whole process was repeated until the performance of this NLP system was improved. The k-fold cross validation was used for the NLP model validation. The NLP pipeline version with the best performance was validated by the test dataset. Therefore, the words associated with LBP were extracted from the X-ray, CT, and MRI reports, using the patterned regular expressions. The extracted words are expressed in a matrix. The source code of the NLP system can be found at https://github.com/YJK96/Low-Back-Pain-NLP. git, which was accessed on 28 November 2022.

#### *2.4. Statistical Analysis*

The chi-square test was used to compare ratios, and the t-test was used to compare quantitative variables. The performance of the NLP system was measured using recall, precision, accuracy, and the F1 score.


*Number o f true positives* + *Number o f true negatives*

*Number o f true positives* <sup>+</sup> *Number o f true negatives* <sup>+</sup> *Number o f f alse positives* <sup>+</sup> *Number o f f alse netgatives* (1)

The F1 score is defined as <sup>2</sup>×*recall*×*precision prcision*+*recall* and is used to estimate performance. This score is the harmonic mean of precision and recall, and it is considered more useful than accuracy because of the prevalence of class imbalance in text classification [23].

#### **3. Results**

#### *3.1. Dataset Characteristics*

The dataset (*n* = 900) was divided by the training datasets (*n* = 720) and the test datasets (*n* = 120). The median age, standard deviation of patients, and the proportion of male or female patients for each dataset were calculated.

In the dataset evaluation, the training and test dataset was classified by 23 radiologic findings. The prevalence of radiologic findings ranged from 1.4% (disc sequestration) to 56.3% (disc bulge) in the training set and from 0% (disc herniation, facet hypertrophy) to 53.9% (disc bulge) in the test set. In both the training and test sets, disc bulge was the most common finding (56.3% and 53.9%, respectively). Spondylosis was the second-most common finding (39.7% and 42.8%, respectively), followed by listhesis as the third-most common finding (36% and 33.3%, respectively) (Table 2). Between the training and test datasets, 4/23 radiologic findings (facet hypertrophy, lateral recess stenosis, any stenosis, disc sequestration) show a statistically significant difference. (*p*-value < 0.05) (Figure 2).

**Table 2.** Characteristics of training and testing datasets for the development of the natural language processing pipeline.


This table shows the number of diseases diagnosed in 720 training set and 180 test set reports.

#### *3.2. NLP System Performances*

To evaluate the best performances of the NLP system, the NLP model was trained and datasets tested by calculation of the recall, precision, accuracy, and the F1 score. The raw data of the NLP system is available at https://github.com/YJK96/Low-Back-Pain-NLP.git, which was accessed on 28 November 2022.

In the test datasets (*n* = 120) of the rule-based NLP model, all 23 radiologic findings had scores of more than 0.9 for recall, precision, accuracy, and the F1 score. These four scores were 1.0 for 10 radiologic findings—listhesis, annular fissure, disc bulge, disc extrusion, disc protrusion, endplate edema or type 1 Modic change, lateral recess stenosis, Schmorl's node, osteophyte, and any stenosis. The lowest F1 score was 0.9802 for facet hypertrophy.

The highest recall score was 1.0 for 14 radiologic findings—listhesis, fracture, annular fissure, disc bulge, disc extrusion, disc herniation, disc protrusion, endplate edema or type 1 Modic, facet hypertrophy, lateral recess stenosis, spondylolysis, Schmorl's node, osteophyte and any stenosis. The lowest recall score was 0.9840 for central stenosis.

**Figure 2.** Frequency of the radiologic findings between the training and test datasets. (\*\*) indicates a statistically significant difference (*p*-value < 0.05).

The highest precision score was 1.0 for 16 radiologic findings—listhesis, scoliosis, spondylosis, annular fissure, disc bulge, disc extrusion, disc protrusion, endplate edema or type 1 Modic, central stenosis, nerve root displaced/compressed, lateral recess stenosis, Schmorl's node, osteophyte, disc space narrowing, any stenosis, and disc sequestration. The lowest precision score was 0.9611 for facet hypertrophy. The highest accuracy score was 1.0 for 10 radiologic findings—listhesis, annular fissure, disc bulge, disc extrusion, disc protrusion, endplate edema or type 1 Modic, lateral recess stenosis, Schmorl's node, osteophyte, and any stenosis. The lowest accuracy score was 0.9611 for facet hypertrophy (Table 3).

The F1 score for the potentially clinically important 7 radiologic findings is as follows: 1.0 for disc extrusion, 1.0 for endplate edema or Type 1 Modic, 1.0 for any stenosis, 0.9919 for central stenosis, 0.9882 for foraminal stenosis, 0.9969 for nerve root displaced/compressed, and 1.0 for lateral recess stenosis.


**Table 3.** Performance of the natural language processing pipeline in the test dataset (N = 180).

#### **4. Discussion**

The purpose of study was the development of the NLP system for extracting the clinical findings, which are associated with LPB on X-ray, CT and MRI radiologic reports. For NLP system modeling, the training datasets (*n* = 720; 240 X-ray, 240 CT. and 240 MRI) and the test datasets (*n* = 120; 40 X-ray, 40 CT, and 40 MRI) were extracted from 18,640 radiologic reports. Between the training and test datasets, 19/23 radiologic findings show no statistically significant differences.

In this study, we developed a rule-based NLP pipeline for identifying 23 radiologic findings in X-ray, CT, and MRI reports. In detail, regular expression and negation detection were used for modeling. In assessing the performance of this NLP system, there are four parameters for validation—recall, precision, accuracy, and F1 score. This NLP system presented accuracy more than 0.9611. To correct the effects of dataset bias, the NLP system was evaluated using the F1 score. The F1 score is the most important thing among the parameters for NLP system evaluation. The F1 score is an indicator that reflects recall and precision. In this NLP system, all radiologic findings had an F1 score of 0.9802 or higher. Recall means the sensitivity of the statistics. The developed NLP system had a sensitivity of at least 0.984. Compared to the previous rule-based NLP system [16], the sensitivity of the rule-based NLP pipeline was improved.

From the previous study, the potentially clinically important radiologic findings were defined—Disc extrusion, endplate edema or Type 1 Modic, any stenosis, central stenosis, foraminal stenosis, nerve root displaced/compressed, and lateral recess stenosis [16]. In this study, the 7 potentially clinically important radiologic findings had an F1 score of 0.9882 to 1.0. Moreover, the developed NLP system presented a higher sensitivity of potentially clinically important radiological findings than did the previous rule-based NLP system [16]. In this study, a unique NLP system was developed for extracting the clinical findings associated with LBP pain from the radiologic reports of CHA University Bundang CHA Medical Center.

The imaging findings associated with LBP are not explicitly coded in the medical database that is part of the electronic health record. An NLP system automatically identifies such findings from free-text radiology reports, reducing the burden of manual extraction [16]. In this study, a rule-based NLP system was developed and validated for identifying radiologic findings related to LBP in radiology reports. Radiology reports contain a substantial amount of content within electronic health records and can be used for communication and documentation of the imaging. The same radiologic findings can be expressed in different words in radiology reports, and many radiologic findings remain unstructured in medical databases. A trained person can perform the information extraction task, but manual extraction is labor-intensive, costly, and time-consuming. An NLP system can be used to select useful information from unstructured text in electronic medical records, reducing the burden of manual extraction [13].

Numerous NLP systems now exhibit similar accuracy to humans and have already been used for radiologic reports. Knirsch et al. compared an NLP system with specialist review for chest reports [24]. The chest reports of patients with a positive culture and suspected to have tuberculosis were identified. The focus was on identifying six keywords in the chest reports, and a consensus of 89–92% was gained. Chapman et al., developed an NLP system for disease chronicity, certainty, presence, and examination technical quality [25]. The NLP system showed high sensitivity (86–98%) and high specificity (89–93%), but did not exhibit such high results for chronicity (60% and 99%). Jujjavarapu et al. compared different preprocessing and featurization of the NLP pipeline for classifying radiologic findings associated with LBP from X-ray and MRI radiologic reports. In that study, if the NLP pipeline was developed in the same system (e.g., healthcare institution, EMR system), N-grams was the preferred NLP method. If NLP pipelines were developed in multiple systems, document embeddings were considered the best NLP method [17]. Travis Caton et al., investigated the relationship between the radiologic findings of degenerative spinal stenosis on lumbar MRIs (LMRI) and patient characteristics (e.g., age, sex) using rule-based NLP systems. The NLP system identified the patterns of lumbar spine degeneration [19]. Caton et al. investigated lumbar spine MRIs (LMIR) to identify a performance metric for measuring global severity of lumbar degenerative disease (LSDD) [20]. Huhdanpaa et al. developed a rule-based NLP system to identify the patients with Type 1 Modic endplate changes from lumbar MRI reports [18]. Lewandrowski et al. developed deep learning neural network modes to identify features from MRI digital imaging and communications in medicine (DICOM) datasets to produce automatic MRI reports [26]. Galbusera et al. developed BERT-based NLP system to generate annotations for radiographic images from lumbar X-ray reports [21]. In many previous studies, there were no NLP systems for classification or generation annotations from CT radiologic reports.

Some previous studies used clinical notes or electronic health records to assist in clinical diagnosis or treatment decisions for LBP patients. Miotto et al. developed a convolutional neural network model for classifying acute LBP episodes from free-text clinical notes [27]. Walsh et al. developed a support vector machine model for identifying the axial spondyloarthritis (SpA) concepts from electro medical records [28]. Additioanlly, Walsh et al. developed an NLP system for identifying the patients with axial SpA from clinical charts [29]. Zhao et al. developed the NLP algorithms to classify the axial SpA patients from the free text data of the electronic health records [30]. Particularly, the NLP system was applied to identify the electronic medical record or operation notes during surgery or post-operation. Ehresman et al. developed the NLP system for identifying incidental durotomies from intra-operative electronic health records [31]. Karahade et al. developed the NLP system for identifying incidental durotomies from the free-text operation notes [32]. Additionally, a few previous studies identified the complications of spinal surgery from operative notes or electronic medical records. Karhade et al. developed the machine learning algorithms for the prediction of intra-operative vascular injury (VI) and the NLP systems for identifying VI from free-text operative notes [33]. Karhade et al. developed the NLP algorithms for identifying the post-operative wound infection that

requires reoperation after lumbar discectomy from free-text operation notes [34]. Karhade et al. investigated the NLP algorithms for predicting the 90-day unplanned readmission of the lumbar spine fusion patients from free-text notes during the hospitalization period [35]. Dantes et al. developed an NLP system based on IDEAL-X to identify venous thromboembolism from electronic medical records [36]. The free-text type clinical note, operation notes, and electronic health records could be used for developing NLP algorithms for identifying or predicting specific spinal diseases or complications of spinal surgery (Table 4) [12].

Recently, NLP system research on clinical reports written using official language has been conducted globally. Emilien et al. developed the NLP system for learning contextual embeddings from free-text form clinical records in France [15]. Kim Y et al. developed an NLP system for identifying a Korean medical corpus using BERT models [37]. Dahl et al. developed an NLP system for classifying Norwegian pediatric CT radiology reports. A bidirectional recurrent neural network model, a convolutional neural network model, and a support vector machine model were used for training and testing the NLP system [38]. Fink et al. developed an NLP system for identifying the oncologic outcomes from structured oncology reports created in the German language [39].

In previous research, machine learning- or deep learning-based NLP system predicted multiple findings using the same datasets. Moreover, the machine learning-based NLP system was more scalable than the rule-based model. However, machine learning- or deep learning-based NLP systems required larger datasets and a higher set-up cost than rule-based NLP systems [16].

In this study, the NLP system is developed using the rule-based approach. Clinicians determined and set the regular expressions using "re." The various patterns of the free-text form clinical reports could be defined as regular expressions. Because rules are set by the user, and the rule-based method is easier to debug and has higher precision that what can be expected using the machine learning-based method. Several limitations of the previous studies were overcome. First, this NLP system was trained with relatively smaller datasets (*n* = 720) than those used in previous studies. Second, the NLP system showed improved sensitivity (from 0.984 to 1), which has been pointed out as low in previous rule-based NLP systems [16]. Third, the utilization of the NLP system was extended to the CT radiology reports for extracting words associated with LBP.

There are several limitations of this study. First, the radiologic findings were collected from one medical center. The pipeline may not properly process reports with different styles of reporting or expressions. In future works, numerous radiology reports with different styles of reporting should be obtained, as well as datasets from multiple institutions. Second, only the rule-based method was used for developing the NLP model. In future work, other featurization methods (e.g., N-grams, document embeddings) should be used and compared to each other. Despite these limitations, the feasibility of the NLP system to identify lumbar imaging results associated with LBP from sampled radiographic reports was confirmed.

In the future, the NLP system will be used to identify lumbar imaging results among subjects without LBP, to predict whether LBP will become chronic, and to investigate the accuracy of predictions of chronic LBP persistence by investigating these predictions along with other test results, such as hematological tests. Additionally, the transformer models (i.e., BERT model) can be used for further studies. The patient datasets with LBP and without LPB will be learned in the BERT model. Through this method, we are expected to discover the sentences or words with a high correlation with LBP. More keywords, in addition to the 23 radiologic findings used in the current study, will be found.


**Table 4.** Classification of previous studies to develop NLP systems regarding LBP pain, spinal disease, and complications of spinal surgery.

#### **5. Conclusions**

A rule-based NLP system was developed to identify 23 findings associated with LBP from X-ray, CT, and MRI radiology reports sampled from CHA University Bundang CHA Medical Center. This rule-based NLP system presented good performance. In particular, the potentially clinically important seven radiologic findings exhibited high F1 scores.

Through research, an NLP system was developed to identify and extract the necessary terms from the free-text form clinical data used in medical institutions. Through this study, an NLP system for extracting words associated with LBP from not only X-ray and MRI reports, but also CT reports, was developed.

The utilization of this NLP system will be expanded to extract data needed for medical research or data that becomes a marker for diagnosis or treatment from free-text form clinical data of the medical institution. This NLP system has potential for use in clinical diagnosis, treatment decisions, and large-scale research.

In future studies, the criteria of radiologic findings will be subdivided according to a lesion location. Additionally, results of each patient's complete blood count (CBC) test and blood chemistry test should be used for the prediction of chronic LBP persistency, along with the radiologic findings from the rule-based NLP system. Moreover, the correlation between the degree of pain and the radiologic findings from the NLP system will be investigated.

**Author Contributions:** Conceptualization and methodology; H.-W.H. and I.H., writing; Y.K., S.B.K., G.S. (Gyuseon Song) and C.S. (Chanyoung Song); data acquisition; Y.K., S.B.K. and G.S. (Gyuseon Song). All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the Korean Health Technology Research and Development Project, the Ministry for Health and Welfare Affairs (HR16C0002 and HH21C0027), and the Institute of Information and Communications Technology Planning and Evaluation (IITP) grant funded by the Korea government (MSIT) (No. 2019-0-00224, AIM: AI based Next-Generation Security Information Event Management Methodology for Cognitive Intelligence and Secure-Open Framework).

**Institutional Review Board Statement:** This study was approved by the Institutional Review Board of the Bundang CHA Medical Center (IRB No. 2021-04-212).

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are openly available in [repository name e.g., FigShare] at [doi], reference number [reference number].

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


#### MDPI

St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

Academic Open Access Publishing

www.mdpi.com ISBN 978-3-0365-8199-6