Glucocorticoids Significantly Influence the Transcriptome of Bone Microvascular Endothelial Cells of Human Femoral Head.

BACKGROUND: Appropriate expression and regulation of the transcriptome, which mainly comprise of mRNAs and lncRNAs, are important for all biological and cellular processes including the physiological activities of bone microvascular endothelial cells (BMECs). Through an intricate intracellular signaling systems, the transcriptome regulates the pharmacological response of the cells. Although studies have elucidated the impact of glucocorticoids (GCs) cell-specific gene expression signatures, it remains necessary to comprehensively characterize the impact of lncRNAs to transcriptional changes.
METHODS: BMECs were divided into two groups. One was treated with GCs and the other left untreated as a paired control. Differential expression was analyzed with GeneSpring software V12.0 (Agilent, Santa Clara, CA, USA) and hierarchical clustering was conducted using Cluster 3.0 software. The Gene Ontology (GO) analysis was performed with Molecular Annotation System provided by CapitalBio Corporation.
RESULTS: Our results highlight the involvement of genes implicated in development, differentiation and apoptosis following GC stimulation. Elucidation of differential gene expression emphasizes the importance of regulatory gene networks induced by GCs. We identified 73 up-regulated and 166 down-regulated long noncoding RNAs, the expression of 107 of which significantly correlated with 172 mRNAs induced by hydrocortisone.
CONCLUSIONS: Transcriptome analysis of BMECs from human samples was performed to identify specific gene networks induced by GCs. Our results identified complex RNA crosstalk underlying the pathogenesis of steroid-induced necrosis of femoral head.

INTRODUCTION

Osteonecrosis of the femoral head (ONFH), of which avascular necrosis comprises the principal part, is a progressive and disabling condition in young and middle-aged patients. Typically, patients are asymptomatic in early stages while eventually requiring a total hip replacement (THR). Morbid vascularity (i.e., thrombophilia) in cancellous bones is etiological determinants of this disease. Vascular endothelial growth factor (VEGF) is a critical mediator of vasculogenesis and angiogenesis. It has been reported that erythropoietin protects the femoral head against glucocorticoids (GCs)-induced osteonecrosis in part by stimulating the expression of VEGF.[1] In support of these findings, VEGF receptor antibody can block angiogenesis and induce ONFH in animal models.[2] Being the most potent and effective therapy for immune-related diseases, GCs have been a mainstay of therapeutic recourse for most inflammatory disorders.[3] Among patients younger than 40 years of age, administration of GCs is the most prominent causative factor (about 60%) for ONFH.[4] Besides thrombophilia and other genetic propensities, GCs have been proven to be crucial in early stages of steroid-induced necrosis of femoral head (SINFH). Because their treatment results in reduced blood flow and hypercoagulability, GCs damage the bone microvascular endothelial cells (BMECs) of the femoral head.[5] Among several pathogenic extrapolations for SINFH, the vascular hypothesis-in that local microvascular dysfunction leads to decreased blood flow in the femoral head has become the most widely supported. Based on the specific microstructures of the femoral head, investigators have proposed both direct and indirect effects of GCs on BMECs. Those effects that GCs particularly have on BMECs would destroy the prop of osseous microcirculation, hereby leading to SINFH.[6]

In the current study, we assumed that GCs significantly impact the transcriptome of BMECs during the early stage of SINFH. We simulated SINFH of early stage by GCs-treatment in cultured BMECs and investigated their transcriptomes by microarray analyses, expecting to provide clues for remedying SINFH.

METHODS

Ethics statement and patients

The protocol described herein were approved by the ethics reviewing council of China-Japan Friendship Hospital, which abides by the Declaration of Helsinki on Ethical principles for medical research involving human subjects (Institutional Clinical Trials Register Number is 2014-34, Beijing 100029, China). Written informed consent was obtained from all participants. All patients were enrolled at China-Japan Friendship Hospital between October 1, 2013 and February 28, 2014, with diagnosis classified as traumatic transcervical fracture. Any traumatic fracture complicated with other etiological mechanisms (e.g., osteonecrosis, arthritis or osteoarthritis) or secondary fracture related to precipitating factors, such as osteoporosis, was excluded.

