Microarray Data Sets

Previously published mouse embryonic fibroblast (MEF) cell microarray datasets were used in our analyses [4]. Briefly, the cells were synchronized by either serum starvation or hydroxyurea treatment. We used the data sets obtained from the serum starved cells for the analysis of re-entry into G1 from G0 (0, 6, 12, 15, 18, 21, 24 hours after serum starvation), and those from the hydroxyurea exposed fibroblasts for the G1-S analysis (0, 3, 6, 9, 12, 15, 18 hours after the treatment). The detailed methods used to obtain these microarray data have been previously described [4].

Selection of the Subset Database

The original gene expression data, comprising about 6437 genes, were screened for genes that showed at least a 2.0-fold change (up- or down-regulation) using GenMAPP [24]. The distribution and frequency of the fold changes (relative to the time 0) at each time point were analyzed by MAPFinder 1.0 beta, an accessory tool of GenMAPP, to identify the optimal biological maps. From this collection of maps, we selected those related to cell cycle processes that had a “z” score greater than 1.95 (the z score represents the difference between the observed number of genes meeting the criteria and the expected number of genes meeting the criteria in each map based on gene ontology). As detailed in Table 1, 10 maps were selected based on gene ontology (denoted MAPP) and the relationship to the cell cycle. A subset of 39 genes was chosen from among the MAPP maps selected (Table 2). The abbreviated names of the genes that were analyzed in this report are presented according to the displays listed in GenMAPP.

Mathematical Models

We applied the expression-associated network modeling method previously developed by Yamanaka et al. [21] to the fold-change data from the gene-expression data sets. This method falls under the general area of Bayesian networks, with a likelihood-based selection algorithm used to identify the most promising networks. In general, if X₁, X₂….X_p represents the data obtained for p genes, N denotes a network, and θ denotes parameters in that network, with the likelihood given by:

1. The choice of the best network would be the one with the largest value of the posterior density at the chosen network topology; that is
   find
   
   (2)  Nˆ=arg maxNfN|DN|D
   where
   
   (3)  fN|DN|D∝fNN⋅fD|ND|N
2. The Bayesian network used in this analysis had the following assumptions:
   1. 
      (4)  fNN⇒uniform distribution
   2. 
      (5)  fD|θ,ND|θ,N=∏j=1p∏i=1nfXj|N,θjXji|paXji,N,θj
   3. 
      (6)  fθ|Nθ|N=∏j=1pfθj|Nθj|N

where pa (X_ji) is the collection of genes that link to the j^th gene in the network (a pathway).

with these assumptions,

Thus, it is possible to focus on each gene rather than the whole network and still obtain a global optimum. To quantify rates in the gene-expression network, we used the Bayesian methods developed by Toyoshiba et al. [22, 23]. Supposing that X_i (i=1,2,3…p) represents the natural log of the relative ratio, the functional relationships between the genes could be characterized using the log-linear model below:

where I_ji is an indicator function (-1, 0, 1) characterizing the effect from G_j to G_i, T represents a matrix having I_ij as the (i,j) element, and β_j.=[β_j1, β_j2β_j3β_j4... β_jp] is the vector in which each β_ji is the magnitude by which one unit of gene X_j will affect the expression levels of gene X_i. Thus, if Iji is not equal to 0, Pa(Xi) contains X_j.

If f (X_i| T, θ ) is defined as the distribution of gene expression in the given model, then the likelihood is written as

where θ represents the parameter vector in the model.

By Bayes’ theorem, the prior distribution is given by

The posterior distributions f_θ|X,T were evaluated using the MCMC method. In our analyses, f_X|T,θ was assumed to be normal, with a mean defined by equation (8) and a random variance whose prior distribution was assumed to be uniform with 0 as the lower bound and twice the maximum STD for each gene distribution. The prior distribution for θ and f_θ was assumed to be lognormal with a mean of 0 and a variance of 1.0.

The MCMC analysis was applied as described in [23] and [22]. A typical MCMC run was 100,000 samples with the first 20% of the samples discarded to “burn in” the algorithm. Some runs were much longer depending on convergence and stabilization of the resulting posterior distributions.

The model described in this section is an analysis tool and is not intended to characterize the mechanisms by which the different genes are linked. Instead, it is intended to find the most prominent linkages between cells to provide hypotheses that can be further explored and later modeled mechanistically.

Visualization of Gene Networks and Clustering Analysis

We used a MATLAB script newly developed by Parham et al. (unpublished data) to generate transcriptional regulatory networks using MATLAB version 6.5 (The MathWorks, Inc., Natick, MA).

Establishment of Mouse CDK7 Recombinant Retrovirus

