Genome-wide analysis of differentially expressed long noncoding RNAs induced by low shear stress in human umbilical vein endothelial cells

Abnormal shear stress plays a critical role in development of atherosclerosis via modulating complex signal pathways. Long noncoding RNAs (lncRNAs) function directly as RNAs. However, the differential expression of lncRNAs and their net-work regulating the cellular response to the change of shear stress remains to be unclear. Human umbilical vein endothelial cells (HUVEC) were exposure to lower shear stress (2 dynes/cm2) for 120 minutes. Differential expression (defined as a fold change ≥2.0, P<0.05) of lncRNAs and mRNAs was obtained by microarray and analyzed by Boxand Scatter-plot, Heat Map and Hierarchical Clustering. Gene ontology and pathway analyses were performed thereafter. Of 149 lncRNAs differentially expressed between sheared and normal emdothelial cells, 64 lncRNAs were upregulated and 85 lncRNAs were downregulated. 109 mRNAs were differentially expressed (upregulated and 30 downregulated). We found 84 matched lncRNAmRNA pairs for 62 differentially expressed lncRNAs and 35 differentially expressed mRNAs. 9 pathways through which low shear stress acted on the regulation of cellular function were identified. According to pathway analysis, two up-stream genes (TLR4 and BDNF) were upregulated. Our study provides the expression profiles of lncRNA and mRNA in sheared cells (after 120 min), and reports that TLR4 and BDNF pathways plays crucial role in low shear stress-induced cellular injury, among 84 paired differentially expressed lncRNA-mRNA. Abbreviations: lncRNA: long non-coding RNA; EC: endothelial cell; HUVEC: Human umbilical vein endothelial cells; GO: gene ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; TLR 4: toll-like receptor 4; LSS: low shear stress Introduction Vascular endothelial cells (ECs) form the inner lining of blood vessel wall, are directly exposed to blood flow pattern, and serve important homeostatic functions in response to various chemical and mechanical stimuli [1,2]. Steady laminar flow (10-20 dyn/cm2) is mandatory to maintain ECs functional by releasing many athero-protective cytokines [3]. In contrast, lower (<10 dyn/cm2) or disturbed blood flow is atheroprone evidenced as atherosclerosis more frequently seen in the inner curvature and bifurcated area of coronary artery tree [4]. Previous studies have reported the differentially expressed genes when Ecs were exposed to lower shear stress or disturbed flow [57], with subsequently activation of oxidation, apoptosis, cell death and pro-inflammatory pathways [8,9]. However, the mechanisms regulating those activated genes induced by low shear stress still remains understudied. Recently, genome-wide analysis have unmasked the extensive transcription of large RNA transcripts that do not code for proteins, termed long noncoding RNAs (lncRNAs), and provide an important new perspective on the centrality of RNA in gene regulation [10]. These lncRNAs function directly as RNAs [11], recruiting chromatin modifiers to mediate transcriptional changes in processes ranging from X-inactivation to imprinting. Unfortunately, there is a lack of information about differential expression of lncRNAs and their network regulating the cellular response to the change of shear stress. Accordingly, the current study aimed to provide the genome-wide analysis data to show the differentially expression of lnc RNAs/mRNA and the pathways by which low shear stress acts on the endothelial cells. Materials and methods Cell culture Human umbilical vein endothelial cells (HUVEC) were cultured at 37°C in a humidified incubator with 5% CO2 and maintained in RPMI1640 culture medium supplemented with 10% fetal bovine serum. When grown to 90% confluence, cells were trypsinized, harvested, resuspended and seeded to a 0.1% gelatin-coated glass. After adherence, monolayer cells grown on the glass were imposed with LSS for different time lengths. Flow study A parallel-plate flow chamber that exerts continuous flow was made by sandwiching a silicon gasket between two stainless steel plates with a cover slip sink in the base plate. The chamber and all parts of the circuit tubes were sterilize before the placement of a glass (50 mm × 30 mm) with monolayer cells into the flow chamber. The parameter settings were as follows: flow chamber height, 0.56 mm and pump rate, Correspondence to: Shao-Liang Chen, Nanjing First Hospital, Nanjing Medical University, 68 Changle Road, Nanjing 210006, China, Tel: +86-25-52208048, Fax: +86-25-52208048; E-mail: chmengx@126.com


