Developmental Expression of Morphoregulatory Genes in the Mouse Embryo: An Analytical Approach Using a Novel Technology

BIOCHEMICAL AND MOLECULAR MEDICINE 60, 8 1 - 9 1 (1997) ARTICLE NO. MM972576

Developmental Expression of Morphoregulatory Genes in the Mouse Embryo: An Analytical Approach Using a Novel Technology J. C. CRAIO,* J. H. EBERWINE,t J. A. CALVIN,~: B. WLODARCZYK,* G. D. BENNETT,* AND R. H. FINNELL* *Department of Veterinary Anatomy and Public Health and $.Department of Statistics, Texas A&M University, College Station, Texas 77843; and tDepartment of Pharmacology, University of Pennsylvania School of Medicine, Philadelphia, Pennsylvania 19107 Received J a n u a r y 9, 1997

critical regulatory events directing early mammalian development. This has largely been facilitated by the advent of new experimental methodologies such as in situ hybridization and immunohistochemical technologies, which have enabled investigators to examine the expression of selected genes and their products in embryonic tissues. While such progress has elucidated the temporal and spatial pattern of expression for a rapidly expanding number of individual genes, it has not addressed how the interactions among genes provide regulatory control over the events involved in embryogenesis. This is because morphogenesis is not regulated by individual genes, working in isolation, but is instead the result of m a n y genes functioning in concert. For example, an event such as neural tube closure (NTC), which involves cellular proliferation, migration, and adhesion, along with the acquisition of new cellular phenotypes (1-3), is far more likely to be controlled by the combined effects of several genes, rather t h a n isolated single major gene effects. Therefore, it is far more informative to study the developmental regulation of gene groups, in order to gain a better understanding of their coordinated interactions during morphogenesis. Using in situ transcription together with amplified antisense RNA (aRNA) technologies, it is possible to monitor changes in gene expression within anatomically restricted regions of the developing mouse embryo. This approach has provided an enormous amount of novel data on the coordinate expression of a large number of developmentally regulated genes during murine NTC. In order to rationally approach the interpretation of this new information, we have attempted to develop a statistical framework for the

The molecular techniques of in situ transcription and antisense RNA amplification (IST/aRNA) have allowed for the monitoring of coordinate chnnges in the expression of multiple genes simultaneously. However, the analysis of their concurrent behavior during murine embryogenesis has been problematic. Studies involving the investigation of temporal and spatial gene expression during embryogenesis have focused solely on the nnAlysis of isolated, single gene events. Such an approach has failed to provide an integrative picture of genetic control over the varied and complicated cellular processes governing embryogenesis. In order to interpret the enormous amount of gene expression data generated by these procedures, we have attempted to develop an analytical framework by employing the statistical concepts of principal components analysis (PCA). For the current study, we performed IST/aRNA on neural tubes dissected from the highly inbred LM/Bc murine strain collected during four gestational time periods. A subset of these genes, representing a partial signAllng pathway in the developing neuroepithelium, was then subjected to PCA. Here, we report that PCA highlighted the trnn.~criptional interplay among the genes p53, wee-l, Tgffl-2, and bcl-2 such that the combined reciprocal regulation of their gene products is suggestive of a predominant proliferative state for the developing neuroepithelium. The application of PCA to the gene expression data has elucidated previously unknown interrelationships among cell cycle genes, growth, and transcription factors on a transcriptional level during critical stages of neurulation. The information gleaned from this analysis, while not definitive, suggests distinct hypotheses to guide future research. ~ 1997AcademicPress Recent advances by developmental geneticists have significantly altered perceptions concerning 81

1077-3150/97 $25.00 Copyright © 1997 by Academic Press All rights of reproduction in any form reserved.

• 82

CRAIG

analysis of gene expression data using principal components analysis (PCA). PCA is a statistical technique designed both to condense information and to create a new, more interpretable structure in the data. Thus, it is useful in finding hidden relationships in data and reducing the amount of information to a manageable level. We utilized PCA to explore the gene expression data in an effort to elucidate biologically meaningful combinations of gene expression behavior within each of the gestational time periods under study. In addition to a better understanding of each time period, our goal was to compare the results from Gestational Day (GD) 9:0, representing neural tube closure, to the other three post-NTC time periods in an effort to find combinations of the genes which were capable of distinguishing NTC from post-NTC embryos. This multivariate statistical approach was utilized to examine the coordinate changes in gene expression in 7 (out of an original panel of 40) genes that represent a partial signaling pathway in the developing neuroepithelium (Table 1). This gene subset was chosen because it may reflect regulation of critical cellular events during neural tube closure. The expression patterns were examined in murine embryos collected during the period of GD 9:0 (Day 9 plus 0 h postconception) and compared to those in embryos collected at GDs 9:12, 10:0, and 10:12. The hybridization intensity of aRNA probes generated from the neural tube tissue to cDNA clones of selected candidate genes immobilized on reverse Northern slot blots allowed the relative abundance of the mRNAs of interest to be quantified and statistically analyzed. Fluctuations in embryonic mRNA population levels during this period were monitored, and relationships among the selected genes with respect to developmental time during early neural tube development were assessed. Ultimately, these comparisons will allow for the quantification of covariant molecular aspects of normal neural tube development, for the purpose of providing vital information on the molecular pathogenesis of neural tube defects (NTDs). This approach was previously utilized in murine embryos to describe coordinate changes in gene expression in response to teratogenic insult (4). However, the present study provides the first composite examination of coordinate changes in mRNA abundance over time. MATERIALS AND METHODS

Embryo collection and morphological staging. The highly inbred LM/Bc/Fnn mice were maintained

ET

AL.

