Abstract

   The goal of the present study is to identify the differential
   expression of circular RNA (circRNA), miRNA, and piwi-interacting RNA
   (piRNA) after lineage commitment towards osteo- and chondrogenesis of
   human bone marrow mesenchymal stromal cells (hMSCs). The cells were
   maintained for 7 days in either osteogenic or chondrogenic medium. RNA
   sequencing was performed to assess the expression of miRNA and piRNA,
   while RNA hybridization arrays were used to identify which circRNA were
   differentially expressed. qPCR validation of a selection of targets for
   both osteogenic and chondrogenic differentiation was carried out. The
   differential expression of several circRNA, miRNA, and piRNA was
   identified and validated. The expression of total and circular isoforms
   of FKBP5 was upregulated both in osteo- and chondrogenesis and it was
   influenced by the presence of dexamethasone. ZEB1, FADS2, and SMYD3
   were also identified as regulated in differentiation and/or by
   dexamethasone. In conclusion, we have identified a set of different
   non-coding RNAs that are differentially regulated in early osteogenic
   and chondrogenic differentiation, paving the way for further
   investigation to understand how dexamethasone controls the expression
   of those genes and what their function is in MSC differentiation.

   Keywords: human mesenchymal stem cells, RNA sequencing, circular RNA,
   miRNA, piwi-interacting RNA, dexamethasone, differentiation,
   transcriptome, non-coding RNA

1. Introduction

   Bone marrow-derived mesenchymal stromal cells (MSCs) potentially
   represent a good source of cells for bone and cartilage tissue
   engineering. MSCs are widely studied and used in several clinical
   trials, with >80 active trials (clinicaltrials.gov, September 2019) and
   at least 20 for the treatment of musculoskeletal diseases. Despite all
   the effort, there is still little evidence that the use of MSCs is
   effective in musculoskeletal tissue repair. The reasons are manifold,
   one being the lack of control of a cell’s final differentiation. For
   example, the formation of mechanically inferior fibrous cartilage,
   instead of hyaline, is one the most common outcomes for articular
   cartilage tissue engineering [[34]1]. Moreover, the extensive ex vivo
   manipulation required to obtain large quantities of cells may alter
   their phenotype and response [[35]2]. Thus, there is a strong need for
   a better understanding of the underlying molecular mechanisms in order
   to achieve an efficient differentiation. Besides the genes coding for
   lineage-specific proteins, such as transcription factors, surface
   receptors, or matrix proteins, cellular processes comprise a variety of
   different non-coding RNAs (ncRNA) which have regulatory roles and
   actively participate in lineage specification. The study of ncRNA is
   therefore of utmost importance to understand the molecular mechanisms
   of differentiation and how to tune their regulation to be able to build
   and replace tissues.

   The most studied type of ncRNA are microRNAs (miRNA), and their role in
   stem cell maintenance and differentiation has been investigated
   [[36]3,[37]4,[38]5]. A related class of small ncRNA is represented by
   piwi-interacting RNA (piRNA). piRNA associates with the Piwi subfamily
   of Argonaute proteins, protecting genome integrity of germline cells by
   silencing transposable elements through the formation of RNA-induced
   silencing complexes with Piwi proteins, also known as piRISC
   [[39]6,[40]7]. Accordingly, their expression has been identified in
   testis [[41]8], ovary [[42]9], and epididymis [[43]10]. Nonetheless,
   there is a growing body of evidence that piRNAs also possess important
   regulatory functions in somatic and cancer cells
   [[44]7,[45]11,[46]12,[47]13,[48]14,[49]15,[50]16,[51]17,[52]18].

   Circular RNA (circRNA) represents a relatively newly recognized class
   of long ncRNA. Even though their existence has been long acknowledged,
   circular forms of RNA were thought to merely represent splicing errors
   with no biological function. Through the advancement of next generation
   sequencing, which has allowed the rise in RNA sequencing from diverse
   samples, circRNAs are now recognized to have distinct biogenesis and to
   regulate gene expression and biological processes through different
   mechanisms [[53]19,[54]20], with some miRNA-sponging circRNAs
   identified [[55]21,[56]22,[57]23,[58]24,[59]25,[60]26]. Moreover, since
   circRNAs derive from precursor mRNA, they are believed to influence the
   transcription and/or translation of the linear, protein-coding form
   [[61]20].

   In any case, the role of circRNA in different biological processes is
   becoming more and more evident. Up to now, a few publications have
   reported the expression and possible functions of circRNA in MSCs from
   different organisms, tissue sources, and experimental settings
   [[62]27,[63]28,[64]29,[65]30,[66]31,[67]32]. circRNA were
   differentially expressed in osteogenic differentiation of rat bone
   marrow MSCs engineered with estrogen receptor β-targeting shRNA
   [[68]33]. In human cells, the involvement of circRNAs in osteogenic
   differentiation of periodontal ligament stem cells
   [[69]34,[70]35,[71]36] or maxillary sinus membrane stem cells [[72]37]
   was reported. One study focused on the osteogenic differentiation of
   bone marrow stem cells [[73]38]. A map of circRNA expression in
   clinically relevant tissues, including MSCs and differentiated cells,
   was also suggested [[74]39] but, until now, there is a paucity of data
   regarding differential circRNA expression during osteogenic and
   chondrogenic differentiation of human bone marrow MSCs.

   The aim of the present study is to identify differential expression of
   circRNA, miRNA, and piRNA after lineage commitment towards osteo- and
   chondrogenesis. Based on previous work highlighting the importance of
   early events during differentiation, we investigated after day 7 of
   initiation of differentiation induction. Thus, human MSCs were
   maintained for 7 days in either osteogenic or chondrogenic medium with
   relative controls, and RNA sequencing or RNA hybridization arrays were
   performed. We then identified and validated several differentially
   regulated circRNA, miRNA, and piRNA, and revealed for some genes the
   influence of dexamethasone.

