|Home | About | Journals | Submit | Contact Us | Français|
Nonmalignant human mammary epithelial cells (HMEC) seeded in laminin-rich extracellular matrix (lrECM) form polarized acini and, in doing so, transit from a disorganized proliferating state to an organized growth-arrested state. We hypothesized that the gene expression pattern of organized and growth-arrested HMECs would share similarities with breast tumors with good prognoses. Using Affymetrix HG-U133A microarrays, we analyzed the expression of 22,283 gene transcripts in 184 (finite life span) and HMT3522 S1 (immortal nonmalignant) HMECs on successive days after seeding in a lrECM assay. Both HMECs underwent growth arrest in G0–G1 and differentiated into polarized acini between days 5 and 7. We identified gene expression changes with the same temporal pattern in both lines and examined the expression of these genes in a previously published panel of microarray data for 295 breast cancer samples. We show that genes that are significantly lower in the organized, growth-arrested HMEC than in their proliferating counterparts can be used to classify breast cancer patients into poor and good prognosis groups with high accuracy. This study represents a novel unsupervised approach to identifying breast cancer markers that may be of use clinically.
Loss of growth control and disruption of tissue architecture are among the earliest hallmarks of cancer. We hypothesized that the gene expression changes that occur during the organization and growth arrest of cultured mammary acini would share similarities with breast tumors that had a good prognosis. To test this hypothesis, we examined temporal changes in gene expression in two nonmalignant human mammary epithelial cells (HMEC) grown in a three-dimensional laminin-rich extracellular matrix (lrECM) assay (1). In this model, single nonmalignant breast epithelial cells form polarized growth-arrested multicellular structures that resemble acini over a period of several days. This transition from the unpolarized actively dividing state to a polarized nondividing state is the reverse of what happens in the early stages of tumorigenesis.
We used two nonmalignant HMECs, a nonimmortalized HMEC strain, 184 (2, 3), and a spontaneously immortalized cell line, HMT3522 S1 (4), and monitored the changes in gene and protein expression as the cells formed acinar structures in three-dimensional lrECM cultures. During the process of self-organization and withdrawal from the cell cycle, we noted a progressive hypo-phosphorylation of the retinoblastoma (Rb) gene and an induction of the cyclin-dependent kinase (CDK) inhibitor p27kip1. Analysis of the changes in gene expression in both cell lines allowed the identification of sets of commonly regulated genes.
Although established breast cancer prognostic markers, such as tumor size, grade, nodal, and hormone receptor status, are useful in predicting survival in large populations (5–7), there is a pressing need to develop better prognostic signatures to predict recurrence and overall survival. A particular benefit would be the identification of patients with good prognoses, whose tumors are highly unlikely to recur and who nevertheless are being treated with cytotoxic chemotherapies (8). The advent of gene expression technologies has greatly aided the identification of molecular signatures with value for tumor classification and prognosis prediction (9–14). van de Vijver et al. have developed a 70-gene signature that effectively stratifies patients into good and poor prognosis groups (15, 16). Paik et al. (17) have proposed a 21-gene signature with which to calculate a “recurrence score” that predicts the likelihood of recurrence in estrogen receptor (ER)-positive lymph node–negative patients. In each of these studies, the predictive signatures have been derived by using a training set of patients of known outcome followed by testing these signatures in a validation set of patients. In contrast, our approach has been to identify directly genes, the expression of which changes as cultured mammary epithelial cells transition from a disorganized to an organized state, and to then test the proof of principle and the possible use of these genes as prognostic markers in a validation set of patients.
Finite life span 184 HMECs were obtained from reduction mammoplasty tissue and grown in a serum-free MCDB 170 medium (mammary epithelial growth medium, Clonetics Division of BioWhittaker, Walkersville, MD) as described previously (2, 18). HMT3522 S1 mammary epithelial cells were cultured in a H14 medium (DMEM/F12 containing 250 ng/mL insulin, 10 μg/mL transferrin, 2.6 ng/mL sodium selenite, 10−10 mol/L estradiol, 1.4 × 10−6 mol/L hydrocortisone, 10 ng/mL epidermal growth factor, and 5 μg/mL prolactin). The cells were cultured in a three-dimensional lrECM (Matrigel, BD Biosciences, Franklin Lakes, NJ) as described (19). The colonies were isolated from the Matrigel in ice-cold PBS/5 mmol/L EDTA after 3, 5, 7, and 10 days. For Western blot analysis, the colonies were lysed in 150 mmol/L NaCl, 1% NP40, and 50 mmol/L Tris (pH 8.0). RNA was isolated using the RNeasy kit (Qiagen, Valencia, CA) according to the manufacturer’s instructions.
Acinar structures were fixed on glass slides in methanol-acetone (1:1) at −20°C for 10 minutes and air dried. A primary block was done in an immunofluorescence buffer (130 mmol/L NaCl, 7 mmol/L Na2HPO4, 3.5 mmol/L NaH2PO4, 7.7 mmol/L NaN3, 0.1% bovine serum albumin, 0.2% Triton X-100, 0.05% Tween 20) plus 10% goat serum for 1 hour at room temperature. A secondary block was done in immunofluorescence buffer plus 10% goat serum plus 20 μg/mL goat anti-mouse F(ab′)2 fragment for 30 minutes. The primary antibody rat anti-α6 integrin (Chemicon International, Temecula, CA) was diluted 1:50 in the latter blocking buffer and incubated overnight (15–18 hours) at 4°C. Secondary antibodies conjugated with fluorescent dye FITC or Texas red were diluted 1:100 in blocking buffer and incubated 1 hour at room temperature. Immunofluorescent images were acquired using an inverted microscope equipped with a digital camera and SPOT software. Confocal analysis was done using a Zeiss 410 confocal microscopy system (Zeiss, Thornwood, NY). The images presented are representative of two or more independent experiments. All the images were converted to TIFF format and arranged using Adobe Photoshop 5.0.
To analyze the DNA content, acini were suspended in 300 μL 0.25% trypsin and incubated at 37°C for exactly 10 minutes. The dispersed cells were washed thrice in PBS and fixed in 40% ethanol at 4°C overnight. They were then incubated with 500 μg/mL RNase A in PBS at 37°C for 30 minutes and stained with 69 μmol/L propidium iodide solution in PBS at room temperature for 30 minutes. The DNA content was determined by flow cytometry using a FACScan (Becton Dickinson, Hialeah, FL), and the data were analyzed with CellQuest software (Becton Dickinson). Statistics for the reduction in S-phase kinetics were based on the means of repeated measurements: S1, day 3, n = 4; S1, day 5, n = 5; S1, day 7, n = 6; 184, day 3, n = 2; 184, day 5, n = 1; and 184, day 7, n = 2.
Whole-cell lysates (20–100 μg protein) were separated by SDS-PAGE and transferred electrophoretically to polyvinylidene difluoride membranes (Millipore, Billerica, MA). After blocking in 10% nonfat dry milk for 90 minutes at room temperature, they were incubated with the following antibodies at 1:500 to1:1,000 dilution in 2% nonfat dry milk for 2 hours at room temperature or overnight at 4°C: rabbit anti-phosphorylated Rb at S807/811 and S795 (Cell Signaling Technology, Beverly, MA), rabbit anti-phosphorylated Rb at S612 and T821 (Biosource, Camarillo, CA), and mouse anti-p27kip1 (BD Transduction Laboratories, Franklin Lakes, NJ). The blots were then incubated with horseradish peroxidase–conjugated sheep anti-mouse IgG or anti-rabbit IgG at 1:2,000 in 2% nonfat dry milk for 90 minutes at room temperature. The blots were developed using SuperSignal West Pico Chemilumiscent Substrate (Pierce, Rockford, IL).
A database consisting of the microarray profiles of 295 human breast tumors with the associated clinical data (15) was obtained from Rosetta Inpharmatics (http://www.rii.com/publications/2002/nejm.html). For survival analysis of the 19 individual marker genes, the patients were stratified into quartiles for expression of each marker, and the survival curves were computed using the method of Kaplan-Meier. Statistical significance was determined using the log-rank test. Statistical analyses were done using GraphPad Prism. For survival analysis of the set of 249 marker genes, the patients were stratified into two groups using GeneSpring software by hierarchical cluster analysis with a distance metric of the expression pattern of all 249 genes. Kaplan-Meier survival curves, log-rank statistics, and the estimated hazard ratio for these two groups were computed using the Excel add-in EcStat.
Cell samples were harvested in duplicate at three time points, 3, 5, and 7 days, after seeding in lrECM. Purified total cellular RNA was biotin labeled and hybridized to human oligonucleotide microarrays (Affymetrix HG-U133A; Affymetrix, Santa Clara, CA) as described previously (20). Experiments with Affymetrix-present P-call rates of >30% were included in the analysis. Signal values from each of the 22,283 probe sets were calculated by means of robust multiarray analysis (21) using Bioconductor (http://www.bioconductor.org) in the R computing environment. The signal values were inverse log2 transformed and then imported into GeneSpring software (Silicon Genetics, Palo Alto, CA), and each array was normalized to its median signal intensity. The genes were normalized to the mean of the 3-day time point for each cell type independently. For method 1 of selecting significantly differential, temporally coregulated genes, significantly up-regulated genes in each cell specimen were identified by first selecting the genes induced at least 1.5-fold in at least one of the six conditions and then doing an ANOVA analysis as a function of time. Variances were calculated using the cross-gene error model (GeneSpring), Ps cutoff 0.05 (multiple testing correction: Benjamini-Hochberg false discovery rate). About 5% of the identified genes in each set would be expected to pass this restriction by chance. Significantly down-regulated genes were identified in the same manner after normalizing to the 7-day time point. Genes that were up-regulated or down-regulated early in each cell line were selected from the significantly up-regulated or down-regulated gene lists. The early genes were defined as those with a mean expression at 5 days of at least 50% of their mean expression at 7 days. For method 2 of selecting significantly differential, temporally coregulated genes, we first selected all genes that were at least 1.5-fold differential (up or down) in at least one of the four samples from the two later time points (days 5 and 7) for either cell specimen. We then did an ANOVA as a function of time. Variances were not assumed equal (Welch ANOVA), Ps cutoff 0.05 (multiple testing correction: Benjamini-Hochberg false discovery rate). From this list of genes that were significantly differential in either cell line, we then identified genes that were up-regulated or down-regulated early (mean expression at 5 days of at least 50% of their mean expression at 7 days) or late in each cell line. We then identified those genes that were coordinately regulated in both cell lines.
To identify consistent changes in gene and regulatory protein expression levels, we characterized two independently derived nonmalignant HMECs: one finite life span strain (184; ref. 2) and one spontaneously immortalized line (HMT3522 S1; ref. 4). Both cells formed acinus-like structures with similar morphology and basal polarity when cultured from single cells in lrECM (Fig. 1A and B; see also ref. 1). We did temporal studies to determine when, and in what order, changes in critical cell cycle regulatory molecules occurred in nonmalignant HMECs in three-dimensional lrECM cultures. Flow cytometric analyses indicated that, after undergoing a limited number of cell divisions, the majority of nonmalignant S1 cells accumulated in the G0–G1 phase of the cell cycle by day 7 (Fig. 1C). The product of the Rb susceptibility gene, a central player in the G1-S transition, is inactivated by phosphorylation, allowing cell cycle progression. Rb inactivation occurs through the sequential actions of cyclin D/CDK4 and CDK6 and cyclin E-CDK2 complexes (22). We analyzed cell cycle regulators known to affect G1 checkpoints in a time-dependent manner in S1 cells in three-dimensional cultures. Phosphorylation of several sites on Rb was found to gradually decrease between days 5 and 10 in these nonmalignant cells, consistent with the growth-suppressive role of the hypophosphorylated form (Fig. 1D). Cyclins D1, E, and A, as well as their binding partners, CDK4, CDK6, and CDK2, also decreased during this period (data not shown). In contrast, the CDK inhibitor p27kip1 increased between days 5 and 10 (Fig. 1D). The pronounced down-modulation of Rb phosphorylation and the elevation of p27kip1 protein levels were changes observed in both the S1 and 184 cells (data not shown). Although other studies have used different HMECs (MCF10A) to determine how cells can escape normal growth control (23–26), the mechanisms by which mammary cells actually initiate and maintain growth arrest during the process of acini formation in the context of three-dimensional lrECM had remained to be determined.
To probe systematically the molecular changes that accompany acinus formation, we analyzed the expression profiles of 22,283 transcripts using Affymetrix HG-U133A microarrays. Microarray experiments were done with biological duplicates using RNA samples harvested from S1 and 184 cells after 3, 5, and 7 days of culture in lrECM. Comparison of the decrease in the percentages of cells in S phase showed that growth arrest occurred with kinetics that were significantly correlated (>95% confidence level) in both cell types (correlation coefficient, 0.89; Fig. 2A). Further, the S-phase decrease was significant (P = 0.05, ANOVA). Hence, we reasoned that the gene expression changes important for these processes would follow a common temporal pattern in both cell lines and that changes that were cell type-specific could be disregarded.
We first identified genes that showed at least 1.5-fold changes during the time course in the individual cell specimens (P < 0.05, ANOVA; within this window, 363 genes were up-regulated and 117 genes down-regulated in 184 cells and 234 genes were up-regulated and 351 genes down-regulated in S1 cells). We then divided these lists into genes, the expression of which changed ‘early’ by our definition (between days 3 and 5) or ‘late’ (between days 5 and 7) in S1 and 184. Finally, we identified the genes from each temporal group that were common to both cell types (Fig. 2B). A total of 60 genes with common temporal patterns were identified, including 21 genes that were up-regulated early, 11 genes that were up-regulated late, 6 genes that were down-regulated early, and 22 genes that were down-regulated late (Fig. 2C and D; Supplementary Data). The magnitude of the expression changes of the 22 down-regulated genes in HMEC ranged from 2.5-fold (ACTB) to 5.4-fold (TRIP13).
To relate the process of acinar development in three-dimensional lrECM cell cultures to the changes that occur in breast cancer, we examined the expression levels of the differentially regulated genes identified using our model using previously published microarray data for a panel of 295 breast cancer samples from the fresh-frozen tissue bank of the Netherlands Cancer Institute (Amsterdam, the Netherlands), including 151 lymph node–negative disease and 144 lymph node–positive disease patient samples (15). Fifty-five of the 60 genes selected in our three-dimensional culture analysis were included on these microarrays. We looked at 5- and 10-year survival data and applied Student’s t test to determine how many of the genes modulated in three-dimensional cultures showed survival-associated expression changes. Student’s t tests were done to determine whether the difference in the expression level of a given gene in two groups (e.g., patients who survived 5 years versus patients who did not) was large enough that it was not likely to be due to chance. The numbers and percentages of genes exhibiting significantly different expression in the tumors of patients with differential survival (P < 0.05) were tabulated for (a) all the genes represented on the microarrays, (b) genes selected based on differential expression during the three-dimensional lrECM time course, or (c) randomly generated gene lists (Table 1). The percentage of genes with survival-associated expression changes was highest for those genes down-regulated late (between days 5 and 7) in the time course. The percentage for this gene list exceeded those of the unfiltered list of all 25,773 genes represented on the arrays, 5 random gene lists, and all other three-dimensional lrECM gene lists. The list of genes that were down-modulated late in the lrECM time course showed a marked enrichment in genes, the expression level of which correlated to 5-year (68%) and 10-year (53%) survival. The levels of the majority of the late down-regulated genes were higher in patients who died within 5 or 10 years (ACTB, VRK1, ODC1, CKS2, FLJ10036, FLJ10540, FOXM1, RRM2, TRIP13, CDKN3, STK6, FLJ10517, TUBG1, ACTN1, TNFRSF6B, and EPH2), whereas the levels of three genes (DUSP4, HBP17, and EIF4A1) were lower in these patients. The magnitude of the expression changes of the down-regulated genes in the 250 tumor samples ranged from 1.5-fold (FLJ10036) to 9.3-fold (HPV17).
We identified 22 genes in all that were down-regulated in both HMECs between days 5 and 7 of lrECM culture. Of these, 19 genes were represented in the published data set for the 295 patient tumor samples. We stratified the 295 tumors into quartiles based on the relative expression level of each of the genes in the selected set and further analyzed the relationship of the expression level of each individual gene to survival (Fig. 3). The resulting Kaplan-Meier curves showed that gene expression levels correlated significantly with outcome for 14 of the 19 selected markers. For 13 of the 14 markers, gene expression was lower in tumors from patients with better outcomes, whereas in one case (DUSP4) gene expression was lower in tumors from patients with poorer outcomes.
To test whether expression levels of the 19 selected marker genes correlated with lymph node status, we calculated Pearson product-moment correlation coefficients. Although no genes were correlated at the 95% confidence level, expression levels of 4 of the 19 genes showed a trend toward a correlation with lymph node number or status (80% confidence level), including DUSP4, HBP17, TNFRSF6B, and TUBG1.
Collective gene signatures have the potential to discriminate among clinical end points more accurately than markers used individually. Hence, we tested the ability of our set of 19 genes to classify breast cancer patients into prognostic groups. We used hierarchical cluster analysis to separate the patients into groups and then determined the overall 10-year survival rates for these groups (Fig. 4). The cluster analysis separated the patients into five groups, three of which had tumors that expressed comparatively lower levels of most of the 19 genes and two of which expressed higher levels. The 10-year survival rates for these five groups were 95%, 84%, 67%, 61%, and 54%, respectively.
To test whether other sets of genes down-regulated late in the lrECM time course identified by using other selection strategies would also include useful breast cancer markers, we applied a second selection strategy (Fig. 5A) and tested the ability of the resulting gene set to predict breast cancer prognosis. This second method was less restrictive than the first and resulted in the identification of 287 genes that were significantly down-regulated late in the three-dimensional time course of both HMEC specimens (see Supplementary Data for complete gene list and gene expression information). Seventeen of the 22 genes selected using method 1 were also included in the 287 genes selected using method 2. We tested the ability of these genes to predict breast cancer prognosis by using hierarchical cluster analysis in the same set of previously published microarray data from 295 breast cancers (15). A large majority, 249 of the 287 genes, was included on these microarrays. Of the 17 overlapping (methods 1 and 2) genes on the Affymetrix chips, 15 were present in the Rossetta data set: ACTB, ACTN1, CDK3, CKS2, DUSP4, EPHA2, HBP17, FOXM1, ODC1, RRM2, STK6, TNFRSF6B, TRIP13, TUBG1, and VRK1.
Hierarchical cluster analysis using the 249-gene signature classified the samples into two groups of approximately equal numbers of tumors (Fig. 5B). Overall 10-year survival rates were 90% (138 of 154) for the good prognosis group and 59% (83 of 141) for the poor prognosis group. To assess the significance of these predictions and take into account patients that could not be followed the entire length of the study, we did a Kaplan-Meier analysis. The results show that the 249-gene profile was highly informative in identifying patients with poor outcome (P = 2.7 × 10−10, log-rank test; Fig. 5C). The estimated hazard ratio for poor outcome (failure to survive) in the group with the poor prognosis signature compared with the good prognosis signature was 4.7 [95% confidence interval (95% CI), 2.8–7.9].
Three-dimensional lrECM cultures permit nonmalignant cells to exhibit self-organizing properties. Such cultures provide models that allow the study of processes that are aberrant in breast cancer (1, 19, 27, 28). To understand how breast epithelial cells transition from a proliferating, unorganized state to a resting, organized state and to relate this process to the opposing changes that occur in clinical breast cancer, we did genome-wide gene expression profiling for two nonmalignant HMECs in three-dimensional cultures and used published data from a panel of 295 breast cancer samples. We show that genes that are down-regulated as HMEC transition from a proliferative to an acinar phenotype can be used collectively as signatures to predict clinical breast cancer prognosis.
Our approach represents a new way to identify genome-wide cancer prognostic markers. We based all the marker selection steps on HMEC cultured in three-dimensional lrECM. The three-dimensional model system provided a means to focus on epithelial cells themselves and a defined, highly relevant biological process, the formation of breast acini. Whereas the stroma is absent, the three-dimensional lrECM assays seem to substitute for myoepithelial cells and other signals that are needed to form organized acini (29). Differentially expressed genes identified using this model system are likely to be functionally linked to the transformation relevant process. Further, we have applied an unsupervised method (hierarchical cluster analysis) to classify the patient samples using selected markers. Hence, neither our method of marker selection nor our sample classification method relies on any clinical information.
Gene expression profiling of tumors using DNA microarrays is a promising method for predicting prognosis and treatment response in cancer patients (17, 30–32). Two studies have recently used genome-wide microarray analysis to identify gene signatures that predict prognosis in breast cancer. The profiles studied by both groups of researchers were reported to be more powerful predictors of the outcome of disease than standard systems based on clinical and histologic criteria. The study by van’t Veer et al. (16) used the supervised classification of a primary data set of 78 tumor samples to identify a 70-gene signature that divided the samples into classes of poor and good prognosis. Validation of the classifier by a second overlapping set of 295 tumors showed that it accurately predicted 10-year survival in breast cancer patients (15). Overall, the 10-year survival rates were 54.6 ± 4.4% and 94.5 ± 2.6% for the poor and good prognosis groups, respectively. Their estimated hazard ratio for poor outcome (distant metastases) in the group with a poor prognosis signature compared with the group with the good prognosis signature was 5.1 (95% CI, 2.9–9.0; P < 0.001). In a similar study, Wang et al. (33) also used supervised classification in a training set of 115 breast tumor samples. They identified a 76-gene signature, which was verified in an independent set of 171 breast tumors. Their estimated hazard ratio for poor outcome (distant metastases within 5 years) in the group with the poor prognosis signature compared with the group with the good prognosis signature was 5.7 (95% CI, 2.6–12.4). Our 249-gene signature predicted 10-year survival rates of 59% and 90% for poor and good prognosis groups, respectively. The estimated hazard ratio for poor outcome (failure to survive) in the group with the poor prognosis signature compared with the good prognosis signature was 4.7 (95% CI, 2.8–7.9).
Our 249-gene signature overlapped by 11 genes with the 70-gene signature of van’t Veer et al. (including NUSAP1, UCHL5, RAMP, DC13, PRC1, COL4A2, PITRM1, CENPA, MELK, KNTC2, and MCM6; ref. 16) and by 7 genes with the 76-gene signature of Wang et al. (including MTB, POLQ, SUPT16H, FEN1, DUSP4, PLK1, and SMC4L1; ref. 33). The van’t Veer and Wang signatures overlapped by a single gene (cyclin E2). Our 19-gene signature had no genes in common with either of the previously published signatures.
Our 19-gene signature included several genes encoding proteins with roles in the cell cycle and in cell division. The cell cycle genes were previously identified as important markers of prognosis in ER-positive younger patients (34). In this earlier study, 50 cell cycle–related genes divided 83 ER-positive younger patients into two groups of good versus poor prognosis. The overall 10-year survival rates were 46% and 96% for the poor and good prognosis groups, respectively. Similarly, we found that a core group of predominantly cell cycle and mitotic organizing center genes (CDKN3, RRM3, FLJ10540, FOXM1, STK6, TRIP13, EIF4A1, FLJ10036, VRK1, TUBG1, CKS2, and FLJ10517) made a strong contribution to stratifying tumors into good versus poor prognostic groups. One gene from our 19-gene signature, STK6, which encodes Aurora-A (35, 36), was also included in the 50-gene signature of Dai et al. (34), whereas 15 genes from our 249-gene signature were included in the Dai et al. signature (including BM039, DKFZp762E1312, LOC51203, LOC51659, ID-GAP, KNSL6, PRC1, STK6, CDC45L, SNRPA1, H2AFZ, CENPA, CDC6, BIRC5, and BLM). In addition to cell cycle genes, our prognostic genes also encoded products with other functions, including genes involved in cytoskeletal regulation (ACTB and ACTN1; ref. 37), cell survival (TNFRSF6B; refs. 38, 39), polyamine biosynthesis (ODC1; ref. 40), and cell-cell interactions (EPHA2; ref. 41). The genes in these additional functional groups were important in subdividing the patients into subgroups with differing survival rates.
In conclusion, we report that the gene expression changes that commonly occur in nonmalignant HMEC grown in three-dimensional lrECM cultures provide gene expression signatures that effectively stratify patients into prognostic groups according to overall survival rates. Our 249-gene signature achieved a hazard ratio of 4.7, which is comparable with hazard ratios achieved by large-scale supervised breast cancer microarray studies. Our results underscore the relevance of three-dimensional lrECM cultures for studies of malignant transformation and suggest potentially valuable new biomarkers for further clinical evaluation.
Grant support: Grant DE-AC03 SF0098 (M.V. Fournier, P. Yaswen, and M.J. Bissell); U.S. Department of Energy, Office of Biological and Environmental Research Distinguished Fellowship Award, National Cancer Institute grant 2 R01 CA064786-09, and U.S. Department of Defense Breast Cancer Research Program Innovator Award grant BC012005 (M.J. Bissell); National Institute of Allergy and Infectious Diseases grant U19-AI057319 (K.J. Martin, K. Xhaja, and I. Bosch); and Postdoctoral Training Fellowship grant DOD BCRP DAMD17-00-1-0224 (P.A. Kenny).
We thank Martha Stampfer for kindly providing the 184 HMEC and her and Saira Mian for critical reading of the article and other members of our laboratories for fruitful discussions.