Genes involved in TGFβ1-driven epithelial-mesenchymal transition of renal epithelial cells are topologically related in the human interactome map

Abstract

Background

Understanding how mesenchymal cells arise from epithelial cells could have a strong impact in unveiling mechanisms of epithelial cell plasticity underlying kidney regeneration and repair.

In primary human tubular epithelial cells (HUTEC) under different TGFβ1 concentrations we had observed epithelial-to-mesenchymal transition (EMT) but not epithelial-myofibroblast transdifferentiation. We hypothesized that the process triggered by TGFβ1 could be a dedifferentiation event. The purpose of this study is to comprehensively delineate genetic programs associated with TGFβ1-driven EMT in our in vitro model using gene expression profile on large-scale oligonucleotide microarrays.

Results

In HUTEC under TGFβ1 stimulus, 977 genes were found differentially expressed. Thirty genes were identified whose expression depended directly on TGFβ1 concentration. By mapping the differentially expressed genes in the Human Interactome Map using Cytoscape software, we identified a single scale-free network consisting of 2630 interacting proteins and containing 449 differentially expressed proteins. We identified 27 hub proteins in the interactome with more than 29 edges incident on them and encoded by differentially expressed genes. The Gene Ontology analysis showed an excess of up-regulated proteins involved in biological processes, such as "morphogenesis", "cell fate determination" and "regulation of development", and the most up-regulated genes belonged to these categories. In addition, 267 genes were mapped to the KEGG pathways and 14 pathways with more than nine differentially expressed genes were identified. In our model, Smad signaling was not the TGFβ1 action effector; instead, the engagement of RAS/MAPK signaling pathway seems mainly to regulate genes involved in the cell cycle and proliferation/apoptosis.

Conclusion

Our present findings support the hypothesis that context-dependent EMT generated in our model by TGFβ1 might be the outcome of a dedifferentiation. In fact: 1) the principal biological categories involved in the process concern morphogenesis and development; 2) the most up-regulated genes belong to these categories; and, finally, 3) some intracellular pathways are involved, whose engagement during kidney development and nephrogenesis is well known. These long-term effects of TGFβ1 in HUTEC involve genes that are highly interconnected, thereby generating a scale-free network that we named the "TGFβ1 interactome", whose hubs represent proteins that may have a crucial role for HUTEC in response to TGFβ1.

Background

Epithelial-to-mesenchymal transition (EMT) of renal tubular cells is a fundamental sign of epithelial cell plasticity in physiological processes such as regeneration and wound healing, but it also characterizes pathological conditions such as fibrosis and carcinogenesis.

The adult mammalian renal tubular epithelium exists in a relatively quiescent to slowly replicating state, but has great potential for regenerative morphogenesis following severe ischemic or toxic injury [1]. Dedifferentiation, i.e. the acquisition of mesenchymal markers such as vimentin and N-cadherin, seems to represent a crucial step in the recovery of tubular integrity and precedes the reconstitution of a well-differentiated morphology. In the adult kidney, however, the tubular cells' acquisition of a mesenchymal phenotype represents one of the crucial steps towards transdifferentiation into myofibroblasts, the effector cells of tubulo-interstitial fibrosis [2].

Transforming growth factor β1 (TGFβ1) is a key modulator of EMT in a variety of epithelial cells, but is also capable of inducing the myofibroblast phenotype, i.e. the acquisition of alpha smooth muscle actin (αSMA) microfilaments in fibroblasts during wound healing, in mesangial cells in culture and in renal tubular cells [3].

TGFβ1-induced EMT appears to depend primarily on intact Smad signaling. To date, Smad proteins are the only TGFβ1 receptor substrates with a demonstrated ability to propagate signals [4]. It is now becoming evident, however, that EMT is not a uniform process. Its role and features clearly differ, depending on the physiological context and type of epithelia (developmental EMT, oncogenic EMT, non-oncogenic EMT) [5].

Using primary human tubular epithelial cells (HUTEC), we demonstrated that chronic exposure to TGFβ1 prompted morphological, molecular and biochemical changes towards a mesenchymal phenotype, but this gave rise to no de novo expression of αSMA gene or myofibroblast phenotype [6]. We hypothesized that the process triggered by TGFβ1 in our model is a dedifferentiation event that may be part of the vital plasticity of renal tubular cells.

Our results prompted us to further characterize this EMT process. Since microarray technology powerfully monitors gene expression and has led to the discovery of pathways regulating complex biological processes, we explored the molecular mechanisms underlying this transition using this approach. A global view of the EMT process was obtained identifying the Gene Ontology (GO) classes enriched by differentially expressed genes and analyzing KEGG pathways involved in signal transduction. To obtain an overview of their topological properties, we also mapped differentially expressed proteins in the human interactome map using Cytoscape software. This analysis enabled us to establish that about 50% of the genes up- and down-regulated by TGFβ1 were strongly interconnected and formed a large network that we named the "TGFβ1 interactome".

