Transcriptomic profile of adverse neurodevelopmental outcomes after neonatal encephalopathy

A rapid and early diagnostic test to identify the encephalopathic babies at risk of adverse outcome may accelerate the development of neuroprotectants. We examined if a whole blood transcriptomic signature measured soon after birth, predicts adverse neurodevelopmental outcome eighteen months after neonatal encephalopathy. We performed next generation sequencing on whole blood ribonucleic acid obtained within six hours of birth from the first 47 encephalopathic babies recruited to the Hypothermia for Encephalopathy in Low and middle-income countries (HELIX) trial. Two infants with blood culture positive sepsis were excluded, and the data from remaining 45 were analysed. A total of 855 genes were significantly differentially expressed between the good and adverse outcome groups, of which RGS1 and SMC4 were the most significant. Biological pathway analysis adjusted for gender, trial randomisation allocation (cooling therapy versus usual care) and estimated blood leukocyte proportions revealed over-representation of genes from pathways related to melatonin and polo-like kinase in babies with adverse outcome. These preliminary data suggest that transcriptomic profiling may be a promising tool for rapid risk stratification in neonatal encephalopathy. It may provide insights into biological mechanisms and identify novel therapeutic targets for neuroprotection.

Scientific RepoRtS | (2020) 10:13100 | https://doi.org/10.1038/s41598-020-70131-w www.nature.com/scientificreports/ inflammation, infection, metabolic and genetic causes, although these may be clinically indistinguishable 5 . These difficulties have hindered the development of novel neuroprotective strategies. A better understanding of molecular mechanisms of hypoxia-induced brain injury and of the different underlying pathologies has translational potential that may lead to future discoveries and novel neuroprotective therapies. The use of gene expression for disease stratification and for elucidation of underlying disease mechanisms has previously been shown in sepsis, paediatric tuberculosis and Kawasaki Disease [6][7][8] . More recently, we have reported that babies with NE have a unique gene expression profile when compared with healthy controls and septic babies, and have an upregulation of the hypoxia inducible transcription factor 1α (HIF1α) 9 . Here, we examined if babies who developed adverse neurodevelopmental outcome after NE, had a unique host blood gene expression profile at birth.

Results
We analysed the blood samples from the first 47 neonates recruited to a randomised controlled trial of whole body cooling versus usual care-HELIX (Hypothermia for Encephalopathy in Low and Middle-Income Countries) trial 10 . Two neonates with positive blood culture sepsis (Klebsiella pneumoniae and Pseudomonas aeruginosa) were excluded and therefore 45 infants were included in the analysis. The adverse outcome group was defined as death or moderate or severe disability (motor disability and Bayley III). Neurodevelopmental outcome data, based on Bayley assessments carried-out between 18 to 22 months of age, were available in all infants, of which 23 (51%) had adverse outcome and 22 (49%) had good outcome. No significant demographic differences were found between NE infants with adverse and good outcome ( Table 1).
As the trial is ongoing, authors were masked to the treatment (therapeutic hypothermia) allocation, which was provided as A or B. There was no statistically significant difference in outcome between treatment groups A and B, and the clinical severity of NE in two groups was balanced at the time of randomisation allocation ( Supplementary Fig. 1). Differential gene expression analyses comparing infants with adverse and good outcome revealed 855 significant differentially expressed genes with a false discovery rate (FDR) value < 0.05 and an absolute log2 fold change ≥ 0.4, after adjusting for the randomisation allocation and gender (Fig. 1A). Of these genes, 523 (61%) were over-expressed and 332 (38%) under-expressed in NE infants with adverse outcome as shown in the volcano plot.
Finally, we re-analysed the data after adjusting for the randomisation allocation, gender and blood cell type proportions (neutrophils and lymphocyte cell fractions) and identified 660 significant differentially expressed genes with an FDR value < 0.05 and an absolute log2 fold change ≥ 0.4 ( Supplementary Fig. 2). Of these genes, 389 (59%) were over-expressed and 271 (41%) under-expressed. The most significant gene remained RGS1 (FDR value < 0.001, log2 fold change 2.28) together with SMC4 (FDR value < 0.001, log 2 fold change 0.63). We then applied pathway enrichment analysis (Ingenuity Pathway Analysis-www.Ingen uity.com) to all the differentially expressed genes (FDR < 0.01 and absolute log 2 fold change ≥ 0.4). The most statistically significant pathways based on the FDR and ratio values were Mitotic Roles of Polo-Like Kinase (FDR = 0.002, ratio = 0.15) and Melatonin (FDR = 0.003, ratio = 0.11) (Fig. 1C).  CaMK-IV in the cytosol has a proapoptotic effect and is responsible for hypoxic neural cell death both through activation of MAP kinases signalling in the cytosol and through phosphorylation of cAMP response element-binding protein (CREB) in the nucleus, which enhances the expression of pro-apoptotic proteins. Melatonin binds to its plasma membrane receptor MTNR1, to calmodulin and to nuclear receptor retinoidrelated receptor alpha (RORA) increasing its expression. RORA is also considered a downstream target of HIF-1α and its levels have been found upregulated in the cellular response to hypoxia. MTNR1A and MTNR1B activation increases PKC activity through activation of Gαq, which stimulates the PLC signalling cascade and leads to inhibition of Ca 2+ /calmodulin dependent protein kinase (CAMK). Both MTNR1A and MTNR1B activation by melatonin inhibits cAMP formation. Furthermore, activation of MTNR1B decreases the expression of the glucose transporter GLUT4, which in turn decreases glucose uptake. www.nature.com/scientificreports/ We assessed the relationship of genes of interest with different blood cell types (neutrophils and lymphocytes) 11 and found that RGS1, SMC4, as well as the differentially expressed genes of the melatonin pathway were most positively correlated with the lymphocytes fraction (Supplementary Table 1).