2. Materials and Methods

2.1. Cell Isolation and Culture

   MSCs were isolated from human bone marrow, as previously described
   [[75]40]. The samples were obtained with full ethical approval (Bern
   Req-2016-00141). Cells from a total of 16 donors were used (11 M/5 F;
   age mean 60 years; age range 33–80). Cells were cultured maintaining an
   initial cell density of 3 × 10^3 cells/cm^2 and grown until passage 2
   (p2) in Minimum Essential Medium Eagle-alpha modification (α-MEM,
   Gibco, Thermo Fisher, Zürich, Switzerland) with the addition of 10%
   MSC-qualified FBS (Pan-Biotech, Aidenbach, Germany), 100 U/mL
   penicillin, 100 μg/mL streptomycin (Gibco) and 5 ng/mL basic fibroblast
   growth factor (bFGF, Fitzgerald Industries International, Acton, MA,
   USA). Cultures were kept in a 37 °C/5% CO[2] humidified atmosphere and
   the medium was refreshed every second day.

   For initial RNA sequencing and hybridization arrays, MSCs from three
   donors were used. Cells were profiled for standard MSC markers with
   flow cytometry, as previously described [[76]41]. For further
   validation, MSCs from 6 additional donors for osteogenic
   differentiation and chondrogenic differentiation were isolated and
   differentiated following the same procedure (final n = 9). At p3, cells
   were induced to either osteo- or chondrogenic differentiation or
   cultured under control conditions.

   For osteogenic differentiation, cells were seeded at a density of 1.5 ×
   10^4 cells/cm^2 in 6-well plates. After overnight cell attachment, the
   medium was changed: control (CRL monolayer) samples were cultured in
   Dulbecco’s Modified Eagle Medium (DMEM) 1 g/L glucose (Gibco), 10%
   heat-inactivated FBS (Gibco), 100 U/mL penicillin and 100 μg/mL
   streptomycin (Gibco). Osteogenic differentiation was induced by the
   addition of 50 µg/mL ascorbic acid 2-phosphate, 5 mM β-glycerol
   phosphate, and 10 nM dexamethasone (all from Sigma-Aldrich) to the
   control medium. Samples were prepared in triplicate for each donor and
   condition.

   For chondrogenic differentiation, cells were seeded into pellet
   culture, using 2 × 10^5 cells/pellet. Cells were directly seeded into
   either control (CRL pellet) or chondrogenic medium. Control medium was
   composed of DMEM 4.5 g/L glucose, 100 U/mL penicillin, 100 μg/mL
   streptomycin, 1% Corning ITS+ Premix Universal Culture Supplement, and
   1% non-essential amino acids (NEAA). Chondrogenic medium was prepared
   from control medium with the addition of 50 µg/mL ascorbic acid
   2-phosphate, 100 nM dexamethasone, and 10 ng/mL TGF-β1 (Fitzgerald
   Industries). The medium was refreshed every second day to keep stable
   levels of differentiation factors in the media and samples were
   collected at day 7 for RNA isolation.

   Moreover, cultures were maintained up to day 21 for assessing final
   differentiation outcomes (Alizarin Red staining for osteogenic
   differentiation and Safranin-O/Fast Green staining on cryosections for
   chondrogenic differentiation). The RNA samples derived from those three
   donors were not used for sequencing but for validation purposes only. A
   summary of donor cohorts used for all experiments is reported in
   [77]Table S1.