Results

At genome-wide level, we investigated the expression profile changes occurring in the EMT of primary HUTEC under chronic TGFβ1 treatment. Our in vitro model of human renal EMT has been described in detail elsewhere [6]. By immunocytochemistry, we demonstrated that, in addition to the front-end to back-end cell morphology, TGFβ1 can induce a markedly dose-dependent up-regulated expression of mesenchymal markers, including collagen III, with a parallel drop in the expression of epithelial markers such as E-cadherin and cytokeratin. In our in vitro model, however, TGFβ1 stimulation did not induce conversion to the myofibroblast phenotype [6].

Cultures treated for 4 days with 5–10–50 ng/ml of TGFβ1 were considered for the microarray studies. We refer to these experiments as "coll-5 ng", "coll-10 ng", "coll-50 ng", "plast-5 ng", "plast-10 ng" and "plast 50 ng". To ascertain the role of the substrate, the expression profile of cells cultured on plastic plates was also compared with that of cells grown on collagen I in basal conditions; we refer to this experiment as "coll-vs-plast".

Each of the seven experiments was performed in duplicate, differentially expressed genes were obtained independently for each comparison and results are given in the Additional file 1. Of ~21000 unique genes of the human oligo chip, we found 977 differentially expressed, represented by 993 spots (Additional file 1).

Since our experiment had a typical dose-response design, the expected effect was a dose-related increase in differentially expressed genes, as shown in Figure 1. Only the up-regulated genes increased considerably, however, going from 5 to 50 ng/ml of TGFβ1, while the down-regulated genes remained much the same in this range. This result did not depend on the cell-growing substrate, though it was more prominent for collagen I. Probably the higher the TGFβ1 concentration, the greater the number of genes exceeding the significance threshold.

Figure 1

Number of differentially expressed genes in each experimental condition. Cells were grown on plastic (plast) or collagen I (coll) plates and exposed to different TGFβ1 concentrations (5–10–50 ng/ml). The last experiment (plast vs. coll) is the comparison between cells grown on the two different substrates and not treated with TGFβ1.

To establish whether the differentially expressed genes in each experiment were coordinately regulated, we compared expression ratios using cluster analysis (Figure 2). We identified four clusters of coordinately regulated genes that were common to the six experiments, demonstrating that the expression profile is similar for all the TGFβ1 concentrations whatever the substrate used. This finding also means that no illegitimate induction took place and seems to support the suggestion that more differentially expressed genes are obtained with the higher dosages because more genes exceed the threshold, not because there are more "new" genes.

Figure 2

Cluster analysis of the differentially expressed genes using TMEV software. The 993 genes in Additional file 1 are clustered on the basis of the log₂ expression ratio. Up-regulated and down-regulated genes are colored in red and green, respectively. Columns refer to the seven experiments performed using different TGFβ1 concentrations (5–10–50 ng/ml) and substrates (collagen type I and plastic). Note that genes are grouped in large blocks, in which they tend to be up- or down-regulated in all the first six experiments.

Given these results, a global list of 919 differentially expressed genes, excluding those identified in the "coll-vs-plast" experiment, was considered for the GO, KEGG and other analyses.

We performed a statistical analysis using the SAM software to assign a significance level to the genes whose expression depended on the TGFβ1 dosage. We identified 25 genes related directly and five genes related inversely to the TGFβ1 dosage with a q value of zero, indicating a strong correlation between dosages and expression levels. The complete list of the genes is given in Table 1. To avoid false negatives, we applied very stringent criteria in the statistical analysis, selecting only genes with a q value of zero – hence their small number. Five of these genes were also identified by Zavadil et al. [7] in time-dependent microarray experiments on human keratinocytes. THBS1, WEE1, PLS3 and ITGAV expression was increased in both our own dose-dependent experiments and in Zavadil's time-course experiments, while MAL expression was decreased, suggesting the existence of a common regulatory mechanism regardless of the cell type.

Table 1 Dose response genes identified using SAM software.

The large number of genes regulated by TGFβ1 makes it difficult to understand the role of this growth factor in the EMT process using a gene-by-gene study. We opted for an analysis designed to comprehensively delineate the genetic programs associated with EMT in response to TGFβ1.

The biological significance of differentially expressed genes was explored using two databases as resources, i.e. the GO and the Kyoto Encyclopedia of Genes and Genomes (KEGG) maps (see Material and Methods).