on a 12-h light/dark cycle in the Laboratory Animal Research and Resource Center at Texas A&M University. Pathogen-free, virgin females, 50 to 70 days of age, were mated overnight with experienced males, and the dams were checked for the presence of a vaginal plug the following morning. The start of gestation was set at 10 PM of the previous evening, the midpoint of the dark cycle (5). Embryos were collected at selected time points between Gestational Days 9:0 and 10:12. At the assigned hour, pregnant dams were killed by cervical dislocation, the abdomen was opened, and the uterine contents were removed. The locations of all viable embryos and resorption sites were recorded. Using watchmaker's forceps, the embryos were dissected free of the decidual capsule and its chorion and amnion while in cold PBS, under a Wild M8 dissecting microscope (Heerbrugg, Switzerland). The embryos were grossly examined morphologically and were classified as to the basis of their stage of NTC using previously described standardized staging criteria (6). These morphological data have been reported elsewhere (7).

Removal of neural tube from NTC-stage embryos. With the aid of watchmaker's forceps, the neural tube proper was dissected away from supporting paraxial mesodermal tissue under the dissecting microscope (8, 9). Once dissected, it was examined to ensure that only the intact neural tube tissue had been collected and placed in a hybridization buffer containing 5 mM DTT, 100 U RNasin (Promega, Madison, WI), and 0.1% digitonin. Following a brief pulse with a sonic dismembranator, an additional 50 U of RNasin was added to the hybridization buffer and the tissue was frozen at -80°C until further processed for in situ transcription and amplified antisense RNA studies.

In situ transcription and aRNA amplification procedures. These procedures are described elsewhere (10). Briefly, an RNA/DNA hybrid molecule was made from the population of mRNAs present in the isolated neural tube samples with the addition of avian myeloblast reverse transcriptase (Seikagaku America Inc., Bethesda, MD), and an oligo(dT)-T7 oligonucleotide primer that binds to the poly(A) ÷ tail of the mRNAs. The single-stranded DNA produced in this reaction was isolated and made into doublestranded cDNA by hairpin loop formation. Radiolabeled, amplified, antisense RNA was generated in an enzymatic reaction driven by T7 RNA polymerase (Epicentre Technologies, Madison, WI) in the presence of [32p]CTP. The aRNA, representing the entire population of mRNAs from the neural tube tissue,

83

GENE EXPRESSION DATA ANALYSIS TABLE 1 Candidate Genes Examined for Evidence of Developmental Regulation i n t h e M o u s e E m b r y o b e t w e e n G e s t a t i o n a l D a y s 9:0 a n d 10:12 Gene symbol

Gene name

Classification

c-los wee-1 p53 bcl-2 Tgf/3-1 Tgf/~-2 TgfB-3

Cellular fos "wee" Protein 53 B-cell leukemia type 2 Transforming growth factor beta 1 Transforming growth factor beta 2 Transforming growth factor-beta 3

Transcription factor Cell cycle gene Oncogene/cell cycle Oncogene/cell cycle Growth factor Growth factor Growth factor

was allowed to hybridize to the "reverse" Northern blots in order to quantitatively determine the patterns of gene expression.

"Reverse" Northern blots and genetic expression profiles. Equimolar concentrations of 40 cDNA clones for the candidate genes were immobilized on nylon membranes (Zetaprobe, Bio-Rad, Richmond, CA) using a Bio-Rad slotting apparatus. Seven of these clones (Table 1) were selected on the basis of their known important functions in early murine morphogenesis. Each blot was prehybridized for 30 min in buffer (7% SDS (w/v), 0.12 M Na2HPO4 (pH 7.2), 0.25 M NaC1, and 50% formamide) at 42°C. The heat-denatured aRNA was applied to the hybridization buffer and allowed to hybridize for 24 h. Subsequently, the slot blots were washed down to 0.1× SSC containing 0.1% SDS at 42°C, dried, wrapped in plastic wrap, and placed in an Ambis 101 twodimensional radioanalytics imaging detector (Scanalytics, Billerica, MA) which directly measured the radioactivity (cpm) of each slot on the reverse Northern blots. The individual signals were normalized to cyclophilin gene expression. The selection of cyclophilin as the normalizing cDNA simply enabled us to make comparisons between different blots, as the individual hybridization intensities of each cDNA on each blot can be expressed as a ratio of its expression to that of cyclophilin. The hybridization pattern obtained from each embryonic neural tube-derived aRNA probe was referred to as its "genetic expression profile" (Fig. 1). For each of the time points examined, no fewer than 14 expression profiles were generated, with at least four different females providing at least three embryos each.

Basis for selection of candidate genes. The seven cDNA clones selected for this particular analysis encode three cell cycle genes (wee-l, p53, and bcl-2), three growth factor genes (transforming growth factor beta 1, 2, and 3 (Tgfi3-1, Tgf~-2, and Tgf~-3)) and