Discussion
Our study shows that infants who will go on to develop adverse outcome after NE have a specific gene expression profile, when compared to those who have good outcome. Our data also suggest that polo like kinase and melatonin pathways play a key role in those babies with NE who later develop adverse outcomes. Melatonin is a powerful free radical scavenger, which controls cell death, free radical production, ATP consumption and integrity of the electron transport chain 12,13 . Polo-like kinases instead, are activated by hypoxia/reoxygenation insult and result in a decreasing survival and increasing cell apoptosis 14 .
Despite the widespread use of therapeutic hypothermia as standard treatment in NE, up to 55% of infants with moderate/severe NE are reported to have adverse long-term outcome, in high-income countries 15 . Therefore, there is need to identify adjunct neuroprotective strategies to improve the outcome of infants with moderate or severe NE. Establishing the biological function of differentially expressed genes associated with adverse outcome after neonatal encephalopathy may enable better understanding of brain injury, permitting the identification of new therapeutic targets.
Different pre-clinical studies suggest the potential of melatonin as a neuroprotectant 16 . Melatonin counteracts the effects of hypoxia firstly through the inhibition of ERK/MAK signalling with down regulation of c-AMP response element binding (CREB) protein, secondly through the inhibition Ca 2+ /calmodulin-dependent protein kinases and finally through the upregulation of RORα, which causes a decreased expression of pro-apoptotic proteins and anti-inflammatory effects [17][18][19] . RORA is also considered a HIF-1 target gene and its levels have been found upregulated in the cellular response to hypoxia (Fig. 1C) 20 .
Animal evidence has shown that after hypoxia, there is an increased nuclear Ca 2+ influx and nuclear Ca 2+ / calmodulin kinase IV (CaMKIV) activity, which enhances the expression of different transcription factors responsible for programmed cell death through the upregulation of CREB protein 21,22 . This latter is responsible for the neuronal apoptosis following cerebral hypoxia and is a crucial point of action of therapeutic hypothermia neuroprotective effects [21][22][23] . A recent study has underlined the importance of this cascade, showing that its synergic downregulation together with therapeutic hypothermia, potentiates neuroprotection with reduced cell death and improved neuropathology in a piglet hypoxia model 24 . Of note, CaMKIV levels have been found to increase with the increase in the degree of cerebral tissue hypoxia assessed by cerebral tissue high energy phosphates 25 . Therefore, the present upregulation of CaMKIV in our study, may reflect that infants with adverse outcome were those exposed to more severe degrees of cerebral hypoxia.
The polo-like kinase family are cell cycle process regulators, which are important mediators of the cellular responses to hypoxia and reactive oxygen species exposure 26 . Hypoxia-induced HIF-1α expression is strongly associated with a significant down-regulation of polo-like kinases expression in HeLa cells leading to increased cell survival and proliferation 27 . HIF-1 in fact, has as an important role under hypoxic conditions in regulating cell proliferation and metabolism. Taken together, our findings stress the leading role of HIF-1α in NE, which is involved in both the melatonin and polo-like kinase cascades [27][28][29][30] and interacts directly with both RGS1 and SMC4, which are downstream targets of HIF-1α 30,31 . These results are also consistent with our previous data, comparing healthy infants and NE, which showed an upregulation of HIF1A and MALAT1 in NE 9 . MALAT1 inhibits HIF-1α ubiquitination and in this way, enhances its activity and increases anaerobic glycolysis.
In our study a marked differential expression was found in RGS1 and SMC4. SMC4 is responsible for regulating chromosome organization and dynamics 32 . Recent studies demonstrated also that SMC4 is associated with tumor progression and invasion and the expression of SMC4 is positively correlated with HIF-1a 30 . RGS1 protein is the most abundant RGS protein in the microglia and is a key gene of the immunomodulatory response to neuroinflammation as well as a key target of different immunological and neurodegenerative diseases such as multiple sclerosis, inflammatory bowel and Parkinson's disease 33,34 . Recent evidence suggests that RGS1 protein is involved in the hypoxia-dependent inflammatory response 35,36 , through activation of the AKT signalling pathway, which is switched-on after chronic moderate hypoxia. A leading role of neuroinflammation is now increasingly recognized in neonatal brain damage following NE 37,38 . Therefore, the high expression levels of this gene around birth in babies with adverse outcome later in life, may reflect the increased hypoxia-dependent neuroinflammation in babies with an adverse outcome.
When we assessed the gene expression signature of different blood cell types, we found that the top differentially expressed genes were expressed the most in lymphocytes. This may reflect the immunomodulatory effect of HIF signal, which has been shown to increase the potency of regulatory T cells by affecting the expression of their transcriptional activator Foxp3 39 .
Our study has several limitations. Firstly, the sample size of the study was small and these findings need to be validated in a larger cohort. Secondly, we examined gene expression at only one timepoint, within six hours after birth. Gene expression profiles are likely to change over time in the post-natal period and therefore our findings may differ if gene expression profiles are examined beyond this time window. Finally, this study was conducted in low-middle income countries and therefore, caution is needed before translating these findings to different settings.
Our research was exploratory without a predefined hypothesis on hypoxic-ischemic injury induced gene expression changes and was intended to examine the early host gene expression profile associated with later adverse outcomes. Our preliminary data suggest that babies with adverse outcome after NE had a characteristic gene expression profile that can be measured in whole blood soon after birth, and that selected transcripts correlated better with final clinical outcome than immediate clinical assessment of NE severity. If these findings

