Hemoglobin electrophoresis and thalassemia gene detection
The hemoglobin electrophoretic patterns of the six cord blood samples used for RNA sequencing (RNA-Seq) analysis are shown in Supplementary Figure S1-2. The average HbF levels of the full-term female (FF) group and premature female (PF) group statistically differed (84.5% ± 8.1% and 57.8% ± 8.1%, respectively; P = 0.016). The gene detection results showed that none of the six newborns carried any of the three α-globin gene deletion mutations (--SEA, -α 3.7, or -α 4.2), three non-deletion mutations (Hb CS, Hb QS, or Hb WS), or 17 β-globin gene mutations most commonly seen in the Chinese population (Supplementary Figure S3).
Isolation Of Nrbcs
NRBCs were successfully isolated from each cord blood sample using magnetic bead technology based on the CD71 and CD45 cell surface markers. All isolated cells were analyzed by flow cytometry to determine purity. The average purity of all samples was > 92% (Fig. 1a, b).
The isolated cells were stained with Wright’s stain and observed under a microscope. Most of the cells in the field of vision were NRBCs; the NRBCs were scattered, and mostly comprised orthochromatic normoblasts with rare polychromatic erythroblasts, basophilic erythroblasts, and proerythroblasts (Fig. 1c). The enucleated red blood cells were further identified as reticulocytes by Brilliant cresyl blue staining (Fig. 1d). NRBCs and reticulocytes accounted for > 99% of the sorted cells.
Differentially Expressed Small Rnas
More than 23 million clean 50-bp single-end reads were generated for each sample. A Poisson distribution method was used to screen the differentially expressed small RNAs between samples, and NOIseq software was used to screen the differentially expressed small RNAs between the two groups of newborns. We detected 1,767 differentially expressed miRNAs between the FF and the PF group, including 1,033 miRNAs that were upregulated and 734 miRNAs that were downregulated in the FF group compared with the PF group (Fig. 2a).
With a fold change threshold of ≥ 2 and a deviation probability value of ≥ 0.6 as our screening criteria, we identified 79 differentially expressed miRNAs in the FF group compared with the PF group, including 53 that were upregulated and 26 that were downregulated. The volcano map (Fig. 2b) and hierarchical clustering heat map (Fig. 3) were constructed to illustrate the distribution of the differentially expressed miRNAs; some of the significantly differentially expressed miRNAs are marked on the volcano map. Several members of the let-7 miRNA family— let-7a-3p, let-7b-3p, let-7b-5p, let-7f-5p, and let-7 g-3p—were more highly expressed in the FF group, whereas miR-675-3p and miR-675-5p were expressed at significantly lower levels (Table 1).
Table 1
Differential expression of the let-7 family and miR-675
GENE ID | Expression in FF | Fold change | Probability |
hsa-let-7a-3p | Upregulated | 7.58 | 0.85 |
hsa-let-7b-3p | Upregulated | 22.41 | 0.95 |
hsa-let-7b-5p | Upregulated | 10.49 | 0.88 |
hsa-let-7d-3p | Upregulated | 2.51 | 0.56 |
hsa-let-7f-5p | Upregulated | 6.58 | 0.83 |
hsa-let-7f-1-3p | Upregulated | 1.24 | 0.27 |
hsa-let-7 g-3p | Upregulated | 14.30 | 0.94 |
hsa-let-7i-3p | Upregulated | 2.63 | 0.57 |
hsa-miR-675-3p | Downregulated | 38.12 | 0.88 |
hsa-miR-675-5p | Downregulated | 26.03 | 0.82 |
GO enrichment analysis was performed for the target genes of the miRNAs that were differentially expressed between the two groups. We identified 49 significant enrichment terms at the biological process level, including histone modification, covalent chromatin modification, negative regulation of protein phosphorylation, positive regulation of neuron projection development, cell growth, transforming growth factor beta receptor signaling pathway, protein methylation, stem cell differentiation, and Wnt signaling pathway. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was also performed for the target genes of the differentially expressed miRNAs; it revealed enrichment for involvement in two pathways (Q ≤ 0.05): morphine addiction and nicotine addiction (Supplementary Figure S5).
Differentially Expressed Mrna And Lncrna Analysis
A total of 114–164 million clean 150-bp paired-end reads were generated per sample. Between the FF and PF groups, there were 12,793 differentially expressed mRNAs (11,464 known mRNAs and 1,329 novel mRNAs), and 9,121 differentially expressed lncRNAs (8,417 known lncRNAs and 704 novel lncRNAs). The analysis parameters to determine differential expression with the DEGseq software were set to a fold change of ≥ 2 and a P-value of ≤ 0.001. We found 7,166 differentially expressed mRNAs (6,549 known mRNAs and 617 novel mRNAs) and 3,243 differentially expressed lncRNAs (3,009 known lncRNAs and 234 novel lncRNAs) that met these criteria (Supplementary Figure S6). Volcano maps illustrating the distribution of the differentially expressed mRNAs and lncRNAs are shown in Fig. 4.
Table 2
Differential expression of globin genes and related genes
GENE ID | Expression in FF | Fold change (FF/PF) | P-value | Q-value |
HBB | Upregulated | 2.78 | 0 | 0 |
HBD | Upregulated | 21.53 | 0 | 0 |
HBE1 | Downregulated | 26.45 | 0 | 0 |
HBG2 | Downregulated | 1.56 | 0 | 0 |
BCL11A | Upregulated | 11.26 | 7.53E-29 | 1.19E-28 |
ZBTB7A | Upregulated | 2.75 | 2.18E-117 | 1.21E-116 |
KLF1 | Upregulated | 1.43 | 9.32E-273 | 1.28E-271 |
SOX6 | Upregulated | 1.19 | 9.83E-16 | 1.03E-15 |
GATA1 | Upregulated | 1.36 | 1.70E-33 | 3.01E-33 |
HBS1L | Upregulated | 1.41 | 9.57E-81 | 3.76E-80 |
MYB | Upregulated | 4.81 | 2.01E-35 | 3.74E-35 |
NR2C2 | Upregulated | 2.09 | 1.21E-31 | 2.05E-31 |
LIN28B | Downregulated | 5.83 | 1.35E-02 | 5.73E-03 |
HMGA2 | Upregulated | 4.42 | 2.35E-04 | 1.27E-04 |
HMGB2 | Upregulated | 11.79 | 0 | 0 |
SIRT6 | Upregulated | 3.50 | 2.76E-16 | 2.97E-16 |
RUNX1 | Upregulated | 2.25 | 1.01E-03 | 5.09E-04 |
RUNX1-T1 | Upregulated | 4.28 | 8.39E-07 | 5.62E-07 |
JAZF1 | Upregulated | 3.07 | 3.86E-252 | 4.80E-251 |
CA1 | Upregulated | 9.78 | 0 | 0 |
Among the differentially expressed transcripts, HBB (encoding β-globin) expression was 2.78-fold higher, HBG2 (encoding Gγ- globin) expression was 1.56-fold lower, HBG1 (encoding Aγ- globin) expression was 1.47-fold lower, and HBD (encoding δ-globin) expression was 21.52-fold higher in the FF samples than in the PF samples (Table 2). The genes for the transcription factors BCL11A, ZBTB7A, KLF1, and SOX6, which are known to be involved in γ-globin repression, were 11.26-, 2.75-, 1.43-, and 1.19-fold higher, respectively, in the FF compared with the PF samples. Expression of LIN28B, which is associated with high HbF levels [11], was 5.83-fold lower in the FF samples. In addition, we found that H19 was expressed at much lower levels (85.28-fold) in the FF samples (Table 3).
Table 3
Differential expression of lncRNAs
GENE ID | Expression in FF | Fold change (FF/PF) | P-value | Q-value |
DANCR | Upregulated | 3.87 | 3.33E-24 | 4.64E-24 |
H19 | Downregulated | 85.28 | 0 | 0 |
RN7SK | Upregulated | 3.69 | 0 | 0 |
SBF2-AS1 | Downregulated | 11.87 | 8.66E-159 | 6.54E-158 |
GO enrichment analysis of the 6,549 known differentially expressed mRNAs, with a corrected P-value of ≤ 0.05 as the threshold, yielded 967 terms for enriched biological processes (Fig. 5a), including DNA replication, DNA conformation changes, and RNA splicing. KEGG pathway enrichment analysis revealed that the differentially expressed mRNAs were enriched in 57 pathways (Q ≤ 0.05) related to the spliceosome, RNA transport, cell cycle, systemic lupus erythematosus (SLE), the Fanconi anemia pathway, DNA replication, and other critical regulatory processes (Fig. 5b).
We also performed GO enrichment analysis of the target genes of the differentially expressed lncRNAs. For biological processes, 457 GO terms were enriched, such as DNA conformation changes, RNA splicing, DNA replication, and RNA transport (Fig. 6a). KEGG pathway enrichment analysis of the differentially expressed lncRNA target genes revealed enrichment of 11 pathways (Q ≤ 0.05), including the spliceosome, ubiquitin mediated proteolysis, SLE, the mRNA surveillance pathway, and chronic myeloid leukemia (Fig. 6b).
Quantitative Polymerase Chain Reaction Analysis
We collected 75 umbilical cord blood samples from which we isolated NRBCs to verify the differential expression of seven mRNAs and lncRNAs by quantitative polymerase chain reaction (qPCR). Compared with the preterm group, the full-term group had higher mRNA expression levels of BCL11A, HBB, and HMGB2 (3.24-, 3.52-, and 1.59-fold, respectively). The expression of H19 and LIN28B was 69.98- and 2.38-fold lower in the full-term group. The results were consistent with those of the RNA-SEq. We observed no significant differences in the expression of DANCR and SIRT6. In addition, we found that BCL11A expression significantly differed based on biological sex. Cord blood from preterm male neonates had 4.25-fold higher expression than cord blood from preterm female neonates, and cord blood from full-term male neonates had 2.64-fold higher expression than cord blood from full-term female neonates (Fig. 7).