one transcription factor gene (cfos) (Table 1). The basis for selection of these cDNA clones was their documented interactive effects on epithelial cell differentiation and proliferation through the cell cycle. The selected cell cycle genes encode checkpoint proteins that operate at the G2 (wee-l) or the G1 (p53 and bcl-2) cell cycle phases and are involved in various aspects of assessment and regulation of DNA replication and repair (11). wee-1 encodes a protein kinase that negatively regulates entry into mitosis by inactivation of the maturation-promoting factor (12). Overexpression of wee-1 inhibits the cellular transition from G2 to mitosis in the cell cycle, ultimately leading to apoptotic cell death (12,13). At the G1 cell cycle phase, it has been reported that the protein product of bcl-2 mediates the growth-inhibiting and apoptotic effects of p53 protein product (14). Furthermore, the expression of bcl-2 is induced by the Tgf~s in the CNS (15). All three of these cell cycle genes are reportedly important in murine development and are expressed in a number of different tissues, especially in epithelium-like cells, p53 is expressed in the developing murine neural tube and has important implications in orchestrating neural tube closure, via cell cycle regulation (16), while bcl-2 is normally highly expressed in fetal CNS tissues, as well as in rapidly proliferating epithelial cells (17,18). Tgfl3 is a member of a superfamily of proteins that has important implications in growth and development via a variety of functions, including the modulation of cell cycle progression (19). The Tgf~ isoforms, Tgf~-l, Tgf~-2, and Tgf/~-3, are known to possess similar structural as well as biological features (20). They have all been implicated as potent negative or positive growth regulators, depending on the cell type. In adult and embryonic epithelial cells, these proteins act as growth inhibitors, although the exact mechanism is not well understood (20). Several

84

CRAIG ET AT_,.

.....J

I Hoxc8

Pax-3

I cfos

Wnt-1 -

--

cqm

Q

m N-Cam

CREB

m CRBP-1

n TGFD3

Ern~l

i

D

TGFa

m

NFL

I~1

Igfll

BDNF

=--~m,m-.-~Ngf-R

Emx-2

GAPD

g CRBP-2 m W

pro3

TGFbI

m

TGFb2

Q

m

wee- 1

mtrk

Gc-R

m

NT3

NGF

bFGF

g

Shh

cph

b-myo

bcl-2

ng

FBP

FIG. 1. A representative "reverse Northern blot" illustrating the genetic expression profile for 40 candidate cDNAs, plus the cyclophilin standard cph from an individual LM/Bc/Fnn embryo collected at Gestational Day 9:0.

studies have demonstrated a negative role for Tgf~ signaling in cell cycle regulation, specifically at the late G1 phase, through inhibition ofcyclin-dependent protein kinase 4 (CDK4) and cyclin A expression in epithelial cell lines (21, 22). In addition, it is known that Tgf~-I suppresses the expression of several early gene products in the epithelium, including the proliferin gene cfos, most likely through its inhibitory effects on CDK2. The net result is to block the progression of cells to the S phase of the cell cycle. Tgf~-l, Tgf~-2, and Tgf~-3-induced down-regulation of CDK4 also inhibits CDK2 activation, thus providing a link between the two latter isoforms to the prevention of cell cycle progression via inhibition of cfos (21). These results underscore a potential functional link between the Tgf~s and other CDK inhibitors, such as p21 and p53, the latter of which has been postulated to positively regulate TgfB gene expression in the epithelium (23,24). Taken together, the interassociations of these seven genes may constitute a partial signaling pathway in the developing neuroepithelium that regulates the transition between cellular proliferation and differentiation. Such a pathway may have important implications for neuroepithelial cell function during neural plate folding and neural fold fusion.

Statistical analyses. The data for the statistical analysis were generated from LM/Bc/Fnn embryos collected at GDs 9:0, 9:12, 10:0, and 10:12. Based on results described in previous work (7), logs of the individual gene expression measurements were used in all statistical procedures described below. That is,

ln(p53) was analyzed, rather than p53 directly. This was motivated by the nonconstant variance observed in the original measurements. All statistical analysis of these data was performed using SAS (Statistical Analysis System) and completed in three stages. The first of these stages involved univariate analysis of individual candidate genes.

Univariate analysis. In this first stage, each gene was studied to determine if its mean expression differed among the four time periods. Statistical significance was set a t P < 0.05 (25). Due to the embryo sampling scheme, which involved the collection of up to four embryos from each dam, the model for the analysis was more complicated than the onefactor ANOVA model, which would have been appropriate had the embryos all come from different dams. To adjust for any potential correlation among embryos coming from the same dam, a nested two-factor mixed model was used (see 26 for details). Estimates of the mean expressions for each GD were found using least squares means (LSMEANs). The LSMEANs were estimates of the mean response which were adjusted for the correlation in the data and any inequality in the sample sizes. For a description of the least squares estimation formula see Neter et al. (27), and for a s u m m a r y of the LSMEANs option see SAS (28). Mean separation comparisons based on the four LSMEANs was performed using Bonferroni adjusted significance levels. The SAS procedure MIXED was used in this stage. While the above analysis was described in terms

GENE EXPRESSION DATA ANALYSIS

of individual genes, this approach could also be performed on any predefined combination of the gene expression values. Based on previously published results, and an inspection of the individual gene expression results, two simple ratios of gene expression (p53:bcl-2 and Tgf~-2:bcl-2) were also analyzed using the techniques described above. Principal components analysis. PCA is a statistical technique which looks at a group of data with the goal of creating new variables which are (i) recombinations of the original data and (ii) capable of providing the same information as the original data, but which are also more insightful when trying to assess interdependence among the original measurements. This technique starts with a measure of the interdependence of the original data, the covariance matrix, and creates a set of new variables with a sample covariance matrix that indicates their independence from each other. This is accomplished by finding linear or additive combinations of the original measurements based on the eigenvectors of the covariance matrix. The goal is to find such a transformation which will yield biological, as well as statistical, sense. For a discussion of PCA see Rao (29), Rohlf and Bookstein (30), or Johnson and Wichern (31). As a technical footnote, PCA is usually performed using observations which are independent. This is not our situation, but due to the extensive controlled breeding program for this strain of mouse, it is reasonable, in this exploratory setting, to set aside the embryo interdependence when looking for target gene relationships. A separate PCA was performed for each GD. Since we started with expression data on seven genes, the PCA results consisted of seven new variables (principal components, PCs), which are linear combinations of the original measurements and, for convenience, sorted based on their variability. Thus, the first new variable, called PC1, is the combination of the original measurements which is the most variably expressed, whereas PC7 is the combination which is the least variably expressed. Thus, PC6 and PC7 are much more consistent in their respective values, among all embryos within the selected time period. As such, they make good candidates for describing interdependencies among the original measurements which are consistently expressed for all embryos within a GD. For each GD, PC1 through PC7 was computed, and the combination of genes going into each PC was assessed for biological viability. The SAS procedure PRINCOMP was used for these calculations. Follow-up analysis. Once a specific PC was deemed interesting, the interdependence was sim-