Materials and methods
Study design and subjects, samples. The HELIX trial was a multi-country randomised controlled trial of whole-body cooling versus usual care in NE and recruited 408 babies from seven tertiary neonatal units in India, Sri Lanka and Bangladesh between Aug 2016 and Feb 2019 10 . The study was approved by the research ethics committee of Imperial College London, and recruiting hospital sites, and informed parental consent was obtained prior to recruitment. All methods were carried out in accordance with relevant guidelines and regulations. All babies had a structured neurological examination based on the Eunice Kennedy National Institute of Child health and Human Development (NICHD) system by a certified examiner within six hours of birth and the babies were classified as having moderate or severe encephalopathy 40,41 . We collected 0.5 ml of blood (venous or arterial) from the recruited infants within six hours of birth and prior to initiation of cooling therapy, if randomised to the cooling arm. The blood was gently mixed with 1.4 mL RNA stabilizing solution (PreAnalytiX BD/QIAgen) and subsequently stored in a − 80 °C freezer until analysis. The blood samples from the first 47 infants recruited from three Indian sites were included in this work. We excluded babies with culture positive sepsis from the analysis. A flow-chart of the study is in supplementary Fig. 1.
Between 18 to 22 months of age, a certified examiner performed neurodevelopmental evaluation using Bayley scales of Infant development Version III. The primary outcome was death or moderate or severe disability. We defined severe disability as any of the following: a Bayley III cognitive score < 70; a Gross Motor Function Classification System (GMFCS) level of 3 to 5; blindness; or hearing loss (inability to understand commands despite amplification) 41 . We defined moderate disability as a Bayley III cognitive score of 70 to 84 and either a GMFCS level of 2, seizure disorder or a hearing deficit requiring amplification to understand commands 41 .
RNA extraction, alignment and next generation sequencing. We extracted total RNA from whole blood according to the manufacturer's instructions and removed ribosomal RNA and globin mRNA from 4 microg of total RNA. RNA extracts were sequenced using Illumina HiSeq2500 to generate 30 M, 2 × 100 bp reads/sample (Medgenome Labs Ltd, Bangalore, India). Based on the quality report of fastq files, high quality sequence (Q ≥ 30) was retained for further analysis and the low-quality sequence reads were excluded from the analysis. Adapter trimming was performed using fastq-mcf (version-1.04.676) and cutadapt (version 1.8dev). Contamination removal was performed using Bowtie2 (version-2.2.4) and the paired end reads were aligned to the reference human genome (GRCh37/hg19) using "STAR-v2.4.1d". Coverage analysis for genes was performed with "HTSeq-v0.6.1" (reference gtf: GRCh37 version 75) and "bedtools-v2.17.0". The aligned reads were used for estimating expression of the genes and transcripts using the cufflinks program.

Data analysis
Data analysis was performed using R (version 4.0.0). For differential expression analysis we used DESeq2 (v1.17.18) 42 . We corrected all the p values for multiple testing using the Benjamini-Hochberg FDR method to control for type I error. The comparison of babies with adverse and good outcomes was performed first unadjusted and then adjusted for the randomisation allocation (masked as ' A' and 'B'), gender and neutrophils and total lymphocyte cell counts, estimated by using CIBERSORT (deconvolution) 11,43 . The differentially expressed genes were subjected to pathway analysis (Ingenuity Pathway Analysis) and Fisher's Exact was used to assess the significance of the association between the differentially significant genes identified and the canonical pathways. Additional details of analysis are provided in the Supplementary Methods.