Primary cell culture and phenotyping

The current study was performed in Beijing Key Lab for Immune-Mediated Inflammatory Diseases, which is affiliated to China-Japan Friendship Hospital, China. To harvest BMECs for microarray analysis, cancellous bone of femoral head was obtained from patients who underwent routine THR surgery. The crunched cancellous bones were promptly shipped to the cell culture facility, where they were digested overnight with 0.2% of type one collagenase and subsequent trypsinized of 20 min. After centrifugation, BMECs were collected and plated in 100 mm gelatin-precoated culture dishes, cultured at 37°C. Complete culture medium was M199 medium (Thermo Fisher Scientific, Waltham, MA, USA) supplemented with 2 mmol/L L-glutamine, 100 μg/ml streptomycin and 100 IU/ml penicillin (Sigma-Aldrich), 20% fetal bovine serum, and recombinant human VEGF 165 (rhVEGF 165, Sino Biological Inc., Beijing, China). The medium was changed after 24 h to remove nonadherent cells. Once grown to 80% confluency, cells were detached by trypsinization and passaged to two duplicates of cultures with totally 106 cells per culture. Once the BMECs in each culture had grown to 80% confluence, one of the duplicate cultures was treated with hydrocortisone (Tianjin Kingyork Group Co. Ltd., Tianjin, China) to final concentration of 0.1 mg/ml for 24 h (i.e., GCs-treated experimental sample). The paired-control culture, which was from the same patient, was left untreated. In total, eight pairs of duplicate cultures from eight different patients were investigated in this study.

For analysis of cell surface markers, the cells were grown on sterile coverslips in a basal culture medium for 48 h, then washed and fixed in 4% paraformaldehyde following permeabilization with 0.1% Triton X-100. After blocking with 5% bovine serum albumin, cells were incubated overnight at 4°C with specific primary antibody, which was one of CD31 (PECAM-1), CD54 (ICAM-1), CD106 (VCAM-1) and CD144 (VE-cadherin). After washing off the primary antibody, the fluorescein isothiocyanate conjugated secondary antibody and hoechst 33342 were added and incubated for 1 h at room temperature. The slides were mounted using anti-fading mounting medium and digital slides were acquired using a fluorescence microscope IX71 (Olympus, Tokyo, Japan). For negative controls, an isotype IgG was used in place of a primary antibody.

Microarrays

Total RNA was extracted from cells using the Trizol reagent and purified with mirVana miRNA Isolation Kit (Ambion, Austin, TX, USA). The purity and concentration of RNA were determined from OD260/280 readings using spectrophotometer (NanoDrop ND-1000). RNA integrity was determined by capillary electrophoresis using the RNA 6000 Nano Lab-on-a-Chip kit and the Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). Only RNA extracts with RNA integrity number values >6 underwent further analysis. Higher yields of cDNA were labeled with a fluorescent dye (Cy5 and Cy3-dCTP) using CapitalBio cRNA amplification and labeling kit (CapitalBio, Beijing, China). The labeled cRNAs from lncRNAs and mRNAs were purified and hybridized to Agilent Human lncRNA + mRNA Array V3.0 (Agilent). The labeled and purified microRNAs were hybridized with Agilent Human microRNA Microarray Release 19.0 according to the manufacturer's instructions. Images were scanned with the Agilent microarray scanner, gridded, and analyzed using Agilent Feature Extraction software version 10.10 (Agilent). The raw data were summarized and normalized using the GeneSpring software V12.0 (Agilent).

Data analysis