Among the categories of biological processes (BP), six featured a large number of down-regulated genes and 20 contained a large number of up-regulated genes (Figure 3B). Among the latter, we found proteins involved in cell proliferation, the cell cycle and cell growth regulation, the cytoskeleton and protein transport. Among the cell component (CC) categories (Figure 3A), there was an over-representation of the up-regulated proteins involved in extracellular matrix modification and cell-cell adhesion and of proteins localized in subcellular compartments of the exocytotic pathway, confirming the importance of these pathways in the EMT reprogramming process. The differentially under-expressed genes belong to the CC class of the apical and lateral plasma membranes, indicating the loss of the typical epithelial morphology and functions of cells undergoing EMT.

Figure 3

Gene Ontology annotation of differentially expressed genes. Up-regulated (black bars) and down-regulated (gray bars) genes shown in Additional file 1 were analyzed using GoMiner software. In (A) and (B), genes are classified respectively on the basis of Biological Process and Cellular Component. Only classes with a p value lower than 0.01 (A) and 0.05 (B) are reported.

Up-regulated genes are particularly over-represented in the BP categories such as "morphogenesis" and "cellular morphogenesis", as well as those belonging to developmental processes such as "cell fate determination", "myogenesis", "regulation of development" and "muscle development", indicating that an embryological program could be awakened.

All the genes containing an EntrezGene ID (886) were mapped to the KEGG database to define the signaling pathways and cellular structures that can participate in the EMT process. Only 267 were classified by the software, nevertheless we identified 14 pathways with more than nine differentially expressed genes that were downloaded and saved (Additional file 2). Up- and down-regulated genes were colored red, orange, blue and light blue (see Material and Methods). A visual inspection of these pathways is useful in clarifying the specific (activator or inhibitor) role of the differentially expressed genes, obtaining more detailed information than from the GO analysis alone. In the "apoptotic" pathway, for instance (Additional file 2), both activators (which are up-regulated) and repressors (which are down-regulated) were identified, providing a clear picture of the cell death process.

The signal transduction pathways involved in EMT were better detailed by KEGG analysis, while other processes such as cell adhesion or cytoskeleton remodeling were also identified using the GO. Seven important intracellular pathways were involved in the EMT process in our model, i.e. TGFβ1/SMAD, MAPK, WNT, JAK/STAT, calcium signaling pathways, the cell cycle and apoptosis.

TGFβ1 is known to activate a large group of genes through both Smad-dependent and Smad-independent pathways [8, 9].

We hypothesize that the genes we identified may be closely interconnected in the human interactome, generating a sort of "TGFβ1 interactome". We identified 449 differentially expressed proteins that interact either directly or with undifferentially expressed proteins. The connected component that we highlighted contains 2630 proteins (nodes) and 4183 edges with an averaged connectivity degree of 3.15. The large number of proteins belonging to the cluster identified makes the result noteworthy and indicates that most of the genes regulated by TGFβ1 at transcriptional level during EMT encode proteins that work in a highly coordinated manner.

An analysis performed with the tYNA software [10] confirms that the distribution of the connectivity degree in the network follows the power law characteristic of scale-free networks (data not shown). These networks are typical of cell systems and the majority of nodes have few links in such an arrangement, but a few nodes have numerous links, thereby ensuring that the system is fully connected. A few highly-interconnected nodes thus act as hubs that shape the network's overall operation, but may also represent sites of system vulnerability [11]. We selected the proteins with a connectivity degree – i.e. the number of edges incident on them – higher than 29, calling these proteins "hubs". We identified 27 hubs. The complete list of the hubs is given in Table 2. All these proteins were differentially expressed and three of them, PXN, SMAD3 and TGFBR1, interacted with more than 12 differentially expressed proteins. If we consider these 27 hubs and the proteins interacting directly with them, we identify a sort of central core of the TGFβ1 interactome (network) consisting of 1235 proteins (Figure 4), 115 of which were differentially expressed. Using tYNA software, moreover, we found that the hub proteins had the "shortest path length" and were consequently closely connected within the network – in fact, 24 of them interacted directly.

Table 2 List of the TGFβ1 interactome hubs.

Figure 4

Cluster of interacting proteins. The TGFβ1 network consists of 27 hub proteins with more than 29 edges incident on them and proteins that interact directly with them identified in the "TGFβ1 interactome". The hubs are clearly visible in the network from their high degree of connectivity. Proteins encoded by up-regulated and down-regulated genes are 115 and colored respectively in red and green using a color gradient reflecting the mean log₂ ratio given in Additional file 1, column N. This TGFβ1 core consists of 1235 proteins.

