Repetitive in vivo manual loading of the spine elicits cellular responses in porcine annuli fibrosi

Back pain and intervertebral disc degeneration are prevalent, costly, and widely treated by manual therapies, yet the underlying causes of these diseases are indeterminate as are the scientific bases for such treatments. The present studies characterize the effects of repetitive in vivo manual loads on porcine intervertebral disc cell metabolism using RNA deep sequencing. A single session of repetitive manual loading applied to the lumbar spine induced both up- and down-regulation of a variety of genes transcribed by cells in the ventral annuli fibrosi. The effect of manual therapy at the level of loading was greater than at a level distant to the applied load. Gene ontology and molecular pathway analyses categorized biological, molecular, and cellular functions influenced by repetitive manual loading, with over-representation of membrane, transmembrane, and pericellular activities. Weighted Gene Co-expression Network Analysis discerned enrichment in genes in pathways of inflammation and skeletogenesis. The present studies support previous findings of intervertebral disc cell mechanotransduction, and are the first to report comprehensively on the repertoire of gene targets influenced by mechanical loads associated with manual therapy interventions. The present study defines the cellular response of repeated, low-amplitude loads on normal healthy annuli fibrosi and lays the foundation for future work defining how healthy and diseased intervertebral discs respond to single or low-frequency manual loads typical of those applied clinically.


Introduction
Persistent low back pain is a global problem responsible for more years-lived-with-disability than any other condition [1,2]. Hence, the health, societal, and economic burdens associated structural integrity of the disc-in responding to loads typical of those in clinical practice when applied repeatedly by in vivo spinal manipulation of a porcine model using a discoverybased, broad-spectrum gene expression analysis. Before it is possible to credibly interpret a "clinical dose" of manual therapy (i.e., the episodic application of short-term loads), it is necessary to define positive and negative controls to frame the interpretation of clinical loading paradigms. The present study seeks to identify any gene candidates of manual therapy-induced mechanotransduction by evaluating differences between the gene expression of IVD cells in positive and negative controls. In this initial study, a porcine model is used to define any changes in gene expression, as indices of cell responsiveness, in healthy normal annuli fibrosi of positive controls (treated with multiple applications of manual therapy over several hours) and sham-treated (non-loaded) negative controls. These initial proof-of-principle studies intend to serve as a foundation for future gene expression analyses of normal and diseased annuli fibrosi loaded using clinical doses of manual therapy of the spine.

Results
The mean (±SD) yield of total RNA was 0.089 ± 0.035 (μg/mg wet weight) for the central ventral wedge of the L 3-4 annuli fibrosi. There was no significant difference in RNA yield among intervertebral discs between Hyperloaded and Control specimens (p = 0.49).

RNA deep sequencing
After quality control (see Methods), 17.7 million raw reads per sample were aligned on average, revealing 13,402 unique expressed transcripts; of these, 11,370 were annotated sufficiently for DAVID analysis. Following filtering for transcripts, unsupervised hierarchical cluster analysis (JMP Genomics) detected two distinct groups that corresponded to Non-loaded Controls and Hyperloaded groups as evidenced clearly on the heatmap of raw data, which is presented with a sorted dendrogram of transcripts of similar expression profiles (Fig 1). Using a cutoff of p<0.01, annuli fibrosi tissues treated with Hyperloading, compared to Non-loaded Controls, were calculated to have 348 mRNAs significantly up-regulated, 430 mRNAs significantly down-regulated, with 10590 mRNAs identified as "background," (i.e., transcribed mRNAs with levels p�0.01). A Principal Component Plot reveals close clustering of Hyperloaded samples implying that any effects of loading duration or frequency are indistinguishable (Fig 2).  Table 1 lists, by fold-change in transcript expression, an assortment of 20 differentially expressed transcripts selected by magnitude of change and relevance to IVD pathobiology. The full lists of up-and down-regulated genes are provided online as S1 and S2 Tables.

RT-qPCR analysis of select RNA sequencing targets
RT-qPCR results are listed as fold-change in mRNA copy number in Non-loaded Control and Hyperloaded samples and are listed in comparison to RNA sequencing fold-change (Table 2). Note the relatively diminished response to Hyperloading in the L 1-2 versus L 3-4 discs.

