Identifying gene network patterns and associated cellular immune responses in children with or without nut allergy.

BACKGROUND: Although evidence suggests that the immune system plays a key role in the pathophysiology of nut allergy, the precise immunological mechanisms of nut allergy have not been systematically investigated. The aim of the present study was to identify gene network patterns and associated cellular immune responses in children with or without nut allergy.
METHODS: Transcriptome profiling of whole blood cells was compared between children with and without nut allergy. Three genes were selected to be validated on a larger cohort of samples (n = 86) by reverse transcription-polymerase chain reactions (RT-qPCR). The composition of immune cells was inferred from the transcriptomic data using the CIBERSORTx algorithm. A co-expression network was constructed employing weighted gene co-expression network analysis (WGCNA) on the top 5000 most variable transcripts. The modules were interrogated with pathway analysis tools (InnateDB) and correlated with clinical phenotypes and cellular immune responses.
RESULTS: Proportions of neutrophils were positively correlated and CD4+ T-cells and regulatory T-cells (Tregs) were negatively correlated with modules of nut allergy. We also identified 2 upregulated genes, namely Interferon Induced With Helicase C Domain 1 (IFIH1), DNA damage-regulated autophagy modulator 1 (DRAM1) and a downregulated gene Zinc Finger Protein 512B (ZNF512B) as hub genes for nut allergy. Further pathway analysis showed enrichment of type 1 interferon signalling in nut allergy.
CONCLUSIONS: Our findings suggest that upregulation of type 1 interferon signalling and neutrophil responses and downregulation of CD4+ T-cells and Tregs are features of the pathogenesis of nut allergy.

Introduction

Today, food allergy affects up to 11% of infants and 3.8% of children, with progressively increasing prevalence over the last few decades. Children with food allergy are at risk for potentially life-threatening allergic reactions including breathing difficulties, swollen tongue, throat tightness and wheezing, particularly when exposed to nut allergens. Moreover, nut allergy tends to be persistent over the full duration of a person's lifetime and can cause anaphylactic reactions., With the exception of treatment of peanut allergy, immunotherapy for desensitizing in food allergy is still not recommended for clinical practice. Hence, the best management for food allergies is to strictly avoid specific food allergens. As such, it is crucially essential to elucidate the fundamental mechanisms of food allergy to help develop effective prevention and treatment for this condition in children.

Food allergy is mechanistically characterized by the development of an overactive immune response to an otherwise harmless allergen, resulting in a T helper 2 (Th2) polarized cytokine response to the allergen and the production of an Immunoglobulin E (IgE) antibody response. Upon recognition of food allergens by antigen presenting cells, mainly dendritic cells (DCs), naïve T-cells are instructed to differentiate into allergen-specific Th2 cells, which are characterized by the expression of type-2 cytokines including interleukin (IL)-4, IL-5, and IL-13. These cytokines then mediate immune responses to food allergens by supporting B-cell proliferation,8, 9, 10 promoting IgE isotype switching, and inducing the activation of mast cells and basophils.11, 12, 13, 14 Given that immune responses to food allergens are complex and involve multiple cell populations, we reasoned that an unbiased systems biology approach could reveal a unique level of insight into the underlying immunological mechanisms. RNA sequencing has been widely used for investigating the pathogenesis of complex human diseases, because it enables the systematic study of the molecular states that underpin pathogenic states. RNA sequencing data can be analysed with systems biology methods such as weighted gene co-expression network analysis (WGCNA), to elucidate the global architecture of the gene expression program and unmask systems-level properties of the biological systems under study. In the current study, we applied an RNA sequencing approach coupled with cell deconvolution and weighted gene co-expression network analysis to identify gene network patterns and associated cellular immune responses in children with or without nut allergy, in order to better understand immunological mechanisms of nut allergy.

Methods