85

plified by eliminating any genes which had minimal impact on the PC (eigenvector values < 10.21 ). This new, reduced combination of the original measurements was then constructed as a multigene ratio for all embryos, regardless of GD. To graphically assess the change in the PC from one GD to the next, a plot of the log of the numerator versus the log of the denominator was produced. If the PC distinguishes among the time periods, then the grouping for one GD will look quite different from the rest, in both the location of the points and the pattern of the points within a GD. Simple descriptive statistics were calculated for this ratio and a test of equal variances was performed using Hartley's F ~ test (32) for each GD. To ensure that the exploratory process did not misstate the temporal relationship due to the lack of adjustment for correlation due to our sampling, the new PC was also assessed using the same univariate procedures described for the first stage. RESULTS

Univariate analysis. Analysis of variance, using the MIXED procedure in SAS (Statistical Analysis System), was performed to generate and compare the singular mean expression of the seven genes of interest, as well as selected pairwise, simple ratio relationships across gestational time. An exhaustive univariate comparison using these data was performed previously (7). In the current study, the mean expression levels of all three cell cycle checkpoint genes, as well as cfos and Tgf/~-3, were significantly increased from GD 9:0 to 9:12 (Fig. 2; P < 0.05). Of these, only the mean expression level of wee-1 was significantly altered over all four gestational time periods, increasing to 15 times its original expression level by GD 10:12 (Fig. 2; P < 0.05). The mean expression level of Tgf/~-2 decreased significantly at GD 10:12, compared to that in GD 9:0 embryos (Fig. 2; P < 0.05). Two of the cell cycle genes continued to change over time, with significant decreases seen in the mean expression levels of p53 and bcl-2 from GD 9:12, compared to GDs 10:0 and 10:12 (Fig. 2; P < 0.05). Similar decreases were noted in the mean expression levels of Tgf/~-2 and Tgf;3-3 (Fig. 2; P < 0.05). Finally, the mean level of cfos expression significantly decreased, while that of Tgf/~-I significantly increased from GD 9:12 to 10:0 (Fig. 2; P < 0.05), though these levels remained somewhat negligible. The temporal relationships of two simple ratios, p53:bcl-2 and Tgfl~-2:bcl-2, were compared using the same univariate approach described for the above

86

CRAIG ET AL. TABLE 2 Least S q u a r e Mean Response of Select Univariate Ratio Comparisons

A

i

1.2

JL

c-fos

•

+x~-~

1.0

-

wee-1

0.8

•

p53

2

o.s

~"

0.4

~

0.2

.I =

,

t

9:0

9:12

10:0

10:12

p53:bcl-2 Tgfl3-2:bcl-2

-0.081 -0.069*'**

0.070* -0.240

-0.086 -0.219

-0.236 -0.256

* P < 0.05 from GD 10:12. * * P < 0.05 from GD 9:12.

t

-0.2

!

i

i

GD 9:0

GD 9:12

GD 10:0

i

GD

10:12

b. "

0.4" I C w Q

:E

Gene ratios m~

o.o

0.6

Gestational day

0.2

°

0.0

"

~

*

--

Tgfb-1

•

Tgtb-2

1

Tgfb-3

Principal components analysis and multivariate ratio comparisons. Due to the fact that it is infeasi-

$ ,=J

-0.2

-0.4 GD 9:0

sons for this ratio were found to be nonsignificant (P > 0.05). In contrast, the temporal comparison of Tgfl3-2 to bcl-2 across time revealed significant differences between several of the time periods, specifically between GDs 9:0 and 9:12, and GDs 9:0 and 10:12 (Table 2; P < 0.05). While these individual gene expression averages and prespecified simple ratios provide some fundamental information, it is difficult to divine any more complicated relationships from these results. This is due to the fact that these results are based on averages of groups of embryos and not on interrelationships within an embryo. Thus, it is possible that important gene interactions may go unnoticed without a mechanism for looking at within-embryo gene relationships.

GD 9:12

GD 10:0

GD 10:12

FIG. 2. The least square m e a n response of the seven candidate genes at four gestational time points. Expression values are reported

as the mean hybridization intensities relative to cyclophilin and were log transformeck (a) Expression values for cell cycle and transcription factor genes. Co) Expression values for growth factor genes. *Indicative of significant differences from gestational day (GD) 9:0 (iv < 0.05). tIndJcative of significant differences from GD 9:12 (P < 0.05).

analysis. Table 2 indicates a significant difference (P < 0.05) for the p53:bcl-2 ratio between GDs 9:12 and 10:12 only, while the other temporal compari-