The hubs could be categorized in different ways, e.g. on the basis of a specific GO term being over-represented in the protein group formed by the hub and the interacting proteins. Some nodes, such as SMAD3, CDKN1A and CCND1, interacting with proteins that belong to the same cellular compartment, could be defined as "party hubs" [12]. Other hubs, such as TGFBR1 and PXN, interact with proteins that have a different cellular localization and could be defined as "date hubs". Moreover, TGFBR1 and PXN tend to interact with proteins that act as bridges with other hubs, thus becoming more interconnected than other proteins. As shown by Han et al. [12] in yeast, party and date hubs may have different functions. In particular, date hubs seem to take part in a wide range of integrated connections required for the global organization of biological modules in the whole proteome network.

To validate the microarray data, RT/PCR was performed for six selected genes belonging to the most relevant GO categories identified by our analysis. The RNA samples used for RT/PCR analysis were those used for the microarrays. The up-regulation of tenascin (TNC), fibronectin 1 (FN1), matrix metalloproteinase 2 (MMP2) and connective tissue growth factor (CTGF) and the down-regulation of SMAD3 and collagen IV were confirmed by the RT/PCR experiments. All these genes proved to be similarly regulated by TGFβ1 in all three stimulation experiments.

Discussion

We used the global gene expression profile approach to identify context-dependent markers of the EMT obtained from the long-term TGFβ1 treatment of HUTEC in primary culture. Based on our previous data, we had speculated that the context-dependent EMT process we obtained was a dedifferentiating event. One of the aims of the present study was to further substantiate this hypothesis.

Various studies have shown that genes with a similar expression pattern frequently display common functions and form networks of interacting proteins [13]. Assuming that the genes identified in our experiments belong to the TGFβ1-regulated pathways, we searched for interactions between the proteins encoded by the differentially expressed genes given in Additional file 1. We reasoned that microarray analysis might identify only a part of the complex TGFβ1 network, due to other effects such as post-transcriptional regulation, so we used protein-protein interaction data to identify proteins interacting with those encoded by differentially expressed genes. We obtained a single connected component (interactome) consisting of 2630 proteins and containing 449 differentially expressed proteins that interact directly or with undifferentially expressed proteins. This analysis is extremely useful not only for detecting the network of interacting proteins that respond to TGFβ1, but also for identifying the network hubs, i.e. proteins with a high degree of connectivity that could have a crucial role in response to TGFβ1. We identified 27 hubs with more than 29 edges incident on them and encoded by genes found differentially expressed in our experiments. Of the three hubs identified as having more than 12 interactions with differentially expressed genes, SMAD3 was classified as a party hub while TGFBR1 and PXN were considered as date hubs. As shown by Han et al. [12] in yeast, party and date hubs might have markedly different global properties in the interactome network. In accordance with this model, TGFBR1 and PXN, which interact with proteins that have a different cellular localization, might represent global or "high level" connectors between different biological modules and SMAD3, which functions within a module, works at a "lower level" in the organization of the proteome.

Among hubs, thrombospondin 1 (THBS1) – the principal activator of the TGFβ1 peptide [14] – emerges because it links a wide range of matrix proteins mediating their interaction with cell-surface receptors. Its central role is confirmed by the 35 links connecting this protein to the other proteins in the interactome map. It could be considered as a date hub in the TGFβ1 interactome, according to the definition given by Han et al. [12]. In Figure 5, the topological connections of THBS1 are magnified (A) and the differentially expressed genes interacting with it are indicated (B). Moreover, THBS1 is the gene that has the most significant correlation with TGFβ1 dosage (Table 1).

Figure 5

Topological connections of the date hub THBS1. Thrombospondin 1 (THBS1) links directly 35 proteins in the TGFβ1 interactome: A) the topological connections are shown; B) the proteins encoded by differentially expressed genes are indicated and colored respectively in red and green using a color gradient reflecting the mean log₂ ratio reported in the square.

THBS1, PLS3 and ITGAV are known to have a pivotal role, both in cytoskeleton remodeling and in de novo extracellular matrix synthesis, two of the most relevant processes in EMT. The importance of these processes is highlighted by the considerable number of differentially expressed genes in the "extracellular matrix" and "cytoskeleton" GO categories (Figure 2).

In the "extracellular matrix (EMC)" category, the different regulation of the various collagen types emphasizes the transition from an epithelial-like to a mesenchymal ECM. The up-regulation of TNC and FN1, one of the hubs of the TGFβ1 interactome, confirmed the switch to a mesenchymal ECM. The laminin α chain was found differentially down-regulated in our model, as in other EMT models [5], while the β1 chain is up-regulated, as occurs in the precondensed mesenchyma [15], suggesting that an embryonal mesenchymal ECM is being synthesized. This suggestion is reinforced by the types of integrins up-regulated during EMT. In fact, the "ECM receptor interaction" of the KEGG pathways (Additional file 2) reveals not only the up-regulation of integrins generating the network for mesenchymal cell adhesion (ITGA5 and ITGAV) [16], but also – and above all – the strong over-expression of ITGA11, which constitutes the receptor for the interstitial collagens involved in cell migration and collagen reorganization on mesenchymal cells during development [17].