To identify differentially expressed genes, we used threshold values of ≥2 absolute fold change and a Benjamini-Hochberg corrected P ≤ 0.05. The data were Log2 transformed and median centered by genes using the Adjust Data function of Cluster 3.0 software (open source). These genes were classified according to the GO analysis provided by Molecular Annotation System 3.0 (http://bioinfo.capitalbio.com/mas3/). Signaling pathway analysis of these genes was performed with KEGG PATHWAY Database (http://www.kegg.jp/kegg/pathway.html).

The lncRNAs were classified according to whether the corresponding transcript was mapped to the opposite strand of an exon, to an intron, or within an intergenic region of a protein-coding gene. The differentially expressed lncRNAs and mRNAs with Pearson correlation coefficients no. <0.99 were chosen to draw the co-expression network. The cis-acting lncRNA prediction was carried out based on their tight correlation (Pearson's correlation coefficient minimum of 0.99) to a group of differentially expressed protein-coding genes. The trans-prediction was conducted using blat tools (Standalone BLAT v. 35 × 1 fast sequence search command line tool) to compare the full sequence of the lncRNA with the 3’-untranslated region (3’ UTR) of its co-expression mRNAs. The target prediction were performed by their tight correlation (Pearson's correlation coefficient minimum of − 0.7) to a group of differentially expressed protein-coding genes. The correlated genes were intersected with an assembly of target genes predicted by anyone of 9 online programs, including miRanda (http://www.microrna.org/microrna/home.do), miRDB (http://mirdb.org/miRDB/), DIANAmT, etc.

Quantitative real-time polymerase chain reaction

To validate the microarray results, quantitative real-time polymerase chain reaction (qRT-PCR) was performed using the Lightcycler-Faststart DNA master SYBR green I PCR kit (Roche, Basel, Switzerland) according to the manufacturer's specifications. The comparative threshold cycle (CT) method was used for the calculation of amplification fold. To verify the microRNA expression profile, stem-loop microRNA qRT-PCR was performed using the small RNA U6 for an internal control.[7] The relative quantity of each microRNA in each sample, normalized to U6 RNA and relative to the expression in control samples, was calculated using the comparative CT (ΔΔCT) method, with the equation of RQ = 2−ΔΔCT, where ΔΔCT = (CTmicroRNA − CTU6) Exp − (CTmicroRNA − CTU6) Control

RESULTS

Expression profiles of mRNAs and lncRNAs in glucocorticoids-treated bone microvascular endothelial cells

Consistent with previous observations, a greater percentage of lncRNAs was detected above background compared to encoding ones.[8] We found that 78% of lncRNAs interrogated on the microarray exhibited expression above background and that 0.9% of these were significantly differentially expressed between GCs-treated and untreated BMECs. By contrast, 63% of protein-coding mRNA transcripts on the microarray were expressed above background, and 2.2% of these were significantly differentially expressed [Table 1]. In addition, statistical analysis showed that lncRNAs expressed above background were scattered on all chromosomes. The ratio (expressed probes/total probes) expressed from nuclear chromosome was lower than the mitochondrial genome, although few of them were differentially expressed [Figure 1a]. When we evaluated the expression levels of lncRNAs in paired samples (GCs-treated to untreated ratio) by log fold-change, 166 lncRNAs were found to be significantly down-regulated, and 73 lncRNAs were significantly up-regulated in the GCs-treated BMECs. There were more than twice as many lncRNAs down-regulated as there were lncRNAs up-regulated [Figure 1b]. Conversely, differentially up-regulated mRNAs (336 out of 518) were more common than significantly down-regulated mRNAs (182 out of 518). In line with our expectations, some annotated lncRNAs were differentially expressed in GCs-treated BMECs, such as NAV2-IT1, NAV2-AS5, NAV2-AS1, NTM-IT2, and ARHGEF19-AS1. We analyzed the genomic location – which partly reflects regulatory mechanisms - of the differentially expressed lncRNAs. We performed this analysis based on the lncRNA orientation to local protein-coding genes by categorizing their associated protein-coding genes as small nucleolar, divergent, intronic, antisense and intergenic according to our modified definition [Figure 1c]. All the records have been submitted to Gene Expression Omnibus (GEO) of National Center for Biotechnology Information, have been approved, and assigned GEO accession numbers referred as GSE60333 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE60333). This SuperSeries record provides access to two linked SubSeries that are GSE60332 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE60332) and GSE60093 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE60093), along with Platform GPL19072 for Agilent-052909 CBC_lncRNAmRNA_V3 (Probe Name version), which is initially submitted by us (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL19072).

Table 1

Summary of microarray analysis results

Probe class	Total (n)	Expressed above background (n,%)	Differentially expressed* (n,%)
lncRNAs	34,235	26,646 (78)	239 (0.9)
mRNAs	37,585	23,833 (63)	518 (2.2)
miRNAs	2027	368 (18)	5 (1.4)
Combined	73,847	50,847 (69)	762 (1.5)

