Comparison of gene expression profiles between human erythroid cells derived from fetal liver and adult peripheral blood

Isolation of primary fetal and adult CD34⁺ cells

Human tissue samples were obtained from donors with written informed consent for research use in accordance with the Declaration of Helsinki and approved by Committee on human rights related to research involving human subjects (Faculty of Medicine Ramathibodi Hospital, Mahidol University, Bangkok, Thailand (IRB no. MURA2013/363)). Adult and fetal CD34⁺ cells were isolated from adult healthy donor peripheral blood and fetal liver, respectively. Fetal livers were surgically obtained from five fetuses therapeutically aborted after 13–21 weeks of amenorrhoea. Fetal liver tissue was dissected within 12 hours after collection. CD34⁺ cells derived from fetal liver were prepared using a procedure modified from Bility et al. (2012). Briefly, fetal liver was manually minced into small pieces and homogenized by automated tissue dissociation (gentleMACS Dissociator, Miltenyi Biotech Inc., Auburn, CA, USA) in digestion buffer containing DNaseI grade II, from bovine pancreas (Sigma-Aldrich, St. Louis, MO, USA) and collagenase type I (GIBCO, Grand Island, NY, USA). A single-cell suspension was obtained using a 70 µm cell strainer. Mononuclear cells were isolated from the fetal liver cell suspension and adult healthy peripheral blood by centrifugation on a density gradient of 1.077 g/mL (Lymphoprep, Axis-Shield PoC AS, Oslo, Norway) and subsequently selected for CD34⁺ cells using the positive immune-magnetic selection method (CD34 MicroBead Kit, Miltenyi Biotech Inc.) according to the manufacturer’s recommendations.

Erythroid differentiation of CD34⁺ cells

Primary CD34⁺ cells were cultured by a two-phase liquid culture system. Cells were cultured for 4 days in phase I medium consisting of Iscove’s Modified Dulbecco’s Medium (IMDM; GIBCO, Grand Island, NY, USA) supplemented with 20% of fetal bovine serum (FBS; Sigma-Aldrich), 300 ng/mL holo-transferrin (holo-TF; PromoCell, Heidelberg, Germany), 25 ng/mL interleukin-3 (IL-3; Cell Signaling Technologies, Beverly, MA, USA), 50 ng/mL human stem cell factor (SCF; Cell Signaling Technologies), and 2 units/mL human recombinant erythropoietin (EPO; CILAG GmbH, Zug, Switzerland). After 4 days, non-adherent cells were collected and reseeded in phase II medium consisting of IMDM supplemented with 20% FBS, 300 ng/mL holo-TF, and 5 units/mL EPO. The culture was maintained under an atmosphere of 5% CO₂ at 37 °C.

Erythroblast differentiation was monitored throughout culture by flow cytometry analysis using a FACSCalibur flow cytometer (BD Biosciences, San Jose, CA, USA), in which cells were immuno-stained with phycoerythrin-conjugated anti-transferrin receptor (CD71-PE; BD Biosciences Pharmingen, San Diego, CA, USA) and allophycocyanin-conjugated anti-glycophorin A (GPA-APC; BioLegend, San Diego, CA, USA) antibodies. In addition, cell maturation was analyzed by light microscopy of May-Grünwald-Giemsa-stained cytospin preparations. To perform RNA analysis, the culture was sampled at day 8 (4 d in phase I and 4 d in phase II). The culture was continued in phase II medium for another 6 d. At the end of the culturing period at day 14 (4 d in phase I and 10 d in phase II), cells were harvested for analysis of hemoglobin content by high performance liquid chromatography according to the manufacturer’s recommendations (HPLC, VARIANT II hemoglobin testing system; Bio-Rad Laboratories, Hercules, CA, USA).

RNA isolation and microarray hybridization