The "cytoskeleton" category is characterized by a large number of up-regulated genes, some of which have a crucial role in cytoskeleton remodeling. It is worth noting that four hubs of the TGFβ1 network, i.e. filamin A, alpha (FLNA), alpha 1 actin (ACTA1), paxillin (PXN) and beta actin (ACTB), belong to this category.

"Morphogenesis" is the GO category with the largest number of up-regulated genes, as shown clearly in Figure 3. The significance score assigned by the GOMiner software is not very high, however, because this category is one of the most represented in the microarray chip. It should nonetheless be noted that some of the genes belonging to this category are among the top ten up-regulated differentially expressed genes (Additional file 1). Stimulated by TGFβ1, HUTEC seem to reactivate the developmental processes: this may point to a sort of stemness of tubular cells that enables them to dedifferentiate when stimulated (kidney repair and maintenance), but also to reawaken a silenced embryogenetic program. Sox 11, GADD45B, N-cadherin (CDH2), Activin A (INHBA), CTGF, FGF1/5, Angiopoietin (ANGPTL4), natriuretic peptide precursor B (NPPB), calcitonin receptor (CALCR) and caldesmon 1 (CALD1) are the most over-expressed genes belonging to the "morphogenesis" category. That tubular cells have the ability to reactivate an embryogenetic program has been demonstrated very recently by Kitamura et al. [18] and Maeshima et al. [19], who identified and isolated cells of the S2 segment with a potential tubulogenic ability and a capacity for integration in the developing kidney. Our results are in line with their findings.

CD133+ cells have been found to have stem cell potential in the adult kidney [20]. Hypothesizing that CD133 up-regulation should occur if a stemness property is activated by TGFβ1, we looked specifically for prominin 1 (CD133) mRNA activation in our EMT model. Surprisingly, we found CD133 down-regulated, indicating that CD133 is expressed in control conditions. This finding supports the view taken by Florek et al. [21], who showed that prominin 1 transcript and the alphahE2 epitope immunoreactivity of CD133 (obtained using a novel antibody instead of AC133) occur in several adult tissues and in the proximal tubular cells of the adult kidney in particular.

Signaling members of the Wnt (WNT5B) and FGF (1/5) families and transcription factors such as Sox 11, known to have a crucial role in nephrogenesis and cell fate determination during kidney development [22, 23], were found up-regulated. The involvement of the Wnt pathway was also confirmed by KEGG analysis.

The non-canonical Wnt signaling pathway centered on WNT5/Ca²⁺ [24, 25] seems to be activated through Nemo-like kinase (NLK) in our EMT model and to antagonize the canonical beta-catenin Wnt signaling (Additional file 2). In fact, we observed the down-regulation of both wingless-type MMTV integration site family, member 2B (WNT2B) and catenin, beta 1 (CTNNB1) and the up-regulation of both WNT5B and the calcium signaling pathway (see KEGG pathways). This effect seems to be reinforced by the up-regulation of dapper, antagonist of beta-catenin, homolog 1 (DACT1), which is a known beta-catenin Wnt signal inhibitor. It has been reported that inhibiting the beta-catenin system strongly inhibited TGFβ1-induced αSMA expression in tubular cells [26]. Our data confirm this finding, since no αSMA expression was triggered by TGFβ1 in our EMT model [6].

The non-canonical Wnt pathway, which includes planar cell polarity – an important process in embryonal axis development involving cytoskeletal polarity, as well as in the calcium pathway regulating cell adhesion [24]-, was thus found up-regulated in our EMT model, reinforcing the idea that an embryological program is awakened. Very recently, Osafune et al. [25] reported that the capacity of renal progenitor cells of the metanephric mesenchyme to form colonies in vitro and undergo mesenchymal-epithelial transition is positively regulated by planar cell polarity pathways downstream from Wnt.

Although the canonical Wnt signaling was repressed, the final effector of the pathway – and one of the most important – Cyclin D1 (CCND1) was up-regulated, whereas Cyclin B2 (CCNB2) – which is assumed to bind to TGFβ R2 and thus play a key part in TGFβ-mediated cell cycle control – was down-regulated.

Apoptosis and EMT are two distinct and opposite signal modules for TGFβ1-downstream effects. There is growing evidence that SMAD3 is an important signaling anchor for the apoptotic network for TGFβ1 too. In particular, the loss of SMAD3 function due to a decrease in its expression might be a requirement for epithelial cells to survive in the presence of prolonged TGFβ1 stimulation [27]. This was also confirmed in our EMT model: visual inspection of the TGFβ-SMAD KEGG pathway (Additional file 2) reveals what we demonstrated previously using RT/PCR analysis [6], i.e. that Smad signaling was down-regulated. Despite the number of up-regulated inducers (TGFβ1, INHBA), the key transducers are all down-regulated, as are the Id genes (the effectors of cell differentiation). The ID2B gene was specifically down-regulated. Mad expression and ID2 down-regulation are important events in the TGFβ1 cytostatic program in epithelial cells and ID2 suppression by TGFβ1 is essential for EMT to occur [28, 29].