*Significant differential expression was defined as probes with P≤0.05 and absolute fold-change ≥2.

Figure 1

Overall profiles of the glucocorticoid versus control transcriptomes of bone microvascular endothelial cells (BMECs). (a) Relative distribution of expressed lncRNAs (expressed probes/total probes) derived from each chromosome; (b) Relative chromosomal distribution of differentially regulated (up- or down-regulated probe number/total probes number) lncRNAs; (c) Annotation of genomic context of differentially expressed lncRNAs; (d) Validation of microarray results by quantitative real-time polymerase chain reaction.

To validate the microarray data, we investigated seven randomly selected RNAs (2 microRNAs, 2 lncRNAs and 3 mRNAs) which were identified to be differentially expressed after GCs administration using quantitative real-time PCR (qRT-PCR). The primers specific to these RNAs were designed based on their sequences. U6 served as an internal control for microRNA and GAPDH served as an internal control for lncRNA and mRNA detection, respectively. The qRT-PCR results showed that expression levels of ENSG00000259007.1 (probe name, p4493), XLOC_011117 (probe name, p19376), peroxisome proliferator-activated receptor gamma, BMP binding endothelial regulator and SMAD family member 3 were concordant with microarray data with a strong correlation [Figure 1d].

The microRNA transcriptome in bone microvascular endothelial cells after glucocorticoids treatment

As for microRNAs, the percentages for above-background expression and the differential expression were 18% and 1.4%, respectively [Table 1]. Hierarchical clustering analysis for the 16 sets of paired data (two data each patient) was summarized in a dendrogram [Figure 2a]. Eight sets of data from each group were combined for significance analysis of microarrays. After filtering out microRNAs that were undetectable, 368 microRNAs were reliably detected in control and GCs-treated groups. Of them, 5 microRNAs were up-regulated by more than 2.5-fold, while 11 microRNAs were down-regulated by more than 2.5-fold after 24-h treatment. Only 5 microRNAs among them were consistently up/down-regulated at least 2.5-fold at all 8 samples after treatment with corrected P ≤ 0.05 [Figure 2b]. Except for hsa-miR-339-5p that was up-regulated, all the others were down-regulated [Table 2].

Figure 2

Bioinformatic profiles of the differentially expressed microRNAs. (a) Hierarchical clustered heat maps showing the Log2 transformed expression values for differentially expressed microRNAs between glucocorticoids (GCs)-treated and untreated bone microvascular endothelial cells of eight patients; (b) Scatter and volcano plots of microRNAs with significant differential expression induced by GCs-treatment.

Table 2

Summary of microRNAs with significantly differentially expressed genes in GCs-treated versus control BMEC cells*

Systematic name	P (corrected)	FC (absolute)	Regulation	Active-sequence	Chromosome	Mirbase accession	Strand
hsa-miR-339-5p	0.47040328	2.565216	Up	CGTGAGCTCCTGGA	chr7	MIMAT0000764	+
hsa-miR-100-3p	0.36541834	5.2831144	Down	CATACCTATAGATACAAGCTT	chr11	MIMAT0004512	+
hsa-miR-222-5p	0.26147932	6.893519	Down	AGGATCTACACTGGCTA	chrX	MIMAT0004569	+
hsa-miR-23b-5p	0.32848376	4.3506756	Down	AAATCAGCATGCCAGGAACC	chr9	MIMAT0004587	−
hsa-miR-933	0.41049543	3.6403227	Down	GGGAGAGGTCTCCCT	chr2	MIMAT0004976	+

*Significant differential expression was defined as P≤0.05 and absolute fold-change ≥2.5. BMEC: Bone microvascular endothelial cell; GCs: Glucocorticoids; FC: Fold-change.

Correlation among differentially expressed microRNAs and lncRNAs with protein-coding genes