ble to envision all possible relationships among the genes of interest by looking at each gene individually, PCA was used to provide insight into more complicated aspects of gene expression. From the seven original measurements, seven PCs were computed for each of the four GDs using the PRINCOMP procedure in SAS. For the sake of brevity, only those coefficients for the GD 9:0 PCs, along with their respective percentage of original variability, are given in Table 3. Note that PC1 is extremely variable when compared to PC6 or PC7. Of these new variables, PC6 was deemed to be the most biologically interesting. It contains only 0.7% of the original variability and thus is expressed at a very consistent level within all embryos at GD 9:0. Inspecting the eigenvector weights associated with each gene, we decided that a simplification involving onlyp53, Tgf~-2, wee1, and bcl-2 adequately represented this PC. This subset of genes is used in any further discussion of PC6. Thus, for the plot of In(numerator) (natural log numerator) versus In(denominator), p53, Tgf~2, and wee-1 are in the numerator and bcl-2 is in the

87

GENE EXPRESSION DATA ANALYSIS TABLE 3 PC Weights for the PCA Performed Genes

on GD 9:0 Embryos

PC1

PC2

PC3

PC4

PC5

PC6

PC7

bcl-2 p53 weeol cfos Tgf~-I TgfB-2 Tgf~-3

0.394 0.411 -0.367 0.267 -0.375 0.403 0.409

0.218 0.074 0.402 0.779 0.407 0.112 -0.021

0.560 -0.037 0.630 -0.477 -0.420 0.135 0.204

0.062 0.029 -0.213 -0.285 0.418 0.718 -0.449

0.030 0.449 -0.170 -0.190 0.696 -0.284 0.408

-0.635 0.450 0.474 -0.002 -0.116 0.320 0.124

-0.276 -0.617 -0.087 -0.038 0.139 0.329 0.638

Variance:

80.4%

12.8%

2.75%

1.59%

1.2%

0.7%

0.52%

Note. The PCA was based on the expression values of the three cell cycle checkpoint, one transcription factor, and three growth factor genes. The proportion of variance explained by each PC is given. Boldface values on the selected PC6 signify heavily weighted genes.

denominator (Fig. 3). The small variability of this combination can be interpreted by observing that as the numerator genes increase so too does the denominator, thereby keeping the ratio constant. In Fig. 3, note that the embryos for GD 9:0 are much more closely clustered around their line, indicating the consistency of their expression of this gene combination. Also note that the cluster for GD 9:0 is noticeably different than the clusters for the other time 3

• A 0 []

GD GD GD GD

9:0 9:12 10:0 10:12

periods, indicating a temporal change in the gene relationships. The mean and variance for this gene combination were computed for each GD and are displayed in Fig. 4. The means and variances vary across gestational time. Using Hartley's F~,~ test for constant variance, we can reject constant variance (P < 0.05). Thus, this gene combination is much less predictable at later time periods. A follow-up test for equality of means also rejects at P < 0.05, with the mean for the GD 9:0 embryos being significantly dis-

A 1.25 -

?

O

*

2

[]

OA il

A

t t

1-

t

T

±

0.75 -

|

0.5 o

E

0.25

_o m

-1

i tS3 O

~. 9

