Next Article in Journal
Pathophysiological Involvement of Mast Cells and the Lipid Mediators in Pulmonary Vascular Remodeling
Next Article in Special Issue
Spleen Transcriptome Profiling Reveals Divergent Immune Responses to LPS and Poly (I:C) Challenge in the Yellow Drum (Nibea albiflora)
Previous Article in Journal
Liposomes in Cancer Therapy: How Did We Start and Where Are We Now
Previous Article in Special Issue
Transcriptome and DNA Methylome Analysis of Two Contrasting Rice Genotypes under Salt Stress during Germination
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
IJMS
Volume 24
Issue 7
10.3390/ijms24076611
ijms-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Analysis of RNA Polyadenylation in Healthy and Osteoarthritic Human Articular Cartilage

by
Phaedra Winstanley-Zarach
1,
Gregor Rot
2,3,
Shweta Kuba
1,4,
Aibek Smagul
1,
Mandy J. Peffers
1 and
Simon R. Tew
1,*
1
Centre for Integrated Research into Musculoskeletal Ageing (CIMA), Department of Musculoskeletal and Ageing Science, Institute of Life Course and Medical Sciences, University of Liverpool, Liverpool L7 8TX, UK
2
Institute of Molecular Life Sciences, Winterthurerstrasse 190, 8057 Zurich, Switzerland
3
Swiss Institute of Bioinformatics, Amphipôle, Quartier UNIL-Sorge, 1015 Lausanne, Switzerland
4
School of Health and Life Sciences, National Horizons Centre, Teesside University, Darlington DL1 1HG, UK
*
Author to whom correspondence should be addressed.
Int. J. Mol. Sci. 2023, 24(7), 6611; https://doi.org/10.3390/ijms24076611
Submission received: 13 January 2023 / Revised: 17 March 2023 / Accepted: 27 March 2023 / Published: 1 April 2023
(This article belongs to the Special Issue Molecular Mechanisms of mRNA Transcriptional Regulation)
Browse Figures
Review Reports Versions Notes

Abstract

:
Polyadenylation (polyA) defines the 3′ boundary of a transcript’s genetic information. Its position can vary and alternative polyadenylation (APA) transcripts can exist for a gene. This causes variance in 3′ regulatory domains and can affect coding sequence if intronic events occur. The distribution of polyA sites on articular chondrocyte transcripts has not been studied so we aimed to define their transcriptome-wide location in age-matched healthy and osteoarthritic knee articular cartilage. Total RNA was isolated from frozen tissue samples and analysed using the QuantSeq-Reverse 3′ RNA sequencing approach, where each read runs 3′ to 5′ from within the polyA tail into the transcript and contains a distinct polyA site. Differential expression of transcripts was significant altered between healthy and osteoarthritic samples with enrichment for functionalities that were strongly associated with joint pathology. Subsequent examination of polyA site data allowed us to define the extent of site usage across all the samples. When comparing healthy and osteoarthritic samples, we found that differential use of polyadenylation sites was modest. However, in the genes affected, there was potential for the APA to have functional relevance. We have characterised the polyadenylation landscape of human knee articular chondrocytes and conclude that osteoarthritis does not elicit a widespread change in their polyadenylation site usage. This finding differentiates knee osteoarthritis from pathologies such as cancer where APA is more commonly observed.

1. Introduction

Osteoarthritis is a multifactorial disease in which initiation and progression is driven by a combination of genetic, environmental and mechanical factors. The functionality of the chondrocytes within the articular cartilage of the joint is significantly altered by the late stages of the disease, and this is associated with altered gene-expression patterns that have been examined in a range of earlier transcriptomic studies [1].
The processes regulating gene-expression changes during osteoarthritis are a source of ongoing study. Genetic-association studies have identified numerous polymorphisms that are linked to osteoarthritis incidence [2], which can affect the processes that control transcriptional regulation [3]. Mechanisms of transcriptional regulation of known effectors of articular cartilage degradation, such as extracellular proteases, have also been determined [4]. Transcriptomic analysis of osteoarthritic cartilage in a variety of species has demonstrated that widespread changes in gene expression are induced by the disease [5].
A growing number of studies have implicated post-transcriptional gene-regulation in the disruption of transcript abundance during osteoarthritis. Mouse genetic models have identified significant disruption to skeletal development when regulators of mRNA decay or translation such as microRNAs and RNA binding proteins are disrupted [6,7,8], indicating the importance of these factors in chondrocyte function. Our previous studies have demonstrated that mRNA decay rates are affected in osteoarthritic human chondrocytes, with an increased number of short-lived mRNAs evident in osteoarthritic cells [9].
One aspect of post-transcriptional gene-regulation that has received little attention in cartilage biology is the alternative polyadenylation (APA) of transcripts. It is well established that a high proportion of human genes display evidence of APA, leading to functional and structural heterogeneity in the transcript pool [10]. By using different polyadenylation sites, transcripts from the same gene can be produced which have altered 3′ untranslated region sequence lengths. This can result in the inclusion or exclusion of protein and small RNA interacting regions, generating pools of transcripts that code for the same protein but have different regulatory responses [11]. Functionally, altered 3′UTR length induced by APA influences the decay, localisation and translation of mRNAs. Furthermore, the use of polyadenylation sites upstream of the 3′UTR, often contained within intronic sequences, can lead to transcripts where the protein-coding region is altered resulting in the production of truncated proteins with altered functionality [12,13].
Evidence now suggests that modulation of the transcriptome by APA can be dynamic, is driven by differentiation state and can be tissue-specific [14]. Furthermore, links between APA and pathologies are beginning to emerge. Global shortening of 3′UTRs because of APA has been demonstrated to be a feature of cancer pathology [15,16] whilst impaired APA has been implicated in the development of alcohol-induced fatty liver disease [17]. To date however, there has been no study into how chondrocytes utilise this method of transcript mRNA regulation or whether it is affected by degenerative diseases such as osteoarthritis.
In this study we used the QuantSeq-Reverse 3′ end sequencing approach, which produces reads that run 3′ to 5′ from the polyA tail into the 3′ transcript sequence [18], to examine RNA purified from healthy or osteoarthritic human knee articular chondrocytes. By using age-matched human samples these data allow complementation of existing studies by analysis of how transcript levels change in the disease. Furthermore, these data also allow us to define how chondrocytes polyadenylate their transcripts and how this is affected by disease. We determined that chondrocytes, like many other cell types, make widespread use of APA to generate diversity in their transcript pool. We also found that disease-associated changes in APA are observed at a relatively low level but that there is potential to produce significant alterations to both gene function and regulatory kinetics where it occurs.

2. Results

2.1. Analysis of Transcript Abundance in Control and Osteoarthritis Articular Cartilage Using QuantSeq Reverse

We used the QuantSeq-Reverse approach to create 3′ RNA sequence libraries that provided transcriptome-wide read coverage across human chondrocyte polyadenylation sites. Total RNA isolated directly from snap-frozen tissue was used for this and mean sample age for healthy tissues was 75 years with standard deviation (SD) of ± 6 (n = 6) and for osteoarthritic tissues mean age was 70 years with SD± 5 (n = 6). Analysis with an unpaired t-test indicated no statistical difference between the ages of the groups (p = 0.16). Sequencing resulted in median sample reads of 4.4 × 107 with the highest number of reads in a sample being 6.0 × 107 and the lowest 3.7 × 107 reads. We began our analysis by examining differential expression of steady-state transcript levels between healthy and osteoarthritic samples (Figure 1).
Principal component analysis of the data demonstrated that the first principal component (PC) was responsible for 40% of the experimental variance and the second principal component, 13% of the variance. There was a clear separation between the healthy and osteoarthritic samples across PC1 (Figure 1A). Analysis of differential expression revealed that osteoarthritic tissue had 561 transcripts that were differentially upregulated and 630 transcripts that were differentially downregulated (where log2fc > 1 and adjusted p-value was < 0.05) (Figure 1B and Table S1). There was overlap with OA-regulated genes observed in other studies. For instance, when comparing genes regulated in a study of non-OA vs. OA human tissue [5], we found that 423 (31%) of the genes we identified as being significantly regulated were also regulated in the other dataset with strong correlation in the direction of the fold changes (Figure 1C). Functional annotation enrichment analysis revealed a strong association of regulated genes with osteoarthritis-related characteristics such as extracellular structure/matrix organisation, and skeletal development in addition to associations with cell motility and migration (Figure 1D). Examination of differentially regulated genes using RcisTarget allowed us to identify enriched transcriptional regulatory sequences associated with OA (Figure 1E). This identified that genes with promoter sequences containing SMAD3 and IKZF1 binding motifs were over-represented in genes that were regulated during osteoarthritis. Regulated genes also consisted of a significant number of long non-coding RNAs that were differentially regulated in osteoarthritis (Figure 1F). We noted that very few small regulatory RNAs such as microRNAs (miRNAs) or small nucleolar RNAs (snoRNAs) were represented in the data set. When examining data from across both control and osteoarthritic samples, only 4 miRNAs (miR125b-1, miR320d-1, miR3671, miR7851) and 11 snoRNAs (SNORA11F, SNORA12, SNORA59B, SNORA72, SNORA79B, SNORD104, SNORD13, SNORD13E, SNORD38C, SNORD38-1, SNORD38-2) were detected across more than 50% of the samples analysed. No significant effect of osteoarthritis was observed on the levels of the small RNAs detected with SNORD13 which exhibited a log2 fold decrease of 2.0 in osteoarthritic samples with an adjusted p-value of 0.093, being the only one to exhibit a trend toward regulation by disease.