This study was approved by the Curtin Human Research Ethics Committee (Curtin HREC HRE2016-0178) and Child and Adolescent Health Service (CAHS) Human Research Ethics Committee (CAHS HREC 2016046 EP) and conducted in accordance with the National Health and Medical Research Council National Statement on Ethical Conduct in Human Research. Written informed consent was obtained from all parents on behalf of participants. Our study included a total of 92 individuals (84 children with a nut allergy and 8 children without a nut allergy). Children (aged 1–16 years old) with suspected or diagnosed allergies to tree nuts and peanut were recruited from the outpatient clinic of the Immunology Department. Children that were older than 16 years were excluded. Diagnoses of nut allergy were determined by an immunologist on the day of recruitment based on the clinical outcomes of the patients. Nut allergy was defined by immediate symptoms (1–2 h) after nut ingestion combined with either failed oral food challenge or with positive skin prick test wheal diameter≥3 mm to nuts). In contrast, a negative skin prick test and/or a passed oral food challenge for any food allergens defined the subject as children without a nut allergy. Venous blood was collected into PAXgene RNA blood tubes (PreAnalytiX, Qiagen, Hilden, Germany) by a trained phlebotomist.

RNA extraction

Total RNA samples were extracted using PAXgene Blood RNA Kit according to the manufacturer's instructions (Qiagen). Briefly, samples were incubated at room temperature for at least 2 h before RNA extraction. After washing with RNase-free water, the pellet was resuspended in 350 μl resuspension buffer and incubated with 300 μl binding buffer and 40 μl proteinase K for 10 min at 55 °C. The lysate was transferred into a PAXgene shredder spin column and centrifuged for 3 min at 14 000 g. The flow-through fraction was mixed with 350 μl ethanol and transferred to a PAXgene RNA spin column. After washing the column with washing buffer 1, samples were incubated with 80 μl DNase 1 of (RNFD) incubation mix at room temperature for 15 min. The sample was washed several times with washing buffer before RNA was eluted with 40 μl elution buffer and stored at −80 °C. Total RNA concentration and purity were assessed by determining the A260/280 and A260/230 ratios, respectively (NanoDrop).

Library preparation, RNA sequencing and quality control (QC)

The library preparation, sequencing read, QC and read alignment of 30 samples were performed at the Australian Genome Research Facility (AGRF). Briefly, Ribosomal RNA was removed using the Ribo-Zero Gold rRNA removal kit prior to preparation of the cDNA libraries. Library was prepared using TruSeq stranded mRNA Library Prep Kit (Illumina) as per the manufacturer's instructions. The process included purification of mRNA purification via oligo (dT) beads, fragmentation of mRNA with divalent cations and heat, as well as the synthesis of first strand cDNA and 2nd strand cDNA. cDNA libraries were prepared by DNA fragment end repair, 3′ adenylation of DNA fragments, sequence adaptor ligation and amplification of library via PCR. In total, 30 cDNA libraries in 2 separate batches (13 samples and 17 samples) were constructed for sequencing. Samples were then sequenced on Illumina HiSeq 2500 platform with a 100 bp single end read.

Post run processing, including demultiplexing and generation of Fastq files was performed the Illumina bcl2fastq 2.20.0.422 pipeline. The quality of Fastq data were assessed for quality check using FastQC. The data were also screened for the presence of any Illumina adapter/overrepresented sequences and cross-species contamination.

