Aberrantly Expressed lncRNAs in Primary Varicose Great Saphenous Veins

Long non-coding RNAs (lncRNAs) are key regulatory molecules involved in a variety of biological processes and human diseases. However, the pathological effects of lncRNAs on primary varicose great saphenous veins (GSVs) remain unclear. The purpose of the present study was to identify aberrantly expressed lncRNAs involved in the prevalence of GSV varicosities and predict their potential functions. Using microarray with 33,045 lncRNA and 30,215 mRNA probes, 557 lncRNAs and 980 mRNAs that differed significantly in expression between the varicose great saphenous veins and control veins were identified in six pairs of samples. These lncRNAs were sub-grouped and mRNAs expressed at different levels were clustered into several pathways with six focused on metabolic pathways. Quantitative real-time PCR replication of nine lncRNAs was performed in 32 subjects, validating six lncRNAs (AF119885, AK021444, NR_027830, G36810, NR_027927, uc.345-). A coding-non-coding gene co-expression network revealed that four of these six lncRNAs may be correlated with 11 mRNAs and pathway analysis revealed that they may be correlated with another 8 mRNAs associated with metabolic pathways. In conclusion, aberrantly expressed lncRNAs for GSV varicosities were here systematically screened and validated and their functions were predicted. These findings provide novel insight into the physiology of lncRNAs and the pathogenesis of varicose veins for further investigation. These aberrantly expressed lncRNAs may serve as new therapeutic targets for varicose veins. The Human Ethnics Committee of Shanghai East Hospital, Tongji University School of Medicine approved the study (NO.: 2011-DF-53).


Introduction
Varicose veins, whose manifestations can vary from leg edema to chronic, disabling venous ulceration, affect around 25% of the adult population and can lead to considerable morbidity and congestion of health service resources [1]. Great saphenous veins (GSVs) or saphenofemoral junction accounts for about 70% of varicose veins [1][2][3]. Many risk factors, such as increased age, female gender, diet rich in processed foods, obesity, work requiring large amounts of standing, tight undergarments, cigarette smoking, hypertension and diabetes mellitus, contribute to the prevalence of GSV varicosity [1,3]. However, the molecular mechanisms underlying the connection between these factors that lead to the prevalence of varicose veins remain obscure.
Recent progress addressing the molecular mechanism of GSV varicosities includes venous valvular dysfunction-causing reflux [4], leukocyte diapedesis and local inflammation, smooth muscle cell apoptosis and proliferation, endothelial cell injury, increased vein wall tension, vein wall dilation, extracellular matrix degra-dation and tissue remodeling [3,5,6]. Many molecules from many different pathways are involved in the pathological processes associated with varicose veins. These molecules include hypoxiainducible factor 1-alpha (HIF-1alpha) in the hypoxia pathway [7], the Janus kinase/signal transducers and activators of transcription (JAK-STAT) and nuclear factor kappaB (NF-kappaB) in the inflammatory pathway [8], poly ADP ribose polymerase (PARP) and bax in the apoptotic pathway [9], adhesion molecules and cytokines such as intercellular adhesion molecule 1 (ICAM-1), interleukin-1 alpha (IL-1alpha), and tumor necrosis factor alpha (TNF-alpha) [10]. However, these molecules account for only some of the tangled mechanisms, and many more remain to be identified. In addition, molecules such as long non-coding RNAs (lncRNAs) and micro-RNAs whose concentrations are correlated with those of their targets at the mRNA level, post-transcriptional level, or protein level must be systematically screened and validated.
LncRNAs are a class of transcripts whose lengths exceed 200 nt [11]. They are found throughout the genome. LncRNAs play key regulatory roles in regulating transcription in both cis form and antisense form, localization of proteins, and organizational frameworks of sub-cellular structures. LncRNAs also posttranscriptionally control mRNAs by affecting splicing, editing, translation, and degradation. In addition, many lncRNAs are processed into small RNAs or, modulate how other RNAs are processed. It is becoming increasingly clear that lncRNAs function in several different ways and play key roles in many intracellular regulatory processes [12,13]. The abnormal regulation of lncRNAs is involved in several human diseases, such as cancer [14][15][16], Alzheimer's disease [17], spinocerebellar ataxia (SCA) [18] and cardiovascular disease [19][20][21][22][23][24][25]. However, the relationship between lncRNAs and varicose veins is still unclear.
The present study was designed for screening and validation of lncRNAs at the micro-array level and examination of their relationship with the prevalence of varicose veins mRNAs and mRNA pathways associated with the aberrantly expressed lncRNAs, were identified at the whole micro-array level. These mRNAs and pathways may be relevant to the incidence of GSV varicosity.