The central role of SMAD3 is also demonstrated by its position in the TGFβ1 network: it is one of the hubs, most likely a date hub since it works within a single module (Figure 6), at a low level of network organization. This might explain why apoptosis seems to be induced in our model, despite SMAD2 and SMAD3 down-regulation. In fact, a visual inspection of the KEGG apoptosis pathway clearly shows the up-regulation of caspase 3 (CASP3), a known inducer of cell death, and the under-expression of BCL2 and BIRC3 (IAP gene in the KEGG map), which counteract this action.

Figure 6

Topological connections of the party hub SMAD3. SMAD3 links directly 119 proteins in the TGFβ1 interactome: A) the topological connections are shown; B) the proteins encoded by differentially expressed genes are indicated and colored respectively in red and green using a color gradient reflecting the mean log₂ ratio reported in the square.

On the other hand, the up-regulation of several genes implicated in the cell cycle pathway (Additional file 2), such as CCND1, GADD45, YWHAG (tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein), indicates that cells are entering the cell cycle. In this pathway, however, the up-regulation of wee1 tyrosine kinase (WEE1), one of the genes strictly regulated by TGFβ1, seems to indicate a sort of actual control of cell proliferation, so both apoptosis and cell cycle entry seem to be strictly controlled in the EMT process. Indeed, neither apoptosis nor proliferation were seen in our model by immunocytochemistry [6]. It is tempting to speculate that a concerted and strictly controlled action between signals for cell death and proliferation might be taking place in cells after long-term TGFβ1 exposure, miming the developmental process in which morphogens such as TGFβ1 act not as positive regulators of cell differentiation but as key regulators of cell survival [30]. It is worth noting that key proteins of the Wnt signaling and apoptosis/cell cycle control pathways (i.e. CTNNB1, CCND1, CASP3, BCL2) were found to constitute some of the hub proteins of the TGFβ1 network.

Finally, KEGG analysis of microarray data highlighted that RAS/MAPK signaling was the principal downstream effector of chronic TGFβ1 stimulation in our EMT model, confirming the suggestions advanced by other authors [5, 8], i.e. that both Smad-dependent and Smad-independent signaling cascades are activated by TGFβ1 and that they regulate mesenchymal transition in a context- and cell-dependent manner. The MAPK signaling pathway has an important role in connecting the signal triggered by TGFβ1 to important downstream processes such as apoptosis/proliferation and the Wnt cascade [5]. Our results confirm reports from other authors on the role of this signaling in other EMT processes [31–33].

Conclusion

Understanding how mesenchymal cells arise from epithelial cells could have a strong impact in unveiling the mechanism behind fibrosis and cancer progression. Moreover, it might reveal mechanisms of epithelial cell plasticity underlying kidney regeneration and repair.

In the kidney, tissue regeneration and repair occur through three, not mutually exclusive, cellular and molecular mechanisms: differentiation of the somatic stem cells, recruitment of circulating stem cells and, more importantly, proliferation/dedifferentiation of mature cells. Dedifferentiation seems to represent a critical step for the recovery of tubular integrity and precedes the reconstitution of a well-differentiated morphology. Understanding the cellular and molecular events involved in renal tubule regeneration is indispensable to design cell-based and other therapeutic strategies in order to potentiate this innate capacity. EMT is now considered a part of tubular cell plasticity.

The purpose of our study was to substantiate our initial hypothesis that the process of EMT induced by TGFβ1 chronic exposure in HUTEC is a dedifferentiating process. Our earlier results suggested that this might be so and our present findings support that impression. In fact: 1) the principal functional category involved in the EMT process concerns morphogenesis and development; 2) the most up-regulated genes belong to this category; and, finally, 3) some intracellular pathways are involved, whose engagement during kidney development and nephrogenesis is well known. Our results thus suggest that the recapitulation of embryological programs is an integral part of the EMT process on long-term exposure to TGFβ1 and that tubular cells may have the capacity, under appropriate environmental cues, to redifferentiate not only back to an epithelial type but also towards another cell type, i.e. myofibroblasts or endothelial cells. Judging from our data, moreover, TGFβ1 seems (in the context of our experiments) to act as a morphogen regulating cell survival by means of strictly balanced signals for cell death and proliferation.