The cleaned sequence reads were then aligned against the Homo sapiens genome (Build version HG38) by the Tophat aligner (v2.1.1). The transcripts were assembled with the Stringtie toolv1.3.3 (http://ccb.jhu.edu/software/stringtie/) utilising the reads alignment with hg38 and reference annotation based assembly option (RABT) generating assembly for known and potentially novel transcripts. The reads corresponding to each gene were summarized using the featureCounts v1.5.3 utility of the subread package (http://subread.sourceforge.net/).

Gene expression data of RNA sequencing

The raw gene read counts were converted into counts per million by using edgeR package of R. Genes with less than 20 counts in the sample with the lowest sequencing depth were excluded to reduce potential noise. Counts were then converted to log2-counts-per-million (logCPM) with precision weights, based on voom method using Limma package of R prior to statistical analysis. Empirical Bayes statistics were applied to compute model statistics and calculate log-fold change (log2FC) of differential expressed genes. The cut-off criteria for significant analysis was: |log2FC|>0.3 and p < 0.05. Age, sex and batch effect were considered as potential confounders.

Reverse transcription-polymerase chain reactions (RT-qPCR)

RNA was reverse transcribed to cDNA using the QuantiTect Reverse Transcription Kit (Qiagen, Germany). One μg of the total RNA from each sample was mixed with gDNA wipeout buffer and RNase-free water to prepare annealing mix, followed by incubation at 42 °C for 2 min. The annealing mix was stored on ice until reverse transcription master mix was added. The final mixture was incubated for 30 min at 42 °C and 3 min at 95 °C to inactivate reverse transcription reaction before storing at −80 °C.

RNA-specific primers for reverse transcription-polymerase chain reactions (RT-qPCR) were designed on the specific region of the genes (Supplemental Table 1) (Supplemental Figure 1). Quantitative PCR was performed on a ViiA7 ™ Real-Time PCR System (Thermo Fisher Scientific). qPCR thermal cycling was set as follows: initial denaturation at 95 °C for 5 min, followed by 40 cycles at 95 °C for 10 s and 55 °C for 30 s and a final extension at 95 °C for 15 s, 60 °C for 1min and 95 °C for 15 s 18s rRNA was measured as the housekeeping gene for qPCR. The fold change in the expression of targeted genes was computed using the ΔΔCt method.

Quantification of immune cells

The RNA-Seq data was then used to estimate the cellular composition of the whole blood samples using CIBERSORTx. CIBERSORTx is a deconvolution algorithm, which can estimate the abundance of different cell types from RNA-Seq data based on cell type specific reference gene expression profiles. CIBERSORTx provides 22 types of functionally defined human immune cells (LM22 signature matrix) as a reference, and we focused on proportions of 12 human immune cells (B cells, CD8+ T-cells, CD4+ T-cells, follicular helper T-cells, Tregs, gamma delta T-cells, NK cells, monocytes, dendritic cells, mast cells, eosinophils and neutrophils). The cell type proportions in children with and without nut allergy were compared employing an independent t-test, and were correlated with coexpression modules derived from WGCNA.

Gene co-expression network construction

Genes with low variation from RNA sequencing data were then filtered out using a coefficient of variation cutoff (CV > 0.05). This resulted in the identification of 9462 stably expressed genes in all the samples. Top 5000 genes with high coefficient of variation values were selected for additional analysis.

The selected genes were then used to construct a co-expression network by using the WGCNA. First, the goodSamplesGenes function was used to filter genes with many mixing values. Then, hierarchical clustering of samples was analysed using the hclust function. The soft-thresholding power was then calculated by using a scale-free topology and the soft-thresholding power was determined based on a scale-free topology index (R2) of 0.85 (Supplemental Figure 2). According to Zhang and Horvath, the gene distribution conformed to a scale-free network if R2 value (model fitting index) is close to 1. Next, the selected soft-thresholding power of 8 was used to calculate the adjacency matrix (correlation strength) between the corresponding pair of genes. The adjacency matrix was then transformed into a Topological Overlap Matrix (TOM) to measure the connectivity of network of the genes and to remove any spurious association. TOM was later performed average linkage hierarchical clustering in order to identify modules of highly co-expressed genes. Network modules were subsequently identified using a dynamic tree cut algorithm with a minimum cluster size of 30 and merge cut height of 0.25 (default). The minimum size of the modules was set to 10 to ensure that small as well as large modules are detected. Subsequently, highly co-occurred genes were merged into modules and these modules later assigned to different colours for visualization.

Module trait relationship construction

Module eigengenes were used to perform principal component analysis of the expression matrix from each gene module. The correlation between module eigengenes (ME), proportions of immune cells, phenotype as well as demographics traits such as age, gender, and batch effect were calculated using Pearson correlation coefficient. Gene modules, which have p-value <0.05, were identified to have significant correlations with nut allergy and these modules were selected for further analysis.

Hub genes selection and visualization

Next, an intramodular analysis was performed to determine the hub genes in the selected modules by summing the connection strengths with other module genes. Hub genes were defined based on the standard cut off of module membership (MM) > 0.8 and gene significance (GS) > 0.5. Module membership (MM) is defined as the correlation between expression profile of a gene and each module eigengene. MM takes values between 0 and 1 and tells “how well a gene belongs to a module”; hub genes have an MM value closer to 1. Gene significance (GS) is defined as the correlation between gene expression and the outcome, nut allergy. The correlation between a module and nut allergy status is in fact a correlation between the module eigengene and nut allergy status. Hub genes of the significant modules were then visualized using Cytoscape v3.8.0.

Functional enrichment analysis

The biological function and associated pathways enriched in each module was characterized by employing gene ontology (GO) and using InnateDB.

Result

Study population

A total of 92 whole blood samples (30 whole blood samples for RNA sequencing and 86 whole blood samples for RT-qPCR) were obtained from children with nut allergy and children without nut allergy (Supplemental Table 2). The 2 study cohorts (RNA sequencing and RT-qPCR cohorts) were well-balanced (p > 0.05) for all the major demographic characteristics. Around sixty three percent of the subjects were boys, with the median age for RNA sequencing cohort and RT-qPCR cohort of 10.5 years and 9.3 years, respectively.

Gene expression profiling of whole blood in children with or without nut allergy

We utilized RNA sequencing to compare gene expression patterns of whole blood samples of children with nut allergy (n = 23) to whole blood samples of children without nut allergy (n = 7). RNA sequencing produced a total of 520 million reads with an average of 17 million mapped reads per sample. A total of 12 523 genes were detected in both samples with and without nut allergy. Differential expression analysis adjusted for age, sex, and batch effect identified 184 upregulated genes (log2FC > 0.3) and 490 downregulated genes (log2FC < −0.3) (Fig. 1). However, it is noteworthy that these results were not significant after multiple corrections.

Fig. 1

The volcano plot for differentially expressed genes. The horizontal axis represents the log2 fold change, and the vertical axis represents the -log10 (pvalue). Red denotes upregulated genes with log2FC > 0.3 and p-value<0.05 while blue denotes downregulated genes with log2FC  <  −0.3 and p-value<0.05. The top 15 differentially expressed genes are labelled.

Validation of RNA sequencing results with RT-qPCR

We selected 3 genes, which were G Protein Subunit Beta 4 (GNB4), Golgi Associated Kinase 1B (GASK1B) and LysM Domain Containing 2 (LYSMD2), for validation by RT-qPCR on a larger cohort of samples according to the following criteria: 1) dysregulated genes based on RNA-seq; 2) high abundance based on logCPM of RNA-seq data; 3) these genes are particularly interesting, associating with infection and inflammation based on literature. We observed the results of RT-qPCR were significantly correlated with RNA sequencing (p < 0.01) (Supplemental Figure 3). These results suggested that our RNA-Seq data was reliable.