Patients and Tissue Samples
Thirty-two samples of human primary GSVs were retrieved from 32 patients (14 men, 18 women) who were undergoing GSVs varicose vein excision in Shanghai East Hospital, Tongji University School of Medicine, China. The diagnosis of primary varicose GSVs was based on the clinical signs and duplex ultrasound scanning. All patients were characterized as having primary varicosities, and patients with secondary varicosities were excluded. None of the participants had any history of deep venous thrombosis, superficial thrombotic phlebitis, post-thrombotic syndrome, Klippel-Trenaunay syndrome, May-Thurner syndrome or any other venous disease. The clinical, etiological, anatomical and pathological elements classification system (CEAP) was used to classify chronic lower-extremity venous disease (CVD) in this case [26,27]. The patients' clinical signs placed all of them in classes 4-6, 30 of them were in class 4. Preoperative lowerextremity venous duplex ultrasound scanning assessment was performed on all patients, and examinations of both the superficial and the deep venous systems examination were conducted. All patients exhibited reflux in the GSVs. Patient demographics and clinical risk factors are given in Table 1.
Paired tissues were used to evaluate the differences in expression level between varicose veins (VVs) and adjacent normal segments of saphenous veins (NVs). According to visual inspection and pathological examination, VV was the obvious varicose vein, and NV was the adjacent normal vein. NV samples were collected about 3-4 cm from the VV area. The tissues were snap-frozen into liquid nitrogen immediately after resection for later RNA extraction. Six pairs of samples were taken from six patients from 32 patients were used for lncRNAs expression microarray and all 32 samples were used for quantitative real-time PCR (Q-RT PCR) validation. Written informed consent was obtained from all participants. The study was approved by the Human Ethnics Committee of Shanghai East Hospital, Tongji University School of Medicine (NO.: 2011-DF-53).

RNA Extraction and RNA Quantity Control
Total RNA was extracted from 32 pairs of samples which had been snap frozen using TRIzol reagent (Invitrogen, Carlsbad, CA, U.S.) according to the manufacturer's protocol. The amount and quality of RNA were measured using NanoDrop ND-1000 and RNA integrity was assessed by standard denaturing agarose gel electrophoresis.