Functional analyses of RNAseq transcripts
DAVID Functional Annotation analysis enumerated transcript count and the percentage of differentially expressed transcripts, of which the top 20 (of 237) are listed in Table 3 along with  their calculated Enrichment scores. PANTHER Functional Analyses plot the percentage of listed background, up-, and down-Regulated genes for various Biological Processes, Molecular Functions, and Cellular Components (Fig 3) using Fisher's Exact test with the FDR multiple-test correction (n = 348 upregulated; n = 430 down-regulated, and n = 10590 unchanged genes). The numbers of up-regulated and down-regulated genes were highly correlated for each of these analyses (>0.95) and both were highly correlated (>0.92) to the number of background genes (i.e., the number of transcribed genes in each category). Compared to background expression, the relative overand under-representation of differentially expressed transcripts in Gene Ontogeny (GO) classes of various biological and cellular processes is given in Table 4  that strongly associated with the trait Treatment (Non-loading, Hyperloading- Fig 5). Based on calculated enrichments for Gene Ontology, the top 10 ranked transcripts for each module are listed in Table 5. The blue module is enriched with genes that participate principally in inflammatory processes; the turquoise module is enriched with genes that participate principally in skeletogenesis.

Discussion
Spinal diseases, particularly low back pain, are commonly treated with some form of physical intervention. Even if the underlying mechanisms of these physical interventions are unknown and potentially complex, clinical reports suggest such interventions may influence low back pain in some, though not all, people [21]. Hence, it is unsurprising that extensive clinical and basic research has investigated the influence of physical, i.e., mechanical, loads on the spine pathophysiology. While pathology of the spine, including the intervertebral disc, is commonly documented in patients with a history of back pain, the exact source of back pain remains indeterminate, and it is unclear if manual therapies have any influence on IVD metabolism that might promote repair, preserve health, or diminish pain. This initial study defines changes from baseline gene expression in healthy discs in response to manual loads of magnitudes relevant to clinical treatment. As with all experiments, these studies have several assumptions and limitations as well as certain advantages. Although the use of skeletally immature, quadrupedal pigs as a model of human IVD biology is an obvious limitation, the size of the pigs in this study enables both precise sampling of IVD tissues and the relevant application of spinal manipulation by an experienced clinician, which makes this model highly advantageous for the purposes of this study. The goal here was to define a positive control for applying compressive manual loading to normal healthy IVD, which leaves future studies to determine, whether or not, any molecular changes might occur after a typical, clinically applied spinal manipulation in healthy and diseased IVDs. Nevertheless, the mRNA changes documented here indicate that IVD cells have a distinct, near-to-short-term (within 4 hours) stimulation of gene transcription in response to applied, repetitive loads, which is in accord with the studies of MacLean et al. [20] and others. It must also be acknowledged that although an equal magnitude of load was applied to all Hyperloaded pigs, two subgroups received different frequencies of load, which would be an insufficient number to discern any graded effect of load. Nevertheless, in the present study, all pigs were entered as individuals into an unsupervised analysis, which statistically determined that they belonged in two distinct groups that corresponded to the binary variable of load, i.e., hyperloaded versus non-loaded individuals. Principal component analysis supports a distinct treatment effect of loading, but no appreciable differences of load history. Lastly, it is noteworthy that the present report pertains only the ventral segment of annuli fibrosi of healthy spines, and based on previously published work in humans [22,23] and bovines [24], it seems likely that different parts of the intervertebral disc might receive different types and magnitudes of loads and would likely have different cellular responses based on region and health of the disc. While it is well known that the skeleton responds reliably to repeated loads with biological adaptations (e.g., muscle tone and skeletal density) over long-time scales (weeks-to-months), short-term loads activate cellular metabolism as evidenced by changes in nucleic acid transcription [20], which can be detected and evaluated very precisely by RNA sequencing and RT-qPCR. In the present study, it is noteworthy that the yield of RNA extracted from annuli fibrosi is similar to previous reports of other dense connective tissues such as ligament and tendon [25]. And, of the reverse-transcribed mRNA, 6.8% (778 of 11,368 total transcripts) of ventral annulus transcripts were determined (by unsupervised clustering) to be differentially expressed between Non-loaded Controls and Hyperloaded groups. That transcripts were both up-regulated and down-regulated suggests that a systematic bias is unlikely to account for such differential expression. Moreover, the confirmation by RT-qPCR of transcript changes that generally mimic those of RNA sequencing ( Table 2, R 2 = 0.75), and a "dose-responsive" effect of loading (L 3-4 > L 1-2 ; p< 0.05) supports a genuine biological response to manual loading. The discrete quadrants readily visible on the RNA sequencing heatmap (Fig 1) and the segregation on the PCA plot (Fig 2) highlight the distinct differential expression of RNA transcripts between the Hyperloaded and Non-loaded Control groups, supporting an unbiased and unexpectedly clear treatment effect of a single session of repeated manual loading.
The selection of up-regulated transcripts in Table 1 reveal a preponderance of chemokines and enzymes with catabolic actions on extracellular matrix, while the down-regulated transcripts are notable for diminished expression of the (anabolic) structural (extracellular) matrix components aggrecan and collagen XI. If persistent, such changes infer net overall matrix  catabolism, which is opposite to the intended goal of tissue repair and rebuilding, however, initial changes such as these are consistent with tissue inflammation, which necessarily precedes tissue repair. In particular, the strong upregulation of four C-C motif chemokines (CCL2, CCL4, CCL8, CCL24), a collection of specialized secondary mediators of inflammation capable of responding to primary inflammatory mediators such as IL-1beta, implies that annuli cells are preloaded and prepared to respond to mechanical loading with classical inflammation mediators.
The distinct fortification of immune and inflammatory genes identified in the blue module by WGCNA is reinforcing evidence that inflammatory mediators are "first responders" to repetitive mechanical loading, whereas the fortification of genes related to skeletogenesis reinforces evidence for recapitulation of skeletal development as is typical of skeletal and connective tissue responses to injury and inflammation.
Although it is unclear exactly how manual therapy might activate annulus cells, the overrepresentation of membrane and transmembrane transcripts (Tables 3 and 4) is consistent with various models of cellular mechanotransduction that involve outside-in signaling [26] as well as membrane-associated signaling of immune-recognition and inflammation. In particular, there are noteworthy differentially regulated transcripts encoding for integrins and cytoskeleton (e.g., integrin subunit 4), ion pores (e.g., aquaporin 9), and various membrane receptors (e.g., TLR-2) in response to Hyperloading.
Structural molecules notwithstanding, there are signs of anabolic signaling. For example, there is a noteworthy increase in insulin-like growth factor 1 (IGF-1) and a curious concomitant decrease in the IGF-1-antagonist cartilage intermediate layer protein (CILP). Should such changes result in a net increase in the IGF biosynthesis it could be viewed as a hopeful outcome of manual loading as IGF-1 reportedly has a number of "positive" biological effects, including cell proliferation and matrix synthesis, on intervertebral disc cell metabolism [27]. Even while the overall balance appears to favour catabolism over anabolism, it should be recognized that the targets selectively listed in Table 1 are but a small percentage (20/778 = 2.5%) of the total number of significantly changed transcripts, many of which have unknown functions, hence, it seems fair only to broadly conclude that the cellular response is rich and complex.
MacLean and colleagues used RT-PCT to evaluate a small, select set of mRNA changes in rat (caudal segment) annuli fibrosi (dynamically loaded for 1.   [20]. Nevertheless, baseline (time zero) mRNA expression of all these mRNAs was not significantly different compared to endogenous non-loaded controls, but became elevated in all mRNAs except MMP-2 when sampled after 8, 24, or 72 hours after loading [20]. Hence, the methods and system of MacLean et al. were relatively insensitive to changes at baseline, yet clearly defined long-lasting changes in gene expression well after the application of mechanical stimulation. In view of MacLean's findings, the present studies support early initial changes in gene expression that are likely to lead to downstream changes that ultimately could lead to adaptive changes of tissue structure and function. With respect to discogenic pain, the noteworthy decrease in the expression of galanin receptor 3, which binds the nociception-inhibiting neuropeptide galanin [28,29] and reportedly has anti-inflammatory activity [30]. These findings provoke speculation that mechanical load might somehow be involved in intervertebral disc neuroinflammation and nociception, which receives some support with the detection of inflammatory pathways in WGCNA analysis.
As with any large dataset, it is difficult to avoid over-interpreting such a long list of differentially regulated transcripts, so it is important to recognize that these studies are a first attempt to define potential targets and pathways that have changed in response IVD loading. Nevertheless, the present studies demonstrate that repetitive manual loading of the living, intact multielement organ system of the lumbar spine transduces a detectable biological signal, which further reinforces the potential role of mechanobiology in spinal pathobiology, repair, and therapy. Notwithstanding the distinct differential expression of IVD transcripts induced by hyperloading reported here, additional studies are needed to discover whether such changes are robust, repeatable, and lead to functionally significant biological responses to clinical applications of manual therapies.

Materials and methods
The overall experimental design is outlined in Fig 6. Gene expression of annuli fibrosi exposed to in vivo repeated manual loading (as applied in routine practice) or sham manual loading was evaluated using discovery-based RNA deep sequencing, and a selection of several up-regulated and down-regulated gene targets where chosen for validation by RT-qPCR.

Animals
Approval for this experiment was provided by the Animal Care and Use Committee at the University of Alberta (573/07/12). Ten domestic Duroc cross (Large White × Landrace) pigs of mixed sex (6 females; 4 males) were included in these studies. Animals ranged in age from 3-4 mos. old and had body mass of 48 kg ± 4 kg (mean ± SD), a size that facilitated the simulation of manual loading in humans. Animals were cared for according to the guidelines of the Canadian Council of Animal Care under the supervision of a veterinarian. Anesthesia was induced using an initial intramuscular injection of azaperone, glycopyrrolate, and a mixture of ketamine-acepromazine-xlyazine, and was maintained for four hours on a mixture of inhaled isoflurane (~1.5%) and oxygen (1%). At the completion of these acute studies animals were euthanized by intravenous injection of barbiturate (Euthanyl, 150 mg/kg).

Manual spinal therapy
Six animals (four females; two males) were assigned randomly to a "Non-loaded Control" group, and four animals (two females; two males) were assigned to a "Hyperloaded" group. Manual therapy of the spine typical of clinical practice was applied by a trained chiropractor (GK) as described previously [31]. The Non-loaded Controls were anesthetized but did not receive any manual loading. Hyperloaded animals received spinal manipulative therapy applied bilaterally to the L3 vertebra once every 1 or 2 minutes for 4 hours; two pigs received loading once a minute (total manipulations = 240), two pigs were loaded once every two minutes, for four hours (total manipulations = 120). The magnitude of manual spinal loading in pigs was calibrated to 400N, the load measured in the human lumbar spine undergoing clinically simulated manual therapy [32]. Briefly, a thin, flexible 10×10 array of 1 cm 2 pressure sensors (sensitivity 0.15%, 120 Hz sampling rate [Pressure Profile Systems, Los Angeles]) was inserted between the hand of the manual therapist (GK) and the animal, which monitored real-time loads throughout the experiment. For these experiments, the repetitive frequencies and total number of applications used for these pigs greatly exceed that delivered clinically with the express purpose of inducing an unequivocal positive loading response, hence, we termed this treatment "hyperloading."

Tissue sampling and processing
At the completion of these acute studies, animals were euthanized and the lumbar spine was dissected to expose the intervertebral discs. The annuli fibrosi from the central ventral (anterior) sector of the L 3-4 intervertebral discs were dissected and removed (Fig 6). Based on the mid-dorsal application of load, this mid-ventral segment of annulus was chosen for analysis as it was expected to experience the maximum resultant load, thus eliciting the maximal mechanobiological response in the cells of the ventral annuli fibrosi. Annulus segments were cut into small pieces (<30 mg), placed into microfuge tubes, weighed on RNAse-free aluminum foil on a microbalance, snap-frozen in liquid nitrogen, placed immediately into RNAse-free containers, and stored at -80˚C. Annuli fibrosi was also sampled from the L 1-2 IVD based on the assumption that this level experienced less load than the L 3-4 disc at the point of loading. Frozen annuli fibrosi tissue samples were weighed and pulverized into a fine powder using a freezer mill (MikroDismembrator; 2600 rpm for 45 seconds) as described previously [33]. Frozen powdered annulus tissue was immersed immediately in 300 μl RLT buffer containing 1% beta mercaptoethanol for extraction and purification of total RNA using the RNeasy Fibrous Tissue Mini Kit (Qiagen Cat No. 74704). Briefly, RLT-tissue lysate was digested with proteinase K at 55˚C for 10 min, centrifuged 3 min at 10,000 g, and the supernatant transferred and mixed by pipette with a 0.5 volume of absolute ethanol, which was then pipetted over an RNeasy mini-spun column (Qiagen Cat. No. 74104) and centrifuged for 15 sec at >10,000g at 20-25˚C. The column was treated with DNAse I at 22˚C for 15 min, washed twice with RW1 Buffer and centrifuged, then eluted with 25 μL for 5 minutes before centrifugation.

RNA sequence analysis and bioinformatics
Twelve samples from annulus level L 3-4 were used for mRNA (cDNA) sequencing (see below): six samples from six Control animals and six samples from four Hyperloaded animals. The two extra Hyperloaded sample replicates from annulus level L 3-4 were powdered, isolated, and purified independently to assess repeatability.

Generation of cDNA libraries and sequencing
The TruSeq RNA sample Prep v2 LS protocol (Illumina, San Diego, CA, U.S.A.) was used for preparation of mRNA libraries. Messenger RNA was purified from total RNA samples using poly-T oligo-attached magnetic beads followed by mRNA fragmentation for first-and secondstrand cDNA synthesis. The overhangs, which resulted from the fragmentation of doublestranded (ds) cDNA, were repaired to form blunt ends and a single adenosine was added to the 3 0 ends of the blunt fragments to prevent them from ligating to one another during the adapter ligation reaction. Multiple indexing adapters were ligated to the ends of ds cDNA to prepare them for hybridization on a flow cell followed by a PCR amplification step. The libraries were quantified using the qPCR technique and analyzed on a Bioanalyzer 2100 (Agilent Technologies) using a DNA-specific chip. Base calling and demultiplexing of transcriptome sequencing reads were performed using the Consensus Assessment of Sequence and Variance (CASAVA) v 1.6 and Novobarcode software (http://www.novocraft.com/).

Quality control and alignment of reads
Reads were mapped to the porcine genome Sscrofa11.1 (Ensembl, http://www.ensembl.org/) using JMP Genomics 7.1 (SAS, Cary, NC, USA), allowing two mismatches per read. Total counts and transcripts-per-million (TPM) values were generated for all transcripts. Only transcripts detected with at least 10 raw reads in all biological replicates for control and/or treated samples were considered present and included in further analyses.

Data analysis
All data analysis was carried out using TPM values. Quality control, unsupervised clustering, grouped correlation analysis, and Principal Component Analysis (PCA) were performed with JMP Genomics 7.1. Unsupervised cluster analysis is an assumption-free approach, hence, each pig was entered individually (i.e., not as treatment and control groups). One-way analysis of variance ANOVA was used to determine differentially expressed genes. Transcripts displaying a fold-change �2 and a p-value of � 0.01 were considered differentially expressed. Unsupervised hierarchical cluster analysis was done using the minimum variance method [34], which minimizes within-cluster variance and assigns mutually exclusive subsets of transcriptome profiles from all samples. These analyses independently segregated groups into Non-loaded Controls and Hyperloaded samples (Fig 1). Subsequent pathway analyses were used using these two distinct groups. Principal component analysis (PCA) was performed to distinguish treatment groups. PANTHER 9.0 (Protein Analysis Through Evolutionary Relationships) Classification System (http://www.pantherdb.org/) and the Database for Annotation, Visualization and Integrated Discovery (DAVID 6.8; https://david.ncifcrf.gov/) [35] were used to classify differentially expressed genes according to their gene ontology (GO) and to statistically determine over-or under-representation of categories (with Bonferroni correction for multiple testing). All remaining (unchanged) transcripts present in Non-loaded Control and the Hyperloaded samples were entered into the analysis as "background." Resulting data were supplemented with additional information from Entrez Gene (www.ncbi.nlm.nih.gov/ entrez/query.fcgi?db=gene) and from the literature (http://www.ncbi.nlm.nih.gov/ PubMed/). PANTHER was used to identify biological process, molecular function, and cell processes of selected individual genes that were up-and down-regulated in Hyperloaded discs versus Non-loaded Controls.
A systems-level view of genes differentially expressed in Hyperloaded discs was carried out using the WGCNA package in R [36,37]  Briefly, pairwise correlations and a power function are calculated to develop a weighted coexpression network of differentially expressed genes, which segregate into discreet clusters termed 'modules.' Pig genes in enriched modules were converted to human orthologs before they could be evaluated (Enrichment P determined by Fisher's Exact test) for functional significance in Gene Ontology (Bioconductor v. 3.12) pathways.

Quantitative RT-PCR analysis of select RNA sequencing targets
Total RNA was transcribed into cDNA using the Stratagene first-strand reverse transcription kit (Stratagene Cat#200420) according to the manufacturer. PCR primers were designed using Primer-BLAST [38] and qPCR amplification of template cDNA was performed in triplicate in a real-time thermocycler (Bio-Rad iCycler) using Sybr Green detection system (iQ™ SYBR 1 Green Supermix Bio Rad Cat#170-8880). qPCR targets were quantified using the "Fit Point Method" by iCycler Bio-Rad software (2 -ΔΔCT) and normalized to the expression levels of the housekeeping gene glyceraldehyde phosphate dehydrogenase (GAP) mRNA. Reaction specificity was ascertained by performing melt-curve analysis at the end of the amplification protocol, and the efficiency of the PCR reaction was evaluated from a dilution series of template (1:1, 1:5, 1:10, and 1:100) using the R 2 value of the linear regression of Ct values for GAP, CCL8, and CCL2 (respective R 2 values were 0.79, 0.94, and 0.76).