We constructed a coexpression network on the top 5000 most variable genes as described in the methods. Through WGCNA, we identified 12 modules of co-expressed genes, and each module comprised between 71 and 1807 genes (Table 1). Among the highly variable genes, only 536 genes (10%) were not assigned to any module, and these genes were clustered into the grey module as per default. Pathways analysis of the modules with InnateDb demonstrated that the modules were significantly enriched for coherent biological functions (Table 1).

Table 1

Modules and associated biological processes.

Module colours	No. of genes	Associated biological processes	Adjusted p-value
Black	190	malonyl-CoA biosynthetic process	0.020
Blue	643	innate immune response	<0.001
Brown	384	negative regulation of dendritic cell differentiation	0.032
Green	266	positive regulation of humoral immune response	0.122
Greenyellow	90	chronic inflammatory response to antigenic stimulus	0.003
Magenta	132	establishment of T-cell polarity	0.037
Pink	187	gene expression	0.001
Purple	116	regulation of cytokine secretion	0.015
Red	229	IMP biosynthetic process	0.174
Tan	71	type 1 interferon signalling pathway	<0.001
Turquoise	1807	transcription, DNA-templated	<0.001
Yellow	349	protein import into peroxisome matrix, translocation	0.211

The module eigengenes were further compared between children with and without nut allergy using module trait association analysis to identify the nut allergy-associated modules. Four modules out of 12 modules were identified to be significantly associated with nut allergy (Fig. 2), which included tan module (r = 0.43, p = 0.03), purple module (r = 0.4, p = 0.04), green module (r = −0.48, p = 0.01) and grey module (r = −0.56, p = 0.003). However, we did not try to find a hub gene for grey module since genes in this module were not co-expressed. It was also observed that greenyellow module was marginally associated with nut allergy (r = 0.4, p = 0.05) and there was no hub gene identified in this module. The upregulated modules (tan and purple modules), were negatively correlated to CD4+ T-cells and positively correlated to neutrophils. In addition, the purple module was also negatively correlated to Tregs. In contrast, the downregulated module (green module) was positively correlated to CD4+ T-cells and Tregs as well as negatively correlated to neutrophils. However, two modules that were not significantly associated with nut allergy (red and greenyellow modules) were also observed to be strongly correlated with neutrophils.

Fig. 2