Finally, our findings are the first to show that genes involved in the TGFβ1-driven EMT process are highly interconnected and topologically related in the human interactome map. They generate a single scale-free network whose hub proteins were found differentially expressed, pointing to a crucial role for them in the EMT process. The main role of one of them, thrombospondin 1, emerges from its high degree of regulation by TGFβ1 and from the 35 links connecting this protein to the other proteins in the interactome map.

Methods

Cell cultures

Primary HUTEC cultures were established as explained elsewhere [6]. Cells at passage 1 were used for TGFβ1 stimulation experiments, designed to monitor the effect of TGFβ1 at both phenotypic and molecular levels simultaneously [6]. Cells were seeded at subconfluence and incubated at 37°C in a 5% CO₂ atmosphere for 24 hrs under quiescent conditions (1% FCS in RPM1 1640) in 6-well plastic or collagen I-coated plates for RNA extraction. Cells were cultured for 4 days in the presence of 1, 5, 10, 50 ng/ml human TGFβ1 (Prepro Tech EC, London, UK). Stimulation experiments were conducted in triplicate and morphological, immunocytochemical and molecular analyses were performed. Control conditions were represented by cells maintained for 4 days in 1% serum without TGFβ1.

RNA extraction and quality control

Total RNA was extracted using RNAzolB (BIOTEX, Houston, USA) according to the protocol.

The Agilent 2100 Bioanalyzer (Agilent Technologies, PaloAlto, USA) was used to assess RNA integrity. RNA was quantified with the UV/VIS Spectrometer (Lambda 2S PerkinElmer).

Microarray platform

Microarray expression was analyzed using the Operon 70 mer oligos collection (Human Version 2.0) containing 21,329 oligonucleotides spotted in duplicate (GPL2136136 record in the GEO database) on MICROMAX glass slides – SuperChip I provided by PerkinElmer Life Sciences Inc. (Boston, USA).

Oligos were printed using the Biorobotics Microgrid II spotter; spots were spaced 115 μm and microarrays consisted of 48 subarrays.

RNA amplification and labeling

1 μg of total RNA was amplified using the Amino Allyl MessageAmp™ aRNA Kit (Ambion, Austin, USA). Five μg of amplified RNA (aRNA) were labeled with Cy3/Cy5 fluorophores using CyDye Post Labeling Reactive Dyes (Amersham Biosciences, Piscataway, USA).

Following purification, dye incorporation was quantified by spectrophotometric analysis.

Hybridization

Approximately 2.5 μg of aRNA labeled with 100 pmoles of fluorophore were used for each hybridization.

Labeled aRNA was precipitated using NH₄Ac and EtOH following standard protocols and resuspended in hybridization buffer (5× SSC; 0.1% SDS; 40 ng/μl SS-DNA; 25% formamide).

Microarray slides were pre-hybridized in GeneMachines™ chambers (Genomic Solutions, Cambs., UK) for 2 h at 48°C with 70 μl of pre-hybridization buffer (5× SSC; 0.1% SDS; 40 ng/μl SS-DNA; 1× Denhardt's solution) using a coverslip.

Slides were washed with water and dried with compressed air.

Hybridization was carried out using Hybridization Station ArrayBooster™ (Advalytix, Brunnthal, Germany) at 48°C for 12 h.

Microarrays were washed with 1× SSC plus 0.2% SDS for 4 minutes, 0.1× SSC plus 0.2% SDS for 4 minutes, twice with 0.2× SSC for 4 minutes, and twice with 0.1× SSC for 3 minutes.

Microarray scanning and image analysis

Microarrays were scanned using ScanArray Lite (PerkinElmer™). Images were analyzed using ScanArray Express (PerkinElmer™).

Statistical analysis

The following procedure was used to remove spots with a low fluorescence intensity or high variability between replicates.

1) "Intensity-dependent calculation of standard Z-score"

spots with a median fluorescence pixel intensity below 700 (calculated considering negative control intensity) on both Cy3 and Cy5 channels were filtered out; those with a median fluorescence pixel intensity of zero or less in only one channel were set to 100 to prevent their elimination during normalization. Files were saved in "tav" format to make them suitable for reading with MIDAS software [34] and normalization was done using the LOWESS method.

Two different procedures (named A and B) were used to eliminate outliers, as follows.

(A) Based on the method suggested by Yang et al. [35], we calculated R1 = (CH1 intensity/CH2 intensity) and R2 = (CH1 intensity/CH2 intensity) values for the two replicates of the same gene on the microarray and log₂(R1/R2). We indicated the two replicates of the spot as R1 and R2. Then we calculated the mean and SD (standard deviation) for the log₂(R1/R2) values of all microarray spots. Those with a log₂ ratio higher than |3 SD| were rejected due to replicate inconsistency.

The geometric mean for the two replicates of the remaining genes was calculated and the output files were saved in "tav" format.

