Comprehensive Bioinformation Analysis of the MRNA Profile of Fascin Knockdown in Esophageal Squamous Cell Carcinoma

Fascin is an actin-bundling protein that induces cell membrane protrusions and increases motility of normal and transformed epithelial cells (Otto, 1994; Edwards et al., 1995). Fascin protein distributes in structures containing actin bundles including filopodia and stress fibers, causing F-actin to aggregate side-byside into bundles (Yamashiro-Matsumura et al., 1986). Accumulated evidences showed overexpression of fascin was found in several human epithelial cancers, such as colonic (Jawhari et al., 2003), gastric (Hashimoto et al., 2004), breast (Grothey et al., 2000; Yoder et al., 2005), skin (Goncharuk et al., 2002), urothelial (Tong et al., 2005). Fascin immunoreactivity in these tumors usually correlates with high-grade, extensive invasion, distant metastasis, or poor prognosis. Our previous studies also illustrated that the overexpression of fascin might be correlated


Introduction
Fascin is an actin-bundling protein that induces cell membrane protrusions and increases motility of normal and transformed epithelial cells (Otto, 1994;Edwards et al., 1995).Fascin protein distributes in structures containing actin bundles including filopodia and stress fibers, causing F-actin to aggregate side-byside into bundles (Yamashiro-Matsumura et al., 1986).Accumulated evidences showed overexpression of fascin was found in several human epithelial cancers, such as colonic (Jawhari et al., 2003), gastric (Hashimoto et al., 2004), breast (Grothey et al., 2000;Yoder et al., 2005), skin (Goncharuk et al., 2002), urothelial (Tong et al., 2005).Fascin immunoreactivity in these tumors usually correlates with high-grade, extensive invasion, distant metastasis, or poor prognosis.Our previous studies also illustrated that the overexpression of fascin might be correlated

Comprehensive Bioinformation Analysis of the MRNA Profile of Fascin Knockdown in Esophageal Squamous Cell Carcinoma
Bing-Li Wu 1& , Lie-Wei Luo 2& , Chun-Quan Li 3 , Jian-Jun Xie1 , Ze-Peng Du 4 , Jian-Yi Wu 1 , Pi-Xian Zhang 1 , Li-Yan Xu 5 *, En-Min Li 1 * with the malignant transformation of normal esophageal squamous cells to cancer cells (Rong et al., 2004).Moreover, upregulation of fascin was markedly correlated with cell proliferation and lymph node metastasis, which implicated the potential of fascin as a novel biomarker for esophageal squamous cell carcinoma (ESCC) (Zhang et al., 2006).Later, we have stably knockdowned fascin by RNAi in ESCC cell line EC109.Downregulation of fascin resulted in a suppression of cell proliferation, a decrease in cell invasiveness as well as tumor growth in vivo.The effect of fascin on cell invasiveness correlated with the activation of matrix metalloproteases such as MMP-2 and MMP-9.These results suggested that fascin might play crucial roles in regulating neoplasm progression of ESCC.To better understand the molecular mechanism of fascin, the Affymetrix GeneChip Human genome U133 plus 2.0 arrays has been applied to analyze the mRNA expression profile in fascin knockdown cells.However, the biological information of this mRNA expression profile has not been fully mining and described in our previous reports (Xie et al., 2010).
To comprehensive understand the mRNA profile and illustrate the functional roles of fascin to a greater extent, we analyzed the mRNA profile with multiple bioinformatic methods in this study.