Module-trait associations. Each row corresponds to a module eigengene (ME) while each column corresponds to either phenotype (NA: nut allergy) or demographic traits such as age and gender. Each cell contains the corresponding correlation coefficient (display at the top of the cell) and corresponding p-values for each module (display at the bottom of the cells within parentheses). Blue and red colours of the spectrum on the right denote low and high correlation, respectively.

The proportions of the 12 human immune cells in children with or without nut allergy are shown in Fig. 3. The cellular composition was dominated by neutrophils, Monocytes, CD8+ T-cells, CD4+ T-cells, B cells and Tregs. Among all the cell types, CD4+ T-cells and macrophages were observed to have a significant difference between children with and without nut allergy.

Fig. 3

Differential immune cell type expression was observed between children with and without nut allergy. Blue colour represents children without nut allergy, while red colour represents children with nut allergy.

Significant modules and associated biological process

The tan module was revealed to have the highest positive correlation with nut allergy. In the tan module, 71 genes were identified to be correlated with nut allergy and these genes were found to be involved in type 1 interferon signalling pathway. Particularly, Interferon Induced With Helicase C Domain 1 (IFIH1) was identified as a hub gene of the tan module (GeneSignificance>0.5 and Module Membership>0.8, Fig. 4).

Fig. 4

Co-expression network of top 30 genes in 3 distinct modules (purple, tan and green), hub gene is indicated with yellow triangle shapes and other highly correlated genes are in round shapes and colour coded according to the module colour.

The purple module was also positively correlated with nut allergy. In the purple module, 116 genes were identified to be correlated with nut allergy and these genes were found to be involved in the regulation of cytokine secretion. Particularly, DNA damage-regulated autophagy modulator 1 (DRAM1) was identified as a hub gene within the purple module (GeneSignificance>0.5 and Module Membership>0.8, Fig. 4).

In contrast, the green module was found to have a negative correlation with nut allergy. In the green module, 266 genes were identified to be correlated with nut allergy and these genes were found to be involved in the positive regulation of humoral immune response. Particularly, Zinc Finger Protein 512B (ZNF512B) was identified as a hub gene of the green module (GeneSignificance>0.5 and Module Membership>0.8, Fig. 4).

Discussion

Food allergy can cause severe, life-threatening anaphylactic reactions. Here, we profiled the transcriptome of immune cells to provide a window into the regulation of immune function in the context of food allergy and anaphylaxis. The application of WGCNA identifies the global connectivity structure of the gene expression program, and unmasks systems-level features of the underlying biology., Herein, we analysed gene expression profiles obtained from whole blood transcriptome of children with and without nut allergy using WGCNA to identify molecular and cellular immune features associated with nut allergy. In the present study, we identified upregulation of the type I interferon production (tan) and cytokine production (purple) modules and downregulation of the humoral immune responses (green) module in children with nut allergy. These changes were positively correlated with proportions of neutrophils and negatively correlated with propotions of CD4 T cell/Treg. The upregulated modules were characterized by hub genes, namely IFIH1, DRAM1, which were involved in type 1 interferon signalling pathway and the regulation of cytokine secretion while the downregulated module was characterized by a hub gene, namely ZNF512B, which was involved in the positive regulation of humoral immune response. Together, our data suggest that the pathogenesis of nut allergy is associated with the upregulation of type 1 interferon and neutrophil responses, and downregulation of CD4 T cell/Treg responses. Consistent with our studies, previous studies have also identified the enrichment of type 1 interferons in subjects with allergic diseases compared to control., The exposure to allergens activates toll like receptors, which in turn mediates the production of type 1 interferons, inducing phosphorylation of STAT1 and inducing expression of pro-inflammatory cytokines.

A limitation of gene expression profiles obtained from whole blood transcriptome is that the data are potentially confounded by variations in cellular composition. The integration of the co-expression network and cell deconvolution approaches allowed us to link co-expression patterns within each module with specific immune cells. CD4+ T-cells and Tregs are already the focus of investigation in studies of food allergy.28, 29, 30 Mouse models of food allergy demonstrated a critical role for CD4+Foxp3+ Treg cells in suppressing food allergy. These studies observed a significant increase in the percentage of CD4+Foxp3+ Treg cells in ovalbumin sensitized mice with mucosal tolerance induction compared to the intolerant group.31, 32, 33 In addition, CD4+Foxp3+ Treg cells were found to inhibit the activation of dendritic cells, mast cells, basophils, and eosinophils, suppress the production of allergen-specific IgE, inhibit Th1, Th2, and Th17 migration patterns and effector functions as well as promote the secretion of IgG4. Consistent with these studies, we also observed a downregulation of CD4+ T-cell and Treg responses in the pathogenesis of nut allergy.