A full length cDNA fragment of mouse cyclin-dependent kinase (Cdk)7 (NIH Mammalian Gene bank accession number: NM_009874) was obtained by RT-PCR from the total RNA extracts of 13.5 day mouse embryo using a previously described method [25]. A hemagglutinin (HA) protein tag sequence was then introduced at the carboxyl terminus of this mouse Cdk7 cDNA using a tailed PCR method. The Cdk7 cDNA was next subcloned into the EcoRV site of pBluescript SKII+ (Stratagene, La Jolla, CA) by blunt end ligation, and the resulting constructs were validated using a cycle sequencing reaction in an ABI 310 genetic analyzer (Applied Biosystems, Foster City, CA). The subcloned Cdk7 cDNA fragment was then transferred into the multiple cloning site of an LXIN retrovirus vector (Clontech, Mountain View, CA). Both empty LXIN vector and LXIN vector harboring the mouse Cdk7 cDNA were introduced into PT67 retrovirus packaging cells (Clontech) using Fugene6 (Roche, Basel, Switzerland). Infected cells were then selected with 1mg/ml G418 (Invitrogen, Carlsbad, CA ) in the growth media for one week.

Measurement of the Retrovirus Titers in the Producer Cells

Conditioned medium from the producer cells was diluted 1:10 and 1:50 with DMEM containing 10% calf serum, and then used for the infection of NIH3T3 cells to measure the titer of the synthesized retrovirus. NIH3T3 cells were grown in media with diluted retroviruses for two days under the same conditions that are described below for mouse embryonic fibroblasts. The infected NIH3T3 cells were diluted 1:100 and 1:1000, and then selected with 1mg/ml G418 for one week. Retrovirus titers of the original conditioned medium were calculated based on the number of colonies demonstrating G418 resistance.

Preparation of MEF Cells that Expresses the Recombinant Mouse CDK7

Mouse embryonic fibroblasts (MEF) were prepared using a previously described method [25]. Second passage primary fibroblasts at a 70% confluency were infected with conditioned medium containing PT67 producer cells at a 1:2 dilution 1:2 with basal MEF medium for 2 days in the presence of 1µg/ml polybrene (Sigma-Aldrich, St. Louis, MO). When the infected MEF cells reached confluence, they were diluted 1:5 as above and selected with 200µg/ml G418 for one week. Control experiments confirmed that non-infected MEF cells did not survive in the presence of 200µg/ml G418 (data not shown). Infected cells selected with G418 were subjected to lysis and protein extraction for western blot analyses.

Western Blot Analyses of MEF Cells Exogenously Expressing Mouse CDK7

Total proteins were isolated from MEF cells infected with either control or Cdk7 recombinant retrovirus using a standard methodology [25]. Heat-denatured proteins were separated by 10% SDS-PAGE and the proteins in the gel were transferred to polyvinylidene difluoride (PVDF) membranes (Immobilon P, Millipore, Billerica, MA). After blocking with 1% non-fat dry milk-Tris buffered saline and 0.1% Tween 20 (TBST), the membranes were probed with anti-HA (High affinity HA 3F10, 1:5,000 dilution, Roche), anti-α-CDK7 (sc-723, 1:5,000 dilution, Santa Cruz, CA), anti-EGFR (kindly provided by Dr. DiAugustine, RP) and anti-N-MYC1 (sc-791, 1:1,000 dilution, Santa Cruz, CA) and anti-c-FOS (sc-52, 1:1,000 dilution, Santa Cruz, CA) antibodies. Blots were then incubated with horseradish peroxidase (HRP)-conjugated rabbit anti-rat IgG (A5795, 1:5,000 dilution, Roche) or donkey anti-rabbit IgG (1:5,000 dilution, GE Healthcare Bioscience) secondary antibodies, respectively. Immunoreactive proteins were detected by enhanced chemiluminescence (P90720, Millipore). Signal intensities from the western blots were detected with X-ray film and quantified using NIH3T3 image software.

Strategy and Analysis of Gene Network Structures

Our experimental strategy is illustrated in Fig. (1) and consists of three steps: selection of datasets, visualization and analysis by mathematical modeling, and prediction of biological function through the analysis of transcripts. Genome-wide expression data can provide information linking diverse genes and may be useful as a classification tool to identify alterations in biological processes linked to disease. In contrast, carefully designed analyses of a limited gene group associated with a specific biological process can be used to quantify the dynamics of a gene regulatory network. The genes associated with cell cycle regulation are an obvious target for this type of analysis and are the focus of our current study.

Fig. (1) Strategy used to identify, analyze, and validate regulatory gene networks.