2.2. RNA Isolation

   Total RNA was isolated from day 7 control (monolayer), day 7 control
   (pellet), day 7 osteogenic differentiation (OSTEO), and day 7
   chondrogenic differentiation (CHONDRO). For control (monolayer) and
   OSTEO samples, RNA was isolated from 1 well in triplicate for each
   donor. For control (pellet) and CHONDRO samples, 4 pellets were pooled
   in one sample and three samples were collected for each condition.
   Total RNA was isolated using standard TRIreagent (Molecular Research
   Center Inc., Cincinnati, OH, USA) extraction with
   1-Bromo-3-chloropropane (Sigma-Aldrich, St. Louis, MO, USA). After
   phase separation, RNA was precipitated from the aqueous phase with the
   addition of 2-propanol (Sigma-Aldrich) with an overnight incubation at
   −20 °C to improve the recovery of small RNA species. Total RNA
   concentration was measured with NanoDrop 1000 (Thermo Fisher) and
   purity was assessed by evaluation of the A260/280 and A260/230 ratios.

2.3. RNA Sequencing and Hybridization

   Two replicates for each group were selected for sequencing (n = 24
   samples). Two µg of total RNA for each sample was collected and dried
   in RNAstable^TM tubes (Biomatrica, San Diego, CA, USA) following the
   manufacturer’s instruction. As previously described [[78]42], samples
   were sent to ArrayStar (Rockville, MD, USA) for library generation and
   sequencing and RNA microarray hybridization.

   Libraries for RNAseq were denatured as single-stranded DNA molecules,
   captured on Illumina flow cells, amplified in situ as clusters, and
   sequenced for 51 cycles on an Illumina NextSeq500 system per the
   manufacturer’s instructions.

   For miRNA, the Solexa CHASTITY quantity filtered reads were harvested
   after sequencing as clean reads. The adaptor sequences were trimmed
   with cutadapt and the adaptor-trimmed reads (≥ 15 nt) were left.
   miRDeep2 software was used to quantify known miRNA and predict the
   novel miRNAs. The CPM value for miRNA and the differentially expressed
   miRNA were calculated and filtered with R package edgeR. Fold change
   (cutoff 1.5), p-value (0.05), and CPM (1 mean in one group) were used
   for filtering differentially expressed miRNAs. Hierarchical clustering
   was performed. miRNA target prediction was performed by targetscan
   ([79]www.targetscan.org) and miRDB ([80]www.mirdb.org), then the GO and
   KEGG pathway analyses were performed based on the top 10 differentially
   expressed miRNAs.

   For piRNA, the quality of sequencing was examined by FastQC software
   ([81]https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and
   the trimmed reads (that passed Illumina quality filter, trimmed
   3’-adaptor bases by cutadapt) were aligned to piRBase
   ([82]www.pirbase.org), a manually curated database for human piRNA
   [[83]43,[84]44], using NovoAlign software (v2.07.11, Novocraft
   Technologies Sdn Bhd, Selangor, Malaysia). The maximum mismatch of ≤2
   reads were kept. The expression profiling and differential expression
   of piRNAs were calculated based on normalized TPM. Hierarchical
   clustering, scatter plots, and classification analysis were performed
   with the differentially expressed piRNAs in R or Perl environment for
   statistical computing and graphics. piRNA annotation and references