Global Expression Profiling Identifies a Novel Hyaluronan Synthases 2 Gene in the Pathogenesis of Lower Extremity Varicose Veins

Lower extremities varicose veins (VV) are among the most easily recognized venous abnormalities. The genetic mechanism of VV is largely unknown. In this study, we sought to explore the global expressional change of VV and identify novel genes that might play a role in VV. We used next-generation ribonucleic acid (RNA) sequence (RNA seq) technology to study the global messenger RNA expressional change in the venous samples of five diseased and five control patients. We identified several differentially expressed genes, which were further confirmed by conventional reverse transcription polymerase chain reaction (RT-PCR). Using these significant genes we performed in silico pathway analyses and found distinct transcriptional networks, such as angiogenesis, cell adhesion, vascular injury, and carbohydrate metabolisms that might be involved in the mechanism of VV. Among these significant genes, we also found hyaluronan synthases 2 gene (HAS2) played a pivotal role and governed all these pathways. We further confirmed that HAS2 expression was decreased in the venous samples of patients with VV. Finally, we used a zebrafish model with fluorescence emitting vasculature and red blood cells to see the morphological changes of the venous system and blood flow. We found that HAS2 knockdown in zebrafish resulted in dilated venous structural with static venous flow. HAS2 may modulate the transcriptional networks of angiogenesis, cell adhesion, vascular injury, and carbohydrate metabolisms in venous tissues and downregulation of HAS2 may underlie the mechanism of VV.