The first step in our approach was to select a data subset from a pool of genes associated with various aspects of cell cycle-related processes. The gene choices were based on the gene ontology of the mouse genome using GenMAPP, a computer application designed for the visualization of gene expression data by using maps representing biological pathways. This technique provided a qualitative tool for grouping genes (see Materials and methods). To gather gene expression data associated with the cell cycle, mouse embryo fibroblasts (MEFs) [4] were serum starved or exposed to hydroxyurea to synchronize and control their movement through the cell cycle. At various time points following the release from G0 and cell cycle re-entry, the mRNA expression levels for 6437 genes was measured using a microarray. For 10 cell-cycle related maps (Table 1), 145 genes were measured in the microarray assay. Of these 145 genes, 50 genes met the criteria of at least 2-fold higher or lower levels as shown in Table 1. Since the number of genes analyzed using TAO-gen had to be reduced due to computer processing limitations, 39 out of these 50 genes were finally selected for further analysis based on tissue-specific expression information and their biological significance from published articles after removing overlapping genes.

Two separate maps linking our selected 39 genes to a network were generated using the G0 course data subset (serum starvation) and the G1/S course data subset (hydroxyurea treatment). Although the expression of these genes is dynamic during the cell cycle, the networks were modeled by assuming equilibrium between the genes and by evaluating those using formal statistical methods that quantified any linkages and assessed their significance. Nodal genes (genes that appeared to be linked to a large number of other genes) were positively identified in the network. In the final verification step, the promoter regions of the genes targeted by each nodal gene were analyzed for common transcriptional factor binding sites. Finally, we discuss the roles of the central nodes and the dynamics of the quantified network in relation to the murine cell cycle.

Identification of a Gene Network Based on Expression Profiles

Representative maps using our 39 gene networks were developed separately for the G0 course (Fig. 2A) and the G1/S course (Fig. 3A). The number of linkages in these two networks is summarized in Table 3. Name abbreviations of the genes analyzed in this paper are shown according to their listing in GenMAPP. These networks were developed using Bayesian networks and a mathematical model allowing each of these mRNAs to connect to any other through direct or indirect transcriptional regulation leading to gene expression changes.

Fig. (2) Representative maps and expression graphs of the transcriptional regulatory networks for selected genes associated with cell-cycle control in MEF cells. Shown are (A) the network identified for the G0 course data and also the isolated linkages associated with nodal genes Cdkn2a (B) and Cdk7 (C). Bold lines indicate linkages from Cdkn2a or Cdk7 as a nodal gene. Red arrows indicate linkages associated with upregulation and blue arrows indicate linkages associated with downregulation for any two genes within the network.

Fig. (3) Networks identified for the G1/S course data (A) and the isolated linkages associated with nodal genes Ccna2 (B), Egfr (C), Fgf3 (D) and Trp53 (E). Red arrows indicate linkages associated with upregulation and blue arrows indicate linkages associated with downregulation for any two genes within the network. Bold lines indicate linkages with nodal genes.

The network from the G0 course data subset in which the cells had been serum starved indicates that the cyclin-dependent kinase inhibitor 2A (Cdkn2a) and Cdk7 are central nodes (Fig. 2B, C), whereas E2f1, known to regulate the G0-G1 transition, plays a lesser role. Although the Cdkn2a and Cdk7 gene products and related molecules have been suggested to functions in regulating G1 entry and progression from side supportive data [26, 27], it was not clear until our current findings whether these molecules functioned as central nodes in the gene networks. Cdkn2a and Cdk7 were not classified as a G0 cluster via k-means in the first report of these microarray data [4]. CDKs are known to be key components of the core cell cycle machinery and are inhibited by cyclin-dependent protein kinase inhibitors (CDKNs). CDK7 and CDKN2A are members of the CDK and CDKN families, respectively. Cdkn2a also encodes p16^INK4a, a protein that indirectly regulates the activities of both pRB and p53 through the inhibition of CDK4 and CDK6. The predictive pathway from Cdk7 suggests that CDK7 down-regulates Ccna2, Egfr, Il1a, Mybl2, Nmyc1 and Nras, and up-regulates Etv6, Figf, Pgf and Rbl1 (Fig. 2C). For the time course of the expression levels of Cdk7, Ccna2, Egfr, Mybl2, N-myc and Nras following the release from serum starvation, when the cells enter G1, the expression levels of Cdk7 are reduced, resulting in the elevated expression of Ccna2, Egfr, Mybl2, N-myc and Nras (data not shown).

In the network found for the G1/S start data subset in hydroxyurea treated MEFs, the structure was observed to be more complicated and have no obvious central nodes. In this network, the number of connections from Cdk7 and Cdkn2a to other genes was greatly decreased, whereas the connections from Ccna2, Egfr, Fgf3, Trp53, Nmyc1, Ptn, and Rbl2 were increased (Fig. 3A). These changes suggested that growth factors, such as Egfr, Fgf3, and Ptn, and proliferation regulators, such as Ccna2, Trp53, and Rbl2, have more prominent roles during S phase progression (Fig. 3B-E). From these data, it becomes obvious that the gene networks which regulate the progression of the cell cycle completely differ between the G0-G1 and G1-S transitions.