Cells harvested at day 8 of culture were stained for CD71 and GPA cell surface markers. Erythroid populations comprising cells with high CD71 and GPA expression were FACS sorted using a BD FACSAria™ III cell sorter (BD Biosciences). Additionally, cells were sorted to remove non-viable cells and debris using forward scatter/side scatter gating. Total RNA was isolated using the RNeasy Plus Micro kit (Qiagen, Mountain View, CA, USA). RNA integrity was determined with the Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA); the RNA integrity number for all samples was greater than 8.5. Biotinylated, fragmented cRNA was synthesized from 100 ng of each total RNA sample using the Affymetrix GeneChip^® WT PLUS Reagent kit (Affymetrix, Santa Clara, CA, USA), and hybridized to Affymetrix Human Gene 2.0 ST arrays (HuGene-2_0-st, Affymetrix) following the manufacturer’s protocol. Arrays were then washed and stained using the FS450_0002 fluidics protocol and scanned using an Affymetrix 3000 7G scanner. The scanned images were inspected for hybridization efficiency and CEL files generated from AGCC (GeneChip Command Console Software) were imported into Expression Console (EC) 1.3 software for reporting array quality control metrics, including perfect match mean (PM_Mean), background mean (Bgd_Mean), positive and negative probes (POS vs NEG AUC), bacterial spike controls and Poly-A controls.

Reverse-transcription quantitative PCR (RT-qPCR)

Total RNA samples as described above for microarray experiments were tested by RT-qPCR. First-strand cDNA was synthesized from 1 µg of RNA sample using SuperScript III reverse transcriptase (Invitrogen) following treatment with 1 unit of DNase I (Thermo Fisher Scientific Baltics UAB, Vilnius, Lithuania) according to the manufacturer’s protocol. RT-qPCR was performed using an ABI PRISM 7700 sequence detection system instrument and software (Applied Biosystems, Foster City, CA, USA). RT-qPCR assays of eight globin transcripts were performed using SYBR^® Select Master Mix (Applied Biosystems) with primers and PCR conditions described previously (Tangprasittipap et al., 2015). Absolute globin gene transcript copy numbers were calculated by comparison with standard curves generated from a plasmid DNA encoding each globin template. The transcript levels of other genes identified by microarray as differentially expressed between FL- and AB-derived erythroblasts were measured using the Assays-on-Demand Gene Expression Products (Applied Biosystems) according to the manufacturer’s instructions. The PCR primers are listed in Table S1. Expression data for each transcript were normalized to the reference transcript 18S ribosomal RNA. The relative fold changes were analyzed by the 2^−ΔΔCt method (Livak & Schmittgen, 2001).

Microarray data analysis

Statistical analysis of microarray data was performed using R version 3.4.1 (R Core Team, 2017). The Affymetrix raw data files (.CEL files) were preprocessed using the robust multi-array average (RMA) method (Irizarry et al., 2003) implemented in the oligo package version 1.40.2 for background correction, log-transformation, and quantile normalization. Differential gene expression between FL- and AB-derived erythroblasts was determined using moderated t-statistics (Smyth, 2004) implemented in the Limma package version 3.32.10. The resulting p-values were further adjusted for multiple testing to minimize the number of false positives in the study using the Benjamini & Hochberg False Discovery Rate (FDR) method (Benjamini & Hochberg, 1995). Differentially expressed transcripts were defined as transcripts with an absolute fold change of at least 1.50 and the adjusted p-value less than 0.05. Hierarchical clustering (single linkage) of these transcripts was employed to group genes with similar expression patterns across samples using MeV version 4.9 (Saeed et al., 2003). The microarray data are available from the NCBI Gene Expression Omnibus database under accession number GSE109186. Microarray data of culture day 7 published previously (Xu et al., 2012) were downloaded from the GEO database and analyzed as described above. Meta-analysis of the two studies combined into one dataset was performed by the Rank Product method using the R package RankProd version 3.2.0 (Hong et al., 2006) based on default parameters. Only genes represented on both microarray platforms (Xu et al. study, Affymetrix human genome U133 Plus 2.0 array; our study, Affymetrix Human Gene 2.0 ST Array) were considered for meta-analysis. In each study, for multi-probe genes, the probe with the highest interquartile range (IQR) was selected to represent the gene to be used in meta-analysis as recommended in (Gentleman, 2005; Walsh et al., 2015).

Functional and pathway analyses

The functions and the associated pathways of all differentially expressed genes were analyzed according to information from the Gene Ontology (GO) database (Ashburner et al., 2000; Gene Ontology Consortium, 2015) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (Kanehisa & Goto, 2000). Functional enrichment analysis was performed using the ToppFun feature in the ToppGene Suite (Chen et al., 2009) using the cutoff of p-value = 0.05 with FDR correction. Lists of genes with significantly higher expression in FL- or AB-derived erythroblasts were used as input for ToppFun to identify pathways and GO terms over-represented among genes in each list. The selected pathways were displayed using Pathview (Luo & Brouwer, 2013; Luo et al., 2017). Significant GO terms were summarised using the REVIGO tool (Supek et al., 2011) with default settings.

