Reagents and Cells
Recombinant VEGF was purchased from R&D Systems (Minneapolis, Minnesota). Human umbilical vein endothelial cells (HUVECs) and Normal Human Dermal Fibroblast (NHDF) cells and culture medium were purchased from Cambrex (Walkersville, Maryland). Cells were maintained in EGM with 10% FBS and FGM-2, respectively for HUVECs and NHDFs, according to manufacture's instructions. In vitro angiogenesis co-culture AngioKit and Optimized medium were purchased from TCScellworks (Buckingham, UK).
In vitro angiogenesis assay
In vitro HUVEC cord formation co-culture was performed using the AngioKit commercially available from TCScellworks (Buckingham, UK) adapted to a 24-well format. To test HUVEC and NHDF passage dependence effects, the co-culture was done in house adapted to a 96-well format. Briefly, HUVEC cells were co-cultured onto human dermal fibroblast feeder layers in Optimised medium (TCScellworks, Buckingham, UK) for a total of about 12 days. Treatment started the day when the plates were received; usually three to four days after initial plating. Medium with or without growth factors or compounds was replenished every 2 to 3 days. Cells were fixed with 70% cold ethanol at the end of the culture and processed for anti-CD31/Alexa 488 immunofluorescence. To visualize feeder layer fibroblasts, some cultures were stained using a polyclonal antibody against Vimentin (Abcam, Cambridge, Massachusetts) along with an Alexa 568 conjugated secondary antibody to examine the feeder cell integrity. Nuclear staining was performed using Hoechst staining (Invitrogen Molecular Probe, Carlsbad, California).
Immunofluorescence stained co-culture plates were scanned using ArrayScan VTI (Cellomics, Pittsburgh, Pennsylvania) and the multi-parameter tube formation Bioapplication was used for quantification. Morphological measurements included cord length, width, area, percent connected cord, branching node and segment count, nuclei per cord, angiogenic index (expressed as the percentage of area covered by connected tube), etc. General cytotoxicity was monitored by examination of the nuclear staining of the feed layers.
Attenuation of gene expression by siRNA
HUVECs were trypsinized by 0.05% Trypsin/EDTA (Invitrogen, Carlsbad, CA), and 1 × 106 cells were nucleofected with siMYC (200 pmol each, Ambion, Austin, TX) using a nucleofector device (Amaxa, Gaithersburg, MD) adopting optimized A-034 protocols for HUVECs according to the instructions from the manufacturer. Transfected cells were appropriately diluted and plated in 96-well plates onto the NHDF feeder cells for the initiation of the cord formation assay. Cells were also plated in parallel into 6-well plates and collected in Trizol (Invitrogen, Carlsbad, CA) 48 h or 72 h later for total RNA extraction.
Real Time RT-PCR
Total RNAs were extracted from Trizol and cleaned by RNeasy mini kit (Qiagen, Chatsworth, CA). Genomic DNA contamination was removed by DNA-free kit (Ambion, Austin, TX), and cDNA synthesis was done by using a high capacity cDNA archive kit (Applied Biosystems, Forster City, CA). The cDNAs were used as templates for real-time PCR with a universal PCR master mix (Applied Biosystems, Foster City, CA). Real time PCR was performed on a ABI 7900 HT. The average values were normalized to 18s rRNA. Delta CT method was applied to calculate the relative gene expression level [40].
RNA isolation and microarray experiments
RNA isolation and microarray studies were performed as described previously [41]. Briefly, HUVECs and NHDFs with different passages were cultured to about 80% confluency. RNA was isolated using Trizol™ (Life Technologies, Inc. Invitrogen, Carlsbad, California) according to manufacture's instructions, and followed by cleanup using RNeasy spin columns (Qiagen, Valencia, California). RNA labelling was performed on the Onyx robotic system from MWG-Biotech and Aviso using the standard labelling protocol. Biotin labelled cRNA was fragmented and hybridized to human whole genome U133 plus 2 GeneChip (Affymetrix, Santa Clara, California). Chip processing, imagine capturing, and raw data analysis were performed using the Affymetrix Microarray Suite MAS with default parameters.
Statistical analysis to identify differentially expressed genes
Hybridized microarrays were analyzed by the Affymetrix microarray analysis software microarray suite 5 (MAS 5). The chip signals were normalized using the default method in Affymetrix MAS5, but setting the target at 2%-trimmed mean to 1500 instead of 500. HUVEC were analyzed with four passages. We designated the cell passages 5, 8, 12 and 22 as 1, 2, 3 and 4 respectively. For each probeset, the signal data were fitted to a linear regression model of time to identify consistent gene expression changes over cell passage. Three batches of the NHDF cells were analyzed, including early passage neonatal NHDF, early passage adult NHDF, and late passage adult NHDF cells. To reduce a possible batch to batch effect, the MAS 5 signals for each chip were re-normalized by Local Polynomial Regression Fitting (loess) approach using one neonatal sample as the baseline on a log scale. The re-normalized signal data were fitted to an ANOVA model to identify differential gene expression changes between the different fibroblast cells.
To control for a false positive rate of testing the expression change of thousands of genes simultaneously, the false discovery rate (fdrate, or FDR) was estimated using an algorithm derived from Benjamini and Hochberg [42]. Probesets with a false discovery rate of 0.2 or smaller were considered as significant and followed up for further analysis. We excluded probesets called "absent" in all chips in a data set. These probesets usually have very low expression signals. In addition, we used expression fold change between different treatment conditions to filter out genes to make sure the gene expression changes are biologically meaningful.
PCA analysis was performed in R using the principal component analysis function with standardized data [43]. The first three principal components with > 80% variations were exported into Spotfire (Sommerville, MA) to aid in creating sample scatter plots. Hierarchical clustering analysis of DE genes was done in Spotfire with Euclidean distance by the complete linkage method.
To identify biological processes that were significantly changed with passages, the identified DE genes were analyzed using the DAVID functional annotation system (National Institute of Allergy and Infectious Diseases, NIH) [21]. Significantly changed biological processes were determined by the enrichment test provided by the DAVID system. To control for a false positive rate of testing the expression change of thousands of genes simultaneously, the false discovery rate (fdrate, or FDR) was estimated using an algorithm derived from Benjamini and Hochberg [42].
Construct angiogenesis related interactome and calculate connectivity index
A tentative interactome was constructed using PathwayAssist from the union of the differentially expressed genes identified either from the HUVEC or the NHDF array experiment. PathwayAssist was configured to retrieve five different types of relationships defined in the software, namely, binding, expression, molecular transport, protein modification and regulation. The software introduced non-differentially expressed genes into the interactome in order to establish the shortest connection pathway between any two differentially expressed genes in the initial input gene list. This is an informative function since it could recover potentially non-transcriptional regulatory relationships between two or more genes. We cleaned up this initial interactome by removing genes that were not detected by either array study. In addition, genes that could not connect to any other genes were called orphans and removed from the interactome. Thus, the final version of the interactome we constructed can be described as follows:
Let G be the set of K genes in the interactome I, and R be one of the five relationships: binding, molecular transport, protein modification, inhibition or up-regulation. Thus, the interactome I we built can be formally described as
I = {(g
i
, g
i
) ∈ R | g ∈ G; i, j = 1,2, ... K} (1)
We assumed R is transitive, that is, if gene A regulates gene B, and gene B regulates gene C, it is assumed that gene A regulates gene C. Except for binding among the five types of relationships, molecular transport, protein modification, up-regulation and inhibition are unidirectional.
We wrote a "network walking" algorithm (Figure 6) in a recursive function with which we calculated the number of unique genes, L
ij
, connecting to the ith gene at the jth step along the connection path in the interactome. We define connectivity of the ith gene at the jth step along the connection pathway as
Ci,j= Ci,j-1+ L
ij
/2j-1 (2)
Thus, for the ith gene in the interactome its connectivity at the jth step in its connection pathway is the connectivity at the (j-1)th step plus weighted number of unique genes at the jth step. We assigned 1/2j-1as an arbitrary weight to the number of unique genes connecting to the ith gene at the jth step, that is, the weight exponentially decreases with the number of steps away from the ith gene. Obviously Ci,jis a non-decreasing function. We further define connectivity of the ith gene C
i
as
(3)
where C
i
is the plateau of connectivity of the ith gene and N
i
the number of steps Ci,jtakes to arrive at its limit or the plateau. Since every gene in the interactome is connected, the total number of genes a gene connects to is the same, that is, after N
i
steps the connectivity of the ith gene will arrive at its plateau. However, they differ in the number of steps N
i
or the rate of the arrival at the plateau. Some genes may take more steps than the others to arrive at the platueas, depending upon what and how many genes it connects, especially at early steps. Schematic elucidation of the algorithm is shown in Fig. 6.
Simulation of network perturbation
Let S be a small molecule kinase inhibitor, and Ps, Dbe the level of network perturbation by the small molecule at a hypothetical dose D. Ps, Dis estimated as
(4)
where K is the number of kinases in the interactome I that can be inhibited by the small molecule S, IC 50i,sis the IC50 value of the small molecule S against the ith kinase. Their relative IC50 values against various kinases were measured at Upsate (UBI) or in house. To simulate the dose response curve, we set the theoretical dose D from 0.1 nM to 2 μM. If a term in (4) is more than the maximum number of the connected genes, the term is then set to the plateau. Summed network connectivity was used to calculate the percent perturbation with the plateau to be set at 100% network perturbation. The resulted simulated dose dependent network perturbation was then imported into the program GraphPad Prism to estimate EC50 of network perturbation by the small molecule S by performing Sigmoidal dose-response curve fitting.