Introduction
Vascular endothelial cells (ECs) form the inner lining of blood vessel wall, are directly exposed to blood flow pattern, and serve important homeostatic functions in response to various chemical and mechanical stimuli [1,2]. Steady laminar flow (10-20 dyn/cm 2 ) is mandatory to maintain ECs functional by releasing many athero-protective cytokines [3]. In contrast, lower (<10 dyn/cm 2 ) or disturbed blood flow is atheroprone evidenced as atherosclerosis more frequently seen in the inner curvature and bifurcated area of coronary artery tree [4].
Previous studies have reported the differentially expressed genes when Ecs were exposed to lower shear stress or disturbed flow [5][6][7], with subsequently activation of oxidation, apoptosis, cell death and pro-inflammatory pathways [8,9]. However, the mechanisms regulating those activated genes induced by low shear stress still remains understudied.
Recently, genome-wide analysis have unmasked the extensive transcription of large RNA transcripts that do not code for proteins, termed long noncoding RNAs (lncRNAs), and provide an important new perspective on the centrality of RNA in gene regulation [10]. These lncRNAs function directly as RNAs [11], recruiting chromatin modifiers to mediate transcriptional changes in processes ranging from X-inactivation to imprinting. Unfortunately, there is a lack of information about differential expression of lncRNAs and their network regulating the cellular response to the change of shear stress. Accordingly, the current study aimed to provide the genome-wide analysis data to show the differentially expression of lnc RNAs/mRNA and the pathways by which low shear stress acts on the endothelial cells.

Cell culture
Human umbilical vein endothelial cells (HUVEC) were cultured at 37°C in a humidified incubator with 5% CO 2 and maintained in RPMI-1640 culture medium supplemented with 10% fetal bovine serum. When grown to 90% confluence, cells were trypsinized, harvested, resuspended and seeded to a 0.1% gelatin-coated glass. After adherence, monolayer cells grown on the glass were imposed with LSS for different time lengths.

Flow study
A parallel-plate flow chamber that exerts continuous flow was made by sandwiching a silicon gasket between two stainless steel plates with a cover slip sink in the base plate. The chamber and all parts of the circuit tubes were sterilize before the placement of a glass (50 mm × 30 mm) with monolayer cells into the flow chamber. The parameter settings were as follows: flow chamber height, 0.56 mm and pump rate, at least 3 out of 6 samples have flags in Present or Marginal were chosen for differentially expressed lncRNA or mRNAs screening. Differentially expressed lncRNAs and mRNAs were identified through fold change filtering.

Scatter plot
The Scatter-Plot is a visualization method used for assessing the LncRNA/mRNA expression variation (or reproducibility) between the two compared samples or groups. The values of X and Y axes in the Scatter-Plot were the normalized signal values of the samples or groups (log 2 scaled). The green lines were Fold Change Lines (The default fold change value given is 2.0). The LncRNAs/mRNAs above the top green line and below the bottom green line indicated more than 2.0 fold change of LncRNAs/mRNAs between the two compared samples or groups.

Box-Plot
The Box Plot is a convenient way to quickly visualize the distributions of a dataset. It is commonly used for comparing the distributions of the intensities from all samples. After normalization, the distributions of log 2 -ratios among all tested samples were assessed by this Box Plot.

Differentially expressed LncRNAs screening
To identify significant differentially expressed LncRNAs/mRNAs with statistical significance, we performed a Volcano Plot filtering between the two compared groups. The threshold was Fold Change ≥2.0 and p-value<0.05. Differentially expressed lncRNAs were defined as p value<0.005, FDR<15% and fold-change>2.0.

Heat map and hierarchical clustering
Hierarchical Clustering is one of the most widely used clustering methods for analyzing LncRNA/mRNAs expression data. Cluster analysis arranged samples into groups based on their expression levels, which allowed us to hypothesize the relationships among samples. The dendrogram showed the relationships among LncRNA/mRNA expression patterns of samples.