2.2. Human Articular Chondrocytes Make Widespread Use of Alternative Polyadenylation to Modify Their Transcripts

Making use of the Quantseq-Reverse data property of sequencing 3′–5′ from the polyA tail, we were able to use the dataset to identify chondrocyte polyadenylation site usage. Aligned reads were analysed using expressRNA to identify polyadenylation site usage in each sample. A full list of polyadenylation sites identified by this analysis can be found in Table S2. To understand general chondrocyte polyadenylation site usage, the healthy and osteoarthritic samples were examined together as a complete chondrocyte dataset. We found that 90% of chondrocyte transcripts had <5 polyadenylation sites with 37% having just one isoform (Figure 2A). No transcripts were observed with more than 14 alternatively polyadenylated isoforms in our chondrocyte data set. We performed similar analysis on a human data set obtained from the polyA_DB database version 3.2 (https://exon.apps.wistar.org/polya_db/v3/, accessed on 21 September 2022), which is comprised of pooled polyadenylation site data from the analysis of 107 tissue and cell samples [19]. This polyA_DB data represents a reference of polyadenylation site usage across a wide variety of human cellular transcriptomes. When we compared chondrocyte polyadenylation site usage to that of the reference PolyA_DB data for 10,551 genes common to each dataset, it was evident that chondrocytes utilise a narrow band of the potential polyadenylation sites that have been so far defined in human cells and tissues (Figure 2B). Our chondrocyte data contained 1114 genes that exhibited >5 alternative polyadenylation sites. We performed ontological analysis of these highly alternatively polyadenylated chondrocyte genes to determine whether they were enriched for any specific characteristics. We found that for genes in chondrocytes that utilised this, many polyadenylation sites were strongly enriched for functions involved in mRNA processing and splicing as well as organelle formation and function (Figure 2C).

2.3. Characterisation of Highly Expressed, Multi-Polyadenylated Genes in Human Articular Chondrocytes

To examine the implication of alternative polyadenylation events in chondrocytes, we examined the ten most strongly expressed, multi polyadenylated genes across the healthy and osteoarthritic chondrocyte data (Table 1). We found that these genes were affected in varying ways based upon the polyadenylation site usage identified. Some genes exhibited simple truncations of the 3′UTR when more proximal sites were used (DCN, SPARC, PRELP, CILP, EGR1). Some of these genes (DCN, PRELP, CILP) did not always exhibit a distal polyadenylation site coincident with those present in the genomic annotation, indicating that in chondrocytes, all of their transcripts are truncated relative to these canonical sites. Other genes exhibited alternative polyadenylation that led to significant changes in mRNA structure. For instance, the chondrocyte data showed that two sites are used for CFH, the more proximal of which is a well characterised alternatively spliced variant that leads to a shorter protein-coding sequence. Less well characterised is the chondrocyte site usage for the CLU gene. Here, six sites were identified, three of which were close to existing annotated sites and were utilised in the majority of samples submitted to polyA_DB. Two of the other three sites were identified in a small proportion (<10%) of the samples used by polyA_DB. These three sites occur at the 3′ end of CLU exon 3 and represent around a quarter of all CLU reads in chondrocytes.

2.4. Osteoarthritis Does Not Result in Widespread Changes in Chondrocyte Polyadenylation Patterns

We used the expressRNA tool to determine the differential usage of polyadenylation sites between healthy and osteoarthritic cartilage samples. We found that disease-specific polyadenylation is not a widespread process with just 20 genes identified as being differentially regulated in this way in the diseased samples (where adjusted p-values for both distal and proximal sites relaxed to <0.1). These genes are presented in Table 2 with details of the alternative polyadenylation sites that were affected between healthy and osteoarthritic samples.
As expected, many of these APA events result in a truncated 3′UTR (KMT2A, SMAD5, SLC25A37, TMOD1) however there are examples of other events which would lead to more substantial changes to transcript structure. One of the most strongly regulated polyadenylation events was found for the lncRNA NEAT1. Healthy chondrocytes primarily use a distal NEAT1 polyadenylation site, encoded at chr11:65,426,254, with a lower-level usage of sites 110 bp and 549 bp upstream of this.
Figure 3 presents the read data for three examples of genes where alternative polyadenylation was observed in osteoarthritis. Figure 3A shows gene structure diagrams for each of these genes, OSMR, KMT2A and ANOS1, with differential polyadenylation site location indicated with arrows. For the first example, we found that the use of a previously identified intronic polyadenylation site in the OSMR gene [20] was significantly downregulated in osteoarthritis (Figure 3B) There was also evidence of use of a further OSMR intronic polyadenylation site, 1.3 kb further downstream which was also diminished in use in osteoarthritic cells. For the KMT2A gene, osteoarthritis was associated with the reduced use of a proximal polyadenylation site in the 3′UTR (Figure 3C). Polyadenylation associated with the full-length mRNA for the ANOS1 gene was only observed in osteoarthritic chondrocytes but there was consistent use of a site within intron 11 by control and osteoarthritic cells in all samples (Figure 3D).

2.5. Intronic Polyadenylation of OSMR mRNA Leads to Altered Transcript Decay Rate and can Be Influenced by Inflammatory Cytokine Stimulation

Focussing firstly on examples of genes that exhibited intronic polyadenylation we were able to produce qPCR assays to assess expression of variants of the OSMR and ANOS1 transcripts. We examined read counts for proximal and distal polyadenylation sites for these two genes and found that osteoarthritis led to a reduced use of the proximal, intronic site for OSMR whilst inducing the use of a distal site for ANOS1 (Figure 4A). We next examined the expression and rates of mRNA decay of the isoforms of these two transcripts in osteoarthritic human chondrocytes (Figure 4B).
For the OSMR gene, the use of the intronic polyadenylation site, which results in a truncated protein-coding domain, also led to significant change in post-transcriptional regulation compared to the full-length transcript. The shorter transcript was less stable with a half-life of 9 h whereas the full-length transcript decayed slowly (half-life > 24 h). Interestingly, despite originating from osteoarthritic cartilage, the cultured chondrocytes expressed higher levels of the proximal isoform. For ANOS1, each isoform was very stable (>24 h half-life) and no difference was observed in their decay rates (Figure 4B). We examined whether the inflammatory cytokine IL-1ß, which is known to strongly regulate many osteoarthritis-associated genes, affected the use of the short and long isoforms of each gene in HACs (Figure 4C). We found that there was a significant reduction in the ratio of short to long isoform usage for OSMR over a 6 h exposure to IL-1ß, consistent with the shift in ratio observed between healthy and osteoarthritic tissue. The ratio of ANOS1 isoforms was not significantly affected by the IL-1ß treatment.

2.6. mRNA Decay Is Affected by Polyadenylation Status of the KMT2A Transcript

Several transcripts that were differentially polyadenylated in osteoarthritis were characterised by differences in length of the 3′UTR. To analyse these isoforms with qPCR we were only able to specifically identify longer, distally polyadenylated isoforms because they contain a specific region of 3′UTR sequence. We therefore intended to use ratios of “distal” to “proximal and distal” site usage to allow us to differentiate between APA isoforms and allow us to determine how they affected transcript decay rates.
For many targets, designing discriminatory assays with satisfactory amplification efficiencies was challenging. We were, however, able to identify primer pairs with comparable amplification efficiencies that could be used to analyse of the KMT2A and SLC25A37 isoforms. Each of these genes exhibited significantly reduced use of a proximal polyadenylation site in osteoarthritic samples (Figure 5A). Using qPCR, we found it difficult to differentiate between expression levels of proximal and distal isoforms for the genes KMT2A and SLC25A37 in both HAC and SW1353 cells (Figure 5B). Because these were genes that exhibited lower proximal site usage in OA this indicated that in cultured human chondrocytes the predominant isoform uses the distal polyadenylation site. Unsurprisingly, the results of actinomycin D chase experiments in chondrocytes led to similar results with each primer pair, as the distal isoform was likely to be predominately detected. We still wanted to understand how the APA of these genes might affect mRNA decay so we screened a range of cell types of different tissue origins in an attempt to identify a model system where a higher ratio of short and long isoforms were expressed (Figure 5C). For SLC25A37, we found that a ratio of 1 in qPCR products was generally consistent across all of the cell types we analysed, again indicating mostly proximal site usage. However, for KMT2A we found that two of the cell lines screened, AC10 and MDA-MB-231, exhibited short to long isoform ratios of 2 and above indicating that they were also expressing the proximally polyadenylated form. We also measured APA usage for OSMR and ANOS1 in the cell lines to determine whether there was any common regulation of the osteoarthritic APA sites across different cell types. This did not seem to be the case: ANOS1 chiefly exhibited preponderance to proximal site usage across all of the lines. In contrast, OSMR site usage varied between cell types, with A549 and NTMS cells making less use of the intronic proximal site compared to the other four cell lines. Based upon our cell-line data, we chose to use the AC10 cells to examine whether APA affected mRNA decay rates for KMT2A (Figure 5D). Using our proximal and distal qPCR assays we found that the product of the distal, long-isoform-specific PCR was more stable than the product of the proximal PCR which detects both long and short isoforms. We can infer from this that the shorter KMT2A transcript is less stable than the longer one.

3. Discussion

This is the first study to examine in detail the use of polyadenylation sites by articular chondrocytes and the effect of osteoarthritis on this process. It has confirmed that chondrocytes make widespread use of APA to create heterogeneity in their transcriptomes. However, change in the pattern of APA is not a widespread feature of the transcriptomes of osteoarthritic chondrocytes.
QuantSeq-Reverse leads to 3′ end reads only and is particularly well suited to studying polyadenylation patterns because each read sequences from the polyA tail into the 3′UTR of a transcript. This approach makes the technique less efficient than other 3′ sequencing approaches for transcript expression profiling but does not preclude such analyses from being conducted. Encouragingly, we have found that our dataset allowed us to identify many transcripts with altered abundance between osteoarthritic and age-matched normal cartilage tissue. Around a third of the highly regulated genes that we have observed have been similarly identified in another transcriptional profiling study examining the disease [5].
In line with the nature of osteoarthritis as a disease, ontology analysis of differentially regulated genes showed enrichment for processes associated with the skeleton, joints and connective tissue; areas that are all known to be significantly affected by osteoarthritis. Our use of age-matched human control and diseased samples allowed us to build on previous transcriptomic studies. Currently, our study and that of [5] are the only transcriptomics studies that make use of age-matched healthy control tissue for the study of human osteoarthritis. Other studies have examined transcriptome or epigenome alterations between intact and fibrillated tissue from osteoarthritic joints [21,22] and have observed differences in expression but our data and that described by Soul and colleagues in [5] indicate that this comparison may not be fully reflective of the differences between diseased and healthy tissue. This approach may have value, because control-tissue choice represents a continual obstacle in the study of human osteoarthritis, where appropriately matched, non-diseased human cartilage can be challenging to obtain. However, an acceptance of potential limitations of different approaches is required. Our healthy samples were sourced from a commercial supplier, whilst we collected osteoarthritic tissue from a local hospital following knee arthroplasty procedures. PCA analysis of the transcriptome data showed that there was significantly less variance in the osteoarthritic samples compared with the healthy, a finding that we have observed in other RNA-sequencing studies [23]. This may reflect a greater homogeny of phenotype among late-stage osteoarthritic samples, in line with a growing opinion that differing age related mechanisms can converge into a common end stage form of osteoarthritis [24].
A characteristic property of the QuantSeq-Reverse technique is that sequences originate within the polyA tail and proceed 3′–5′ into the transcript, meaning that every read specifically identifies a transcript polyadenylation site. This has allowed for fresh insights into transcript regulation in chondrocytes at the post-transcriptional level by identifying patterns of chondrocyte polyadenylation site usage as well as determining the extent of alternative polyadenylation that occurred between age-matched healthy and osteoarthritic cartilage samples. We confirmed that chondrocytes make widespread use of alternative polyadenylation to enhance the diversity of their transcript pool although the number of sites used is substantially lower than those identified across multiple human tissue and cell types in the polyA_DB database [25]. It is therefore likely that the chondrocyte employs a specific pattern of polyadenylation usage, some aspects of which contribute to phenotype-specific gene functionality. An example of this would be the clusterin gene where a quarter of the chondrocytes’ reads indicate relatively high use of polyadenylation sites towards the 3′ end of exon 3, which have been previously identified in this region but are not commonly observed in other cell types. Exon 3 terminating transcripts have been annotated as CLU-204 and CLU-214 and are predicted to code a truncated protein, although they would be polyadenylated directly on the protein-coding region with no 3′UTR. The polyadenylation sites identified in the chondrocytes would actually lead to transcripts between 20 to 83 bp shorter again than those annotations. It is currently not understood what the function of these truncated clusterin transcripts is or whether they are even translated.
It was instructive that chondrocytes from osteoarthritic cartilage do not exhibit widespread altered usage of polyadenylation sites, indicating that a systematic change in the processes controlling APA is not a significant feature of the disease. This distinguishes the processes occurring in osteoarthritis from other pathologies. For instance, many cancers exhibit global regulation of alternative polyadenylation events [15,16]. Despite osteoarthritis lacking these widespread changes, a small number of alternatively polyadenylated transcripts emerged from our analysis and potentially have a capacity to contribute to the disease mechanism. One of the most interesting is the oncostatin-M receptor OSMR. The use of an intronic polyadenylation site has already been demonstrated for this gene, which results in a truncated, secreted form of the receptor that can act as an oncostatin-M antagonist [20]. We have shown here that the use of the intronic site is decreased in osteoarthritic cartilage, which is likely to lead to reduced levels of this antagonistic protein isoform and increased oncostatin-M signalling in osteoarthritic tissue. It was interesting to note that stimulation of chondrocytes with IL-1ß led to a decrease in the ratio of polyadenylation usage at the intronic site versus the full-length form, indicating that inflammatory conditions are able to promote a similar reduction in levels of the oncostatin-M antagonistic form. Inflammatory cytokine stimulation can promote many of the catabolic processes that lead to extracellular matrix destruction in osteoarthritis and oncostatin-M has recently been associated with inflammatory subtypes of osteoarthritis [26]. Understanding how loss of oncostatin-M antagonism through altered use of the regulated intronic polyadenylation sites that we have observed might contribute to disease state warrants further investigation.
The gene KMT2A encodes a lysine methyl transferase, and its expression was downregulated overall in osteoarthritic tissue. It has previously been linked to a potential causal role in altered cartilage thickness through genome-wide association studies [27]. It also exhibited an almost complete absence of the use of a proximal polyadenylation site within its 3′UTR in osteoarthritis, which could influence the dynamics of post-transcriptional processing. Because our chondrocyte primary cell models are from osteoarthritic tissue, we were unable to identify whether this alteration in 3′UTR length affected mRNA decay because it appears that the cells only express the longer version. However, following a screen of a variety of cell lines we found that AC10 cells appear to produce both long and shorter isoforms of the KMT2A transcript and that shortening of its 3′UTR leads to a less stable mRNA. This polyadenylation event may therefore play an important role in maintaining KMT2A levels in healthy cartilage and it would be useful to elucidate the mechanisms that lead to its use being suppressed in osteoarthritis.
Long noncoding RNAs have been increasingly studied in relation to osteoarthritis recently and DeSeq2 analysis of the QuantSeq-Reverse data has identified 18 significantly regulated lncRNAs that had not previously been associated with joint pathology. Furthermore, an alternative polyadenylation event was identified for the lncRNA NEAT1, which has been extensively studied in a variety of tissues and has been found to function as microRNA sponge chondrocytes [28]. The functional significance of this 110 bp difference is not understood although the sequence in this region is more strongly conserved than the rest of the 3′ half of most NEAT1 transcript isoforms. A limitation of the QuantSeq-Reverse approach appears to have been a poor coverage of small non-coding RNAs such as microRNAs and snoRNAs. We did not specifically prepare our RNA samples to enrich for small RNAs, which may explain this. These molecules are known to perform important regulatory roles in chondrocytes and can be regulated in osteoarthritis. They are polyadenylated and so our inability to fully characterise them in our samples represents a limitation of the study. Interestingly, among the small RNAs that we were able to detect were potential biomarkers for osteoarthritis, such as miR-320d [29] and SNORD13 [30].
Our study has been able to identify only a low level of APA caused by osteoarthritis but is subject to some other limitations. Firstly, the number of replicates did not allow us to stratify patients by comorbidities and meant that we limited our study to samples from males, meaning that we have not been able to define any gender-specific effects of APA. Therefore, we may have missed some polyadenylation events that are sensitive to these parameters. We may have also not have been able to effectively quantify genes that were expressed at low levels, which may have led to an underestimation of significant APA in osteoarthritis. Further work is also required to better understand the significance of APA events. We have some evidence that mRNA decay can be regulated by APA for KMT2A but have not looked at other parameters such as translation efficiency and subcellular localisation of the transcripts in this study.
Overall, this study has demonstrated that altered polyadenylation site usage is not a widespread characteristic of osteoarthritis, but that a small number of genes are affected, and the implications of the changes could have significant potential to affect their functions. Further understanding of how to control polyadenylation in chondrocytes across the whole transcriptome is therefore unlikely to be required to combat osteoarthritis but influencing specific polyadenylation events, the loss of intronic polyadenylation of OSMR for instance, could still provide important targets for the modulation of this disease.

4. Materials and Methods

4.1. Age-Matched Healthy and Osteoarthritic Articular Cartilage

Osteoarthritic human articular cartilage tissue was obtained from the knee joints of patients undergoing total knee arthroplasty (Clatterbridge Hospital, Wirral, UK) due to late-stage disease. This material was obtained with informed consent and following national research ethics committee approval (IRAS ID 242434). Age-matched, control human articular cartilage tissue was obtained commercially (Articular Engineering, Northbrook, IL, USA). All control samples had been obtained post-mortem (within 24 h of death) from individuals with no record of joint disease whose joints were macroscopically healthy on visual inspection. Six cartilage samples from each group were analysed and all were obtained from male donors. All tissue was harvested from load-bearing regions of the femoral condyle. Samples were snap-frozen and 10 µm thick full-depth sections of the tissue were cut using a cryostat and collected prior to RNA isolation.

4.2. QuantSeq-Reverse Analysis

Total RNA was isolated from tissue cryosections, to ensure consistent use of full depth tissue, using a Trizol extraction (Thermo Fisher Scientific, Horsham, UK) with subsequent purification using the RNeasy spin column (Qiagen, Manchester, UK) purification procedure including on-column DNase digestion [31]. QuantSeq-Reverse library preparation and NextSeq 75 cycle High Output sequencing with V2chemistry, was performed by Lexogen (Vienna, Austria) [18]. For gene-expression analysis, the in-house Lexogen pipeline was used, consisting of read-trimming, alignment to GRCh38 using STAR aligner [32] and then differential gene expression determined using DESeq2 [33]. Correction for multiple hypothesis testing was performed via the Benjamini–Hochberg procedure [34]. Differentially expressed genes, defined as log2 fold change ≥ 1 and ≤−1 and adjusted p-values (≤0.05), were used to distinguish between downregulated and upregulated genes. Raw FastQ files from QuantSeq-Reverse analysis have been deposited in ArrayExpress (https://www.ebi.ac.uk/arrayexpress/) on 1 October 2022 with the accession number E-MTAB-12184.

4.3. Alternative Polyadenylayion Analysis

We constructed a reference polyA database by considering genome alignments of sequence reads from all twelve experiments. To account for internal-priming events (reads mapping to genomic polyA sequences) we filtered out detected polyA sites without an upstream polyadenylation signal (PAS). After obtaining read coverage for each detected polyA site separately for each experiment, we ran DEXSeq analysis comparing six healthy vs. six osteoarthritis samples to identify alternatively polyadenylated (APA) genes. Further polyadenylation analysis details and methodology can be explored at expressRNA [35].

4.4. Gene Ontology and Pathway Analysis

Enriched KEGG [36] and Reactome [37] pathways were identified using the clusterProfiler package [38] and visualised using the pathview package [39]. For gene-set enrichment analysis (GSEA), a pre-ordered list of gene-named log2 fold change values, and for over-representation analysis (ORA), an input list of differentially expressed genes was used to obtain enriched pathways. Upregulated pathways were defined by a normalised enrichment score (NES) > 0 and the downregulated pathways were defined by an NES < 0. Pathways with BH adjusted p-value ≤ 0.05 were chosen as significantly enriched. RcisTarget v1.9.1 [40] was applied on differentially expressed genes to identify overrepresented potential transcription factor regulation. We examined up to 500 bp upstream of the transcription start site, and TF_highConf (TFs, annotated to the motif according to ‘motifAnnot_highConfCat’ of Rcistarget) were considered associated with TFs of interest.

4.5. Human Articular Chondrocyte Isolation

Osteoarthritic tissue samples were obtained as described above. Chondrocytes were isolated from finely diced cartilage tissue by digestion in 0.08% collagenase type II (Worthington Biochemical, Lakewood, NJ, USA) in growth medium (Dulbecco’s modified Eagle’s medium, containing 10% foetal bovine serum, 100 units/mL penicillin and 100 units/mL streptomycin, all from Thermo Fisher Scientific) overnight at 37 °C with agitation. Isolated cells were washed three times in growth medium without collagenase, plated at 20,000 cells/cm2, and grown for 2 passages at 37 °C in 5% CO2 with 1:2 split ratios before use in experiments.

4.6. Cell Culture

For analysis of the ratios of alternatively polyadenylated transcripts the following cell types were screened: SW1353 chondrosarcoma, AC10 ventricular cardiomyocytes, MG63 osteosarcoma, NTM5 trabecular network cells, A549 lung carcinoma and MDA-MB231 breast adenocarcinoma cells. All cell lines were cultured in growth medium at 37 °C in 5% CO2. For mRNA decay analysis, confluent cell monolayer cultures of the cell lines of interest were cultured in standard media supplemented with 1 ng/mL actinomycin for 0, 2, 4, 6 and 8 h before lysis using TRIzol reagent. Purification of cultured cell total RNA was performed by chloroform-based phase separation and subsequent precipitation with isopropanol. RNA was redissolved in RNAse free H2O and stored at −80 °C until required [9]. Cytokine stimulation of human articular chondrocytes was performed on cells at confluence by adding 10 ng/mL interleukin-1β (IL-1β, Peprotech, London, UK) to standard media and stimulating cells for up to 6 h.

4.7. qRT-PCR Analysis

cDNA was synthesized using 1μg total RNA with 200 units Moloney murine leukemia virus reverse transcriptase (Promega, Southampton, UK), primed with 2 μM oligo dT (Thermo Fisher Scientific). qPCR was performed on a Roche Lightcycler 96 using Takyon ROX SYBR MasterMix (Eurogentec, Seraing, Belgium) using transcript-specific primer pairs at a concentration of 300 nM per primer. Primer efficiency for each assay was established as previously described and assays were only used where efficiencies were within 10% of those of the reference gene [41]. Primer sequences for each assay can be found in Table 3. The 2−ΔCt method was used to calculate gene-expression levels [42] with RPS13 used as the reference expression gene [43].

4.8. Statistical Analysis

Where appropriate, analysis of cell culture experiments was conducted using one-way analysis of variance as described in figure legends and a p-value < 0.05 was considered significant. Statistical methodology underpinning DeSeq2 and expressRNA analysis of transcriptomic datasets has been described previously [33,35].

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ijms24076611/s1.

Author Contributions

Conceptualization, S.R.T. and P.W.-Z.; data curation, formal analysis, validation, and visualization, S.R.T., P.W.-Z., G.R., S.K. and M.J.P.; investigation P.W.-Z. and A.S.; writing—original draft preparation, S.R.T.; writing—review and editing, P.W.-Z., G.R., S.K. and M.J.P.; All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Medical Research Council (MRC) and Versus Arthritis as part of the Medical Research Council Versus Arthritis Centre for Integrated Research into Musculoskeletal Ageing (CIMA) [MR/R502182/1]. The MRC Versus Arthritis Centre for Integrated Research into Musculoskeletal Ageing is a collaboration between the Universities of Liverpool, Sheffield and Newcastle. It was also supported by the Medical Research Council through MRC project grant MR/N011333/1.

Institutional Review Board Statement

Human tissue was obtained with informed consent and following national research ethics committee approval (IRAS ID 242434).

Informed Consent Statement

Informed consent was obtained from all subjects involved in the study.

Data Availability Statement

Raw FastQ files from QuantSeq-Reverse analysis have been deposited in ArrayExpress (https://www.ebi.ac.uk/arrayexpress/) on 1 October 2022 with the accession number E-MTAB-12184.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

References

  1. Katsoula, G.; Kreitmaier, P.; Zeggini, E. Insights into the molecular landscape of osteoarthritis in human tissues. Curr. Opin. Rheumatol. 2022, 34, 79–90. [Google Scholar] [CrossRef] [PubMed]
  2. Boer, C.G.; Hatzikotoulas, K.; Southam, L.; Stefansdottir, L.; Zhang, Y.; Coutinho de Almeida, R.; Wu, T.T.; Zheng, J.; Hartley, A.; Teder-Laving, M.; et al. Deciphering osteoarthritis genetics across 826,690 individuals from 9 populations. Cell 2021, 184, 6003–6005. [Google Scholar] [CrossRef] [PubMed]
  3. Reynard, L.N.; Bui, C.; Canty-Laird, E.G.; Young, D.A.; Loughlin, J. Expression of the osteoarthritis-associated gene GDF5 is modulated epigenetically by DNA methylation. Hum. Mol. Genet. 2011, 20, 3450–3460. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Hirata, M.; Kugimiya, F.; Fukai, A.; Saito, T.; Yano, F.; Ikeda, T.; Mabuchi, A.; Sapkota, B.R.; Akune, T.; Nishida, N.; et al. C/EBPbeta and RUNX2 cooperate to degrade cartilage with MMP-13 as the target and HIF-2alpha as the inducer in chondrocytes. Hum. Mol. Genet. 2012, 21, 1111–1123. [Google Scholar] [CrossRef]
  5. Soul, J.; Dunn, S.L.; Anand, S.; Serracino-Inglott, F.; Schwartz, J.M.; Boot-Handford, R.P.; Hardingham, T.E. Stratification of knee osteoarthritis: Two major patient subgroups identified by genome-wide expression analysis of articular cartilage. Ann. Rheum. Dis. 2018, 77, 423. [Google Scholar] [CrossRef] [Green Version]
  6. Katsanou, V.; Milatos, S.; Yiakouvaki, A.; Sgantzis, N.; Kotsoni, A.; Alexiou, M.; Harokopos, V.; Aidinis, V.; Hemberger, M.; Kontoyiannis, D.L. The RNA-binding protein Elavl1/HuR is essential for placental branching morphogenesis and embryonic development. Mol. Cell Biol. 2009, 29, 2762–2776. [Google Scholar] [CrossRef] [Green Version]
  7. Kobayashi, T.; Lu, J.; Cobb, B.S.; Rodda, S.J.; McMahon, A.P.; Schipani, E.; Merkenschlager, M.; Kronenberg, H.M. Dicer-dependent pathways regulate chondrocyte proliferation and differentiation. Proc. Natl. Acad. Sci. USA 2008, 105, 1949–1954. [Google Scholar] [CrossRef] [Green Version]
  8. Son, Y.O.; Kim, H.E.; Choi, W.S.; Chun, C.H.; Chun, J.S. RNA-binding protein ZFP36L1 regulates osteoarthritis by modulating members of the heat shock protein 70 family. Nat. Commun. 2019, 10, 77. [Google Scholar] [CrossRef] [Green Version]
  9. Tew, S.R.; McDermott, B.T.; Fentem, R.B.; Peffers, M.J.; Clegg, P.D. Transcriptome-wide analysis of messenger RNA decay in normal and osteoarthritic human articular chondrocytes. Arthritis Rheumatol. 2014, 66, 3052–3061. [Google Scholar] [CrossRef] [Green Version]
  10. Derti, A.; Garrett-Engele, P.; Macisaac, K.D.; Stevens, R.C.; Sriram, S.; Chen, R.; Rohl, C.A.; Johnson, J.M.; Babak, T. A quantitative atlas of polyadenylation in five mammals. Genome Res. 2012, 22, 1173–1183. [Google Scholar] [CrossRef] [Green Version]
  11. Hoffman, Y.; Bublik, D.R.; Ugalde, A.P.; Elkon, R.; Biniashvili, T.; Agami, R.; Oren, M.; Pilpel, Y. 3’UTR Shortening Potentiates MicroRNA-Based Repression of Pro-differentiation Genes in Proliferating Human Cells. PLoS Genet. 2016, 12, e1005879. [Google Scholar] [CrossRef] [PubMed]
  12. Davis, M.J.; Hanson, K.A.; Clark, F.; Fink, J.L.; Zhang, F.; Kasukawa, T.; Kai, C.; Kawai, J.; Carninci, P.; Hayashizaki, Y.; et al. Differential use of signal peptides and membrane domains is a common occurrence in the protein output of transcriptional units. PLoS Genet. 2006, 2, e46. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  13. Vorlova, S.; Rocco, G.; Lefave, C.V.; Jodelka, F.M.; Hess, K.; Hastings, M.L.; Henke, E.; Cartegni, L. Induction of antagonistic soluble decoy receptor tyrosine kinases by intronic polyA activation. Mol. Cell 2011, 43, 927–939. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Zhou, R.; Xiao, X.; He, P.; Zhao, Y.; Xu, M.; Zheng, X.; Yang, R.; Chen, S.; Zhou, L.; Zhang, D.; et al. SCAPE: A mixture model revealing single-cell polyadenylation diversity and cellular dynamics during cell differentiation and reprogramming. Nucleic. Acids Res. 2022, 50, e66. [Google Scholar] [CrossRef]
  15. Mayr, C.; Bartel, D.P. Widespread shortening of 3’UTRs by alternative cleavage and polyadenylation activates oncogenes in cancer cells. Cell 2009, 138, 673–684. [Google Scholar] [CrossRef] [Green Version]
  16. Xia, Z.; Donehower, L.A.; Cooper, T.A.; Neilson, J.R.; Wheeler, D.A.; Wagner, E.J.; Li, W. Dynamic analyses of alternative polyadenylation from RNA-seq reveal a 3’-UTR landscape across seven tumour types. Nat. Commun. 2014, 5, 5274. [Google Scholar] [CrossRef] [Green Version]
  17. Jobbins, A.M.; Haberman, N.; Artigas, N.; Amourda, C.; Paterson, H.A.B.; Yu, S.; Blackford, S.J.I.; Montoya, A.; Dore, M.; Wang, Y.F.; et al. Dysregulated RNA polyadenylation contributes to metabolic impairment in non-alcoholic fatty liver disease. Nucleic. Acids Res. 2022, 50, 3379–3393. [Google Scholar] [CrossRef]
  18. Moll, P.; Ante, M.; Seitz, A.; Reda, T. QuantSeq 3′ mRNA sequencing for RNA quantification. Nat. Methods 2014, 11, i. [Google Scholar] [CrossRef]
  19. Wang, R.; Zheng, D.; Yehia, G.; Tian, B. A compendium of conserved cleavage and polyadenylation events in mammalian genes. Genome Res. 2018, 28, 1427–1441. [Google Scholar] [CrossRef] [Green Version]
  20. Diveu, C.; Venereau, E.; Froger, J.; Ravon, E.; Grimaud, L.; Rousseau, F.; Chevalier, S.; Gascan, H. Molecular and functional characterization of a soluble form of oncostatin M/interleukin-31 shared receptor. J. Biol. Chem. 2006, 281, 36673–36682. [Google Scholar] [CrossRef] [Green Version]
  21. Liu, Y.; Chang, J.C.; Hon, C.C.; Fukui, N.; Tanaka, N.; Zhang, Z.; Lee, M.T.M.; Minoda, A. Chromatin accessibility landscape of articular knee cartilage reveals aberrant enhancer regulation in osteoarthritis. Sci. Rep. 2018, 8, 15499. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  22. Zhang, Y.; Fukui, N.; Yahata, M.; Katsuragawa, Y.; Tashiro, T.; Ikegawa, S.; Lee, M.T. Identification of DNA methylation changes associated with disease progression in subchondral bone with site-matched cartilage in knee osteoarthritis. Sci. Rep. 2016, 6, 34460. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  23. Peffers, M.J.; Chabronova, A.; Balaskas, P.; Fang, Y.; Dyer, P.; Cremers, A.; Emans, P.J.; Feczko, P.Z.; Caron, M.M.; Welting, T.J.M. SnoRNA signatures in cartilage ageing and osteoarthritis. Sci. Rep. 2020, 10, 10641. [Google Scholar] [CrossRef]
  24. Dell’Isola, A.; Steultjens, M. Classification of patients with knee osteoarthritis in clinical phenotypes: Data from the osteoarthritis initiative. PLoS ONE 2018, 13, e0191045. [Google Scholar] [CrossRef]
  25. Wang, R.; Nambiar, R.; Zheng, D.; Tian, B. PolyA_DB 3 catalogs cleavage and polyadenylation sites identified by deep sequencing in multiple genomes. Nucleic. Acids Res. 2018, 46, D315–D319. [Google Scholar] [CrossRef]
  26. Garcia, J.P.; Utomo, L.; Rudnik-Jansen, I.; Du, J.; Zuithoff, N.P.A.; Krouwels, A.; van Osch, G.J.V.M.; Creemers, L.B. Association between Oncostatin M Expression and Inflammatory Phenotype in Experimental Arthritis Models and Osteoarthritis Patients. Cells 2021, 10, 508. [Google Scholar] [CrossRef]
  27. Castano-Betancourt, M.C.; Evans, D.S.; Ramos, Y.F.; Boer, C.G.; Metrustry, S.; Liu, Y.; den Hollander, W.; van Rooij, J.; Kraus, V.B.; Yau, M.S.; et al. Novel Genetic Variants for Cartilage Thickness and Hip Osteoarthritis. PLoS Genet. 2016, 12, e1006260. [Google Scholar] [CrossRef] [Green Version]
  28. Tu, Y.; Ma, T.; Wen, T.; Yang, T.; Xue, L.; Cai, M.; Wang, F.; Guan, M.; Xue, H. MicroRNA-377-3p alleviates IL-1beta-caused chondrocyte apoptosis and cartilage degradation in osteoarthritis in part by downregulating ITGA6. Biochem. Biophys. Res. Commun. 2020, 523, 46–53. [Google Scholar] [CrossRef] [PubMed]
  29. Ali, S.A.; Espin-Garcia, O.; Wong, A.K.; Potla, P.; Pastrello, C.; McIntyre, M.; Lively, S.; Jurisica, I.; Gandhi, R.; Kapoor, M. Circulating microRNAs differentiate fast-progressing from slow-progressing and non-progressing knee osteoarthritis in the Osteoarthritis Initiative cohort. Ther. Adv. Musculoskelet Dis. 2022, 14, 1759720X221082917. [Google Scholar] [CrossRef]
  30. Castanheira, C.; Balaskas, P.; Falls, C.; Ashraf-Kharaz, Y.; Clegg, P.; Burke, K.; Fang, Y.; Dyer, P.; Welting, T.J.M.; Peffers, M.J. Equine synovial fluid small non-coding RNA signatures in early osteoarthritis. BMC Vet. Res. 2021, 17, 26. [Google Scholar] [CrossRef]
  31. Flannery, C.R.; Hughes, C.E.; Schumacher, B.L.; Tudor, D.; Aydelotte, M.B.; Kuettner, K.E.; Caterson, B. Articular cartilage superficial zone protein (SZP) is homologous to megakaryocyte stimulating factor precursor and Is a multifunctional proteoglycan with potential growth-promoting, cytoprotective, and lubricating properties in cartilage metabolism. Biochem. Biophys. Res. Commun. 1999, 254, 535–541. [Google Scholar] [CrossRef] [PubMed]
  32. Dobin, A.; Davis, C.A.; Schlesinger, F.; Drenkow, J.; Zaleski, C.; Jha, S.; Batut, P.; Chaisson, M.; Gingeras, T.R. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 2013, 29, 15–21. [Google Scholar] [CrossRef] [PubMed]
  33. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef] [Green Version]
  34. Benjamini, Y.; Hochberg, Y. Controlling the False Discovery Rate—A Practical and Powerful Approach to Multiple Testing. J. R. Stat. Soc. B 1995, 57, 289–300. [Google Scholar] [CrossRef]
  35. Rot, G.; Wang, Z.; Huppertz, I.; Modic, M.; Lence, T.; Hallegger, M.; Haberman, N.; Curk, T.; von Mering, C.; Ule, J. High-Resolution RNA Maps Suggest Common Principles of Splicing and Polyadenylation Regulation by TDP-43. Cell Rep. 2017, 19, 1056–1067. [Google Scholar] [CrossRef] [Green Version]
  36. Kanehisa, M.; Furumichi, M.; Tanabe, M.; Sato, Y.; Morishima, K. KEGG: New perspectives on genomes, pathways, diseases and drugs. Nucleic. Acids Res. 2017, 45, D353–D361. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Fabregat, A.; Sidiropoulos, K.; Garapati, P.; Gillespie, M.; Hausmann, K.; Haw, R.; Jassal, B.; Jupe, S.; Korninger, F.; McKay, S.; et al. The Reactome pathway Knowledgebase. Nucleic. Acids Res. 2016, 44, D481–D487. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  38. Yu, G.; Wang, L.G.; Han, Y.; He, Q.Y. clusterProfiler: An R package for comparing biological themes among gene clusters. OMICS 2012, 16, 284–287. [Google Scholar] [CrossRef]
  39. Luo, W.; Brouwer, C. Pathview: An R/Bioconductor package for pathway-based data integration and visualization. Bioinformatics 2013, 29, 1830–1831. [Google Scholar] [CrossRef] [Green Version]
  40. Aibar, S.; Gonzalez-Blas, C.B.; Moerman, T.; Huynh-Thu, V.A.; Imrichova, H.; Hulselmans, G.; Rambow, F.; Marine, J.C.; Geurts, P.; Aerts, J.; et al. SCENIC: Single-cell regulatory network inference and clustering. Nat. Methods 2017, 14, 1083–1086. [Google Scholar] [CrossRef] [Green Version]
  41. Martin, I.; Jakob, M.; Schafer, D.; Dick, W.; Spagnoli, G.; Heberer, M. Quantitative analysis of gene expression in human articular cartilage from normal and osteoarthritic joints. Osteoarthr. Cartil. 2001, 9, 112–118. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Livak, K.J.; Schmittgen, T.D. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef] [PubMed]
  43. de Jonge, H.J.; Fehrmann, R.S.; de Bont, E.S.; Hofstra, R.M.; Gerbens, F.; Kamps, W.A.; de Vries, E.G.; van der Zee, A.G.; te Meerman, G.J.; ter Elst, A. Evidence based selection of housekeeping genes. PLoS ONE 2007, 2, e898. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Ijms 24 06611 g001 550
Figure 1. (A) Principal component analysis following determination of differential gene expression in QuantSeq-Reverse dataset using DeSeq 2. Red points are control samples and blue points are osteoarthritic samples. (B) Volcano plot showing differential expression of genes between datasets, with genes significantly upregulated in osteoarthritis highlighted in red and genes significantly downregulated in osteoarthritis highlighted in blue. (C) Scatterplot showing the relationship between the fold changes of genes significantly affected by osteoarthritis in this study and also significantly identified between intact and degenerate osteoarthritic cartilage in a previous study by Soul et al. [5]. The dashed line and R2 value are the result of linear regression analysis. (D) Gene ontology analysis of genes that were identified as being differentially regulated in osteoarthritic samples. The colour of each circle represents the p-value associated with the enrichment whilst circle size is associated with the number of genes associated with ontological term. (E) Transcription factor-binding motifs enriched in the promoters of genes that are differentially regulated by osteoarthritis. (F) Heatmap illustrating the expression of lncRNAs that are differentially expressed in osteoarthritic samples compared to controls.
Figure 1. (A) Principal component analysis following determination of differential gene expression in QuantSeq-Reverse dataset using DeSeq 2. Red points are control samples and blue points are osteoarthritic samples. (B) Volcano plot showing differential expression of genes between datasets, with genes significantly upregulated in osteoarthritis highlighted in red and genes significantly downregulated in osteoarthritis highlighted in blue. (C) Scatterplot showing the relationship between the fold changes of genes significantly affected by osteoarthritis in this study and also significantly identified between intact and degenerate osteoarthritic cartilage in a previous study by Soul et al. [5]. The dashed line and R2 value are the result of linear regression analysis. (D) Gene ontology analysis of genes that were identified as being differentially regulated in osteoarthritic samples. The colour of each circle represents the p-value associated with the enrichment whilst circle size is associated with the number of genes associated with ontological term. (E) Transcription factor-binding motifs enriched in the promoters of genes that are differentially regulated by osteoarthritis. (F) Heatmap illustrating the expression of lncRNAs that are differentially expressed in osteoarthritic samples compared to controls.
Ijms 24 06611 g001
Ijms 24 06611 g002 550
Figure 2. (A) The number of cartilage genes exhibiting different numbers of polyadenylation sites plotted as a percentage all genes. (B) Counts for the number of genes with different numbers of polyadenylation sites identified across control and osteoarthritis cartilage samples combined (blue circles) and human data from polyA_DB (orange circles). (C) Gene ontology analysis of genes where >5 polyadenylation sites were identified across all the cartilage samples analysed. The colour of each circle represents the p-value associated with the enrichment whilst circle size is associated with the number of genes associated with the ontological term.
Figure 2. (A) The number of cartilage genes exhibiting different numbers of polyadenylation sites plotted as a percentage all genes. (B) Counts for the number of genes with different numbers of polyadenylation sites identified across control and osteoarthritis cartilage samples combined (blue circles) and human data from polyA_DB (orange circles). (C) Gene ontology analysis of genes where >5 polyadenylation sites were identified across all the cartilage samples analysed. The colour of each circle represents the p-value associated with the enrichment whilst circle size is associated with the number of genes associated with the ontological term.
Ijms 24 06611 g002
Ijms 24 06611 g003 550
Figure 3. (A) Gene structure diagrams for OSMR, KMT2A and ANOS1. ENSEMBL IDs for the specific transcripts shown are indicated in parenthesis. Blocks represent exons and lines represent introns. Black shaded exon regions represent the coding domain. Orange arrows indicate the locations of polyadenylation sites that are differentially utilised in osteoarthritis. Note that the diagrams in this panel are not presented to scale. (B–D) Genome browser track example showing QuantSeq Reverse reads for three genes affected by alternative polyadenylation in osteoarthritic chondrocytes. Tracks for the six control (green) and six osteoarthritic (red) samples are shown. (B) Reads containing OSMR polyadenylation sites at intron 7 (left panel), which will lead to a transcript encoding a secreted form of OSMR and at the full-length transcript’s 3′UTR (right panel) which would encode the full plasma membrane receptor transcripts (each track is scaled from 0–3600 reads). (C) Polyadenylation site reads in the 3′ UTR of KMT2A transcripts (each track is scaled from 0–691 reads). (D) Reads containing intronic and full length 3′UTR polyadenylation sites for ANOS1 transcripts (each track is scaled from 0–590 reads). Note that OSMR and KMT2A are encoded on the sense strand and so are presented 5′–3′ left to right in panels B and C. ANOS1 is encoded on the antisense strand and so runs 5′ to 3′ right to left in panel D.
Figure 3. (A) Gene structure diagrams for OSMR, KMT2A and ANOS1. ENSEMBL IDs for the specific transcripts shown are indicated in parenthesis. Blocks represent exons and lines represent introns. Black shaded exon regions represent the coding domain. Orange arrows indicate the locations of polyadenylation sites that are differentially utilised in osteoarthritis. Note that the diagrams in this panel are not presented to scale. (B–D) Genome browser track example showing QuantSeq Reverse reads for three genes affected by alternative polyadenylation in osteoarthritic chondrocytes. Tracks for the six control (green) and six osteoarthritic (red) samples are shown. (B) Reads containing OSMR polyadenylation sites at intron 7 (left panel), which will lead to a transcript encoding a secreted form of OSMR and at the full-length transcript’s 3′UTR (right panel) which would encode the full plasma membrane receptor transcripts (each track is scaled from 0–3600 reads). (C) Polyadenylation site reads in the 3′ UTR of KMT2A transcripts (each track is scaled from 0–691 reads). (D) Reads containing intronic and full length 3′UTR polyadenylation sites for ANOS1 transcripts (each track is scaled from 0–590 reads). Note that OSMR and KMT2A are encoded on the sense strand and so are presented 5′–3′ left to right in panels B and C. ANOS1 is encoded on the antisense strand and so runs 5′ to 3′ right to left in panel D.
Ijms 24 06611 g003
Ijms 24 06611 g004 550
Figure 4. (A) QuantSeq-Reverse read numbers for proximal (intronic) and distal (3′UTR) polyadenylation sites for the OSMR and ANOS1 genes, which were differentially regulated between healthy (control) and osteoarthritic (OA) cartilage samples. OSMR proximal site is at CHR5:38886653, OSMR distal site is at CHR5:38935628, ANOS1 proximal site is at CHRX:8,536,187, ANOS1 distal site is at CHRX:85,28,874. Different colour bars represent read numbers from individual control or OA samples allowing per sample comparison between proximal and distal values in each case. (B) RT-qPCR analysis of mRNA decay rates of OSMR and ANOS1 transcripts that use intronic or canonical 3′UTR polyadenylation sites in human articular chondrocytes treated with actinomycin D. Mean ± standard error shown, n = 3. (C) RT-qPCR analysis of the ratio of intronic to canonical 3′UTR transcript levels of OSMR and ANOS1 in HAC treated with 10 ng/mL interleukin-1ß (IL-1ß). Individual data points are presented (n = 4) whilst line represents mean value and p-values represent result of one-way analysis of variance.
Figure 4. (A) QuantSeq-Reverse read numbers for proximal (intronic) and distal (3′UTR) polyadenylation sites for the OSMR and ANOS1 genes, which were differentially regulated between healthy (control) and osteoarthritic (OA) cartilage samples. OSMR proximal site is at CHR5:38886653, OSMR distal site is at CHR5:38935628, ANOS1 proximal site is at CHRX:8,536,187, ANOS1 distal site is at CHRX:85,28,874. Different colour bars represent read numbers from individual control or OA samples allowing per sample comparison between proximal and distal values in each case. (B) RT-qPCR analysis of mRNA decay rates of OSMR and ANOS1 transcripts that use intronic or canonical 3′UTR polyadenylation sites in human articular chondrocytes treated with actinomycin D. Mean ± standard error shown, n = 3. (C) RT-qPCR analysis of the ratio of intronic to canonical 3′UTR transcript levels of OSMR and ANOS1 in HAC treated with 10 ng/mL interleukin-1ß (IL-1ß). Individual data points are presented (n = 4) whilst line represents mean value and p-values represent result of one-way analysis of variance.
Ijms 24 06611 g004
Ijms 24 06611 g005 550
Figure 5. (A) QuantSeq-Reverse read numbers for proximal and distal polyadenylation sites for the KMT2A and SLC25A37 genes, usage of which leads to different length 3′UTRs and is differentially regulated between control and osteoarthritic (OA) samples. KMT2A proximal site is at CHR11:118523915, KMT2A distal site is at CHR11:118526830, SLC25A37 proximal site is at CHR8:23,572,547, SLC25A37 distal site is at CHR8:235,75,460. Different colour bars represent read numbers from individual control or OA samples allowing per sample comparison between proximal and distal values in each case. (B) RT-qPCR analysis of mRNA decay rates of KMT2A and SLC25A37 using primer pairs that are specific to only longer 3′UTR isoforms (long only) or both proximal and ideally polyadenylated isoforms (both) in human articular chondrocytes treated with actinomycin D. Mean ± standard error shown, n = 3. (C) Heat map showing ratio of expression of intronic to canonical 3′UTR (for OSMR and ANOS1) or long only to both for (KMT2A and SLC25A37) in a panel of different cell lines. (D) RT-qPCR analysis of mRNA decay rates of long only or both polyadenylated isoforms of KMT2A in AC10 cells. Mean ± standard error shown, n = 3.
Figure 5. (A) QuantSeq-Reverse read numbers for proximal and distal polyadenylation sites for the KMT2A and SLC25A37 genes, usage of which leads to different length 3′UTRs and is differentially regulated between control and osteoarthritic (OA) samples. KMT2A proximal site is at CHR11:118523915, KMT2A distal site is at CHR11:118526830, SLC25A37 proximal site is at CHR8:23,572,547, SLC25A37 distal site is at CHR8:235,75,460. Different colour bars represent read numbers from individual control or OA samples allowing per sample comparison between proximal and distal values in each case. (B) RT-qPCR analysis of mRNA decay rates of KMT2A and SLC25A37 using primer pairs that are specific to only longer 3′UTR isoforms (long only) or both proximal and ideally polyadenylated isoforms (both) in human articular chondrocytes treated with actinomycin D. Mean ± standard error shown, n = 3. (C) Heat map showing ratio of expression of intronic to canonical 3′UTR (for OSMR and ANOS1) or long only to both for (KMT2A and SLC25A37) in a panel of different cell lines. (D) RT-qPCR analysis of mRNA decay rates of long only or both polyadenylated isoforms of KMT2A in AC10 cells. Mean ± standard error shown, n = 3.
Ijms 24 06611 g005
Table
Table 1. The 10 most highly expressed alternatively polyadenylated genes in articular chondrocytes regardless of disease state.
Table 1. The 10 most highly expressed alternatively polyadenylated genes in articular chondrocytes regardless of disease state.
Gene Name, (Gene Symbol), Ensembl Gene IDGene Locus (Strand)Sites (Proximal to Distal)Site Distance from Proximal Site (bp)Site Mean Read Count ± SD
Decorin, (DCN), ENSG00000011465chr12:91140483–91182823 (antisense)91145968012,270 ± 8494
911458709819,760 ± 13,585
9114552644251,032 ± 39,570
Noncoding nuclear-enriched abundant transcript 2, (MALAT1), ENSG00000251562chr11:65497687–65506515 (sense)65499339028,517 ± 21,758
6549965932021,804 ± 17,950
65499698359100,707 ± 98,832
6549973639799,063 ± 96,843
6549977643750,554 ± 37,148
654997944555531 ± 4083
6550052211836120 ± 6761
6550055812195446 ± 5383
Complement factor H, (CFH), ENSG00000000971chr11:65497687–65506515 (sense)196701562038,800 ± 68,040
19674750145,9398984 ± 9034
Clusterin, (CLU), ENSG00000120885chr8:27596916–27614699 (antisense)2760902002143 ± 1472
27608995256215 ± 3978
27608957637438 ± 4695
2759799611,02433,471 ± 22,377
2759781311,2078339 ± 9719
2759772811,29221,976 ± 21,562
Osteonectin, (SPARC), ENSG00000113140chr5:151661095–151686974 (antisense)151662435022,006 ± 29,354
1516614499864341 ± 5672
60S ribosomal protein L7a, (RPL7A), ENSG00000280858chr9:133348217–133351425 (Sense)13335141204060 ± 2684
1333514241218,479 ± 12,447
Prolargin, (PRELP), ENSG00000188783chr1:203475805–203491351 (sense)203487161017,465 ± 13,603
203489613245216,419 ± 13,500
Cartilage intermediate layer protein, (CILP), ENSG00000138615chr15:65194759–65211472 (antisense)6519652903135 ± 4161
6519625127813,156 ± 16,572
651960005294965 ± 5868
H3 histone family member 3B, (H3F3B), ENSG00000132475chr17:75776433–75785892 (antisense)7577825002376 ± 2438
757780232272255 ± 2841
7577745679412,733 ± 12,977
Early growth response 1, (EGR1), ENSG00000120738chr5:138465478–138469302 (sense)13846922402641 ±2780
13846924016729 ± 649
1384693017712,484 ± 16,106
Table
Table 2. Genes affected by alternative polyadenylation in osteoarthritic chondrocytes with corrected p-values for distal and proximal sites set to <0.1.
Table 2. Genes affected by alternative polyadenylation in osteoarthritic chondrocytes with corrected p-values for distal and proximal sites set to <0.1.
Gene LocusStrandENSEMBL Gene_idGene SymbolProximal/Distal polyA SitesProximal/
Distal Feature		Proximal log2fcDistal log2fcProximal Adjusted p-ValueDistal Adjusted p-Value
chr11:118482737–118526831+ENSG00000118058KMT2A118523915, 118526830utr3/utr32.06168−1.674040.016250.00558
chr11:10812073–10858795-ENSG00000236287ZBED510852797, 10852747utr3/utr32.29148−2.948750.018120.0418
chr9:97501179–97601742+ENSG00000136842TMOD197600155, 97601252utr3/utr31.51617−1.249270.021170.01554
chr21:15730024–15880068+ENSG00000155313USP2515878905, 15880062utr3/utr3−1.528653.440740.003560.054
chr11:65422773–65445539+ENSG00000245532NEAT165426414, 65426524exon/exon−0.821731.044920.069510.08359
chr19:45349836–45370917-ENSG00000104884ERCC245351594, 45349837utr3/utr3−1.704592.325430.04180.07641
chr11:111540279–111561744+ENSG00000204381LAYN111560953, 111561510utr3/utr3−1.465322.046480.060240.09604
chr5:136132844–136188746+ENSG00000113658SMAD5136179501, 136182732utr3/utr3−3.688031.678260.085760.00919
chr9:129634603–129642168-ENSG00000148331ASB6129636893, 129634604utr3/utr3−1.457831.294240.096040.08311
chr12:51186935–51217707-ENSG00000184271POU6F151189184, 51186936utr3/utr31.70409−1.327880.083590.0527
chr8:23528955–23575462+ENSG00000147454SLC25A3723572547, 23575460utr3/utr30.68296−1.243190.065120.07641
chr2:44316280–44361861-ENSG00000138078PREPL44318763, 44317605utr3/intron−1.549522.273640.016460.04608
chr5:38845857–38945595+ENSG00000145623OSMR38886653, 38935628intron/utr30.7451−1.815060.032280.00233
chr5:132556018–132646348+ENSG00000113522RAD50132616045, 132642473utr3/utr33.07214−5.496430.006830.04908
chr3:48241099–48299252+ENSG00000164048ZNF58948270986, 48273225utr3/intron−0.702270.759390.049080.05682
chr7:77798791–77957502+ENSG00000006576PHTF277940457, 77957223utr3/utr33.12642−2.271330.051890.01812
chr2:85751343–85791382+ENSG00000168874ATOH885775209, 85791376exon/utr33.02271−3.804910.018120.05189
chr6:85449583–85495790+ENSG00000135318NT5E85452154, 85495783intron/utr31.14747−1.703380.049080.0273
chrX:8528873–8732136-ENSG00000011201ANOS18536187, 8528874intron/utr33.04985−3.079180.000080.00008
chr12:19404044–19720800+ENSG00000139154AEBP219441779, 19520222intron/utr33.11654−1.151560.094970.00558
Table
Table 3. Sequences of primers used in qRT-PCR assays.
Table 3. Sequences of primers used in qRT-PCR assays.
Assay NameForward Primer Sequence (5′–3′) Reverse Primer Sequence (5′–3′)
RPS13GGCTTTACCCTATCGACGCAGTCAGATGTCAACTTCAACCAAGTG
OSMR proximalGCAGTAGGTTGTCTGGGTCATGAGCCTGGAACAAACAGCA
OSMR distalTCATTCACAGCGGAGGTGAGGCGTGCATCCATGAGGAGAA
ANOS1 proximalCCAGATCCTGCCTTCCGTAGAGTGCAGCACACAGAATGAC
ANOS1 distalATGCAGATGCCTGGCCATTAGGCTCTGTGGAGTACAGTGA
KMT2A proximal CACCTACAGCGTCTGTCGAACACCGGAGGTGCTAGGAATC
KMT2A distal GGGGTTCCACTAGTGTCTGCTTCAGGACCCATCAGTGCAT
SLC25A37 proximalGGTGCATCTTACCGAGGAGGGTGCATGGTCCCAAGATGGA
SLC25A37 distalGTGTGCTTGTGCGTGTCTACTGTTGCCTTTTCGTTCACCTG
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.

Share and Cite

MDPI and ACS Style

Winstanley-Zarach, P.; Rot, G.; Kuba, S.; Smagul, A.; Peffers, M.J.; Tew, S.R. Analysis of RNA Polyadenylation in Healthy and Osteoarthritic Human Articular Cartilage. Int. J. Mol. Sci. 2023, 24, 6611. https://doi.org/10.3390/ijms24076611

AMA Style

Winstanley-Zarach P, Rot G, Kuba S, Smagul A, Peffers MJ, Tew SR. Analysis of RNA Polyadenylation in Healthy and Osteoarthritic Human Articular Cartilage. International Journal of Molecular Sciences. 2023; 24(7):6611. https://doi.org/10.3390/ijms24076611

Chicago/Turabian Style

Winstanley-Zarach, Phaedra, Gregor Rot, Shweta Kuba, Aibek Smagul, Mandy J. Peffers, and Simon R. Tew. 2023. "Analysis of RNA Polyadenylation in Healthy and Osteoarthritic Human Articular Cartilage" International Journal of Molecular Sciences 24, no. 7: 6611. https://doi.org/10.3390/ijms24076611

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop