miR-26a-5p is a Stable Reference Gene for miRNA Studies in Chondrocytes from Developing Human Cartilage

miRNAs are emerging as key regulators of complex biological systems in several developmental processes. qRT-PCR is a powerful tool to quantitatively assess the profiles and modulation of miRNA expression. In the emerging field of cartilage maturation studies, from precursor to hypertrophic chondrocytes, few data about miRNA regulation are available, and no consensus on the best reference gene (RG) has been reached. This is a crucial pitfall since reliable outcomes depend on proper data normalization. The aim of this work was to identify reliable and stable miRNA RGs, basing the analysis on available high throughput qRT-PCR miRNA data (from the NCBI Gene Expression Omnibus database, GSE49152) obtained from human embryonic cartilage tissues enriched in the precursor, differentiated, and hypertrophic chondrocytes. Four normalization approaches were used, and the stability was quantified by combining BestKeeper, delta-Ct, geNorm, and NormFinder statistical tools. An integrated approach allowed to identify miR-26a-5p as the most stable RG and miR-212-3p as the worst one. RNU44, used in original dataset analysis, performed as second best RG. Applications of different normalization strategies significantly impacted the profiles and modulation of miRNA expression. Herein presented results point out the crucial need of a consensus on data normalization studies aimed at dissecting miRNA role in human cartilage development, to avoid the postulation of unreliable biological conclusions.


Introduction
Gene expression regulation is the biological foundation for the specification of every cell type, tissue, and organ in a multicellular organism. Together with transcriptional regulators, miRNAs have emerged as gene expression repressors directing the post-transcriptional modulation that underlies development [1]. miRNAs are evolutionary conserved, single-stranded, 21-24 nucleotide-long non-coding RNA molecules [2]. They are usually transcribed from DNA sequences first into primary miRNAs, then processed into precursor miRNAs and eventually mature molecules. The main mechanism of action is through direct interaction with 3 untranslated region (3' UTR) of target mRNAs, which induces mRNA degradation resulting in translational repression [3]. In addition, miRNAs may interact with 5 UTR, coding sequence, gene promoters, and under certain conditions, they can also activate translation or regulate transcription [2]. To date, 2300 validated human miRNAs have been reported, and new miRNAs are still being discovered with their roles in gene regulation, starting to be fully deciphered [4]. miRNAs are involved in a variety of biological processes regulating human development in different districts. For some of them, a clear picture has been depicted [5]. As a striking example, for neural development, hundreds of miRNAs are involved in determining the fate of the two major cell types of the nervous system, neurons and glia, at both embryonic and early postnatal stages [6]. An advantage of the neural system is that these processes of neuro-and gliogenesis involve many intermediate cell types that have been exhaustively studied in terms of gene expression and non-coding RNAs regulation [7]. Similarly, in the bone system, the precise roles of many miRNAs have been fully deciphered, giving valuable insights into the treatment of developmental disorders of the skeleton [8]. On the contrary, despite an increasing amount of in vitro data, there is still limited information on the expression and function of specific miRNAs in human cartilage development in vivo.
During embryonic development, mesenchymal progenitor cells differentiate into chondrocytes to form cartilage templates for future bones and cartilage. Proliferating chondrocytes produce extracellular matrix, enriched in type II collagen and aggrecan. Further, they differentiate into hypertrophic chondrocytes that express type X collagen and undergo mineralization to be eventually replaced by mineralized bone, leading to longitudinal bone growth [9]. In this process, different signaling molecules and transcription factors have been shown to regulate the progression of chondrocyte-specific gene expression [10,11]. Regarding miRNAs, few data are available in the in vivo settings, and actual knowledge has predominantly come from studies in mice. let-7 miRNA is required for chondrocyte proliferation [12], whereas miR-140 modulates premature hypertrophic chondrocyte differentiation and delays differentiation of resting chondrocytes to proliferating chondrocytes [13,14]. Nevertheless, data of the human model depicting modulations or differences between cells at different stages of differentiation were missing.
To overcome the limitations of these studies, a pivotal work on miRNA expression patterns was performed by high-throughput qRT-PCR technique within human embryonic (gestational day 54-56) cartilage laser dissected fragments, containing either precursor (PC), differentiated (DC), or hypertrophic (HYP) chondrocytes (Gene Expression Omnibus database, GSE49152) [15]. This report was selected being unique in combining miRNA analysis and specific chondrocyte populations directly involved in regulating cartilage and long bone development. Authors were able to identify differentially expressed miRNAs predicted to regulate growth factors (vascular endothelial growth factor (VEGF), insulin-like growth factor-1 (IGF-1), transforming growth factor beta (TGF-β), bone morphogenetic protein (BMP), fibroblast growth factor (FGF)), Hedgehog and Wnt signaling pathways, and interleukin (IL)-8. Despite this fundamental milestone, only around 50 miRNAs were detectable in all the nine donors, leaving room for further studies aimed either at identifying those candidates that were missed for both technical or low expression reasons, or at studying dysregulation of miRNAs modulating developmental processes involved in skeletal disorders, such as osteoarthritis (OA), chondrodysplasias, or delayed endochondral fracture healing.
For a deep analysis of single or few miRNAs, a robust normalization strategy is a fundamental requisite. Therefore, mining qRT-PCR data of McAlinden dataset, through gold-standard statistical tools (BestKeeper [16], geNorm [17], NormFinder [18], and the comparative delta-Ct method [19]), the aim of this work was to identify stable miRNAs to be used as Reference Genes (RGs) in future studies dissecting these small RNAs involved in cartilage developmental or pathologic processes.

miRNA Selection
Mining qRT-PCR data, it was possible to identify 46 miRNAs with positive amplification values in all the 26 samples (8 PC, 9 DC, and 9 HYP). RNU44, used in McAlinden work as normalizer, was also always scored. The distribution of the quantification cycles (Ct) values of the selected reference genes over the whole sample sets is shown in Appendix A Table A1. miRNAs presented different expression values and variability levels in the three datasets. In PC chondrocytes, miR-26a-5p showed the largest expression (mean Ct = 13.04), while miR-362-3p was the least expressed (mean Ct = 23.80). In DC samples, miR-720 (mean Ct = 13.34) and miR-362-3p (mean Ct = 24.24) were poles apart. In HYP tissues, miR-720 (mean Ct = 14.26) was the most abundant, with miR-362-3p being again low expressed (mean Ct = 25.83). RNU44 was scored with low Ct values (mean Ct = 13.97 for PC, 14.42 for DC, and 15.57 for HYP). In terms of variability, miR-769-5p in PC, miR-202-3p in DC and HYP displayed the lowest standard deviation (SD), 0.51, 1.02, and 1.48, respectively. Eventually, correlation analysis was performed between each dataset. The analysis, shown in Figure 1, pointed out the existence of high correlation (R 2 > 0.8) for PC/DC and DC/HYP when compared to one another, while PC and HYP had lower interrelationship (R 2 = 0.6). This is in agreement with modulation of miRNA expression between precursor and hypertrophic chondrocytes, passing by an intermediate differentiation phase.

Impact of Normalization Strategy on miRNA Profiling
The effects of best/worst normalization strategies on target miRNA profiles were evaluated. The expression levels were computed using either the best (miR-26a-5p) or the worst (miR-212-3p) normalizers and compared with the results obtained with the normalization approach used in McAlinden analysis (RNU44) that also demonstrated its suitability in our ranking. Heat maps in Figure 3A clearly show that all three normalization approaches were able to cluster a group composed of H24243-HYP, H23689-PC, H23814-HYP, H23387-PC, and H23731-DC samples. The hierarchy in this group was maintained for miR-26a-5p and RNU44, while changed using miR-212-3p. Further, with the wrong normalization approach, an apparent and erroneous overexpression of almost all miRNAs in H24243-PC/DC samples emerged. Moreover, despite similar gross outcomes for the two most stable RGs, defined by the identification of those samples clearly different from the others, small rearrangement in both samples and miRNAs dendrograms could be observed, confirming that a most refined strategy allows identifying more specific co-regulation.
To evaluate the crucial impact of the best strategy on the reliable quantification of subtle discrepancies between cartilage regions, a more refined analysis was performed on miR-335-3p levels. miR-335-3p was selected since in McAlinden analysis, performed using RNU44, it was reported as significantly modulated within PC, DC, and HYP samples. In Figure 3B, it clearly appears that in a context of similar outcomes between the two approaches, the use of the most reliable miR-26a-5p does not allow scoring a significant (p-value ≤ 0.05) difference between PC and DC samples.

Impact of Normalization Strategy on miRNA Profiling
The effects of best/worst normalization strategies on target miRNA profiles were evaluated. The expression levels were computed using either the best (miR-26a-5p) or the worst (miR-212-3p) normalizers and compared with the results obtained with the normalization approach used in McAlinden analysis (RNU44) that also demonstrated its suitability in our ranking. Heat maps in Figure 3A clearly show that all three normalization approaches were able to cluster a group composed of H24243-HYP, H23689-PC, H23814-HYP, H23387-PC, and H23731-DC samples. The hierarchy in this group was maintained for miR-26a-5p and RNU44, while changed using miR-212-3p. Further, with the wrong normalization approach, an apparent and erroneous overexpression of almost all miRNAs in H24243-PC/DC samples emerged. Moreover, despite similar gross outcomes for the two most stable RGs, defined by the identification of those samples clearly different from the others, small rearrangement in both samples and miRNAs dendrograms could be observed, confirming that a most refined strategy allows identifying more specific co-regulation.  Values are shown as delta-Ct with respect to normalizing reference genes (RGs). # stands for p-value = 0.1, * for p-value ≤ 0.05, ** for p-value < 0.01, *** for p-value < 0.001, and **** for p-value < 0.0001. (C) Pearson correlation scatter plots for miR-335-3p and miR-26a-5p/RNU44 RGs. x and y-axis indicate Ct values.
To evaluate the crucial impact of the best strategy on the reliable quantification of subtle discrepancies between cartilage regions, a more refined analysis was performed on miR-335-3p levels. miR-335-3p was selected since in McAlinden analysis, performed using RNU44, it was reported as significantly modulated within PC, DC, and HYP samples. In Figure 3B, it clearly appears that in a context of similar outcomes between the two approaches, the use of the most reliable miR-26a-5p does not allow scoring a significant (p-value ≤ 0.05) difference between PC and DC samples. To rule out this variation, correlation analysis for miR-26a-5p, RNU44, and miR-335-3p was performed across all samples. Interestingly, RNU44 and miR-335-3p showed a higher correlation coefficient (R 2 = 0.68 vs. 0.53 for miR-26a-5p/miR-335-3p) ( Figure 3C), leading to more homogenous ratios with reduced standard deviation and eventually higher statistical significance. Notably, the correlation between miR-26a-5p and RNU44 was high (0.82) as expected for RGs with similar stability profiles, again emphasizing that from small differences, major distortions may arise.

Discussion
The present work identified miR-26a-5p as a suitable RG for studies on miRNAs involved in and controlling the developmental process of human cartilage. By a multi-technique quantitative approach, we also demonstrated the good performance of the widely-used normalizer RNU44, in this specific experimental context. The application of different normalization strategies in the exemplary case of miR-335-3p assessment pointed out the criticality of the RG choice to obtain significant information on miRNA modulation in cartilage development.
Normalization strategy is an open question in studies assessing miRNA stability and comparison between samples from different sources or donors. Traditionally, in qRT-PCR studies, the relative quantification method is used, comparing expression levels of target miRNAs with the amount of an endogenous RG. In this context, small nuclear/nucleolar RNAs (snRNAs, e.g., RNU6, RNU44, or RNU48) have been commonly preferred [20]. This approach has some advantages, as well as important pitfalls. In fact, small RNA molecules, such as miRNAs and snRNAs, share similar features, such as stability and size. Further, snRNAs are ubiquitous and abundantly expressed, crucial traits for a reliable RG [21]. Nevertheless, snRNA biogenesis is mechanistically separated from miRNA biogenesis. As an example, RNU6 is not processed by the spliceosome but by the Drosha complex and does not mirror the physicochemical properties of miRNA molecules [22]. Therefore, the selection of the reference RNA on the class of RNAs being investigated is a fundamental issue, suggesting that reference miRNAs would be preferred for miRNAs.
To date, no universal miRNA RGs have been proposed due to high variability between donors, tissues, and developmental stages, making miRNA studies results largely incomparable [23]. Therefore, the selection of RGs that are reliable within the experimental condition under analysis is a pivotal pre-requisite. In this perspective, for large miRNA datasets, the global mean expression value normalization was proposed as highly effective, in terms of both reduction of technical variation and more accurate quantification of biological fluctuations [24]. However, when only a limited number of miRNA molecules are profiled, the selection of specific RGs is mandatory. To answer this need, the identification of invariant miRNAs by algorithms, specifically developed for RG evaluation and selection, resulted in a promising strategy [16][17][18][19].
In the field of developmental biology, such strategy has been successfully applied in a few studies in the plant or animal systems [25][26][27][28][29]. In human studies, a major issue is an availability of starting tissue material, especially when developmental studies are performed at embryonic stages. In the present work, the valuable data regarding miRNA expression in embryonic cartilage fragments, containing either precursor, differentiated, or hypertrophic chondrocytes [15] was mined. RNU44, used by McAlinden group as a normalizer for differentially-expressed miRNAs, was a reliable RG, ranking second when all samples were grouped, and always within the first ten positions in the separated cartilage regions. From our analysis, miR-26a-5p emerged as the most accurate RG, standing first for the analysis of all samples together, and in the top five for disjointed isolates. At present, no reports connecting miR-26a-5p and cartilage development are known, reinforcing its suitability as stable RG, although more focused studies are needed. In a wider context related to cartilage, in human chondrocytes, the downregulation of miR-26a-5p expression by IL-1b through nuclear factor kappa B (NF-kB) enhanced the production of OA-related nitric oxide synthase 2 (iNOS) protein, due to the ability of miR-26a-5p to directly target iNOS mRNA 3'UTR [30]. Therefore, although stable in developing cartilage, miR-26a-5p reliability as chondrocyte miRNA RG in adult tissues under inflammation or OA might be considered carefully, and our suggestion is to verify, through the herein proposed algorithms, a panel of miRNA RGs, if possible selected from studies in the same or similar field.
The marked differences in performance of the tested normalization strategies (miR-26a-5p/RNU44/ miR-212-3p) had a significant impact on miRNAs profiling in the different regions of developing cartilage, possibly precluding accurate biological implications. Although all normalization strategies allowed to clearly distinguish a subgroup of markedly different samples, either best or worst RGs completely changed the arrangement of both sample and miRNA dendrograms ( Figure 3A). Moreover, even using miR-26a-5p, it was not possible to clearly cluster PC, DC, and HYP regions. This may indicate that overall, or at least for the 46 scored miRNAs, no fingerprinting major transcriptional differences are present. In this context of global high similarity, few miRNAs may be crucial to support biological outcomes, and reliable and refined evaluation of small variations gets decisive, as shown for miR-335-3p, in assessing significant modulations.
In conclusion, this study pointed out for the first time the reliability of miR-26a-5p for future studies aimed at dissecting new or low abundant miRNAs regulation in different regions of human developing cartilage. The main limitation of the report is the lack of validation of the proposed analysis on independent samples. Therefore, future evaluation of miRNA expression from developing cartilage should be performed on the path of herein proposed candidates, possibly from multiple samples for each individual to reduce overall variability [31].

Data Retrieval and Ethics Statement
qRT-PCR data can be found in the GEO database [32], with the record GSE49152 [15]. Only miRNAs with reported Ct values in all samples were considered for the analysis (Appendix A Tables A1-A3). For the donors, briefly, human, normal embryonic tissue samples (limbs) at gestational day 54-56 were obtained from nine donors, and laser capture microdissection was performed to obtain tissue from the precursor chondrocyte (PC), differentiated chondrocyte (DC), or hypertrophic chondrocyte (HYP) regions. After RNA extraction, TaqMan ® OpenArray ® technology was used to determine microRNA expression profiles starting from 30 ng total RNA for each sample (Thermo Fisher Scientific, Waltham, MA, USA). TaqMan ® OpenArray ® technology is a fixed-content panel containing 754 validated human TaqMan ® MicroRNA Assays derived from Sanger miRBase release v.14. The panel is specifically designed to provide specificity for only the mature miRNA targets. TaqMan MicroRNA Assays incorporate a target-specific stem-loop reverse transcription primer allowing to work despite the short length of mature miRNAs (~22 nucleotides).
As stated in McAlinden work, the Human Research Protection Office (HRPO) at Washington University in St Louis reviewed the request to work with human embryonic tissue, and the original project was deemed exempt since it did not constitute human subjects research and receiving embryonic tissue from the University of Washington would not involve obtaining data through intervention or interaction with a living individual, and, other than gestational age, no identifying information was provided upon receipt of the tissue.

Assessment of RG Stability
Gene expression stability was evaluated according to four gold-standard statistical approaches: BestKeeper [16], geNorm [17], NormFinder [18], and the comparative delta-Ct method [19]. BestKeeper analysis uses Ct values directly, while geNorm, NormFinder, and delta-Ct method use transformed Ct values of (1 + E) − ∆Ct. The ranking of the RGs according to their stability was generated by each algorithm, and a series of continuous integers starting from 1 was assigned to each RG. The overall performance of the miRNA RGs was evaluated by combining the results of the four approaches through a global ranking obtained as the geometric mean of the rankings given by each analysis [33,34].

Statistical Analyses
Statistical analyses were performed using GraphPad Prism Software version 5 (GraphPad, San Diego, CA, USA). Presence of outliers was scored by Grubbs' test. When two sets of data (PCvsDC, PCvsHYP, DCvsHYP) were compared, as in McAlinden paper, the comparison was performed by using unpaired Student t-test. Significance level was set at p-value ≤ 0.05. Pearson correlation coefficient (R 2 ) was estimated to determine the linear association between samples or miRNA Ct values. The outcome results were interpreted according to the degree of association [35].