Characterization of a Genomic Signature of Pregnancy Identified in the Breast

It is well established that a pregnancy completed to term at a young age reduces the risk of breast cancer later in life. This long-term reduction in risk has been attributed to the early differentiation of breast tissue, which otherwise remains undifferentiated and susceptible to carcinogenic insults (1, 2). Because a first full-term pregnancy (FTP) and ensuing breastfeeding are the most significant physiologic events which transform the breast from an immature to a fully mature organ, Russo and colleagues hypothesized that having completed at least 1 FTP would result in a specific, detectable genomic signature in the breast (3–6). Only 1 study to date has examined gene expression in the healthy breast and it was limited to 64 genes (7).

Identification of a specific genomic fingerprint of pregnancy would open up a broad set of opportunities for understanding, and possibly preventing, breast cancer. We, therefore, undertook to compare the gene expression profiles in breast biopsy specimens of healthy parous and nulliparous volunteers from the general population, using a genome-wide approach. Because we were interested in long-term genomic changes associated with FTP, the study was focused on postmenopausal women.

The study was designed to include 2 phases, a discovery phase (DP) and a validation phase (VP) with a total target sample size of 120 women (40 parous and 20 nulliparous in each phase). Recruitment was conducted without interruption between the 2 phases of the study, using the same source population. The parity distribution was reviewed after every group of 10 eligible volunteers. If the nulliparous to parous ratio differed from 1:2, recruitment was limited to the underrepresented group (usually nulliparous) until the ratio reached the 1:2 target.