The microRNA target prediction was performed for all 5 microRNAs that were significantly influenced by GCs-treatment. We found a total of 55 target genes were predicted to be targeted by these microRNAs. We further expanded our analysis to identify protein-coding genes whose expression inversely correlated to these microRNAs. Some lncRNAs were also implicated in specific diseases by significantly correlating with genes participating in graft-versus-host disease, type 1 diabetes, arthritis, asthma, periodontitis, systemic lupus erythematosus, Graves’ disease, etc., Based on this dataset, we performed the enrichment analysis for GO terms, functionally annotating the genes by a conditional hypergeometric test. Gene functional annotations were classified into three categories: Biological process, molecular function (MF) and cellular component. Because microRNAs and lncRNAs were prone to exert influence on associated protein-coding genes in transcription-related processes, we submitted a list of protein-coding genes which were significantly correlated with differentially expressed noncoding RNAs that is, 172 mRNAs correlated with lncRNAs along with a list of 214 mRNAs correlated with differentially expressed microRNAs, and 518 statistically differentially expressed mRNAs to MAS 3.0 (http://bioinfo.capitalbio.com/mas3/) for GO term enrichment analysis. Some pivotal roles involving both mRNAs and mRNAs were significantly correlated with noncoding RNAs. Differentially expressed mRNAs and their counterparts statistically correlated with noncoding RNAs that were composed of microRNAs and lncRNAs were classified according to their GO annotation terms, respectively [Figure 3a and b]. There were similar portfolios of the highest enriched GO terms that were presented by both mRNAs and mRNAs that significantly correlated with noncoding RNAs.

Figure 3

The Gene Ontology (GO) enrichment analysis provides a controlled vocabulary to describe differentially expressed transcriptomes. These GO terms show a statistically enriched representation. The most significantly enriched terms for mRNAs (a) and the most enriched GO terms for the genes significantly correlated to differentially expressed noncoding RNAs (b) are displayed here. Enrichment Score of the GO equals (−log10 [P value]).

Correlated expression among different transcriptomes

We analyzed correlation for noncoding genes with significantly differential expression following GCs-treatment. For further analyses of protein-coding genes correlated with microRNAs, we constructed a gene interaction network [Figure 4a]. The correlation analyses identified 119 pairs of lncRNAs were expressed in a coordinated manner. A strong correlation coefficient was obtained among some of them. For example, ENST00000553603.1 (probed by p4493) not only correlated with many protein-coding genes, but also had positive interrelated relationships to TCONS_00024496 (probed by p20100) and ENST00000585190.1 (probed by p6957), suggesting the complex regulatory network which p4493 could be implicated [Figure 4b]. Notwithstanding the intricate interrelationship among some lncRNAs, no differentially expressed lncRNAs overlapped with any annotated microRNAs as far as our study is concerned. Through significantly correlating key genes, noncoding RNAs participate in some important processes and functions. In addition, some signaling pathways were conspicuously modulated by noncoding RNAs. Among these pathways, compensatory increased FoxO signaling pathway should be involved in the development of SINFH[9] [Figure 4c].

Figure 4

Co-expression network composed with differentially expressed transcriptomes. (a) Co-expression network composed with differentially expressed microRNAs and their correlated genes (Pearson's correlation coefficient minimum of 0.99). (b) Co-expression network composed with differentially expressed lncRNAs and their correlated genes (Pearson's correlation coefficient minimum of 0.99). (c) KEGG pathway analysis is a functional analysis that maps genes to pathways. The “Foxo Signaling Pathway” shows modulation in apoptotic signatures and may associate with onset of avascular osteonecrosis.

DISCUSSION

MicroRNAs extensively regulate protein synthesis at the posttranscriptional level. They are believed to interfere with their target gene-coding mRNAs by base-pairing and binding to the 3’-UTRs, 5’-UTRs or coding regions.[10] A large number of target genes involved in enriched GO terms were targeted by hsa-miR-100-3p in the present study. For example, by targeting interferon γ receptor 1, hsa-miR-100-3p is believed to be involved in periodontitis. In addition, hsa-miR-339-5p was believed to participate in osteoporosis, periodontitis and rheumatoid arthritis by targeting interleukin 11, while, hsa-miR-933 was related to osteoporosis via targeting cannabinoid receptor 1. Although few reports have proven their involvement in SINFH, these results highlight the need to further elucidate the role of microRNAs in SINFH. In some case, microRNAs can be processed from primary lncRNAs. Complicating this paradigm, however, is that lncRNAs can also be targeted by microRNAs.[11] These noncoding transcripts act as natural microRNA sponges by competing for binding to shared microRNAs.[12] To explore the possibility that some lncRNAs, which are differentially expressed after GCs-treatment, might also act as primary transcripts for microRNAs, we systematically searched for genomic overlap between differentially expressed lncRNAs and known microRNAs. Unfortunately, our study was unable to identify lncRNAs overlapping annotated microRNAs. This result suggested that there might not possibly be any primary precursor or endogenous competitor for microRNAs being significantly influenced by GCs.

The lncRNAs are considerably different from coding RNAs with regards to their pattern of splicing and the narrow distribution of GC content. They typically resided at genomic loci within 10 kb of a protein-coding gene along the genome.[13] The lncRNAs often participated in important, and critical cellular and biological functions regulating gene expression via epigenetic modification or transcriptional co-activation/repression.[14] Numerous studies have demonstrated that noncoding RNAs can regulate the expression of their associated protein-coding genes via a range of systems, including chromatin remodeling, alternative splicing, translational interference and promoter targeting by cis-acting (the same locus) regulatory mechanisms, which is an antisense-mediated regulation by lncRNAs of their neighboring protein-coding counterparts.[15] The high ratio of lncRNAs derived from the mitochondrial genome may be related to the high abundance of mitochondrial lncRNAs in certain tissues.[16] As for their genomic context, antisense transcripts are prevalent throughout the mammalian genome.[17] Among the antisense transcripts that were expressed above background, 51 of them were found to be differentially expressed. Intronic lncRNAs comprise a large class of transcriptional units transcribed from intronic regions of mammalian protein-coding genes in a tissue-specific pattern.[18] Previous studies indicated that intronic lncRNAs could regulate the expression of their host or neighboring genes.[19] We identified 31 intronic transcripts from differentially expressed lncRNAs. Long intergenic noncoding RNAs (lincRNAs) rapidly evolved their sequences despite conserved functions across species.[20] By modulating the expression of their neighboring protein-coding genes or other target genes scattered across the genome via the recruitment of histone-modifying enzymes to chromatin, lincRNAs play important roles in the development of certain diseases.[21] Most lncRNAs within our datasets could fit into this category that is, 126 lincRNAs that approximately made up 58% of the entire differentially expressed lncRNAs.

Except for the housekeeping MFs such as protein binding, the most relevant GO terms were enriched in transcription-related processes and some well-annotated transcriptional regulators belonged to the “cell cycle,” “apoptosis” and “cell proliferation” gene lists, including Gadd45b, Nav1, Fgf5, Lifr, Stat5a, Il11, Cdk10, Il6, Rgs2 and Gas1. It indicated that lncRNAs were likely to perform their function at transcriptional levels by regulating transcription-related genes. In contrast, GO terms related to “response to GC stimulus” or “GC receptor binding” were poorly enriched, which could be caused by the insufficient annotations for GCs-related GO terms. It is noteworthy that the search results for “GC” and “bone” in the GEO DataSets Database are scarce (http://www.ncbi.nlm.nih.gov/gds/?term=(Glucocorticoid+AND+bone)). Although there was few GO term enrichment for GCs-related genes, we still found that 2 of genes had associated GO terms relating to degenerative disorders, which were Fkbp4 and Il6.

In conclusion, we postulated that noncoding RNAs could be involved in microvascular events and related to endothelial dysfunction, hence diffusing apoptosis throughout the femoral head. SINFH impacts such a large population that preventing its attack is of critical significance to GCs users. In addition to GC-induced transcriptional changes, hereditary factors also predispose people to SINFH.[22] Combining the hierarchical clustering data with clinical meta-analysis, we can deduce that hereditary factors were much more impactful than iatrogenic circumstances. We identified significantly differentially expressed transcripts but nonetheless we are still short of a reliable clinical method to control this disorder. The current study outlined the differential expression profiles induced by GCs but didn’t specify any intracellular signaling pathway that is solely responsible for the pathogenesis of SINFH. All the aforementioned propel us to perform further in-depth investigation in order to enhance our understanding of SINFH. Intensive investigations of the signaling pathways screened out in this study should be carried out to verify the pivotal one underlying pathogenesis of SINFH.