The differentially expressed genes
The mRNA expression profile of fascin knockdown by RNAi in EC109 has been uploaded to GEO database (GSE11373) (http://www.ncbi.nlm.nih.gov/geo/), and then subjected to normalization and log transform treatment.The differentially expressed genes (DEGs) were obtained by folds change analysis.

Gene Ontology (GO) enrichment and functional annotation
Analysis of GO annotations and enrichments was performed by applying the Onto-Express (http://vortex.cs.wayne.edu/Projects.html)(Draghici et al., 2003).Onto-Express calculates the expected number of occurrences of each functional category in each cluster.The probability that each functional category was over-represented in a cluster was derived using a binomial model and the P-values were corrected for multiplicity using the Bonferroni method.
Functional Annotation Chart in DAVID (http://david.abcc.ncifcrf.gov/) is able to identify the most relevant (over-represented) biological terms associated with a given gene list (Huang et al., 2009).Currently, it maintains over 40 annotation categories, including GO terms, proteinprotein interactions, protein functional domains, disease associations, pathways, sequence features, homology, gene functional summaries, gene tissue expression, and literature.The terms from Functional Annotation Chart with P < 0.05 were visualized by Enrichment Map plugin of Cytoscape (Merico et al., 2010).

Analysis of KEGG pathway and subpathway
We performed KEGG pathway enrichment analysis by applying the bioconductor SubpathwayMiner package.Except traditional entire KEGG pathway, SubpathwayMiner detects the local structure of an entire pathway (subpathway), which limited to local area of pathway and help to understand the detail information of interested genes (Li et al., 2009).For a given KEGG pathway, the subpathways were obtained by searching all possible paths between start-points and end-points in the adjacency matrix generated by node relationships.SubpathwayMiner adopts the k-clique concept in social network analysis to define subpathway.We set k = 4 in this study, which means that the distance among all genes within the subpathways is no greater than 4.

The promoter sequence patterns of differentially expressed gene
The 2000bp promoter sequences of top 20 most changed downregulated and upregulated DEGs were retrieved from UCSC database (http://genome.ucsc.edu/), respectively.The promoter sequence patterns over-represented or down-represented in these two gene set were discovered by POCO program (http://ekhidna.biocenter.helsinki.fi/poxo/poco/poco).POCO is tool to find over-represented, under-represented and distinctly represented regulatory patterns from either one or two sequence sets (Kankainen et al., 2005).The background organism was set as homo_sapiens_clean and longest pattern length was set as 8 in present study.Subsequently, patterns with significant enrichment were then screened with the JASPAR database to identify recognized transcription factor binding sites (similarity index > 0.70).

GO enrichment and functional annotation analysis
Using 2-fold change as threshold, we totally obtained 276 differentially expressed genes (DEGs) from GSE11373, including 193 upregulated genes and 83 downregulated genes.To understand the functional categories of these DEGs, they firstly were analyzed by DAVID GO enrichment analysis.There were15, 8 and 13 significant GO terms for these DEGs in Biological Process, Cellular Component and Molecular Function, respectively (Figure 1).Not surprisingly, several GO terms associated with cytoskeleton organization were found, including cell adhesion (P=0.015)from Biological Process, actin filament binding (P=0.018) from Molecular    The DEGs were also clustered through Functional Annotation Chart and visualized by Enrichment Map in Cytoscape.As showed in Figure 2, a node represented one functional annotation category.The more significant of the category, the node was bigger.Nodes from the same kind of functional category were showed in the same shape and color background.While the edge width was defined by the overlap coefficient between these categories (overlap coefficient cut-off 0.6).The more same genes in two nodes, the edge width is wider.Except 46 terms from GO categories, the Functional annotation chart results also included other 45 terms from the following annotation categories, 11 INTERPRO, 5 SMART, 14 SP_PIR_KEYWORDS, 14 UP_SEQ_FEATURE and 1 KEGG_PATHWAY.These results would provide more information than the merely GO enrichment.The top two enrichments of SP_PIR_KEYWORDS were cell adhesion and extracellular matrix, which contained 12 DEGs (COL4A3, CNTN5, COL7A1, CTGF, ITGB8, PCDH11X, PVRL2, FBLN7, DSC2, PCDH7, THBS1, CYR61) and 8 DEGs (COL4A3, MAMDC2, COL7A1, CTGF, ITGB8, FBLN7, DCN, FBN2), respectively.These results might explain fascin as an acting-binding protein, its overexpression in tumor cells affects cell-cell adhesion and cell mobility.InterPro is a database that provides functional analysis of protein sequences by classifying them into families and predicting the presence of domains and important sites (Hunter et al., 2012).In the category of INTERPRO, the term IPR006210:EGF-like ranked the first which containing 8 genes (ITGB8, RSPO3, PTGS1, FBLN7, FBN2, NRG1, THBS1, PROS1).It is presumed EGF-like domain containing proteins binds to the cell surface receptors to drive intracellular signal transduction cascades, which eventuate in diverse cell fates including proliferation, differentiation, migration and inhibition of apoptosis (Sanderson et al., 2006).This indicated these DEGs might be involved in the cell proliferation and migration mediated by fascin.Interestingly, hsa04350:TGF-beta signaling pathway, the only term found in KEGG_PATHWAY category, involved 5 genes (NOG, FST, DCN, ID3 and THBS1).SMART (Simple Modular Architecture Research Tool) allows the identification and annotation of genetically mobile domains and the analysis of domain architectures (Letunic et al., 2012).The top overrepresented functional class of SMART was SM00181:EGF, which was similar to the INTERPRO category.

KEGG pathways and subpathways enrichments
To find out how the DEGs from fascin knockdown influence the cell signal transduction, the DEGs were mapped to KEGG pathways by SubpathwayMiner R program package.The DEGs were only enriched in two entire pathways, of which TGF-beta signaling pathway which enriched 5 DEGs and Fructose and mannose metabolism pathway which enriched 3 genes (Table 1).This result was consistent with the KEGG_PATHWAY result from Functional Annotation Chart.By applying SubpathwayMiner package, the DEGs were found significantly enriched in 39 subpathways corresponding to 18 entire KEGG pathwahys (Table 2).Except TGF-beta signaling pathway and Fructose and mannose metabolism pathway, 39 subpathways derived from many other entire pathways were found, which could not detected in traditional KEGG pathway enrichment because of the statistical threshold.As an F-actin biding protein, it is important to reveal how fascin regulate the cytoskeletal structure.To our interest, two subpathways derived from the pathway of Regulation of actin cytoskeleton (hsa04810) have been found.One subpathway path:04810_5 contained three DEGs, the upregulated gene PIP5K1B and two downregulated genes MYL9 and MLCK, while the other subpathway path:04810_8 contained two upregulated genes (PIP5K1B, SCIN) (Figure 3).These results indicated the dyregulation of fascin involved cytoskeletal re-construction.These results indicated the DEGs might involve other pathways, not at the entire level, but also play critical roles in local part of total pathway (subpathway).These results might explain the importance of fascin in the initiation and development of tumor.

Sequence patterns in the promoters of upregulated and downregulated genes
We presumed there are certain sequence patterns in the promoters of DEGs responding to the transcriptional regulation for the co-expression of upregulation and downregulation genes sets after fascin was knockdowned, respectively.POCO is applied to find over-represented and under-represented regulatory patterns between the upregulation and downregulation genes promoter sequences (Kankainen et al., 2005).There were 70 significant sequence patterns that over-represented in the downregulated genes but down-represented in the upregulated genes, of which the top 10 patterns were showed in the Table 3.On the other hand, 67 patterns were found in upregulated genes and simultaneously deficient in downregulated genes, and the top 10 patterns were also presented in the Table 3.All these patterns were presented in all 20 most downregulated or upregulated genes, respectively.Subsequently, the total 137 significant patterns with significant enrichment were then screened  the JASPAR database to identify potential transcription factors.Then four patterns corresponding to four unique transcription factors were found in the downregulated genes set, there were MNB1A, c-ETS, GATA2 and Prrx2 (Figure 4).On the other hand, four unique transcription factors (Arnt-Ahr, ZNF42, Ubx and TCF11-MafG) were found in the upregulated genes set (Figure 4).

Discussion
ESCC is one of the most deadly malignant tumors in the world, with an overall 5-year survival rate less than 20% (He et al., 2005).Fascin is an important mediator of carcinoma invasion and metastasis including ESCC (Rong et al., 2004;Zhang et al., 2006;Xie et al., 2010).
Here we presented a comprehensive analysis of mRNA profile of fascin knockdown in ESCC cell line EC109 by multiple bioinformatic analyses.We obtained 276 DEGs using 2-fold as the threshold.
To better understand the functions of these DEGs, the DEGs were analyzed by Onto-Express GO enrichment tools.The results showed the GO term of cell adhesion was significant, in which the genes including CYR61 and THBS1.This was in agreement with previous results.CYR61 has been reported to participate in the development and progression of various cancers.Overexpression of CYR61 could promote growth, migration, and metastasis of various cancer cells, such as prostate cancer and osteosarcoma (Lin et al., 2005;Sun et al., 2008;Fromigue et al., 2011).CYR61 was downregulated 2.2 folds in the mRNA profiles of fascin RNAi.These results indicated that the reduction of cell proliferation as well as cell invasiveness in fascin knockdown cells might through key molecules such as CYR61.THBS1 was decreased 2 folds after fascin was knocked down in the data that applied in this study.Previous results show that THBS1 is a potent stimulator of prostate tumor cell migration (Firlej et al., 2011).It was reported that loss of THBS1 expression was also significantly associated with distal location, vascular invasion and distant metastasis in gastric cancer and with unfavourable histological grade and shorter overall survival in penile squamous cell carcinoma (Miyamoto et al., 2007;Guerrero et al., 2008).It was evidence that DEGs derived from fascin knockdown were also involved migration and invasionness of tumor cells directly or indirectly.SCIN, NEXN and EPLIN were enriched in the term of actin filament binding.Wang et al. have identified NELIN product was associated with F-actin by immunofluorescence and immunoprecipitation.Stable transfection of NEXN into HeLa cells increased the cell migration by 2.17-fold and the adhesion by 1.67fold, respectively (Wang et al., 2005).EPLIN has at least two actin binding sites.Purified recombinant EPLIN inhibits actin filament depolymerization and cross-links filaments in bundles.EPLIN does not affect the kinetics of spontaneous actin polymerization or elongation at the barbed end, but inhibits branching nucleation of actin filaments by Arp2/3 complex (Maul et al., 2003).
Functional Annotation Chart was applied to class the DEGs based on related gene functions.Annotation categories in Functional Annotation Chart included more than 40 functional categories.The annotation coverage increases the analytic power by allowing investigators to analyze their genes from many different biological aspects in a single space.The results showed the enrichment of genes involved in functions related to cell adhesion, extracellular matrix, These functional terms would be considered as significantly modulated by fascin, which also could be applied to explain the molecular mechanisms of fascin in initiation and development of tumor, especially in the ESCC.
It has proved that subpathway-based approach was more precise and flexible in annotation and identification of disease-risked pathways than entire KEGG pathway enrichment (Li et al., 2012).One of the advantages of subpathway method is it could identify a smaller disease related area in the pathway, but also providing the information of subpathway structure.In this study, our DEGs were found significantly enriched in 39 subpathways corresponding to 18 entire pathwahys.TGFbeta pathway has been confirmed as one critical pathway in the ESCC (Edmiston et al., 2005).The DEGs from fascin RNAi in ESCC also involved this TGF-beta pathway.Except TGF-beta pathway, more pathways related to cell adhesion and actin cytoskeleton were found by the application of subpathway analysis.These subpathway analysis results provided more detail information about the pathways that DEG involved.The DEGs enriched in the subpathways of Regulation of actin cytoskeleton have been reported their importance in the migration or invasionness of tumor cells.PIPKIB, an enzyme that generates phosphatidylinositol-4,5-bisphosphate, spatially regulates integrin endocytosis at adhesion sites to control cell migration fibrosarcoma HT1080 cell line (Edmiston et al., 2005).MYL9 encodes a myosin light chain that may regulate muscle contraction by modulating the ATPase activity of myosin heads (Chao et al., 2010).MLCK, a kinase could phosphorylate the myosin heavy and light chains.Expression levels of MLCK were well correlated with the pulmonary arterial pressure of patients with heart failure.In cultured cardiomyocytes, knockdown of MLCK by specific siRNAs decreased MLC2v phosphorylation and impaired epinephrine-induced activation of sarcomere reassembly (Medjkane et al., 2009).SCIN (scinderin) is a Ca(2+)-dependent actin-severing and -capping protein, might act as a molecular switch in the control of cortical F-actin dynamics during secretion (Seguchi et al., 2007).
To better understand how the significant changed DEGs was transcriptional co-regulated after fascin was knockdown, the promoter sequences of top 20 most changed downregulation and upregulation DEGs were analyzed by POCO program.Dozens of significant sequence patterns were found over-represented and under-represented in the downregulation gene set or upregulation gene set, respectively.These suggested that after fascin was knocked down, the change of signal transduction resulted in distinguishing transcription factors binding to the promoters of the downregulated or upregulated genes, or certain sequence patterns played critical regulation roles, to co-regulate of the significant changed downregulated gene set or upregulated gene set at transcriptional level.Researches have proved the roles of these potential transcription factors in the migration and invasionness of tumor cells.GATA2 is a member of the GATA family of zinc-finger transcription factors.It was reported GATA2 negatively regulated PTEN by preventing nuclear translocation of androgen receptor and by androgen-independent suppression of PTEN transcription in breast cancer to promote tumor invasion and metastases but not initiation (Lejen et al., 2002).c-ETS, also named ETS1, which belongs to the ETS transcription factor family, is a critical regulator in tumor development and metastasis (Wang et al., 2012).Overexpression of ETS1 has been found in ESCC, which was correlated with tumor angiogenesis, lymph node metastasis and patient survival (Kato et al., 2012).The aryl hydrocarbon receptor nuclear translocator (ARNT) belongs to the bHLH-PAS transcription factor family, ARNT binds to DNA in the form of a heterodimer with the aryl hydrocarbon receptor (AhR).The deficiency of ARNT was concurrent with a decline of EGFR and ERK1/2 phosphorylation in human epidermal keratinocyte (Mukherjee et al., 2003).
Consistent with our previous results, fascin transcription may be subjected to the regulation of the EGF/EGFR signaling pathway (Robertson et al., 2012).ZNF42 (also named Myeloid zinc finger 1, MZF1) belongs to the Kruppel family of zinc finger transcription factors, has been found bided to the Axl promoter, transactivated promoter activity, and enhanced Axl-mRNA and protein expression in a dose-dependent manner.ZNF42 induces invasion and in vivo metastasis in colorectal and cervical cancer, at least in part by regulating Axl gene expression (Lu et al., 2010).Transcription factors of TCF11 and MafG form nuclear heterodimer in the regulation of gene transcription mediated by TGF-beta signaling pathway, which is also the pathway that fascin involved (Berg et al., 2007;Mudduluru et al., 2010;Kannan et al., 2012).These results indicated that some migration and metastasis related transcription factors might regulate the DEGs after fascin was knockdown.These transcription factors might be part of fascin-donimated cellular network and play critical roles in the procession or signal transduction mediated by fascin.Though many sequence patterns did not corresponding to known transcription factors, it could not rule out their possibilities and importance in the regulation of DEGs.
In summery, a comprehensive understanding the role of fascin in ESCC was obtained by bioinformatics workflow of mRNA profile after fascin was knockdowned, especially by the means of Functional Annotation Chart and subpathway, which provide more information than traditional enrichment analysis.These analysis methods are well suited for high-throughput data and therefore are proposed as powerful tools for the search novel functional genes and pathways related to interested gene(s).

Figure 1 .
Figure 1.Pie Chart Presentation of GO Enrichment Carried out by Onto-express.Pie slices are proportional to the number of genes, indicated by corresponding GO term and its p-value.Categories shown are significantly overrepresented at p < 0.05

Figure 2 .
Figure 2. Enrichment Map for the Differentially Expressed Genes to Identify Significant Biological Functions (P<0.05).One node represented a significant functional terms from Functional annotation chart.Node size represents enrichment significance.Edges indicate overlap between gene sets, where the thickness indicates the size of the overlap Function, actin cytoskeleton (P=0.001) from Cellular Component.The term of cell adhesion contains 13 DEGs (COL4A3, THBS1, CYR61, CTGF, PCDH7, ITGB8, FBLN7, PPFIBP1, CNTN5, DSC2, COL7A1, PCDH11Xand PVRL2).SCIN, NEXN and EPLIN were enriched in the term of actin filament binding.This indicated fascin might impact the metastasis and invasion of ESCC through the functions of these genes.To our interest, several Transport related terms were discovered, such as organic acid transport, carboxylic acid transport, amine transport.These suggested that cell membrane structure might alter due to the cytoskeletal rearrangement in fascin knockdown.The DEGs were also clustered through Functional Annotation Chart and visualized by Enrichment Map in Cytoscape.As showed in Figure2, a node represented one functional annotation category.The more significant of the category, the node was bigger.Nodes from the same kind of functional category were showed in the same shape and color background.While the edge width was defined by the overlap coefficient between these categories (overlap coefficient cut-off 0.6).The more same genes in two nodes, the edge width is wider.Except 46 terms from GO categories, the Functional annotation chart results also included other 45 terms from the following annotation categories, 11 INTERPRO, 5 SMART, 14 SP_PIR_KEYWORDS, 14 UP_SEQ_FEATURE and 1 KEGG_PATHWAY.These results would provide more information than the merely GO enrichment.The top two enrichments of SP_PIR_KEYWORDS were cell adhesion and extracellular matrix, which contained 12 DEGs (COL4A3, CNTN5, COL7A1, CTGF, ITGB8, PCDH11X, PVRL2, FBLN7, DSC2, PCDH7, THBS1, CYR61) and 8 DEGs (COL4A3, MAMDC2, COL7A1, CTGF, ITGB8, FBLN7, DCN, FBN2), respectively.These results might explain fascin as an acting-binding protein, its overexpression in tumor cells affects cell-cell adhesion Figure 3.The Differentially Expressed Genes Enriched Subpathway from Regulation of Actin Cytoskeleton (path:04810).The two significant subpathways path: 04810_5 and path;04810_8 were indicated in origin entire pathway by red and blue frame, respectively.The structures of path: 04810_5 and path;04810_8 generated by Subpathway package were also showed

Figure 4 .
Figure 4. Motif Logos Were Showed for the Predicted Transcription Factors.The significant sequence patterns were screened with the JASPAR database to identify recognized transcription factors.The potential transcription factors regulating downregulated genes were labeled in blue, while transcription factors regulating downregulated genes were labeled in red

Table 1 . The enriched KEGG pathway of DEGs
# the first number is annotated DEGs in the pathway; the second number is total number of the pathway