RNA Labeling and Microarray Hybridization
Sample labeling and microarray hybridization were performed using a modified version of the Agilent One-Color Microarray-Based Gene Expression Analysis protocol (Agilent Technology). Briefly, rRNA was removed from the total RNA sample, and then the mRNA was purified (mRNA-ONLY TM Eukaryotic mRNA Isolation Kit, Epicentre). Then, each sample was amplified and transcribed into fluorescent cRNA along the entire length of the transcripts without 39 bias using a random priming method. The labeled cRNAs were purified using an RNeasy Mini Kit (Qiagen). The concentration and specific activity of the labeled cRNAs (pmol Cy3/mg cRNA) were measured using NanoDrop ND-1000. Then 1 mg of each labeled cRNA was fragmented by adding 11 ml 10 6 Blocking Agent and 2.2 ml of 25 6 fragmentation buffer. Then the mixture was heated at 60uC for 30 min. Then 55 ml 2 6 GE hybridization buffer was added to dilute the labeled cRNA. Then 100 ml of hybridization solution was dispensed into the gasket slide and placed in the Arraystar Human LncRNA Array v2.0 with 33,045 LncRNAs were collected from the authoritative data sources including RefSeq, UCSC Knowngenes, Ensembl and many related studies. The slides were incubated for 17 h at 65uC in an Agilent Hybridization Oven. The hybridized arrays were washed, fixed, and scanned with using the Agilent DNA Microarray Scanner (part number G2505B). The microarray work was performed by KangChen Bio-tech (Shanghai). The microarray data discussed in this article have been deposited in National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) and are accessible through (GEO) Series accession number GSE51260 (http://www.ncbi.nlm.nih. gov/geo/query/acc.cgi?acc = GSE51260).

Q-RT-PCR
Total RNA was extracted from frozen vein specimens using TRIzol reagent (Invitrogen Life Technologies) and then reverse transcribed using a PrimeScript TM RT Reagent Kit (Takara) according to the manufacturer's instructions. LncRNA expression in VV and paired NV tissues was measured by Q-RT-PCR using Power SYBRH Green PCR Master Mix (Applied Biosystems) on the ABI PRISMH 7900 Sequence Detection System (SDS) instrument. The Q-RT-PCR primers were designed using Primer 3.0 and blasted specifically in NCBI. Then 1 mg of total RNA was converted to cDNA according to the manufacturer's protocol. PCR was performed in a total reaction volume of 10 ml, including 5 ml SYBR Green PCR Master Mix (26), 0.4 ml each of PCR forward and reverse primer (10 mM), 0.5 ml of cDNA, and 3.7 ml of double-distilled water. The quantitative real-time PCR reaction was set at an initial denaturation step of 10 min at 95uC followed by 40 cycles of 95uC for 15 s and, 60uC for 1 min. All experiments were performed in triplicate.

Data Computational Analysis
For lncRNAs and mRNA microarray screening analysis, Agilent Feature Extraction software (version 10.7.3.1) was used to analyze acquired array images. After quantile normalization of the raw data, lncRNAs and mRNAs were chosen for further data analysis. A volcano plot was used to filter lncRNAs/mRNAs that were differentially expressed between two groups and Hierarchical clustering was performed using the GeneSpring GX v12.0 software package (Agilent Technologies). LncRNAs that were significantly differently expressed between all six pairs of VV and NV tissues were selected for Q-RT-PCR validation (P,0.05, $2 fold-change), as were any neighbor mRNAs that were also significantly differently expressed between all six pairs of VV and NV tissues (P,0.05, $2 fold-change).
For Q-RT-PCR validation analysis, all samples were normalized to GAPDH. The mean value in each triplicate was used to calculate relative lncRNAs concentrations (DCt = Ct mean lncRNAs 2Ct mean GAPDH ). Expression fold changes were calculated using 2 2DDCt methods. The differences in the level of lncRNAs expression between VVs and NVs were analyzed using the Student's t-test and SPSS (Version 16.0, SPSS Inc.) with the value of P,0.05 was considered as statistically significant.
Gene Ontology (GO) analysis and pathway analysis were performed to identify differentially expressed mRNA pathways. GO analysis showed significantly up-regulated and down-regulated mRNAs in question to be related to biological processes (BP), cellular components (CC) and molecular functions (MF) (P,0.05) The latest version of the Kyoto Encyclopedia of Genes and Genomes (KEGG) database and GO categories derived from Gene Ontology (www.geneontology.org) were used for pathway analysis, which was performed using the standard enrichment computation method [28]. A coding-non-coding gene co-expression network (CNC network) was drawn using Cytoscape with Pearson coefficient (|r|) .0.99 [29,30]. The analysis work was performed by KangChen Bio-tech (Shanghai P.R. China).

Profiles of the Differently Expressed lncRNAs and mRNAs
Among the 33,045 lncRNAs and 30,215 coding transcripts probed in the microarray, 12,264 lncRNAs and 14,862 mRNAs were detected in all six pairs of samples. There was an average of 2,426 differentially expressed lncRNAs (Table S1) (ranging from 1,442 to 4,049 depending on the subjects) and 3,029 differentially expressed mRNAs (Table S2) (ranging from 1,914 to 5,063). Among them, 557 lncRNAs (302 up-regulated and 255 downregulated) and 980 mRNAs (239 up-regulated and 741 downregulated) were significantly differently expressed ($2 fold) ( Table 2). Among the 20 most significantly differentially expressed lncRNAs (Table 3), BC041954 (fold change< 23.1) was the most significantly up-regulated lncRNA and AK023929 was the most significantly down-regulated (fold change = 4.84). There were more down-regulated lncRNAs than up-regulate lncRNAs.
The GO analytical data of aberrantly expressed mRNAs are shown in Table S3 (for biological processes), Table S4 (for cellular components), and Table S5 (for molecular functions). The top 20 most significantly differentially expressed mRNAs are shown in Table 4, and the top 5 most enriched GO terms, specifically biological processes, cellular components and molecular functions, are shown in Table 5. Pathway analysis showed these 123 significantly differentially expressed mRNAs to be involved in the pathogenesis of GSV varicosities (Table S6). Among them, 26 down-regulated and 2 up-regulated pathways were detected and 7 pathways were found to contain 46 mRNAs (changed $2 fold, enrichment score (2log10 (P -value)) value .2, P,0.01). Among these pathways, six (glycolysis/gluconeogenesis, fatty acid metabolism, tyrosine metabolism, pyrimidine metabolism, peroxisome and maturity-onset diabetes of the young) were focused on metabolic pathways. Differently Expressed lncRNA Subgroups Overall, 48 enhancer-like lncRNAs were found to be differentially expressed (Table S7) using with Gencode annotation [31,32]. Their nearby significantly differentially expressed coding RNAs (distance,300 kb) are shown in Table  S8. The Arraystar Human LncRNA Array v2.0 microarray contains profiling data of all probes in the 4 HOX loci, targeting 407 discrete transcribed regions, lncRNAs and coding transcripts. Among them, 10 lncRNAs and 10 coding transcripts were found to be differentially expressed in the human HOX loci of VVs and those of paired NVs (Table S9). Seventy-eight Rinn's lincRNAs (Table S10), 14 lincRNAs, and 19 nearby coding RNAs found to be significantly differentially expressed (Table S11) [33,34].

Validation of Aberrantly Expressed Candidate lncRNAs
LncRNAs significantly differently expressed between all six pairs of VV and NV tissues with their neighbor mRNAs significantly differently expressed between all six pairs of VV and NV tissues were selected for Q-RT-PCR validation, and 22 lncRNAs were validated. Nine of the 22 lncRNAs whose Q-RT-PCR primers were suitably designed (using Primer 3.0 and blasted specifically in NCBI). Primers are shown in Table S12. These lncRNAs were then replicated with Q-RT-PCR. Of these nine lncRNAs, seven showed the same trends of up-and downregulation as the microarray data and six were statistically significant (P,0.05) (Figure 1), supporting a strong consistency between the Q-RT-PCR results and the microarray data.

Prediction of the Functions of the Validated lncRNAs
To study the relationship between the lncRNAs and mRNAs more visually, the validated six significantly differently expressed lncRNAs were used to establish coding-non-coding gene (CNC) networks. These CNC networks were used to search for correlations between the lncRNAs and mRNAs and to determine the potential function of the lncRNAs. This would increase understanding of lncRNAs biological networks and of the complex pathogenesis of GSV varicosities. Overall 11 significantly aberrantly expressed mRNAs were found to be correlated with four validated lncRNAs. They were used to construct four separate networks with lncRNAs in the center node of the hub. The mRNAs CHAT and TMEM38B were found to be related to lncRNA AF119885; the mRNAs CCNO, EPC2, FAM13C and SHOC2 were found to be related to lncRNA G36810; the mRNAs EMX1 and SMC3 were found to be related to lncRNA NR_027927; and the mRNAs ATXN7, HOXC4,and RTCD1 were found to be related to lncRNA uc.345-. These separate networks are shown in Figure 2 as a whole. Further information (standard deviation (SD), and false discovery rate (FDR)) regarding these 11 mRNAs is shown in Table 6.
The potential functional effects of the six validated candidate lncRNAs were predicted through their correlations with the aforementioned six GSV-varicosities -related mRNA metabolic pathways (glycolysis/gluconeogenesis, fatty acid metabolism, tyrosine metabolism, pyrimidine metabolism, peroxisome and maturity onset diabetes of the young) with a criteria of coexpression gene pairs with |r|.0.9 and P,0.01. Here, eight mRNAs were found to be correlated with specific lncRNAs (uc.345 with NME7; G36810 with POLR3F, IDH1, and ALDH2; AF119885 with ENTPD1, SLC25A17, and PDHA1; NR_027927 with NEUROG3). Table 7 shows the mRNA pathways and pathway IDs.

Discussion
In the present study, the profiles of aberrantly expressed lncRNAs, mRNAs, and lncRNA-related mRNAs and pathways in primary varicose GSVs were evaluated. Six potential pathogenic lncRNAs were identified by Q-RT-PCR. The potential functional effects of these lncRNAs were predicted. This is the first systemic screening and validation of the GSV varicosities related lncRNAs in the genome-wide RNAs expression level.
The aberrantly expressed lncRNA patterns were classified into different subgroups with different analytical methods. Forty-eight enhancer-like lncRNAs were found to be differently expressed using GENCODE annotation, 10 discrete transcribed regions were targeted by 4 HOX loci, and 78 Rinn's lincRNAs were identified using chromatin-state maps. Enhancer-like lncRNAs were associated with decreased expression of their neighboring genes through loss-of-function effects [31]. A total of 407 discrete transcribed regions including exons (101), introns (75), and intergenic transcripts (231) were identified in the 4 HOX loci [35]. Rinn's lincRNAs profiles showed strong purifying selection, clear evolutionary conservation, and possible functions in many different biological processes [33,34]. The aberrantly expressed lncRNAs observed may provide clues to the pathophysiological properties of GSV varicosities. Six validated candidate lncRNAs (AK021444, AF119885, uc.345-, NR_027927, G36810, NR_027830) were identified in GSV varicosities by Q-RT-PCR and their potential functional effects were predicted. AK021444 is a 1611 bp lncRNA with an exon sense-overlapping relationship with POSTN. POSTN encods the protein periostin, which promotes cardiac repair and cardiomyocyte proliferation [36] and induces vascular cell differentiation and migration during the repair of vascular injury [37]. AF119885 is a 1269 bp lncRNA transcribed from the natural   HOXC4 belongs to the homeobox family of genes, which encode a highly conserved family of transcription factors that play an important role in morphogenesis in all multicellular organisms. NR_027927 is a 2249 bp lncRNA with a bidirectional relationship with the gene CAPN7. CAPN7 protein is an intracellular, nonlysosomal cysteine protease whose activity is regulated by calcium influx and oxidative stress [38]. G36810 is a 401 bp lncRNA transcribed from the natural antisense strand of gene TMEM47. This gene encodes a member of the PMP22/EMP/ claudin protein family, which is localized to the ER and the plasma membrane. NR_027830 is a 5928 bp lncRNA with an exon sense-overlapping relationship with the gene BCAP29. BCAP29 may be involved in the anterograde transport of membrane proteins from the endoplasmic reticulum to the Golgi and in CASP8-mediated apoptosis. These data suggest that the six candidate lncRNAs may function in the etiology and pathogenesis of primary varicose GSVs, but further study is required to confirm this.
In this study, bio-informatics strategy was used to show that many lncRNAs expressions are correlated with that of nearby mRNAs expressions (Table S13). This is consistent with one proposed mechanism that lncRNAs are associated with related mRNAs [12]. Four of the six validated lncRNAs were associated with 11 significantly differentially expressed mRNAs. These mRNAs were involved in various biological processes. The biological functions of these genes which transcribed the 11 mRNAs are given in Table S14. These genes are involved in pathways involving the cell cycle (CCNO, SMC3), DNA repair (EPC2), RNA metabolism (RTCD1), neurological function (CHAT, FAM13C, EMX1, ATXN7, and HOXC4) and intracellular signaling (TMEM38B, SHOC2). However, the mechanisms by which they are connected to the designated lncRNAs and so to the prevalence of GSVs remain unknown. The lncRNAs and related gene pathways detected in our study suggest the complicated molecular mechanism of varicose veins.
The significantly aberrantly expressed mRNAs were clustered into 7 pathways that were involved in the pathophysiological properties of GSV varicosity (Figure 3). Of these, six are metabolic pathways. These pathways are involved in lysosomes, peroxisomes, glycolysis, diabetes, and tyrosine metabolism. Lysosomes transport the iron to the cytoplasmic labile iron pool [39] and glycolysis is a central molecule in glucose metabolism [40]. Peroxisomes participate in many crucial metabolic processes such as fatty acid oxidation, biosynthesis of ether lipids and free radical detoxification [41]. Maturity onset diabetes of the young (MODY) is an autosomal dominant form of diabetes involving defective b cells. Another three pathways, fatty acid metabolism, tyrosine metabolism and pyrimidine metabolism, are related to the amino acid metabolism. Taken together, these associations suggest that metabolic processes are involved in the pathogenesis of GSV varicosity.
The strength of this study should be addressed. Paired tissues taken from GSV varicosity patients were used to evaluate differences in expression between tissues from lesions and from adjacent normal segments of the saphenous vein. Unlike studies in which saphenous veins from healthy subjects were served as controls, this study was not subject to confounders associated with differences in genetic backgrounds, environmental exposure, or lifestyle. In addition, subjects with symptoms of valve insufficiency and reflux of GSVs were the focus of the present work. This diminished the heterogeneity of the studied phenotype, and so reduced the likelihood of false positives. The criteria used to select  lncRNA for Q-RT-PCR validation was very robust. The limitations of this study should also be addressed. First, due to the lack of experimental validation, the bio-computationally derived links between lncRNAs and either individual mRNA or mRNA pathways should be viewed as preliminary. Second, varicose veins typically contain inflamed tissue, and inflammation is an important mechanism for the development of varicose veins. In addition, changes in mRNA and lncRNA levels may reflect pathological changes of inflamed tissue [42,43]. However, since we did not strip the NVs of the additional lymphatic tissue and connective tissue that can build up in these samples, we could not distinguish whether or not the enriched lncRNAs and mRNAs are attributed to the additional cell types that are naturally more prevalent in inflamed tissue, which call for further study.
In summary, we systematically screened, validated, and functional predicted the aberrantly expressed lncRNAs primary varicose GSVs. These findings provide novel insights into physiology of lncRNAs and pathophysiological properties of GSV varicosities. These data may be useful in further investigation. The aberrantly expressed lncRNAs may also serve as new therapeutic targets for varicose veins.       Figure 3. Pathways of dysregulated mRNAs with the enrichment scores (2log10 (P-value)).2. The bar plot shows the enrichment scores (2log10 (P-value)) value of the significant enrichment pathways. The white bar shows the pathway in which the up-regulated mRNAs were found to be involved and the blue bars show the pathways in which the down-regulated mRNAs were found to be involved. Pathway analysis involves mapping genes to KEGG pathways. The P-value denotes the significance of the correlation between the pathway and the conditions. Most of the shown here are related to metabolism, which indicates that the varicose veins may be a metabolic disease. doi:10.1371/journal.pone.0086156.g003