|Home | About | Journals | Submit | Contact Us | Français|
Time course and dose-response studies were conducted in HL1-1 cells, a human liver cell line with stem cell–like characteristics, to assess the differential gene expression elicited by 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD) compared with other established models. Cells were treated with 0.001, 0.01, 0.1, 1, 10, or 100nM TCDD or dimethyl sulfoxide vehicle control for 12 h for the dose-response study, or with 10nM TCDD or vehicle for 1, 2, 4, 8, 12, 24, or 48 h for the time course study. Elicited changes were monitored using a human cDNA microarray with 6995 represented genes. Empirical Bayes analysis identified 144 genes differentially expressed at one or more time points following treatment. Most genes exhibited dose-dependent responses including CYP1A1, CYP1B1, ALDH1A3, and SLC7A5 genes. Comparative analysis of HL1-1 differential gene expression to human HepG2 data identified 74 genes with comparable temporal expression profiles including 12 putative primary responses. HL1-1–specific changes were related to lipid metabolism and immune responses, consistent with effects elicited in vivo. Furthermore, comparative analysis of HL1-1 cells with mouse Hepa1c1c7 hepatoma cell lines and C57BL/6 hepatic tissue identified 18 and 32 commonly regulated orthologous genes, respectively, with functions associated with signal transduction, transcriptional regulation, metabolism and transport. Although some common pathways are affected, the results suggest that TCDD elicits species- and model-specific gene expression profiles.
2,3,7,8-Tetrachlorodibenzo-p-dioxin (TCDD) is a ubiquitous environmental contaminant that elicits a broad spectrum of responses including endocrine disruption, immunotoxicity, hepatotoxicity, teratogenesis, and tumor promotion (Denison and Heath-Pagliuso, 1998). It is a multisite carcinogen and liver tumor promoter in rodents (Knerr and Schrenk, 2006), and classified as a human carcinogen by the International Agency for Research on Cancer (IARC) (McGregor et al., 1998). The effects are mediated by the aryl hydrocarbon receptor (AhR) (Okey, 2007), a ligand activated basic helix-loop-helix/PER-ARNT-SIM (bHLH/PAS) transcription factor. Upon ligand binding, the cytosolic AhR undergoes a conformational change, which leads to the dissociation of chaperone proteins and translocation to the nucleus where it heterodimerizes with the aryl hydrocarbon receptor translocator. The heterodimer complex binds to dioxin response elements (DREs) in the regulatory region of target genes to modulate gene expression. AhR-mediated changes in gene expression are believed to account for the toxicity of TCDD and related compounds.
In vivo and in vitro models have been used to investigate the mechanisms of action of TCDD and to assess its potential toxicity to humans and other environmentally relevant species. However, the relevance of data extrapolation from experimental models to target species has been questioned (Olson et al., 2000). In general, minimizing the extrapolation between experimental models and the relevant species is expected to more accurately assess the potential toxicity to the target of concern. Although cell lines model may reflect some in vivo responses, their utility for predicting in vivo toxicity is limited. Nevertheless, human based in vitro models are expected to more accurately indicate potential in vivo human toxicity. For example, human based models are preferred in drug development to investigate the mechanisms of action and assess potential toxicities of drug candidates (Abbott, 2003; MacGregor, 2003; MacGregor et al., 2001; Pritchard et al., 2003). Immortalized cell lines, primary cells and organ slices have been used as drug screening tools (Li, 2001) and to assess potential mechanisms, including toxicant-target interaction and dose/time dependent responses (MacGregor et al., 2001). However, although continuous cell lines are convenient, they are transformed or derived from diseased tissue, which may not accurately reflect normal tissue responses (Abbott, 2003). In contrast, the availability, use, cost, safety, and stability of primary cells and human tissue continues to be a factor.
Human stem cells provide an attractive alternative and potentially unlimited source of normal cells (Davila et al., 2004; Rohwedel et al., 2001). They are defined by their self-renewal and differentiation capabilities, and classified as either embryonic or adult stem cells based on their developmental status. Normal adult human stem cells isolated from liver tissue is an alternative model for toxicity studies, that may more closely mimic human tissue (Carere et al., 2002; Rohwedel et al., 2001; Rolletschek et al., 2004). Stem cells are also amenable to high throughput screening to rank and prioritize chemicals and drug candidates that warrant further investigation or development (Ahuja et al., 2007; Davila et al., 2004; Ek et al., 2007).
HL1-1 cells were derived from adult normal liver tissue. They exhibit high proliferation potential, express stem cell (Oct-4), and liver oval cell markers (AFP, vimentin, and Thy-1), and have the ability to differentiate into hepatocytes (Chang et al., 2004). In this study, HL1-1 cells were used to investigate the time- and dose-dependent gene expression changes elicited by TCDD. HL1-1 gene expression data were also compared with other in vitro (HepG2 and Hepa1c1c7 cells) and in vivo (C57BL/5 hepatic tissue) data sets treated with TCDD using comparable study designs and data analysis approaches. Although some common pathways are affected, the results suggest that TCDD elicits species- and model-specific gene expression profiles.
HL1-1 human liver cell line was derived from a normal healthy liver section obtained during the surgical resection of a male (age 49) with hemangioma. This cell line has been previously shown to possess liver stem cell characteristics (Chang et al., 2004; Tai et al., 2005) including, (1) high self-renewal ability, cumulative population doubling level (cpdl) about 50, (2) deficiency in gap-junctional intercellular communication, (3) expression of Oct 4, and liver stem cell markers, alpha-fetoprotein, vimentin, and Thy-1, and (4) the ability to differentiate into albumin producing cells. The morphology and proliferation potential of HL1-1 cells are shown in Supplementary Figure 1. The finite lifespan of HL1-1 cells and the nontumorigenic nature of HL1-1 cells immortalized by SV40 large T-antigen (no tumor developed in six mice injected with 4 × 106 cells at two sites per mouse, our unpublished results) indicate that HL1-1 cells are normal nontumorigenic cells. The cells at approximately passage 8 with cumulative population doubling level of about 40 were used for this study.
HL1-1 cells (Chang et al., 2004) were maintained in a low calcium (0.09mM) modified MCDB 153 medium (Keratinocyte-SFM, Invitrogen Corporation, Carlsbad, CA) supplemented with 2mM N-acetyl-l-cysteine, 0.2mM l-ascorbic acid 2-phosphate, and 5mM nicotinamide (referred to as K-NAC medium). The medium contains human recombinant epidermal growth factor (5 ng/ml), bovine pituitary extract (50 μg/ml), and 10% fetal bovine serum (Hyclone, Logan, UT). 1 × 106 cells were seeded into vented 25-cm2 cell culture flasks (Corning, Inc., Corning, NY) and incubated under standard conditions (5% CO2, 37°C). For time course studies, cells were treated with 10nM TCDD (S. Safe, Texas A&M University, College Station, TX) or DMSO (Sigma, St Louis, MO) and harvested at 1, 2, 4, 8, 12, 24, or 48 h (Supplementary Fig. 2). For dose-response studies HL1-1 cells were treated with DMSO (vehicle control) or TCDD (0.001, 0.01, 0.1, 1, 10, 100nM) for 12 h (Supplementary Fig. 2). For the cycloheximide (CHX) study, HL1-1 cells were divided into four treatment groups: (1) 10nM TCDD, (2) DMSO as a vehicle control, (3) pretreatment with 10 μg/ml CHX (Sigma-Aldrich, St Louis, MO) for 1 h prior to treatment with DMSO, or (4) pretreated with 10 μg/ml CHX for 1 h prior to treatment with 10nM TCDD, and then harvested at 4 and 12 h (Supplementary Fig. 2).
Cell lysates were prepared in RIPA buffer (1× PBS, 0.1% SDS, 1% IGEPAL CA-630, and 0.5% Na-Deoxycholate) containing a protease inhibitor cocktail (Roche Diagnostics, Mannheim, Germany). Protein concentrations were measured with a modified Lowry assay (Bio-Rad, DC Protein Assay, Hercules, CA). Cell lysates (20 μg) were separated on a 10% SDS-PAGE, transferred onto nitrocellulose membranes (Amersham Biosciences, Inc., Piscataway, NJ), and probed with anti-human AhR (N-19) antibody followed by horseradish peroxidase–conjugated secondary antibodies (Santa Cruz Biotechnology, Inc., Santa Cruz, CA). The blot was imaged by immunochemical staining and fluorescence detection on X-ray film by the enhanced chemiluminescence method (Amersham Biosciences, Inc.).
All treatment groups were harvested in 1.0 ml of Trizol Reagent (Invitrogen) and stored at −80°C. Total RNA was isolated according to the manufacturer's protocol, resuspended in RNA Storage Solution (Ambion, Inc., Austin, TX), and quantified spectrophotometrically (A260). The purity and integrity of each sample was assessed by the A260/A280 ratio and gel electrophoresis.
Gene expression changes were measured on custom human cDNA arrays containing 9684 features representing 6995 unique genes. For the time course study, an independent reference study design was used that included three replicates. Time-matched TCDD treated and vehicle samples were cohybridized at each time point (Supplementary Fig. 3A). Dose-dependent changes in gene expression were analyzed using a spoke design with three replicates at 12 h (Supplementary Fig. 3B). Dye swaps were also performed for both time course and dose-response studies to account for dye biases. For CHX studies, a 2X2 factorial design was used to facilitate appropriate statistical comparisons between all four treatment groups (Supplementary Fig. 3C) (Yang and Speed, 2002). If TCDD-elicited gene response were sustained or enhanced by CHX cotreatment, a direct effect independent of down stream translational activities was inferred and classified as a putative primary response. However, if the gene expression response was attenuated by CHX, it was assumed that the response was dependent on protein production(s), which were blocked by CHX, and therefore classified as putative secondary responses. Candidate genes that passed the first filter were further investigated by comparing CHX treatment alone to vehicle and cotreatment compared with CHX alone in order to exclude unclassified genes whose change in expression was affected by CHX.
PCR amplified cDNAs were robotically spotted onto epoxy-coated slides (Schott-Nexterion, Duryea, PA) by an Omnigrid arrayer (GeneMachines, San Carlos, CA) equipped with Chipmaker 3 pins in a CHP3 printhead (Telechem International, Inc., Sunnyvale, CA) at the Research Technology Support Facility, Michigan State University (http://www.genomics.msu.edu). Selected clones were obtained from EPAMAC, Research Genetics and Van Andel Research Institute. Detailed protocols for processing of microarrays including the labeling of the cDNA probe are available at http://dbzach.fst.msu.edu/interfaces/microarray.html.
Briefly, the Genisphere 900 3DNA Array Detection (Genisphere, Inc., Hatfield, PA) indirect incorporation kit was used to generate cDNA samples for hybridization according to manufacturer's protocol. After 20 h of cDNA hybridization, slides were washed and rehybridized with a Cy3:Cy5 (1:1) dendrimer mixture to indirectly incorporate dyes at the Cy3- and Cy5-dendrimer-tagged cDNA for 16 h. Slides were then scanned at 635 nm (Cy3) and 532 nm (Cy5) using a GenePix 4000B Array Scanner (Molecular Devices, Union City, CA). Images were analyzed for feature and background intensities using GenePix Pro 6.1 (Molecular Devices).
Data were normalized using a semiparametric approach (Eckel et al., 2005). Empirical Bayes analysis was used to calculate posterior probabilities (P1(t) value) of expression change on a per gene and time point or dose group basis using the model-based t-value (Eckel et al., 2004). The data were filtered using a P1(t) and fold change to obtain the most reproducible differentially expressed genes for initial analysis and interpretation. All raw and normalized data were stored in the toxicogenomic information management system (TIMS) dbZach which supports microarray data management, mining, visualization and knowledge management (Burgoon and Zacharewski, 2007; Burgoon et al., 2006). Expression changes that passed the criteria were analyzed by hierarchical clustering (GeneSpring 6.0, Agilent Technologies, Inc., Santa Clara, CA, and Multiexperiment Viewer [MeV] in TM4 software; Saeed et al., 2003) using uncentered Pearson correlation with average linkage. Normalization and empirical Bayes analysis were performed using SAS 9.1 (SAS Institute, Cary, NC) and R 2.0.1 (http://www.r-project.org). Dose-response analysis was performed using Graph Pad Prism 4.0 (GraphPad Software, San Diego, CA). Functional categorization of genes was performed using Database for Annotation, Visualization and Integrated Discovery analysis as a gene ontology tool (Dennis et al., 2003).
Quantitative real-time PCR (QRT-PCR) was performed for a selected number of genes to verify microarray data. Total RNA (1.5 μg) was reverse transcribed by SuperScript II according to the manufacturer's protocol (Invitrogen). cDNA products were then amplified with gene specific primers designed using Primer3 (Rozen and Skaletsky, 2000) and SYBR Green PCR reaction mixture (Applied Biosystems, Foster City, CA) on an Applied Biosystems PRISM 7500 Sequence Detection System. Input copy number was quantified using a standard curve of log copy number versus threshold cycle (Ct). The copy number of each sample was standardized to the geometric mean of β-actin and glyceraldehyde-3-phosphate dehydrogenase to control for differences in RNA loading, quality, and cDNA synthesis. For graphing purposes, the relative gene expression levels were scaled such that the expression level of the time-matched vehicle treated control group was equal to 1. Official names and abbreviations, forward and reverse primer sequences, and product length are listed in Supplementary Table 1.
Intergenic regions 10-kb upstream of RefSeq annotated transcription start sites (TSS) were obtained from the University of California Santa Cruz Genome Browser and deposited into dbZach (Burgoon and Zacharewski 2007). Core DRE sequences (5′-GCGTG-3′) were computationally identified, extended by the flanking seven base pairs and scored using a position weight matrix (PWM) application (Sun et al., 2004). The PWM was developed using experimentally verified DREs; DRE sequences are from genome sequence. Putative DREs are those with a matrix similarity score greater than 0.80 (Boverhof et al., 2006; Kopec et al., 2008; N'Jai et al., 2008).
AhR mRNA was detected in HL1-1 cells, C57BL/6 hepatic tissue, human HepG2 and mouse Hepa1c1c7 cells (Fig. 1A). AhR protein expression was confirmed by Western blot (Fig. 1B). Treatment of HL1-1 cells with TCDD resulted in the dose-dependent induction (> 600-fold) of CYP1A1 mRNA with an EC50 of 8.30nM (Fig. 2), which is less sensitive compared with other human in vitro models (HepG2: 0.68nM; Zeiger et al., 2001), HepG2: 0.2nM and fresh hepatocytes: 0.14nM; Silkworth et al., 2005). There was no indication of TCDD treatment mediated differences with respect to growth or physical appearance. Collectively, these results confirm the functionality and responsiveness of HL1-1 AhR to TCDD.
Temporal and dose-dependent changes in gene expression were determined using custom human cDNA microarrays consisting of 9684 features representing 6995 unique genes. Differentially expressed genes were identified using P1(t) values greater than 0.9999 and |fold change| > 1.5 as the criteria. A total of 113 unique genes (72 induced; 41 repressed) exhibited dose-dependent regulation with a majority exhibiting differential expression between 1 and 100nM TCDD (Fig. 3A). EC50s ranged from 0.18nM for ACSL3 to 37.3nM for MOBP, although induction and repression ranged from 40- to 2-fold for genes CYP1B1 and GALNT1, respectively.
In the time course study, 144 unique genes were differentially regulated by TCDD at one or more time points (Fig. 3B). Hierarchical clustering identified induced and repressed expression with three distinct, temporal clusters of early (2–4 h), mid (8–12 h), and late (24–48 h) responses (Supplementary Fig. 4). Early and mid–time point groups showed vast differences in their expression profiles, illustrating a temporal cascade of responses.
A subset of responsive genes including CYP1B1, ALDH1A3, and SLC7A5 was verified by QRT-PCR (Fig. 4). There was good agreement between the microarray and QRT-PCR data for the temporal expression profiles with comparable EC50 values.
HL1-1 cells were pretreated with CHX to inhibit de novo protein synthesis in order to identify putative primary gene expression responses mediated by the AhR. Following treatment 78 and 203 differentially expressed genes (P1(t) > 0.9999 and |fold change| > 1.5) were identified at 4 and 12 h, respectively (Fig. 5). These genes were classified into putative primary, secondary or unclassified groups. Forty-seven genes at 4 h and 53 genes at 12 h were putatively determined to be primary responses. A total of 78 unique genes were identified as putative primary responses with 22 genes in common at both time points.
Putative primary responses were associated with xenobiotic and lipid metabolism, cell cycle regulation, transcription regulation, transport and signal transduction (Table 1). This included the prototypical xenobiotic metabolism related genes such as cytochrome P450s and aldehyde dehydrogenases which were highly induced. Other primary responses included genes associated with lipid metabolism (APOM, PLD3, ST8SIA1, ACSL3), cell cycle regulation (CDCA5, KANK1, FHIT, CTGF), and transcriptional regulation (MXD3, DEAF1, SERTAD2). The regulation of transcription factors and subsequent changes in gene expression is a hallmark of TCDD action (Reymann and Borlak, 2006). Immune response genes (IL1A, IL1B, CD8A), signal transduction-related genes (MAPK7, PRKCB), and transporters (SLC2A1, SLC7A5, MTCH2) were also identified as putative primary responses. Computational analysis of the regulatory region of these putative primary responses revealed that 71 of the 78 genes had one or more putative DREs with a matrix similarity score greater than > 0.8 (Table 1).
Temporal changes in gene expression elicited by TCDD in HL1-1 cells were compared with intralaboratory time course studies conducted in other model systems. 251 HL1-1 genes were identified as differentially expressed in the time course study using relaxed filtering criteria (P1(t) > 0.999 and |fold change| > 1.3) to include responses on the margins. Using the same relaxed criteria and the same study design, cDNA microarray, and data analysis strategy, 1057 HepG2 genes were identified as differentially expressed at one or more time points following treatment with 10nM TCDD (Dere et al., manuscript in preparation). Only 74 common genes were differentially regulated in HL1-1 and HepG2 cells by TCDD (Table 2, Supplementary Fig. 5A). Of these, 55 exhibited similar temporal expression patterns (38 induced; 17 repressed; Table 2, Supplementary Fig. 5B), and 12 were classified as putative primary responses. Of the 19 genes exhibiting divergent regulation (seven induced in HL1-1 but repressed in HepG2; 12 repressed in HL1-1 but induced in HepG2 (Table 2, Supplementary Fig. 5C), none were putative primary responses in the HL1-1 based on the CHX study.
HL1-1 specific responses were associated with immune and lipid/xenobiotic metabolic processes (Supplementary Fig. 5D), although HepG2-specific responses were involved in the cytoskeleton, calcium signaling, and lipid metabolism. Other functional associations included transport, cell cycle, signal transduction, and transcriptional regulation.
Comparison of HL1-1 differential gene expression to mouse Hepa1c1c7 cell (Dere et al., 2006) and C57BL/6 hepatic tissue (Boverhof et al., 2005) were also examined using relaxed filtering criteria (P1(t) > 0.999 and |fold change| > 1.3) to include those responses approaching the initial criteria (P1(t) > 0.9999, |fold change| > 1.5). All of these studies used comparable study designs, cDNA microarray platforms and data analysis strategies.
A total of 5505 orthologous genes, defined by HomoloGene (http://www.ncbi.nlm.nih.gov/HomoloGene/), were represented across the human and mouse cDNA arrays (Fig. 6A). Comparison of HL1-1 to C57BL/6 liver tissue, and Hepa1c1c7 cells identified only 32 and 18 common differentially expressed genes, respectively (Fig. 6B, Table 3). However, not all common differentially expressed genes exhibited the same expression pattern. For example, of the 18 genes differentially expressed in both HL1-1 cells and Hepa1c1c7 cells, 10 exhibited the same pattern (six induced and four repressed genes), although eight were divergently regulated (six genes induced in HL1-1 cells were repressed in Hepa1c1c7; two genes repressed in HL1-1 were induced in Hepa1c1c7). Similar analyses were conducted between HL1-1 and C57BL/6 hepatic tissue differential gene expression data sets (Fig. 6B, Table 4). Across all four models, only three genes (IRF1, SLC12A7, and ID3) were differentially expressed with only one gene (IRF1) exhibiting the same expression pattern.
Functional annotations of the common TCDD-regulated genes were associated with cell cycle regulation, and development. Several collagenases were differentially expressed in HL1-1 and mouse Hepa1c1c7 hepatoma cells, although some responses (e.g., CDC25B) were divergently regulated suggesting species-specific effects (Table 3) (Boverhof et al., 2006, Sun et al., 2004). Comparison of the functional annotation of differentially expressed genes in HL1-1 cell and hepatic mouse tissue identified transcriptional regulation, signal transduction, metabolism, apoptosis and transport as being commonly regulated by TCDD (Table 4). In vitro and in vivo comparisons have previously reported that some genes were divergently regulated (Dere et al., 2006). For example, the transcription regulation-related genes CEBPZ and ID3, and the signal transduction–related genes DUSP6, ERBB3, and MTMR7 were divergently regulated in HL1-1 cell compared with hepatic mouse tissue.
In this study, AhR-mediated gene expression in HL1-1 human liver stem cells was compared with other in vitro and in vivo TCDD data sets. The normalcy and self-renewal properties of HL1-1 cells provide a unique, in vitro system that may more accurately reflect in vivo human responses to TCDD. As with other models, HL1-1 cells express a functional AhR, as evident by western analysis and the induction of several AhR gene battery members, including CYP1A1 and aldehyde dehydrogenases. Although comparable functional pathways are affected in each model, further examination indicates that different genes within common functions were differentially regulated. Furthermore, there were several examples of orthologs exhibiting divergent regulation (e.g., induced in one model but repressed in another). Interestingly, functional categorization of TCDD-elicited differential gene expression is consistent with reported species-specific responses (Boutros et al., 2004, 2006) suggesting that HL1-1 cells may more accurately reflect in vivo human responses to TCDD.
Human HepG2 cells are a popular model for investigating hepatotoxicity. Comparison of HepG2 and HL1-1 gene expression profiles revealed the coregulation of several putative primary responses including CYP1A1, PRKCB, PPAN, SERTAD2, SLC2A1, and SLC7A5. In contrast, several divergently regulated genes not classified as primary responses were expressed at later time points suggesting that they may be secondary effects of TCDD. In addition, a number of cell specific responses were observed. For example, HL1-1–specific responses included xenobiotic and lipid metabolism and immune responses that correlated with in vivo mice effects (Boverhof et al., 2005), although HepG2 exhibited cytoskeleton and calcium signaling responses related to cell-adhesion, tumor cell motility and tumor promotion (Monteiro et al., 2008; Puga et al., 2000; Tannheimer et al., 1997). This may be due to differences in basal expression levels or other unique cell characteristics (Harris et al., 2004). For example, primary rat hepatocytes cultured on standard collagen had basal gene expression levels more comparable to whole liver rather than rat hepatoma cells (Boess et al., 2003).
Comparisons of HL1-1 cells to other in vitro (i.e., Hepa1c1c7 cells) and in vivo (i.e., hepatic tissue from C57BL/6 mice) models identified several common TCDD responses, as well as model-specific differential gene expression. Although the structure, function, and mechanism of action of the AhR is highly conserved (McGregor et al., 1998), comparative toxicogenomic and computational DRE search studies suggest that TCDD-elicited gene expression profiles may be species-specific. Computational analysis of the regulatory regions of orthologs using a PWM suggests that DREs are not conserved between humans, mice, and rats (Sun et al., 2004). Moreover, in vivo rat versus mouse (Boverhof et al. 2006), and in vitro (HepG2 vs. Hepa1c1c7 vs. H4IIE) (Dere et al., manuscript in preparation) studies indicate that the hepatic gene expression profiles are species specific, despite the conserved induction of xenobiotic metabolizing genes. Moreover, identified putative primary responses differ between species (Dere et al., manuscript in preparation). Nevertheless, TCDD does affect common pathways across species and models, but appears to do so by regulating the expression of nonorthologous genes within those pathways.
Collectively, these studies not only demonstrate the utility of HL1-1 cells but also the limitations of extrapolating from in vitro models to in vivo effects. Although in vitro models have the advantage of reducing the complexity of a tissue response to allow a more focused analysis of the effects of TCDD on a specific cell type, it does not replicate other interactions that may be important in eliciting the toxic responses observed in vivo (Dere et al., 2006). However, in addition to being a normal human cell which may more accurately reflect human responses relative to rodent models, HL1-1 cells are also stem-like which may be novel targets of toxicity (Davila et al., 2004).
National Institute of General Medical Sciences (GM075838); Michigan Agriculture Experiment Station to T.R.Z.
The authors gratefully acknowledge the collaboration of Drs Jin-Lain Tsai and Kung-Kai Kuo of Kaohsiung Medical University, Kaohsiung, Taiwan, with C.C. Chang in developing the HL1-1 cells used in this study. The authors also would like to gratefully thank Michelle Manente, Joshua Kwekel, Anna Kopec, and Heekyong Bae for helpful discussion and critical review of this manuscript.