Verification of the Quantified Network

Further analyses were conducted to determine the statistical significance of the linkages between our identified genes. To find the most prominent linkages between genes of the network from Cdk7 and of the network from Trp53, the G0 course dataset and the G1/S course dataset obtained from MEFs treated with serum starvation or hydroxyurea were used, respectively. This analysis method can predict both the strength of the relationships between genes and the posterior distribution of parameters in the log-linear model [22, 23]. Of the 10 genes associated with Cdk7, 9 had some posterior densities that did not include 0, suggesting very significant associations (Table 4). Only Il1a included 0 in the posterior density, with 18% of the distribution above zero and 82% below. This finding suggested a statistically marginal down-regulation. Fig. (4) illustrated the distribution for a strong down-regulation (Cdk7 → Nras) and for a weak down-regulation (Cdk7 → Il1a). A negative association between Nras and Cdk7 has been reported previously [28], suggesting that the method we employed in our present analyses can extract negative relationships between two genes using simple microarray data.

Fig. (4) (A, B) Frequency histograms approximating the posterior distributions for linkages from Cdk7 to Nras (a statistically significant downregulation) and Cdk7 to Il1a (marginally significant downregulation). Histograms were derived by Bayesian analysis of the gene interaction network shown in Fig. (2C) using 70,000 out of 140,000 Markov-Chain Monte Carlo samples and prior distributions as shown in Table 4.

In the G1/S network, Trp53 suppressed the expression of Cdk7 and Rbl2, and stimulated that of Abl1, Il1a, Nmyc, Elk1, Ret and Thra (Fig. 3E). Trp53 has previously been shown to negatively regulate cyclinD/CDK4, cyclinD/ CDK6, cyclinB/cdk2, and cyclinA/cdk2 through the activation of p21 in normal cells. CyclinD/CDK4/6 on the other hand activates phosphorylated RB (pRb) which leads to the activation of E2F, which in turn negatively regulates p53 through p19^ARF activation and MDM2 suppression [1, 29]. The interaction between p53 and c-Abl is known to play a critical role in the cell growth and G1 arrest response to DNA damage under normal conditions [30]. It has been reported that CDK7 phosphorylates other CDKs, which is an essential step for their activation [31] and that a direct involvement of p53 in triggering growth arrest by its interaction with the CDK activating kinase complex [32]. These reports and our predictive network suggest therefore that CDK7 is essential for mitosis.

Detection of Gene Networks Using a Recombinant Mouse CDK7 Retrovirus System

Our cell cycle network data indicated that CDK7 activation negatively regulates the expression of Egfr and Nmyc1 in MEFs. To validate this observation, we introduced mouse CDK7 into these cells using a recombinant retrovirus system to evaluate negative regulation of CDK7 against EGFR and N-MYC. The titer of the retrovirus obtained from PT67 producer cells was 4.0 X 10⁹ virus copies/ml for the LXIN empty vector and 5.4 X 10⁹ virus copies/ml for the CDK7 recombinant retrovirus. The hemagglutinin (HA) protein tag was added to the carboxyl terminus of recombinant CDK7 so that we could distinguish the recombinant protein from its endogenous counterpart.

As shown in Fig. (5), western blot detection with a HA antibody revealed the expression of recombinant CDK7 protein in infected MEF cells. Increased levels of total CDK7 protein (endogenous plus recombinant CDK7) was also confirmed by immunoblotting with a CDK7 antibody (Fig. 5A). The EGFR, N-MYC1 and c-FOS protein levels detected by western blot were decreased in MEF cells infected with the CDK7-expressing retrovirus when compared with the control cells (Fig. 5A). c-FOS was used as control because there was no direct linkage between CDK7 and c-FOS (see Fig. 2A). The average levels of EGFR and N-MYC1 from three separate experiments are shown in Fig. (5B). Decreased EGFR and N-MYC1 but not c-FOS protein levels indicated that the exogenous introduction of CDK7 negatively influenced their expression. From these results, we concluded that one part of our newly detected cell cycle network had been validated.

Fig. (5) Experimental verification of detected gene network from Cdk7 to EGFR and N-myc1 using a recombinant retrovirus expression system. (A) The protein levels of exogenous Cdk7 (HA), total Cdk7 (Cdk7), EGFR, and N-myc1 were detected by western blotting. Representative blots obtained from two independent samples are shown in the figure. (B) EGFR and N-myc1 protein levels were quantitatively analyzed. Data are the average plus standard deviation of 6 western blots from two independent samples for each group. *, P < 0.05; **, P < 0.01.