Introduction
Varicose veins (VVs) of lower extremities are among the most easily recognized vascular abnormalities with superficial venous tortuosity and enlargement [1]. Clinical manifestations can vary from leg edema to chronic, disabling venous ulceration and it affects around 25% of the adult population [2]. Age, female sex, family history, pregnancy, obesity, and prolonged standing have been shown to be predisposing factors in several epidemiological studies [3][4][5][6][7][8]. The complex pathophysiology of greater saphenous vein (GSV) valvular reflux and altered extracellular matrix proteins with increased proteolytic enzyme have been addressed as potential mechanisms of VV [9][10][11]. Recently, it has been shown that matrix metalloproteinases may play an important role in the remodeling of venous extracellular matrix and mechanism of VV) [8,10,[12][13][14][15]. However, these processes account for only some of the tangled mechanisms, and many more remain to be identified.
Genetic study is a powerful tool to identify the novel mechanism of human diseases. Several candidate genes have been identified as the possible susceptibility genes for VV, such as COL1A2, FOXC2, NOTCH3, TGF-β, thrombomodulin, and VEGF1 [16][17][18][19][20][21][22]. However, whether these genes are true susceptibility genes for VV remains questioned. The study recently provides novel insight into the physiology of lncRNAs and the pathogenesis of varicose veins, but still need to do more for further investigation [2]. Genome-wide association studies (GWASs), a non-candidate-driven approach, investigate the entire genome and can identify unknown or novel disease-associated single nucleotide polymorphisms (SNPs) or genes. [23][24][25][26]. Recently, GWASs have identified several novel genes that may be involved in the mechanism of venous pathogenesis [13,27].
However, it is difficult to correlate the identified significant SNPs by GWAS to the mechanism of disease. Therefore, instead of traditional forward genetic approach, in the present study we used a reverse genetic approach. To this end, we first systematically investigated and compared the global expressional changes in the venous samples of diseased and control patients using the next-generation ribonucleic acid (RNA) sequence (RNA seq) technology. We identified several differentially expressed genes that might be involved in the molecular mechanism of VV. We then used a novel zebrafish model to see the phenotypic changes when these genes were perturbed.

Patient Selection
Five patients with VV were consecutively enrolled. All patients were diagnosed and assessed by a vascular surgeon and underwent duplex ultrasonographic scanning to assess suitability for treatment and entry into the study. Inclusion criteria were an age of 18 years or older, the presence of unilateral or bilateral primary symptomatic varicose veins at CEAP classification grade 3 or higher [28] and reflux of GSV of more than 0.5 s on duplex ultrasonography after calf compression [29]. Patients with current deep-vein thrombosis, acute superficial vein thrombosis, bleeding diathesis, pregnancy, and contraindications to the use of general or regional anesthesia were excluded. This study was approved by the institutional review board of the National Taiwan University Hospital and informed consent was obtained from participating subjects.

Surgical Removal of Venous Samples
GSV ablation using the ClosureFAST (VNUS Medical Technologies Inc, San Jose, CA, USA) with concomitant stab phlebectomy were done. All the procedures were performed in a two-day inpatient setting. The GSV was cannulated at the most distal point of venous reflux, and the catheter tip was positioned 2 cm distal to the saphenofemoral junction. The GSV was ablated with the concomitant phlebectomy [30,31]. The specimen collected from the avulsion phlebectomy were used for analysis.
Five control patients were consecutively enrolled. The controls patients had coronary artery disease and ischemic cardiomyopathy. Patients taking steroids or with conditions that could modify leukocyte activity (e.g., rheumatoid arthritis, vasculitis, or other collagen diseases) were excluded. Their GSV samples were harvested during elective coronary artery bypass to serve as venous grafts and the redundant part was used for study. The samples were put into and stored in liquid nitrogen immediately after they were harvested. GSVs competency and diameter were screened by the pre-operative duplex ultrasonography. The baseline characteristics of patients with VV and controls are shown in Table 1.

Sample Preparation and RNA Sequencing
Total RNA was extracted from the venous samples. The tissue was homogenized with POLYTRON ® PT 2100 (Dispergierund Mischtechnik, Littau, Switzerland), and Trizol solution (Gibco BRL, Grand Island, NY, USA) was added for RNA extraction. The extracted RNA was dissolved in diethyl pyrocarbonate-treated distilled water. Spectrophotometry at 260 and 280 nm was performed to measure the amount and quality of RNA. The total RNA samples were first treated with DNase I to degrade any possible deoxyribonucleic acid (DNA) contamination. Then, the messenger RNA (mRNA) was enriched by using the oligo (dT) magnetic beads (for eukaryotes). Mixed with the fragmentation buffer, the mRNA was fragmented into short fragments (about 200 base pair). Then the first strand of complimentary DNA (cDNA) was synthesized by using random hexamer-primer. Buffer, dNTPs, RNase H, and DNA polymerase I were added to synthesize the second strand. The double strand cDNA was purified with magnetic beads. End reparation and 3 -end single nucleotide A (adenine) addition were then performed. Finally, sequencing adaptors were ligated to the fragments. The fragments were enriched by polymerase chain reaction (PCR) amplification. During the quality control step, Agilent 2100 Bioanaylzer (Agilent, Santa Clara, CA, USA) and ABI StepOnePlus Real-Time PCR System (ABI, Foster City, CA 94404, USA) was used to qualify and quantify the sample library. The library products were ready for sequencing via Illumina HiSeqTM 2000 (Illumina, San Diego, CA, USA).

Bioinformatics Analysis of RNA Sequence Data
High quality reads with less than 10% unknown bases and median quality scores above 20 were extracted from the raw RNA-seq data. Cleaned reads were mapped to the reference genome hg19 using SOAPaligner/SOAP2 (Beijing Genomics Institute, Beijing, China) [33] with no more than two mismatches allowed in the alignment. The expression level for each gene was calculated using the RPKM method [34]. If there was more than one transcript for a gene, the longest transcript was used to calculate expression level and coverage. Function annotation could provide gene ontology (GO) annotation information.

Screening of Differentially Expressed Genes (DEGs)
This analysis included the screening of genes that were differentially expressed among samples, GO functional enrichment analysis and KEGG pathway enrichment analysis for these DEGs. We used false discovery rate (FDR) ≤0.001 and the absolute value of log2 ratio ≥1 as the threshold to judge the significance of gene expression difference.

Pathway and Network Analysis
While so many components are implicated in the mechanism of a complex trait disease, it is difficult to infer which components play a determining role. In this scenario, pathway analysis is a tool used to identify possible biological pathways involving related components. Ingenuity pathway analysis (IPA) is a powerful analysis and pathway-search tool that uncovers the significance of omics data and identifies new targets or candidate biomarkers within the context of biological systems. IPA has been specifically designed for the analysis, integration, and interpretation of data derived from omics experiments. Therefore, in the present study, to derive a possible pathway that may be implicated in the mechanism of AF related thromboembolic stroke, data were analyzed through the use of IPA (QIAGEN Inc., Germantown, MD, USA) [35]. The significance values for the canonical pathways in IPA was calculated using the right-tailed Fisher's exact test. The p-value for a given pathway annotation is calculated by considering (1) the number of DEGs that participate in that pathway and (2) the total number of genes that are known to be associated with that pathway in the selected reference set. Higher the number of DEGs involved, the more likely it is that the association is not due to random chance, and thus more significant the p-value is. Similarly, the larger the total number of genes known to be associated with the pathway, the greater the likelihood that an association is due to random chance, and the p-value accordingly becomes less significant. The p-value identifies statistically significant over-representation of DEGs in a given pathway. The null hypothesis (H 0 ) states.
H 0 = overlap of molecules for a particular biological category is due to chance. If p-value < 0.005, we reject the null hypothesis. This implies the association between a set of DEGs in our experiment and a given pathway is not due to random chance.

Zebrafish Model
Tg(fli1:egfp)xTg(gata1:dsRed) (egfp, enhanced green fluorescent protein), zebrafish with green fluorescence emitting vasculature and red fluorescence emitting red blood cells were obtained from the Zebrafish International Resource Center (Eugene, OR, USA) [36]. Both wild-type (WT; AB strain) and Tg(fli1:egfp)xTg(gata1:dsRed) zebrafish were kept under a 14-h light and 10-h dark photoperiod at 28.5 • C. The procedures for zebrafish culture, embryo collection, and fluorescence observation have been described previously [37]. We followed the standard criteria to design the zebrafish developmental stages [38]. The zebrafish experiments were approved by the Institutional Animal Care and Use Committee of the National Taiwan University College of Medicine. The anesthesia and euthanasia were performed with diluted tricane methanesulfonate (160 mg/L) and benzocaine hydrochloride (250 mg/L).

Statistical Methods
Continuous data were presented as means ± standard deviation (SD) and analyzed by using the Mann-Whitney U test for two-group comparisons. Categorical data were compared with the Fisher's exact test.

Global Expressional Profiling of VV
Global expressional profiling of venous samples from patients with VV and normal venous samples using the RNA seq technique was performed. The expressional profiling of the significantly up or downregulated genes with log2 ratio >7 were listed in Table 2. Among all these differentially expressed genes, many of them encode proteins functioning as extracellular matrix enzymes, such as Kallikrein-related peptidase 5 (KLK5), chitinase-3-like protein 1 (CHI3L1), hyaluronan synthase 2 (HAS2), and carbonic anhydrase 4 (CA4). These proteins may be involved in remodeling of extracellular matrix proteins and play important roles in vascular integrity and scaffolding. The others are genes encoding cytokines such as chemokine (C-X-C motif) ligand 1 (CXCL1), chemokine (C-X-C motif) ligand 8 (CXCL8), colony stimulating factor 3 (CSF3), and interleukin 6 (IL-6), implicating the importance of paracrine factors related to cell growth and inflammation in the mechanism of VV. Because these genes might be implicated in the mechanism of VV, we did a pathway analysis to elucidate the possible translational meaning and significance (Table 3; Figure 1). Interestingly, the results show that decreased angiogenesis, decreased cell-to-cell binding, decreased cell invasion or motility, and enhanced vascular lesion pathways are involved in the mechanism of VV. The results also show that increased carbohydrate accumulation may be involved in the mechanism of VV. Among all these pathways, HAS2 gene downregulation plays a pivotal role to regulate and govern all these pathways (Figure 1).   We then validated several differentially expressed genes identified in the RNA seq by conventional RT-PCR using VV and control venous samples. Among all the up or downregulated genes in VV samples, there were several genes encoding enzymes in the extracellular matrix such as KLK5, CHI3L1, HAS2, and CA4. We have validated the expressional changes by conventional RT-PCR for these genes (Figure 2). The result showed decreased mRNA levels of CHI3L1, HAS2, and CA4 genes and increased mRNA level of KLK5 gene in VV samples, which was similar to the result of RNA seq. We also found another gene Ficolin 1 (FCN1) encoding extracellular space protein the mRNA level of which was also decreased in the VV samples.

Downregulation of Hyaluronan Synthases 2 in Venous Tissues from Patients with VV
Among all these extracellular matrix enzyme genes, we were particularly interested in the HAS2 gene because it played a pivotal role in the pathway analysis, which modulated all the VV-related

Downregulation of Hyaluronan Synthases 2 in Venous Tissues from Patients with VV
Among all these extracellular matrix enzyme genes, we were particularly interested in the HAS2 gene because it played a pivotal role in the pathway analysis, which modulated all the VV-related pathways (Figure 1). The role of HAS in vascular biology has never been addressed before. Because of its possible important role in the mechanism of VV, we repeated the experiment again and confirmed that the mRNA level was significantly lowered in the venous samples from patients with VV (p < 0.0001) (Figure 3). pathways (Figure 1). The role of HAS in vascular biology has never been addressed before. Because of its possible important role in the mechanism of VV, we repeated the experiment again and confirmed that the mRNA level was significantly lowered in the venous samples from patients with VV (p < 0.0001) (Figure 3).

Knockdown of HAS2 Results in Venous Dilation and Blood Flow Stasis in Zebrafish
Tg(fli1:egfp)xTg(gata1:dsRed) zebrafish with green fluorescence emitting vasculature and red fluorescence emitting red blood cells were used to observe the anatomic change of vasculature (green fluorescence) and blood flow (red fluorescence) and to explore the functional role of the identified gene. We first compared the homology of amino acid sequence of HAS2 among species using Clustal W [39]. HAS2 is evolutionally conservative and expressed in both low and high vertebrates, indicating its important biological function ( Figure S1). The similarity of amino acid sequence of

Knockdown of HAS2 Results in Venous Dilation and Blood Flow Stasis in Zebrafish
Tg(fli1:egfp)xTg(gata1:dsRed) zebrafish with green fluorescence emitting vasculature and red fluorescence emitting red blood cells were used to observe the anatomic change of vasculature (green fluorescence) and blood flow (red fluorescence) and to explore the functional role of the identified gene. We first compared the homology of amino acid sequence of HAS2 among species using Clustal W [39]. HAS2 is evolutionally conservative and expressed in both low and high vertebrates, indicating its important biological function ( Figure S1). The similarity of amino acid sequence of HAS2 between zebrafish and humans was up to 81.3% ( Figure S2).   Morpholino experiment was quantitatively performed to further define the role of HAS2 in the mechanism of VV. In fish with mild phenotypes, a local dilation of venous structure in the caudal vein plexus (CVP) was found, which resulted in a protruding bulb in the tail, mimicking human VV phenotype of protruding venous dilatation ( Figure 4A-D). In fish with moderate phenotypes, dilation of veins, such as intersegmental vein (ISV) or sub-intestinal vein (SIV), was found ( Figure 5). Sluggish blood flow in CVP was also noted (Supplementary Movie 1 and Movie 2), mimicking the human VV phenotype of venous blood flow stasis and leg edema. Furthermore, venous branches were also noted in the SIV of the MO-injected fish, which was a consequence of either a dysregulated angiogenetic process or a delay of retraction of leading venous buds during development [40]. It implies that this is a dysregulated angiogenetic process since HAS2 plays a pivotal role in the angiogenesis network, which is one of the key altered pathways in our transcriptome analysis (Figure 1). In those with severe phenotypes, tangled veins were found in the distal tail which resulted in malformation of tail ( Figure 4E Quantitative data after the injection of HAS2 MO in more than 1000 embryos are shown in Table  4. There was a dose-response relationship. Higher dose of MO (375 uM) was associated with a higher percentage of fish with defective vascular phenotypes than lower dose MO (250 uM). In the highest dose group (500 uM), because of the lethality of high-dose MO, most of the embryos were dead (91%). The survived embryos were probably those without or with mild MO effect and thus there was a relatively lower percentage of survived fish with defect vasculature (43%), compared to those with a medium dose of MO injection (72%) (375 uM) (Table 4). Therefore, we confirm that disruption of HAS2 function directly leads to vascular malformation in vertebrates. Quantitative data after the injection of HAS2 MO in more than 1000 embryos are shown in Table 4. There was a dose-response relationship. Higher dose of MO (375 uM) was associated with a higher percentage of fish with defective vascular phenotypes than lower dose MO (250 uM). In the highest dose group (500 uM), because of the lethality of high-dose MO, most of the embryos were dead (91%).
The survived embryos were probably those without or with mild MO effect and thus there was a relatively lower percentage of survived fish with defect vasculature (43%), compared to those with a medium dose of MO injection (72%) (375 uM) (Table 4). Therefore, we confirm that disruption of HAS2 function directly leads to vascular malformation in vertebrates.

Discussion
In the present study, we first used a global expressional profiling to address the molecular mechanism of VV. We further found several important pathways, such as angiogenesis, cell binding, vascular lesion and carbohydrate metabolism possibly involved in the mechanism of VV. We further used a zebrafish model to address the functional role of an identified gene HAS2 in the mechanism of VV.
In the last decade, the pathophysiological theories of VV shifted from purely mechanical factors to hypotheses focusing on complex molecular and histopathologic alterations in the vessel wall and the extracellular matrix, which may be caused by the genetic inheritance [1,18,19]. In our study, we identified the global expressional change of downregulations of CHI3L1, HAS2, and CA4, and upregulation of KLK5 in the varicose venous samples with both RNA seq technique and conventional RT-PCR.
All these genes encode extracellular matrix enzymes, compatible with the concept of histopathologic alterations in the extracellular matrix as a potential mechanism of VV. Extracellular matrix may provide a scaffold for the vascular structure and integrity. When the components of extracellular matrix have been altered (e.g., decomposition of extracellular matrix proteins by excessive enzyme digestion or accumulation of extracellular matrix proteins by decreased enzyme digestion), the scaffold function may deteriorate, resulting in a weakening of vascular wall support and dilatation or protrusion of the vasculature, a pathognomonic picture of VV. Furthermore, after pathway analysis using all the differentially expressed genes, we found decreased angiogenesis, decreased cellular binding and invasion, increased carbohydrate accumulation, and vascular lesion were all involved in the mechanism of VV. All these pathways are novel and have been rarely addressed as the potential mechanism of VV.
Decreased angiogenesis may implicate impaired repair or proliferative activities of vessels and an injury model of VV. Furthermore, an involvement of vascular lesion network also supports the hypothesis of this injury model of VV. Decreased cell-cell binding and invasion may implicate decreased structural integrity of vessels, indicating weakening of vessel wall and possible protrusion of vessel wall, which is the prototype phenotype of VV. Our study also shows that increased carbohydrate accumulation may be involved in the pathogenesis of VV; however, further studies are required to elucidate the detailed mechanism.
As mentioned before, recent studies have also highlighted the pathological changes in the extracellular matrix as the potential mechanism of VV [1,8,10,12,13,15]. Among all these extracellular matrix enzyme genes, HAS2 expression was decreased in VV samples. Hyaluronan or hyaluronic acid (HA) is a high molecular weight unbranched polysaccharide, and is a major constituent of the extracellular matrix and synthesized by hyaluronan synthases (HAS), which are membrane-bound enzymes and use UDP-α-N-acetyl-D-glucosamine and UDP-α-D-glucuronate as substrates to produce the glycosaminoglycan hyaluronan at the cell surface and extrude it through the membrane into the extracellular space. In our zebrafish model, we found disruption of HAS2 function directly led to venous malformation in vertebrates. Although HAS2 has been shown to play a role in vascular morphogenesis [41], its specific role in venous tissues and the pathophysiology of the most important human venous diseases-that is, VV-has never been addressed before. The mechanism by which HAS2 cause VV may be through its pivotal role to modulate and govern multiple pathways that cause VV, such as angiogenesis, cell to cell adhesion, vascular lesion, and carbohydrate metabolism.
There are limitations in the present study. First, there were only five patients in each group and the differences of genetic results we found between VV patients and non-VV patients might be due to something else not related to VV. Second, the genetic changes in the VV samples might be post-transcriptional or happened due to the local environment in the legs, and the zebrafish model only showed a systemic effect. Therefore, HAS2 might not be the gene responsible for VV. Third, we did not have negative scrambled MO control group, positive control group with other MOs targeting to HAS2 and rescue studies. Therefore, we could not rule out a non-specific effect of MO. Finally, some of our control patients had diabetes that could modify leukocyte activity and might affect the genetic results.
However, there is little doubt that the development and progression of VV and venous ulcer disease is multifactorial and includes a moderate-to-strong genetic component [8]. The genetic studies of VV conducted to date are few and have few functional implications. This is the first study to investigate the global gene expressional profiling in determining the candidate gene of VV. In this study, we show the possible role of HAS2 gene in the pathogenesis of VV in both global expression and functional study.