I (~

I 14")

~ o

I tS~

(5

I ~

I

~-. o

bcl-2 FIG. 3. Principal components plot displaying the LM/Bc/Fnn embryos from GDs 9:0, 9:12, 10:0, and 10:12 in the coordinate space defined by PC6. The analysis utilized the log-transformed expression data for the seven candidate genes. The line through the GD 9:0 group represents the linear relationship of the embryonic response to the defined gene ratio.

[

0

~

"

.

,

Gestatlonal day FIG. 4. The LSMEAN expression level values of the log-transformed, multivariate ratio defined by PC6, across gestational time. *Indicates t h a t the means are significantly different from t h a t at GD 9:0 (P < 0.05). tIndicates t h a t the variances are significantly different from t h a t at GD 9:0.

88.

CRAIG ET AL.

tinct from the means at the latter three time periods (Fig. 4). DISCUSSION The current study utilized an exploratory statistical procedure, PCA, in order to examine simultaneous and temporal changes among the expression patterns of seven genes (Table 1) in neuroepithelial tissue dissected from LM/Bc/Fnn murine embryos. This multivariate approach provided a more unified picture of gene expression during NTC than has previously been possible, due to the limitations imposed by univariate or statistical analyses on individual genes or simple ratios. The methodology revealed developmental trends and counteractive associations among these genes at GD 9:0 that were not seen at the later time periods. Though not shown, the PCAs performed separately on the later three time periods yielded PCs that were inconsistent with those generated for the GD 9:0 group. Therefore, the gene associations revealed by the GD 9:0 PCA may require consistent maintenance during this period of active NTC to ensure normal morphogenesis. The univariate analyses were useful in supplying fundamental data on individual gene changes over time and raising intriguing questions regarding potential interactive events. For example, it was interesting to note that the mean expression levels of most of the genes were temporally altered (Fig. 2). In addition, the majority (six of seven genes) were significantly altered (P < 0.05) with respect to GD 9:0, with the exception of Tgffl-2, whose level of expression remained constant until GD 10:0. It was further noted that the mean expression levels of bcl2 and p53 were initially comparable to those of Tgf~2 and Tgf~-3 at GD 9:0 (approximately 0.02). These four genes resumed their parallel expression patterns after GD 9:12, and after the initial lag in Tgf;32 expression. These expression patterns, though uncorrelated in the univariate analysis, are consistent with the report that the protein products of the TgfBs induce bcl-2 gene expression in epithelial cells (15), as well as the report that the protein product of bcl-2 mediates the apoptotic and growth arrest effects of the p53 protein in the neuroepithelium (33). Furthermore, the univariate analysis demonstrated that Tgffl-1 expression is almost nonexistent, while that ofcfos remains relatively high across time. Since the protein product of Tgf/~-I inhibits proliferation in epithelial cells by indirectly inhibiting cfos transcription (34,35), one can speculate that these univariate patterns m a y be indicative of a

proliferative state among these neuroepithelial cells. While the above observations pose intriguing questions, they reveal nothing about the coordinate regulation of these genes. In an effort to understand the gene interactions associated with such a complex process as NTC on the molecular level, we examined the expression data taken as a whole using PCA and seven genes that form an incomplete cell cycling pathway potentially involved in regulating the transition between neuroepithelial differentiation and proliferation to facilitate or aid NTC. This analysis yielded an interesting gene combination that provided a biologically plausible interpretation of neuroepithelial cell activity at GD 9:0. H a d we relied solely on the univariate mean and ratio analyses (Fig. 2; Table 2), it is quite likely that this combination would have been overlooked. Furthermore, PCA revealed some intriguing gene interassociations that broadened our view of the process of neural tube gene regulation and generated novel hypotheses concerning genetic control over the cell cycle in these tissues. Quite often, the goal of a PCA is to examine interrelationships based on the majority of the variability and neglect those PCs that contribute very little variation. However, we were interested in examining uniform gene expression patterning among embryos collected during GD 9:0, a gestational time period representing active NTC. Accordingly, we sought to identify consistent patterns among these embryos that would reveal stable genetic relationships, including those that were potentially sensitive to perturbations by a teratogenic insult. This would provide us with novel insight into the pathogenesis of teratogen-induced NTDs. Therefore, following the PCA, it was determined that the gene combination delineating PC6 represented a negligible percentage of the variability associated with the original data set, and thus a consistent GD 9:0 embryonic response. In addition, PC6 provided the most biologically intriguing scenario in that the combined expression of wee-l, p53, Tgffl-2, and bcl-2 was comprehensibly counterbalanced to describe this period of NTC. PCA also allowed us to observed more complicated ratios which potentially have significantly greater resolution, in terms of describing the behavior of genes during early embryogenesis. As such, the ratio relationship of these four genes can be comprehended to explain apparent contrasts in their regulatory behavior. Specifically, such an interpretation involves contrasting the behavior of the genes in the numerator with that of those genes in the denominator. For example, all of the GD 9:0 embryos responded simi-

GENE EXPRESSION DATA ANALYSIS larly to the ratio of wee-l, p53, Tgf/~-2 to bcl-2, such that if the numerator genes either increased or decreased their combined expression levels, the expression levels of the denominator genes would be altered in parallel, for each embryo. The lack of significant differences among any of the later three time periods for this ratio indicated that these postNTC embryos expressed this gene combination in a similar manner (Fig. 4). That is, the numerator and denominator genes were not expressed in parallel to one another from one embryo to the next at GDs 9:12, 10:0, and 10:12 and that this correlative ratio association did not distinguish these later time periods from each other. Thus, the results of the analysis revealed a dramatic distinction between the GD 9:0 group and the later three time periods, in terms of the balance of this ratio relationship. More importantly, however, it indicated that this ratio relationship m a y not be as vital at the later time periods as it is at GD 9:0. Taken together, these analyses illustrated that gestational age (time) was a significant factor controlling the relative mRNA abundance of this group of genes. During NTC the neuroepithelium is going through a myriad of changes on the molecular, cellular, and tissue levels. There are likely many varied cellular events, such as differentiation, proliferation, adhesion, and migration, that are unified on a molecular level to accomplish successful neural fold elevation and fusion. The cell cycle is doubtless involved in the orchestration of many of these events. For example, it is believed that cells comprising the "hinge points" of the folding neural plate experience changes in their cell cycle length in order to alter (or juxtapose) the position of the cell's nucleus. The new nuclear positioning allows for a favorable cellular arrangement, such that the neural plate can effectively bend and fold to accommodate neural fold fusion (3). In addition, the cell cycle regulates transitions between differentiation and proliferation that occur at various time periods of embryonic development. While both growth and differentiation are likely occurring concurrently during NTC, it is probable that the former is the favored state during the most active period of NTC. This is because the number of cells required to facilitate elevation, folding, and fusion is doubtless higher than that required for later differentiation events that occur after closure has been accomplished. The four genes of interest, wee-l, p53, Tgf/~-2, and bcl-2, are known to operate at various phases of the cell cycle in epithelial tissues and may thereby influence the development of the murine neural tube

89

via cell cycle regulation. The mathematical relationship of these four genes, such that wee-l, p53, and Tgf;3-2 are in the numerator and bcl-2 is in the denominator, can be aligned with their biological relationship to provide meaningful insight into cell cycle regulation in the neuroepithelium. For example, the contrasting regulatory effects of the p53 and bcl-2 gene products are well documented, p53 is a multifunctional protein that has been implicated in modulating normal development, in part by regulating gene transcription and safeguarding cell cycle checkpoints (36). Its main action is to induce growth arrest at the G1 phase of the cell cycle, in order to assess DNA integrity and allow for DNA repair (36,37,38). It has been suggested that the protein product of p53 may have a role in regulating cell cycle arrest and mediating a cellular transition from proliferation to differentiation during NTC (16). The bcl-2 gene product is a potent inhibitor of p53-induced apoptosis and cell cycle arrest that is believed to contribute to the mediation of cell proliferation and differentiation (17,39-41). The growth arrest and apoptotic functions regulated by the p53 and bcl-2 gene products are thought to be activated when the equilibrium between the expression of these genes favors p53 (14,42). Though no gene to gene statistical comparison was performed within a given time period, the univariate analysis indicated a higher mean level of bcl-2 relative to that of p53 at GD 9:0. In addition, the multivariate analysis revealed a counteractive regulatory relationship between these two genes. Taken together, these data argue against the induction of neuroepithelial growth arrest and/ or apoptosis. This observation is further supported by the fact that Tgf;3-2, the gene product of which induces bcl-2 expression (15), is initially expressed at approximately the same level as that of bcl-2, as shown in the univariate analysis (Fig. 2), and is counterbalanced with bcl-2 in the PCA analysis. Furthermore, Tgf/3-2 is reportedly expressed in many embryonic mouse tissues, including epithelial cells and in the prosencephalic region of the developing neural tube (20). The univariate analysis (Fig. 2) also showed a low level of wee-1 expression at GD 9:0. However, the PCA indicates this gene's heavy influence on the patterning of the GD 9:0 embryos, as well as its contrasting behavior to bcl-2. It has been shown that up-regulation of wee-1 results in the induction of apoptosis in certain cells lines (12). Taken together we can speculate that the low level of wee-1 reflects the minimal operation of its protein product at G2 to allow swift progression through the cell cycle, thus complementing the

90.

CRAIG ET AL.

/o, V It ~'~ a-~ ~

Go

.,, /2,. ~'~ checkpointI t (DNAdamage)Tgffl

FIG. 5. Diagram of the primary cell cycle phases and locations of the checkpoints at which the candidate genes are operational. Upward arrows indicate regulatory influence. Encircling arrows indicate direction of neuroepithelial cell progression.

growth-inducing actions of the bcl-2 gene product at G1. Meanwhile, the protein product of bcl-2 is operating at G1 to override p53 growth suppression and allow cycling through this phase. Compositely, the counteractive effects of wee-l, p53, and Tgf~-2 to bcl2 at the gene expression level suggest a cooperative effort to regulate rapid cycling of the neuroepithelial cells, to facilitate the proliferative status required during active NTC (Fig. 5). In closing, this preliminary statistical investigation has revealed new and potentially significant relationships among a subset of seven genes known to govern certain aspects of cell cycle progression. Specifically, it can be inferred that maintenance of the transcriptional counter activities of wee-l, p53, Tgf~-2, and bcl-2 at GD 9:0 may be important to their downstream transcriptional targets or counterparts in the preservation of cell cycle regulatory balance, growth, and differentiation. Although these results are specific to this particular data set, we can tentatively speculate that the combined expression of wee-l, p53, Tgf~-2, and bcl-2 may accurately characterize the GD 9:0 embryos. While these observations do not constitute definitive information on the overall behavior of the neuroepithelium or the morphogenetic processes of neural tube closure, they do suggest distinct hypotheses concerning temporal gene expression patterning in the developing murine neural tube. Investigations involving the assessment of pertinent gene ratio relationships, as well as procurement of additional cDNA clones for analysis of these and other gene interactions, are currently underway in our laboratory. ACKNOWLEDGMENTS This work was supported in part by Grants DE 11303 and ES 07165 from the National Institutes of Health. The authors express

their appreciation to Ms. M. Wheeler and Ms. L. Cook for the well-being of the mice, Ms. M. Davda for her technical expertise, and Ms. A. Miller for her efficient editorial assistance. We thank the following individuals for providing cDNA clones used in these experiments: M. Pisano, G. Wilson, S. Korsmeyer, F. Ruddle, D. Epstein, R. Harland, T. Luikin, A. Brown, S. Miyatani, M. Wrzos, A. McMahon, and J. Eberwine.

REFERENCES 1. Langman J, Welch GW. Excess vitamin A and development of the cerebral cortex. J Comp Neurol 131:15-26, 1967. 2. Bannigan J, Cottell D. Ethanol teratogenicity in mice: An electron microscopic study. Teratology 30:281-290, 1984. 3. SchoenwolfGC, Smith JL. Mechanisms ofneurulation: Traditional viewpoint and recent advances. Development 109:243270, 1990. 4. Finnell RH, Wlodarczyk BC, Craig JC, Piedrahita JA, Bennett GD. Strain dependent alterations in the expression of folate pathway genes following teratogenic exposure to valproic acid in a mouse model. Am J Med Genet, in press, 1997. 5. Snell GD, Fekete E, Hummel KP. The relation of mating ovulation and the estrus smear in the house mouse to the time of day. Anat Rec 76:30-54, 1948. 6. Cole WA, Trasler DG. Gene-teratogen interaction in insulin induced mouse exencephaly. Teratology 22:125-139, 1980. 7. Wlodarczyk B, Bennett GD, Calvin JH, Craig JC, Finnell RH. Arsenic-induced alterations in embryonic transcription factor gene expression: Implication for abnormal neural development. Dev Genet 18:306-315, 1996. 8. Stemple DL, Anderson DJ. Isolation of a stem cell for neurons and glia from the mammalian neural crest. Cell 71:973-985, 1992. 9. Taylor LE, Bennett GD, Finnell RH. Altered gene expression in murine branchial arches following in utero exposure to retinoic acid. J Craniofac Genet Dev Biol 15:13-25, 1995. 10. Eberwine JH, Yeh H, Miyashiro K, Cao Y, Nair S, Finnell RH, Zettel M, Coleman P. Analysis of gene expression in single live neurons. Proc Natl Acad Sci USA 89:3010-3014, 1992. 11. Elledge RM, Lee WH. Life and death by p53. Bioessays 17:923-930, 1995. 12. Davey S, Beach D. RACH2, a novel human gene that complements a fission yeast cell cycle checkpoint mutation. Mol Biol Cell 6:1411-1421, 1995. 13. Igarashi M, Nagata A, Jinno S, Suto K, Okayama H. Weel +like gene in human cells. Nature 353:80-83, 1991. 14. Miyashita T, Krajewski S, Krajewska M, Wang HG, Lin HK, Liebermann DA, Hoffman B, Reed JC. Tumor suppressor p53 is a regulator of bcl-2 and bax gene expression in vitro and in vivo. Oncogene 9:1799-1805, 1994. 15. Prehn JHM, Bindokas VP, Marcuccilli CJ, Krajewski S, Reed JC, Miller RJ. Regulation of neuronal Bcl2 protein expression and calcium homeostasis by transforming growth factor type b confers wide-ranging protection on rat hippocampal neurons. Proc Natl Acad Sci USA 91:12599-12603, 1994. 16. Sah VP, Attardi LD, Mulligan GJ, Williams BO, Bronson RT, Jacks T. A subset of p53-deficient embryos exhibit exencephaly. Nature Genet 10:175-180, 1995.

GENE EXPRESSION DATA ANALYSIS 17. Hockenberry DM. bcl-2, a novel regulator of cell death. Bioessays 17:631-638, 1995. 18. Veis-Novack D, Korsmeyer SJ. Bcl-2 protein expression during murine development. A m J Pathol 145:61-73, 1994. 19. Barnard JA, Lyons RM, Moses HL. The cell biology of transforming growth factor beta. Biochem Biophys Acta 1032:79, 1990. 20. Millan FA, Denhez F, Kondaiah P, Akhurst RJ. Embryonic gene expression patterns ofTGF bl, b2 and b3 suggest different developmental functions in vivo. Development 111:131143, 1991. 21. Ewen ME, Oliver CJ, Sluss HK, Miller SJ, Peeper DS. p53dependent repression of CDK4 translation in TGF-b-induced G1 cell-cycle arrest. Genes Dev 9:204-217, 1995. 22. Geng Y, Weinberg RA. Transforming growth factor beta effects on expression of G1 cylins and cyclin-dependentprotein kinases. Proc Natl Acad Sci USA 90:10315-10319, 1993. 23. Fujiwara T, Mukhopadhyay T, Cai DW, Morris DK, Roth JA, Grimm EA. Retroviral-mediated transduction of p53 gene increases TGF-b expression in a human glioblastoma cell line. Int J Cancer 56:834-839, 1995. 24. Blaydes JP, Schlumberger M, Wynford-Thomas D, Wyllie FS. Interaction between p53 and Tgib-1 in control of epithelial cell proliferation. Oncogene 10:307-317, 1995. 25. Sokal RR, RohlfJF. Biometry. San Francisco, 1981, 2nd ed. 26. Searle SR. Linear Models. New York: Wiley, 1971. 27. Neter J, Wasserman W, Kutner MH. Applied Linear Regression Models. Homewood/Boston: Irwin, 1989, 2nd ed. 28. SAS Institute. SAS/STAT User's Guide Release 6.03. Cary, NC: SAS Institute, 1990. 29. Rao CR. Linear Statistical Inference and Its Applications. New York: Wiley, 1973. 30. RohlfJF, Bookstein FL. Traditional morphometrics. In Proceedings of the Michigan Morphometrics Workshop (Rohlf JF, Bookstein FL, Eds.). Ann Arbor: University of Michigan Museum of Zoology, 1990, pp 77-122.

91

31. Johnson RA, Wichern DW. Principal components. In Applied Multivariate Statistical Analysis (Conmy SR, Ed.). Englewood Cliffs, NJ: Prentice Hall, 1992, pp 356-395. 32. Mason RL, Gunst RF, Hess JL. Statistical Design & Analysis of Experiments. New York: Wiley, 1989. 33. Ryan JJ, Prochownik E, Gottlieb E, Abel CA, Merino IJ, Nunez G, Clarke MF. c-myc and bcl-2 modulate p53 function by altering p53 subcellular trafficking during the cell cycle. Proc Natl Acad Sci USA 91:5878-5882, 1994. 34. Wrana JL, Attisano L, Carcamo J, Zentella A, Doody J, Laiho M, Wang XF, Massagne J. TGFb signals through a heteromeric protein kinase receptor complex. Cell 71:1003, 1992. 35. Chen RH, Ebner R, Derynck R. Cloning of a type I TGFbeta receptor and its effect on TGF-beta binding to the type II receptor. Science 260:1335-1338, 1993. 36. Clarke AR, Purdie CA, Harrison DJ, Morris RG, Bird CC, Hooper ML, Wyllie AH. Thymocyte apoptesis induced byp53dependent and independent pathways. Nature 362:849-852, 1993. 37. Mercer WE, Shields MT, Amin M, Sauve GJ, Appella E, Romano JW, Ullrich SJ. Negative growth regulation in a glioblastema tumor cell line that conditionally expresses human wild-type p53. Proc Natl Acad Sci USA 87:6166-6170, 1990. 38. Diller L, Kassel J, Nelson CE, Gryka MA, Litwak G, Gebhardt M, Bressac B, Ozturk M, Baker SJ, Vogelstein B. p53 functions as a cell cycle control protein in osteosarcomas. Mol Cell Biol 10:572-578, 1990. 39. Reed JC. bcl-2 and the regulation of programmed cell death. J Cell Biol 124:1-6, 1994. 40. Vaux DL, Cory S, Adams JM. Bcl-2 gene promotes haemopoietic cell survival and cooperates with c-myc to immortalize pre-B cells. Nature 335:440-442, 1988. 41. Vaux DL, Weissman IL, Kim SK. Prevention of programmed cell death in Caenorhabditis elegans by human bcl-2. Science 258:1955-1957, 1992. 42. Wang Y, Szekely L, Okan I, Klein G, Wiman G, Wiman KL. Wild type p53-triggered apoptosis is inhibited by bcl-2 in a v-myc-induced T-cell line. Oncogene 8:3427-3431, 1993.