Gene regulatory network construction

To describe the transcriptional regulation of genes co-upregulated either in FL or AB-derived erythroblasts, transcription factors (TFs) that target these genes were identified and the gene regulatory network displaying putative interactions among TFs and their direct target genes was built using iRegulon (Janky et al., 2014). Specifically, regulatory sequences in the vicinity of each gene in the set were analyzed to detect the enriched TF binding motifs or ChIP-seq peaks from databases. Then, these signals were associated with candidate transcription factors and their optimal target sets were determined. The regulatory search regions were set to be 20 kb flanking the transcription start site (from −10 kb to +10 kb) and the enrichment score threshold was 3.0. Other parameters were set as default. The resulting networks were displayed using cytoscape version 3.6 (Cline et al., 2007).

In vitro culture of primary human fetal and adult erythroid progenitors

In vitro differentiated erythroid progenitors were grown from FL- and AB-derived CD34⁺ cells. A greater amplification of erythroid cells was observed from FL than AB cells in our culture system (Fig. 1A). After day 8 of culture, the majority of cells had differentiated into erythroid cells, as shown by the flow cytometric gating of CD71⁺/GPA⁺ cells (Fig. 1B) and morphology (Fig. 1C). The distinctive properties of FL- and AB-derived differentiated erythroblasts with respect to globin expression were assessed by quantification of absolute levels of globin transcripts on day 8 cultured cells (Table 1). The FL-derived erythroblasts demonstrated markedly elevated expression of the fetal γ-globin gene, together with elevated expression of embryonic globin genes (ε- and ζ-globin) compared with AB-derived erythoblasts. In contrast, the levels of adult β-globin and δ-globin transcripts were markedly higher in AB-derived erythroblasts. Transcript levels of the α-, µ-, and θ-globin genes were similar in both cell types. Hemoglobin analysis of day 14 cultured cells (Fig. S1) confirmed the predominant HbF in FL-derived erythroblasts, with a much lower HbF level in AB-derived erythroblasts (FL = 95.4 ± 1.5% vs. AB = 4.4 ± 0.2%). Neither HbZ nor HbE1 embryonic globins were detected by HPLC analysis. These data show that our in vitro culture system is a suitable model for comparing gene expression profiles during human erythroid development at fetal and adult stages.

Comparative gene expression analysis of erythroblasts derived from fetal liver and adult blood

We performed global gene expression profiling to characterize stage-specific gene expression patterns during human erythroid development. Gene expression profiles of FL- and AB-derived erythroblasts were determined from five independent samples as biological replicates by oligonucleotide microarray experiments. 1391 genes (1,414 transcripts) expressed higher in FL compared with AB erythroblasts, whereas 329 genes (340 transcripts) showed the opposite pattern (Fig. 2, Table S2). With the exception of globin genes and genes with unknown function, the neuropilin (NRP) and tolloid (TLL)-like 2 gene (NETO2) is the most significant FL up-regulated gene (adjusted p-value = 1.1E−06, fold change = 9.53) while tenascin XB (TNXB) is the most significant AB up-regulated gene (adjusted p-value = 1E−05, fold change = 2.98). When non-globin genes with known function were ranked by fold change instead, secreted phosphoprotein 1 (SPP1) was the top FL up-regulated gene (adjusted p-value = 2.61E−03, fold change = 14.96) while carbonic anhydrase I (CA1) was the top AB up-regulated gene (adjusted p-value = 1.36E−04, fold change = 10) (Table S2).

The expression patterns of the globin gene family were examined as positive controls to validate the microarray platform for global comparison of FL- and AB-derived erythroblast transcriptomes. The β-globin (HBB) gene was identified as one of the top-ten genes with significantly higher expression in AB compared with FL erythroblasts (fold change = 34.33, adjusted p-value = 4.25E−05) (Table S2), in agreement with RT-qPCR data. In addition, the expression patterns of other globin genes are concordant with RT-qPCR data, including the γ-globin gene (HBG1) which showed a significantly higher expression level in the FL group (fold change = 2.89, adjusted p-value = 6.29E-03). The ζ-globin (HBZ) and ε-globin (HBE1) genes also had significantly higher expression in the FL group (Table S2).

Functional analysis of differentially expressed genes