(B) The "tav" file for each microarray experiment was normalized using MIDAS software and geometric mean values underwent SLICE data analysis, considering only those where |log₂(CH1 intensity/CH2 intensity)|≥1.5 SD.

Each experiment was performed in duplicate using a dye swap procedure and only the genes that independently complied with these filters on both replicates were considered.

2) "SAM analysis"

data were also analyzed using SAM software, but intra-array replicates were not averaged, intensity fluorescence filtering was as described previously and normalization was done using the LOWESS procedure. The four replicates (two intra-array and two deriving from the dye-swap analysis) were t-tested using the SAM software and considering the lowest False Discovery Rate.

In all, 1203 spots were considered at this stage: 589 of them satisfied both statistical procedures 1) and 2) in at least one of the seven experiments (collagen 5–10–50, plastic 5–10–50 and collagen vs. plastic) and 614 spots were identified by the SAM software alone. The geometric mean was calculated from the four replicates of each gene and is given in Additional file 1 (columns G-L). Some genes, especially those identified only by the SAM software and in the "50 ng of TGFβ1" experiments, had a low (<0.65) log₂ ratio and were filtered out. We finally obtained 993 differentially expressed spots, corresponding to 977 genes (some oligos refer to the same gene), given in Additional file 1. This table consists of two parts: spots 1 to 554 are the most significant because they satisfy both statistical procedures 1) and 2); genes from 555 to 993 are less significant because they only passed the SAM analysis.

We also used the "quantitative response" option of the SAM software to check for the presence of genes whose expression is regulated in a dose-dependent manner (Table 1).

Validation of microarray data

RT/PCR-Real Time PCR

six differentially expressed genes belonging to the different GO categories found involved in the EMT process, i.e. TNC, FN1, collagen IV, MMP2, SMAD3 and CTGF, were analyzed. Quantitative comparative RT/PCR and Real Time RT/PCR were performed, as reported elsewhere [6].

Immunocytochemistry

immunocytochemistry was conducted, using antibodies against αSMA, cytokeratin 8–18, vimentin, collagen III, Ki67 and E-cadherin, as described in [6].

Annotation checking procedure

The annotation of the Human Version 2.0 microarray platform (OPERON) was out of date, so we double-checked it using the MatchMiner software [36], loading the GenBank accession numbers and recovering the EntrezGene ID and the gene symbol (only for the genes in Additional file 1). If the GenBank clone had been withdrawn, we performed a BLAST search with the oligo sequence on the UCSC Genome Browser [37] and we updated the GenBank ID. This search was also conducted for the RefSeq sequences with multiple associations. Oligos matching more than one gene were rejected. To verify gene symbols and gene descriptions, we downloaded the "gene_info.gz" file from the NCBI ftp site [38].

Functional analysis of microarray data

Geometric mean expression values of the 993 spots identified were obtained for each of the seven experiments and these data were clustered using TMEV software [34].

Functional categories were identified with the GoMiner software [39] using the Gene Ontology (GO) annotation [40]. We selected a p value threshold of 0.02 for Biological Process categories and 0.05 for the Cellular Component.

To identify KEGG pathways involved in the TGFβ1 response we colored up- or down-regulated genes with the tool available on the KEGG site [41, 42]. Differentially expressed genes that only passed the SAM analysis were colored orange (up-regulated) or light blue (down-regulated); genes that passed both statistical analyses (see Methods) were colored red (up-regulated) or blue (down-regulated).

Mapping the differentially expressed genes in the human interactome

We used data from the literature on protein-protein interactions to verify whether the genes in Additional file 1 are independent or tend to be associated with one another in one or more clusters. We downloaded the protein-protein interactions from three different sources, i.e. (A) the NCBI ftp site [43] (interactions.gz), from which we selected the interactions identified in Homo sapiens and used a Perl script to assign the official symbol to the known proteins; (B) the BioGRID database [44]; and (C) the Biomolecular Interaction Network Database (BIND) [45]. Each data set was loaded separately in the Cytoscape software [46] and three distinct "global" protein-protein interaction networks were generated. For each global network we mapped the genes identified in our microarray experiments (Additional file 1) and selected proteins that interact with them directly, generating a small group of sub-networks. We selected only the three largest sub-networks and they were merged into a single subnetwork using the "union" function in Cytoscape. The network we obtained is composed of 2630 nodes and 4183 edges. We used the "filter option" in Cytoscape to select the nodes interacting with more than 29 proteins. These nodes were called hubs. Gene Ontology analysis was performed on the network using BINGO software [47], while network topological statistics were obtained using tYNA software [10]. The interactome was also deposited in the tYNA database [48] (TGFbeta-network-231106). Microarray experiments were loaded into ArrayExpress with accession number E-MEXP-566 according to the MIAME rules [49].