Gene ontology (GO) analysis
GO analysis is a functional analysis associating differentially expressed genes with GO categories, including biological process, cellular components and molecular function. The GO categories are derived from Gene Ontology (http://www.geneontology.org), which comprise three structured networks of defined terms that describe gene product attributes. The P-value denotes the significance of GO Term enrichment in the differentially expressed gene list. The less the P-value is, the more significant the GO Term is (P-value <= 0.05 was recommended).

Pathway analysis
The predicted target genes above were input into the Database for Annotation. Base on the latest KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www.genome.jp/kegg) database, we performed pathway analysis for differentially expressed mRNAs. This analysis allowed to determine the biological pathway that there is a significant enrichment of differentially expressed genes. Furthermore, we also used the Visualization and Integrated Discovery (DAVID; http://david. abcc.ncifcrf.gov/) and BioCarta (http://www.biocarta.com) to analyze the potential functions of these target genes in the pathways. The P-value denotes the significance of the Pathway. The less the P-value is, 60 times/min. The value of shear stress was adjusted by modulating the flow of after-loading and was automatically calculated using information from the pressure transducers. The shear stress in our flow study was 2 dynes/cm 2 (for 120 min., marked by ABCD 120, EFGH 120 and IJKL 120) and the cyclic fluid used was 1640 medium with 2% fetal bovine serum supplied. Cultured HUVECs without shear stress were considered as control group (marked by ABCD 0, EFGH 0 and IJKL 0 in Figure 1).

RNA extraction
Total cellular RNA was isolated from the cells using an RNeasy Mini Kit (Qiagen, Hilden, Germany) in accordance with the manufacturer's protocol, and then quantified using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA). The RNA integrity of each sample was assessed using standard denaturing agarose gel electrophoresis after RNA extraction and prior to sample labeling.

Labeling and hybridization
Total RNA (200 ng each) from cells was purified from total RNA after removal of rRNA (mRNA-ONLY™ Eukaryotic mRNA Isolation Kit, Epicentre). Then, each sample was amplified and was reversely transcribed into cDNA using an RNA Spike In Kit with one-color (Agilent) in the presence of 0.8 µL of Random Primer and 2 µL of Spike Mix without 3' bias utilizing a random priming method. These cDNA samples were then cleaned and labeled in accordance with the Agilent Gene Expression Analysis protocol using Agilent Quick-Amp Labeling Kit. These labeled cDNA samples were used as probes to hybridize to microarrays for 17 h at 65°C using an Agilent Gene Expression Hybridization Kit in Agilent's SureHyb Hybridization Chambers. The labeled cRNAs were hybridized onto the Human LncRNA Array v3.0 (8 x 60K, Arraystar). After hybridization, the microarrays were washed with an Agilent Wash Buffer kit (Agilent) and scanned with an Agilent microarray scanner G2505C.

lncRNA and mRNA microarray
Arraystar Human LncRNA Microarray V3.0 (Kangcheng Co., Shanghai, China) was designed for the global profiling of human LncRNAs and protein-coding transcripts, with 30,586 LncRNAs and 26,109 coding transcripts being detected. The LncRNAs were carefully constructed using the most highly-respected public transcriptome databases (Refseq, UCSC known genes, Gencode, etc), as well as landmark publications. Each transcript was represented by a specific exon or splice junction probe which can identify individual transcript accurately. Positive probes for housekeeping genes and negative probes were also printed onto the array for hybridization quality control. Each transcript was represented by using 1-5 probes to improve statistical confidence.

Data analysis
The microarrays were scanned at 5 µm/pixel resolutions using an Agilent microarray scanner piloted by GenePix Pro 6.0 software (Axon). Scanned images (TIFF format) were then imported into Agilent Feature Extraction software for grid alignment and expression data analysis. Expression data were normalized by a quantile normalization and the Robust Multichip Average (RMA) algorithm that was included in the Agilent software. Probe-level files and mRNA-level files were generated after normalization. All gene-level files were imported into Agilent GeneSpring GX software (version 11.5.1) for further analysis. After quantile normalization of the raw data, lncRNA or mRNAs that in endothelial cells after treated by low shear stress for 120 min and normal endothelial cells (Figure 1). Using the authoritative data sources, we found that 149 lncRNAs were differentially expressed (fold change ≥2.0, P<0.05) between sheared and normal emdothelial cells. Among them, 64 lncRNAs were upregulated (more than two-fold in sheared vs. the normal cells) and 85 lncRNAs were downregulated (more than two-fold; P<0.05; Table S1). Furthermore, we performed a step-wise genome analysis according to p value. At p value<0.01(more than two-fold), 29 lncRNA were upregulated, with 33 downregulated. When strict criterion (p value<0.005) was used, only 5 lncRNA were upregulated and 8 lncRNA were downregulated (more than two-fold). Finally, when our strict criteria (p value<0.005, FDR<15% and fold-change>2.0) were used, upregulated lncRNAs were C17orf76-AS1, BC079832, and LEF1-AS1; downregulated lncRNA included XLOC-007914 and RP11-473M20.16 (Table S1).

Differentially expressed mRNAs in endothelial cells
In a similar way, we found that 109 mRNAs were differentially expressed (more than two-fold, p<0.05) between sheared and normal endothelial cells (Table S2). Of them, 79 mRNAs were upregulated the more significant the Pathway is (The P-value cut-off is 0.05).

Statistical analyses
All data were expressed as mean ± standard deviation. Statistical analysis was performed using Student's t-test to compare two variables of microarray data. For example, the statistical significance of a microarray result was analyzed by fold change, and a difference with P<0.05 was considered statistically significant. The false discovery rate was also calculated to correct the P value. The threshold value we used to screen differentially expressed lncRNAs and mRNAs was a fold change ≥2.0 (P<0.05). Furthermore, a differential expression of each lncRNA between sheared and the non-sheared endothelial cells was analyzed using Student's t-test with SPSS (Version 16.0 SPSS, Chicago, IL, USA). P<0.05 was considered significant.

Differentially expressed lncRNAs in endothelial cells
To profile differentially expressed lncRNAs in endothelial cells, we performed a genome-wide analysis of lncRNA and mRNA expression and 30 downregulated (p<0.05). When p value was designed as less 0.01, 9 mRNAs were upregulated and 7 downregulated, respectively. Finally, by our stric criterion (p<0.005), 6 mRNAs were upregulated and 5 downregulated. We found 84 matched lncRNA-mRNA pairs for 62 differentially expressed lncRNAs and 35 differentially expressed mRNAs (Table S3).

GO and analyses of net-work
Through GO analysis we found that these under-regulated and upregulated transcripts of lncRNAs were associated with cellular process (ontology: biological process), cell (ontology: cellular component), and binding (ontology: molecular function) (Figure 2 and Table S4). Pathway analysis determined that these lncRNAs may target 23 gene pathways that corresponded to transcripts, i.e., "Intracellular Parts" and "Intracellular" composed of 23 targeted genes (with the recommended P-value <0.05) (Figure 2 and Table S5).
Based on these GO data, we then performed the analysis of net-work of differentially expressed lncRNA ( Figure S1), classified by biological process, cellular components and molecular function, respectively.

Pathway analysis
Subsequently, we performed the enrichment score analysis ( Figure  3) and pathway analysis. We found that 9 pathways through which low shear stress acted on the regulation of cellular function (Table S5). Of these 9 pathways, 2 virus-related (hsa 05162 for measeles, hsa 05161 for hepatitis) and 1 non-small lung cancer-related (has 05223) were excluded from the analysis. Vascular smooth muscle contraction-, adrenergic signaling in cardiomyocytes-and neurotrophin-related signal pathways were also not analyzed from this analysis. Finally, NF-kB, MAPK and Jak-STAT pathways were considered to be as major pathways by which low shear stress regulated endothelial cells, with 9 genes (BDNF, DUSP4, GADD45B, MAP2K5, PLA2G4C, DDX58, GADD45B, TLR4, IL6ST, SPRY1, STAT5A) involved.
According to pathway analysis, two up-stream genes (TLR4 and BDNF) were upregulated (Figure 4 and Table S5). TRL4, one cell membrane located protein, was supposed to be the LSS receptor. As shown in Figure 4, low shear stress significantly upregulated TLR4 and NF-kappa B associated genes. Similarly, low shear stress activated Neurotrophin signaling pathway which also exists in ECs and MAPK pathway both via activation of BDNF. BDNF, a factor that has been reported to be produced by ECs, is a crucial role in the progression of AS.

Discussion
Our study provides the expression profiles of lncRNA and mRNA in sheared cells (after 120 min) and analyzes the possible mechanisms of cellular responses to low shear stress (2 dyn/cm 2 ). The results show that TLR4 and BDNF pathways plays crucial role in low shear stressinduced cellular injury, among 84 paired differentially expressed lncRNA-mRNA.
There is a great interest in the biological function of RNA. According to high-throughput analysis [10,11], regulatory non-coding RNA (ncRNA) can be classified into short (<30 nt in length ), medium (20-300 nt in length) and lncRNA (>200 nt in length). To date, the majority of lncRNA are likely to be functional, but only a relative small proportion has been demonstrated to be biologically relevant, such as regulation of basal transcription mechinery functions in cells [12]. Endothelial cells of different origin express relative high levels of the conserved long noncoding RNAs [13]. Several lncRNAs have been reported to be correlated with atherosclerosis and endothelial cellular dysfunction [14][15][16]. However, there is a lack of differentially expressed lncRNAs in endothelial cells exposed to low shear stress. Here, we for the first reported that three lncRNAs (C17orf76-AS1, BC079832, and LEF1-AS1) were upregulated and two lncRNAs (XLOC-007914 and RP11-473M20.16) were downregulated after endothelial cells were treated by low shear stress (2.0 dyn/cm 2 ) for 120 min.
C17orf76-AS1, also known as LRRC75A-AS1, dynamically conveyed information from the activated epidermal growth factor and extracellular signal-regulated kinase (EGF/ERK) signaling cascade to directly associated proteins and more distant proteins in the network [17,18], which is consistent with our previous finding that low shear stress activated ERK/MAPK signal pathway [9]. LRRC75A-AS1 and BC079832 mainly regulated zinc finger and BTB domain-containing protein 37(ZBTB37) gene which was involved in the regulation of MAPK signals [17][18][19].
Subsequently, we measured the expression of mRNA in endothelial cells after treated by low shear stress. We found that TLR4 and BDNF, two of major membrane genes, were significantly regulated. TLR4 associates with the activation of NF-κB signaling and plays as monitor of NF-κB [22,23]. Numerous studies have demonstrated that TLR4/ NF-κB system participated in the formation and acceleration of AS [24][25][26][27]. Yang et al. identified increased TLR4 expression, activated NF-kB signaling in LSS-treated cell and supressed inflammation in TLR4 mutant [24]. Michelsen also reported that genetic deficiency of TLR4 led to a significant reduction in aortic plaque areas in apo-E deficient mice [27]. Triggering TLR4 pathway leads to the activation of NF-κB signaling and subsequent regulation of inflammatory genes. In this study, we found that LSS significantly promoted the expression of TLR4, DDX58, GADD45B, and these genes were involved in NF-κB signaling. Moreover, our study indicated firstly that TLR4 is the membrane located LSS receptor and initiation of LSS-induced endothelium injury.
We assayed BDNF gene in sheared cells and found an increase and we proposed BDNF contributed to the progression of AS. One decade ago, Nakahashi assayed ECs-produced BDNF under high shear stress (24 dynecs/cm 2 ) and found an decrease [28]. Haplo-deficient BDNF-receptor was associated with reduced AS, smooth muscle cells proliferation and collage accumulation in the lesion [29]. These evidences indicate that BDNF is involved in the development of AS. However, ECs derived BDNF depends on endothelial function rather than the artificial stimulations [29]. Others also reported reduced BDNF in AS patients' serum [30]. Even though the precise role of BDNF is not very clear, it is not advise to deny the association of BDNF with AS, which needs to be further explored.

Conclusions
The present study reported the differential expression of lncRNAs and mRNAs in endothelial cells exposure to low shear stress for 120 minutes. Pathway analysis demonstrated that BDNF and TLR 4 genes played important role in modulating cellular injury induced by low shear stress.

Competing interest
The authors declare that they have no competing interests Author's contribution CSL designed this study and wrote this manuscript. WZM carried out the shear stress experiments, molecular genetic studies, participated in the sequence alignment and drafted the manuscript. HZY carried out the genetic studies. LB participated in the design of the study and performed the statistical analysis. All authors read and approved the final manuscript.