Toppfun functional enrichment analysis was performed on pathways and gene ontology (GO) terms represented among annotated functions of genes with significantly different expression between FL- and AB-derived erythroblasts in order to understand the biological significance of these genes (Table S3). More pathways were significantly associated with genes up-regulated in FL than AB (239 and seven pathway terms from various annotation data sources, respectively) (Table S3). In addition, more GO terms were revealed for FL than AB up-regulated genes, e.g., 457 and 28 enriched GO biological process terms, respectively. Overall, the enriched GO biological process terms of FL up-regulated genes summarized by REVIGO (Supek et al., 2011) indicate gene functions mainly associated with ncRNA metabolism, ribonucleoprotein complex biogenesis, cell cycle, and RNA localization, whereas the enriched terms among AB up-regulated genes indicate gene functions associated with regulation of immune response, transport, and purine deoxyribonucleoside metabolism. The top biological process GO terms were ribonucleoprotein complex biogenesis and one-carbon compound transport for FL and AB up-regulated genes, respectively. The top molecular function GO terms were RNA binding and IgE binding for FL and AB up-regulated genes, respectively. The top cellular component GO terms were ribonucleoprotein complex and side of membrane for FL and AB up-regulated genes, respectively.

Transcriptional regulation networks of differentially expressed genes

Because the switching of globin gene expression between fetal and adult stages is a transcriptionally controlled process, we identified putative transcription factor (TF) binding sites in the vicinity of genes identified by microarray analysis as differentially expressed between FL- and AB-derived erythroblasts (Table 2). Interactions among TFs and their target sites were revealed, in which TF-gene network connections showed largely common target genes among different TFs (Fig. 3). Some of these enriched TFs were also expressed differently between FL and AB, with expression profiles consistent with their target genes, i.e., TFs and their targets were up-regulated in the same cell type. FL up-regulated TFs included MYC, MAX, POLR3A and GABPA, whereas JUND was AB up-regulated. Among the FL up-regulated TFs, MYC was the top regulator with 793 potential targets among FL up-regulated genes.

Meta-analysis with previous microarray study

We reanalyzed the microarray data from Xu et al. (2012), in which a similar comparative study of FL- and AB-derived erythroblasts was performed. We did not include the data reported in Yang et al. (2013) for reanalysis or comparison, since only one sample of each cell type was studied on a different platform (RNA-seq). Moreover, studied cells in the Yang et al. were in a later stage of erythroblast maturation, which could bias the comparison of gene expression profiles. From the Xu et al. data, 420 genes (533 transcripts) were found to be FL up-regulated whereas 582 genes (790 transcripts) were AB up-regulated. This is fewer than reported in the original study (1,057 and 916 up-regulated genes in FL and AB cells, respectively; Table S1 of Xu et al.), suggesting that our analysis method is more conservative. A total of 18,965 genes were identified with features represented on both microarray platforms (Xu et al., Affymetrix human genome U133 Plus 2.0 array; our study, Affymetrix Human Gene 2.0 ST Array). A number of significant genes were identified only in each individual study (237 genes in ours and 50 genes in the Xu et al. study) owing to the presence of different features unique to each microarray platform (Table S2).

From the common genes represented on both microarray platforms, a meta-analysis of the combined dataset from both studies was performed. A total of 599 FL up-regulated and 284 AB up-regulated genes, respectively, were identified by meta-analysis (Table S4). These genes were compared with significant genes from each study dataset analyzed separately. 102 and 44 genes were identified as FL up-regulated or AB up-regulated, respectively in both studies, irrespective of whether the datasets were analyzed separately or combined for meta-analysis (Fig. 4). The combined dataset showed 60 FL up-regulated and 92 AB up-regulated genes in meta-analysis, which were not significant for either dataset analyzed separately (Fig. 4). The finding of 152 genes with significant difference in expression only in meta-analysis indicates that combining datasets gives increased power to detect genes with reproducible changes in expression between the two studies. The heat map representations of meta-analysis significant genes further illustrate the increased power of meta-analysis to detect genes with reproducible patterns (Fig. 5). Despite the increased power to detect reproducible patterns of gene expression by meta-analysis, many genes showed significant changes in expression only in one microarray study dataset analyzed separately, and not the other or the combined dataset in meta-analysis (Fig. 4; Fig. S2).