Our analysis also pointed to a possible contribution of neutrophils to the pathogenesis of nut allergy. However, we also observed a significant association of neutrophils with other modules that were not associated with nut allergy. These contrasting results could be due to the heterogeneity of neutrophils as previous studies have revealed variations in phenotype and functions of neutrophils in the development of allergic diseases. On the one hand, neutrophil Fc gamma receptors, FcγRIIIA and FcγRIV as well as several markers of neutrophil activation, S100A8, S100A9, TLR4, TREM1, S100A9 and interleukin 17A, are shown to induce acute anaphylaxis., On the other hand, neutrophils are revealed to be a source of anti-inflammatory and immunomodulatory cytokines, such as interleukin-10. Moreover, neutrophils can suppress NF-kB activation in the macrophage, which in turn reduces the expression of pro-inflammatory cytokines (tumor necrosis factor, chemokine ligand 8 and interleukin-6), leading to the resolution of inflammation. Therefore our findings require further research on neutrophils-related gene networks with nut allergy.

Our study has several limitations. First, the sample sizes used for RNA sequencing are small (n = 30), but gene expression levels of three selected genes are validated using RT-qPCR in a larger sample (n = 86). Second, it is not known if the observed patterns of gene expression are related to the mechanisms that drive nut allergy or alternatively are reacting to nut allergy due to the cross-sectional nature of our study. Third, our study only focused on nut allergy, and therefore it is not known if the mechanisms we identified are relevant to the pathogenesis of other food allergies. Fourth, flow cytometry-based assays to target multiple cell populations were not feasible due to the volume restriction on blood collection from children. However, we have applied an unbiased deconvolution approach to infer the proportions of 12 human cells in whole blood transcriptome. Fifth, we were unable to deconvolve the co-linearity between the module, the phenotype and the cell counts. Therefore, the module train correlations were all explained by cell compositions. Lastly, our study does not define whether changes in cellular proportions precede or follow transcriptomic changes. Further experiments are needed to identify the cellular origin of the gene expression signals associated with nut allergy.

Our study suggests immunological mechanisms of nut allergy can possibly be explained by changes in cell composition. This is related to type 1 interferon signalling, regulation of cytokine responses and the humoral response. Our findings represent plausible pathways for further mechanistic investigation.

Abbreviations

DCs: Dendritic cells; DRAM1: DNA damage-regulated autophagy modulator 1; GASK1B: Golgi Associated Kinase 1B; GNB4: G Protein Subunit Beta 4; GS: Gene significance; IFIH1: Interferon Induced With Helicase C Domain 1; IgE: Immunoglobulin E; IL: Interleukin; LYSMD2: LysM Domain Containing 2; log2FC: Log-fold change; logCPM: Log2-counts-per-million; MM: Module membership; RT-qPCR: Reverse transcription-polymerase chain reactions; Th: T helper; Tregs: Regulatory T-cells; WGCNA: Weighted gene co-expression network analysis; ZNF512B: Zinc Finger Protein 512B.

Author contributions

K.H.L. prepared the manuscript and figures. K.H.L. carried out the experiments and data analysis with assistance from Y.S., K.Y., G.Z. and A.B. Y.S. designed primers and laboratory protocols used in the study. M.O, J.M., B.J.M. and R.L. helped to supervise the project. All authors conceived the manuscript and revised the manuscript. G.Z. coordinated all aspects of the study.

Ethics approval and consent to participate

This study was approved by the Curtin Human Research Ethics Committee (Curtin HREC HRE2016-0178) and Child and Adolescent Health Service (CAHS) Human Research Ethics Committee (CAHS HREC 2016046 EP) and conducted in accordance with the National Health and Medical Research Council National Statement on Ethical Conduct in Human Research. Written informed consent was obtained from all parents on behalf of participants.

Consent for publication

All authors agreed to publish the work.

Funding

This research was funded by Telethon-Perth Children's Hospital Research Fund (TPCHRF) and the Australian Government Research Training Program (RTP) Stipend Scholarship.

Availability of data

The data that support the findings of this study are available from the corresponding author, G.Z., upon reasonable request.

Declaration of competing interest

The authors declare that they have no competing interests.