To assess within- and between-laboratory reproducibility of gene expression profiles in replicate experiments, we conducted a substudy prior to the start of the discovery and validation phases. Breast tissue samples from 4 subjects were processed, and their gene expression profiles were analyzed at 3 independent laboratories (Fox Chase Cancer Center Breast Cancer Research Laboratory of Dr. Jose Russo; University of Memphis Genomic Laboratory of Dr. Thomas Sutter; and Fox Chase Cancer Center Genomic Laboratory) using identical procedures. Data were preprocessed for each of the 3 laboratories separately using the methods described in the statistical methods section. Supplementary Table S1 shows concordance correlation coefficients (8) and Pearson's correlation coefficients for within- and between-laboratory comparisons. The within-laboratory correlations were very high (Pearson's and concordance correlation coefficients: >97%) and very similar to those reported by others (9). The between-laboratory correlations were also high (Pearson's and concordance correlation coefficients: >92%).

Study subjects were recruited at the Sunderby Hospital in Luleå, Sweden, among women who have had a normal mammogram within the year prior to enrollment. Postmenopausal women (defined as lack of menstrual periods for the previous 12 months) between the ages of 50 and 69 years were approached by a research nurse who explained the study procedures and provided informed consent forms. Volunteers who signed informed consent to participate in the study and to donate biological samples for research were scheduled for an interview.

Women who reported a history of any cancer, the use of any hormonal medications in the 6 months preceding their visit, prior breast biopsy or breast implants were excluded. In the discovery phase, it became apparent that women with fatty breast (transparent mammograms) had to be excluded because of low RNA yield. As a consequence, women with similarly fatty breast were considered ineligible for biopsy in the validation phase of the study. The project was approved by the Regional Ethical Review Board for Northern Sweden at the University of Umeå, Sweden.

The study nurse obtained anthropometrical measures (height and weight) and administered the study questionnaire to eligible and consented women. The collected data included a detailed reproductive history, medical history, first-degree family history of breast cancer, smoking, and use of oral contraceptives (OC), hormone replacement therapy (HRT), and other medications.

An experienced intervention radiologist carried out all breast biopsies with a Bard Monopty (C. R. BARD Inc.) automated core biopsy instrument (14 Gauge, 10-cm long, 22-mm penetration depth) through a single small skin incision after the puncture site had been sterilized and anaesthetized (Xylocain and Adrenalin solution, 10 mg/mL + 5 μg/mL, Astrazeneca). Several (3–5) random biopsies were taken from the upper outer quadrant of 1 breast. One biopsy specimen was placed in 70% ethanol for histopathologic analysis, and the remaining ones were immediately placed in RNAlater (Ambion) solution.

The study pathologist has reviewed all tissues to make sure that research biopsies were free of atypia or cancer using criteria published previously (10). This review resulted in the exclusion of 1 study subject (see, Supplementary Fig. S1).

Prior to sending the samples and data to laboratory at the Fox Chase Cancer Center, Philadelphia, all samples were stripped of any personal identifiers and assigned random numbers. The link between the subject's random number and subject's identifiable information was accessible only by the authorized personnel in Sweden. The laboratory personnel at the Fox Chase Cancer Center were blinded to samples' parity status and other personal information.

Total RNA from the core biopsy samples was isolated using the Qiagen Allprep RNA/DNA Mini Kit, according to the manufacturer's instructions (Qiagen). Total RNA was eluted in a final volume of 60 μL (H₂O) and stored at −80°C until further processing. RNA quantity and quality were assessed by means of the Agilent 2100 Bioanalyzer (Agilent Technologies). The amount of total RNA yielded from the core biopsies ranged from 150 ng to 4 μg depending on the ratio of stroma to epithelial tissue.

The GeneChip Expression 3′-Amplification Two-Cycle cDNA Synthesis Kit was used (Affymetrix). Double-stranded cDNA was synthesized from 100 ng of total RNA. An in vitro transcription reaction was then done to produce biotin-labeled cRNA from the cDNA. The cRNA was fragmented before hybridization. A hybridization cocktail, which included the fragmented target, was prepared. The hybridization cocktail was then hybridized to Affymetrix HG_U133 Plus 2.0 oligonucleotide arrays containing probes to 54,675 transcripts. Standard Affymetrix quality control measures (average background, scale factors, and percent present calls) were applied to assess the quality of RNA samples and their subsequent labeling and hybridization, and chips that did not pass the quality control criteria were rejected. Additionally, graphical criteria based on probe-level model (PLM) analysis were applied.

Raw data from array scans were preprocessed and analyzed using the R language for statistical computing (11) and Bioconductor (12), an open source software for bioinformatics. The data were preprocessed using the Robust Multichip Average method implemented in the Bioconductor package affy that includes background correction, quantile normalization, and summarization of expression values (13–15). Probes for which the proportion of present calls was less than 75%, and the difference in the proportion of present calls between parous and nulliparous women was less than 25% were filtered out. Probes with low coefficient of variation across samples (below first quartile) were also removed. These filtering criteria left 19,028 probes for analysis in the discovery phase and 17,750 probes in the validation phase. The overlap between the 2 sets of probes consisted of 16,002 probes.

The microarray experiments in both phases were conducted in 8 batches. To account for potential between-batch variability, an Empirical Bayes method, implemented in the COMBAT software, developed by (16) and written in R, was used. We also corrected for batch effects in the analysis. Additionally, the quality control duplicate samples were used to evaluate the batch effects and the effectiveness of batch adjustments.

To identify genes differentially expressed between parous and nulliparous samples, we used the following 3 methods: significance analysis of microarrays (SAM; method 1; ref. 17) implemented in the R package samr, surrogate variable analysis (SVA; method 2; ref. 18, 19) implemented in the R package sva, and logistic regression analysis (LRA; method 3).

It has been shown that genetic, environmental, demographic, and technical factors may have substantial effects on gene expression (18–21). In addition to measured variables of interest, there might be sources of signal because of unknown or unmeasured factors. Leek and Storey (18) showed that failing to incorporate these sources of heterogeneity into analysis can result in both spurious and masked associations. They introduced SVA to overcome the problems caused by heterogeneity in gene expression studies and showed that SVA increases the biological accuracy and reproducibility of gene expression studies. SVA uses a residual expression matrix, obtained by removing the effects of the outcome variable (parity status in our study) on expression, to estimate, via singular value decomposition of the residual matrix, the signatures of expression heterogeneity in terms of an orthogonal basis of singular vectors. Statistical procedures are then used to assess the significance of these signatures, to identify the subset of genes driving each signature, and to form surrogate variables (SV) based on the signatures of the corresponding subsets of genes in the original expression data. The resulting SVs are used to adjust the analysis of the associations between genes and parity status. For each gene, an unadjusted P value measuring the significance of that gene as an independent predictor of the outcome variable is calculated using logistic regression that adjusts for SVs (method 2).

We also used LRA (method 3) to identify differentially expressed genes while controlling for the effects of potentially confounding factors that were measured in the study, such as body mass index (BMI), OC use history, HRT, and smoking history. To select a subset of the measured characteristics for inclusion in the LRA for adjustment, we compared them to the significant SVs derived in method 2 using the Spearman's correlation. The top 5 characteristics that were most significantly associated with 1 of the SVs in the discovery phase were selected to be adjusted for in the LRA. These characteristics were BMI, HRT duration, breast density, smoking duration, and OC use history. Other combinations of characteristics significantly associated with SVs were also adjusted for in the logistic regression analyses of parity status. For each gene, an unadjusted P value measuring the significance of that gene as an independent predictor of parity status was calculated using logistic regression adjusting for the selected variables.

The false discovery rate (FDR) approach was used to control for multiple comparisons (17, 22, 23). SAM (method 1) computes a 2-sample t test-like statistic for each gene and uses a permutation procedure to estimate FDRs which are used to select differentially expressed genes. For SVA and logistic regression, the QVALUE method (22) implemented in the R package qvalue was used to adjust P values for multiple comparisons. Genes with an FDR of less than 10% and at least a 1.25-fold change between parous and nulliparous samples were considered statistically significant.

To derive a genomic predictor of pregnancy in the breast, we used 5 classification methods: “neareast shrunken centroids” method implemented in the Bioconductor package pamr (24), support vector machine implemented in the R library e1071 (25), Classification and Regression Trees implemented in the R library rpart (26), boosted classification trees using the AdaBoost algorithm (27) and random forest implemented in the R package randomForest (28). The approximately 500 most significant genes based on a combined FDR and fold-change criterion were used for these analyses except for methods that carry out automatic variable selection (e.g., nearest shrunken centroids classifier). We then used the genomic classifiers identified in the discovery phase to estimate the probability of being parous (parity score) for each woman in the validation phase. The significance of the genomic signature of parity was evaluated using the logistic regression model with FTP as the response variable and the parity score (genomic predictor) along with potential confounding variables, such as BMI, breast density, HRT, and OC duration and smoking history, as independent predictors. Sensitivity and specificity of each genomic classifier derived in the discovery phase were evaluated in the validation phase. The statisticians were blinded to the parity status of the women in the validation phase.

Supplementary Fig. S1 provides a workflow of subject accrual and sample processing. A total of 389 women were interviewed between September 2008 and May 2009. Among these, 134 (34%) were excluded based on the eligibility criteria and 4 (1%) cancelled their interview. This resulted in 251 women (111 nulliparous and 140 parous) included in the study. Two (0.5%) women were excluded later in the study when it was found that one of them had breast cancer and the other one had premenopausal FSH levels. Additionally, 123 women were excluded after RNA extraction because of RNA degradation, absence of epithelial structures or an insufficient amount of RNA. The remaining 126 women (44 nulliparous and 82 parous) were included in the current study. Nineteen microarray chips were rejected based on standard Affymetrix quality control measures (average background, scale factors, and percent present calls) and based on PLM analysis. This left 107 chips for the differential expression analysis of parous versus nulliparous women: 55 in the discovery phase (37 parous and 18 nulliparous) and 52 in the validation phase (30 parous and 22 nulliparous).

Table 1 presents characteristics of parous and nulliparous women in the discovery and validation phases. There were no statistically significant differences between parous and nulliparous women within each study phase or between women in the discovery and validation phases within each parity group.

Using SVA, we identified 11 significant SVs in the discovery and 9 in the validation phase. These variables potentially have a significant effect on gene expression. Two of the SVs in both phases accounted for over 10% of the variation in gene expression (data not shown). Supplementary Tables S2A and SB show Spearman's correlation coefficients of the SVs, including batch. Prior to batch adjustment, the batch variable was significantly associated with SV1 (ρ = 0.56, P = 0.01), indicating that batch adjustment was required. It was no longer significantly associated with any of the SVs after adjustment using the COMBAT method. BMI, HRT duration, OC duration, breast density, and smoking history were significantly associated with the SVs indicating that these factors might impact gene expression. These variables were controlled for in LRA (method 3). Some of the SVs were found to be significantly associated with FTP (e.g., SV1 and SV4 in the discovery phase and SV2, SV4, and SV6 in the validation phase) and were, therefore, excluded from SVA analysis (method 2) because our objective was to identify genes significantly associated with FTP.

Table 2 presents the numbers of statistically significant genes identified using the 3 statistical methods in the discovery and validation phases at an FDR of 10% and with at least a 1.25-fold change. The numbers of differentially expressed genes were much higher in the discovery than in the validation phase. In both phases combined, depending on the statistical method used, between 228 and 288 genes were identified as differentially expressed, whereas 218 genes were identified as differentially expressed by all 3 methods. Upregulated genes were found to be more reproducible than the downregulated genes with 62% to 64% of the significant genes identified in the validation phase being also significant in the discovery phase, compared with 9% to 13% of downregulated genes. SAM and LRA yielded the highest proportion of reproducible genes: 45% of genes significant in the validation phase were also significant in the discovery phase. Using LRA, 305 probe sets identified as significantly differentially expressed in the discovery phase and confirmed in the validation phase. These genes are reported in Supplementary Table S3.

Five genomic predictors of parity were derived in the discovery phase using the 5 classification methods described in the Statistical Methods section, and their significance was evaluated in the validation phase using logistic regression models that adjusted for clinical variables. Table 3 shows the estimated coefficients, standard errors, and P values of the genomic predictors and clinical variables in the logistic regression models applied to the validation phase with FTP as dependent variable. The results indicate that the genomic predictors derived using the discovery phase data remained significant predictors of parity in the validation phase with and without adjustment for other variables (data without adjustment are not shown). Boosted classification trees and nearest shrunken centroids consistently carried out better than the other classification methods. The prediction accuracy of the best classifiers was estimated in the validation phase to be between 65% and 75%.

To our knowledge, this is the first study that seeks to comprehensively characterize in an unselected, population-based group of healthy volunteers the differences in gene expression in breast core biopsy specimens between parous and nulliparous women. Using a discovery validation approach, 274 upregulated and 31 downregulated probe sets were identified, which may constitute the core genomic signature that distinguishes the breast of parous postmenopausal women from that of nulliparous ones. Using supervised learning methods, we derived a genomic signature of parity in breast specimens from the subjects included in the discovery phase of the project and established that it was a significant independent predictor of parity in the subjects in the validation phase.

The genes differentially expressed in parous and nulliparous postmenopausal women are presented in Supplementary Table S3 and involved in regulation of transcription, centrosome organization, RNA splicing, cell-cycle control, adhesion, and differentiation. Among upregulated genes, EZH2 is a member of polycomb group of proteins involved in maintaining transcriptional repression of genes. This gene acts as a tumor suppressor and also functions as a histone methyltransferase (29). NINL, TRAF5, and SFI1 are upregulated in parous breast and involved in centrosome organization and maintaining the microtubule cytoskeleton. CDK3, MCTS1, and SYCP2 are involved in cell-cycle control. PRPF39, LUC7L3, HNRNPA1, HNRNPA2B1, PNN, PABPN1, RBMX, SNRNP200, PRPF4B, and SFPQ are involved in RNA splicing. Among the downregulated genes, CLDN10, CD36, PDZD2, CLSTN2, LAMA4, PCDH9, and SORBS1 play role in cell adhesion. It is of interest that upregulated genes were more prevalent and consistent than downregulated genes in parous compared with nulliparous women. This suggests that parity results mainly in overexpression of genes involved in breast cell differentiation, organization, and tumor suppression as opposed to downregulation of genes that might drive development of cancer. The details and biological significance of the genes and pathways differentially expressed in parous versus nulliparous breast will be discussed in a separate paper.

Recently, Asztalos and colleagues (7) examined gene expression in the normal premenopausal human breast, comparing nulliparous, recently parous (0–2 years since pregnancy), and distantly pregnant (5–10 years after pregnancy) age-matched premenopausal women. They analyzed a customized 64-gene set focusing on the genes involved in inflammation, extracellular matrix remodeling, angiogenesis, and estrogen signaling. They reported that 14 of the 64 selected genes were differentially expressed in parous versus nulliparous breast tissues. Compared with nulliparous breast, parous breast had significant upregulation of genes related to inflammation (CCL21, LBP, SAA1/2, and IGKC) and downregulation of genes involved in angiogenesis (VEGFA) and estrogen signaling (ERα, PGR, and ERBB2; ref. 7).

There was no overlap between the differentially-expressed genes reported here and those reported by Asztalos and colleagues (7), an inconsistency possibly related to differences in menopausal status, study populations, and laboratory methods. Compared with the hypothesis-driven report by Asztalos and colleagues (7), our study took a comprehensive approach and addressed a much broader list of genes. It also focused on older, postmenopausal women because the objective was to detect long-term gene expression changes, that is, gene expression differences between parous and nulliparous that could be observed many years after the first FTP.

The study's strengths include the formal 2-phase approach for analyses of genomic differences between parous and nulliparous breast, with independent discovery and validation phases but identical procedures throughout the study. The comprehensive gene assessment with microarray assays is strength of our study. The study focused on healthy subjects attending a mammography clinic and, therefore, representative of the general population of Sweden, a country where mammography is widely accepted. In addition to using rigorous research procedures, stringent criteria were used to control for multiple comparisons and FDR. The laboratory personnel were blinded to samples parity status. In addition, all data analysts were blinded to the parity status of subjects in the validation phase until the parity status predictions were made and the parity scores were derived for the subjects in the validation phase. The consistency of results across 3 different statistical methods used (SAM, SVA, and LRA) strengthened our confidence in the study results.

There were some limitations as well. The study population was restricted to residents of the northernmost part of Sweden, and all participants were of Swedish or Finnish ethnicity. It was felt that these characteristics would be advantageous to avoid gene expression variations resulting from differences in ethnicity rather than parity. The study results, however, should be confirmed in other populations. In addition, although the vast majority of eligible women accepted rather enthusiastically to participate in the study, some women were excluded from the study based on eligibility criteria (34%) or because of RNA degradation, absence of epithelial structures, or an insufficient amount of RNA (32%). Because pregnancy may affect mammographic density, exclusion of women with low-density mammograms may have resulted in differential selection of women at higher, or lower, risk of breast cancer between parous and nulliparous subjects. This is an issue, though, not easily addressable because examination of gene expression cannot be done unless both epithelial structures and sufficient RNA of good quality are present.

We used a FDR of 10%. However, the true FDR corresponding to our list of significant genes is likely to be much lower than 10% because only probes that passed the FDR of 10% in both the discovery and the validation phases were included in the list and those that passed the FDR of 10% in only 1 of the 2 phases were excluded. Additionally, because we were studying normal, rather than pathologic tissues, we a priori expected modest effect sizes (e.g., fold change of 1.25) for genomic changes associated with pregnancy in healthy postmenopausal women.

In summary, the results provide initial support to the concept that an FTP induces permanent genomic changes in the breast, thus reflecting the well-known permanent phenotypical changes that follow an FTP. Once further confirmed in additional populations with wider ranges of age, ethnicity, and other characteristics, a well-characterized genomic signature of pregnancy could be used as an intermediate marker for instance to assess potential chemopreventive interventions with hormones mimicking the effects of pregnancy.

No potential conflicts of interest were disclosed.

The authors thank the Swedish Norrbotten County Council for their strong support, the women of Norrbotten County for their generous contribution and participation in the project; the staff of the Mammography Department, Sunderby Hospital, Luleå, whose competence and skills made this study possible; Dr. Thomas Sutter for his contributions and Mr. Sören Holmgren for his dedication to the project.

This work was supported by grant 02-2008-034 from the Avon Foundation for Women Breast Cancer Research Program.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.