Although reproducible patterns of gene expression were identified by meta-analysis, the discrepancies between differentially expressed genes identified from our data and those from the Xu et al. study data led us to perform functional enrichment analysis of the combined dataset to determine which pathways/biological processes are reproducible between studies (Table S5). More pathways were significantly associated with genes up-regulated in FL than AB (76 and 29 pathway terms, respectively), although more AB pathways were identified in the combined dataset than from our data analyzed separately (Table S5). Similarly, more GO terms were revealed for FL than AB up-regulated genes, e.g., 268 and 34 enriched GO biological process terms, respectively. Among FL up-regulated genes in the combined dataset, the summary of GO biological process terms was similar to the analysis results of our data (ncRNA metabolism and ribonucleoprotein complex biogenesis) with some additional terms (response to nitrogen compound and glucose catabolism). For AB up-regulated genes, the summarized GO biological process terms from the combined dataset were quite different from our dataset (regulation of cell differentiation, response to oxidative stress) although transport process was in common. The top enriched biological process and molecular function GO terms were the same as our data analyzed separately (ribonucleoprotein complex biogenesis and RNA binding respectively), whereas mitochondrion was the top cellular component GO term. Among AB up-regulated genes in the combined dataset, the top enriched GO terms for each category were different from our data analyzed separately (positive regulation of cell differentiation, ammonium transmembrane transporter activity and nucleosome).

Analysis of TFs targeting meta-analysis significant genes showed some common TFs with the ones identified in significant genes from our microarray data analyzed separately. The common TFs targeting FL up-regulated genes were MYC, MAX, ZBTB33, ZNF143, and POLR3A, whereas the common TFs for AB up-regulated genes were CEBPA, CEBPD, IKZF2, POU2F2, and RELA (Table 2, Table S6 , Fig. S3). Among the FL up-regulated TFs, MYC and its TF partner MAX (Grandori et al., 2000) were the top regulators with 361 and 370 potential target genes with FL up-regulated expression.

Validation of meta-analysis differentially expressed genes by reverse-transcription qPCR

To validate the findings from microarray experiments, the same FL and AB RNA samples as used in microarray were recruited for RT-qPCR. RT-qPCR experiments were performed for thirteen genes with significant differential expression by microarray in our dataset analyzed separately and in meta-analysis of the combined dataset. Five FL up-regulated genes (IGF2BP3, BCAT1, CHD7, NETO2, MYC) and seven AB up-regulated genes (CA1, LGALS3, FAS, MAP3K5, TNXB, JUND, JAZF1) were confirmed (Table 3). Differential expression of the LIN28B gene could not be assessed by RT-qPCR since the expression in adult erythroblasts was lower than the limit of detection (Ct greater than 35).

Common pathways associated with differentially expressed genes

826 and 48 novel genes were identified only from data generated in this study to be FL up- and AB up-regulated, respectively (Fig. 4, Table S2). We explored the potential biological significance of these novel genes by inspection of KEGG pathways. We reasoned that novel genes with functions relevant to erythropoietic transition from the fetal to adult state would be connected in gene function pathways described previously in erythropoiesis. The mitogen-activated protein kinase (MAPK) and the phosphatidyl inositol 3 kinase (PI3K)-Akt pathways are activated in erythropoiesis (Hattangadi et al., 2011). Thirteen differentially expressed genes (CD14, GRB2, MYC, MAX, PRKCA, TGFBR1, MAPKAPK3, MAP3K5, MAP2K2, FAS, VEGFA, DDIT3, FLNA) were identified in the MAPK pathway from meta-analysis together with eleven novel genes from our data (CHUK, CRK, IGF1R, JUND, KRAS, NRAS, MAPK14, NLK, PPP3CA, MAPKAPK5, FGF5) (Fig. 6). Twenty differentially expressed genes (TNXB, COL4A5, SPP1, THBS1, VEGFA, GRB2, HSP90AB1, HSP90B1, ITGA5, ITGAV, MYC, PHLPP2, PIK3R1, PRKCA, PPP2R5B, THEM4, YWHAB, BCL2L1, FOXO3, MAP2K2) were identified in the PI3K-Akt pathway from meta-analysis together with twelve novel genes from our data (CDC37, CDK4, CHUK, CREB3L2, GSK3B, IFNAR1, IGF1R, KRAS, NRAS, RPS6KB1, SGK1, FGF5) (Fig. 7). The MAPK and PI3-Akt pathways overlap with other pathways, including the cell-cycle and cancer. Other differentially expressed genes from meta-analysis together with novel genes identified only from our data mapped to these pathways (Figs. S4 and S5).