Molecular alterations in human milk in simulated maternal nasal mucosal infection with live attenuated influenza vaccination.

Breastfeeding protects against mucosal infections in infants. The underlying mechanisms through which immunity develops in human milk following maternal infection with mucosal pathogens are not well understood. We simulated nasal mucosal influenza infection through live attenuated influenza vaccination (LAIV) and compared immune responses in milk to inactivated influenza vaccination (IIV). Transcriptomic analysis was performed on RNA extracted from human milk cells to evaluate differentially expressed genes and pathways on days 1 and 7 post-vaccination. Both LAIV and IIV vaccines induced influenza-specific IgA that persisted for at least 6 months. Regulation of type I interferon production, toll-like receptor, and pattern recognition receptor signaling pathways were highly upregulated in milk on day 1 following LAIV but not IIV at any time point. Upregulation of innate immunity in human milk may provide timely protection against mucosal infections until antigen-specific immunity develops in the human milk-fed infant.

Introduction

Influenza virus, respiratory syncytial virus (RSV), severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) and other respiratory viral pathogens cause high morbidity and mortality among young infants. Epidemiological studies in low and high income countries show that breastfeeding protects against mucosal infections in infants[1-4]. The immunomodulatory functions of human milk are mediated by pathogen-specific secretory IgA and other immune factors[5,6]. Human milk contains immune cells including monocytes, neutrophils, and lymphocytes (T cells, B cells, and natural killer (NK) cells)[7]. Less is known about innate and cellular defenses in human milk, but they may also play a role[8]. Mothers and infants are usually exposed to the same mucosal pathogens given their close proximity to each other during the infant’s early life. The underlying mechanisms through which pathogen-specific immunity develops in human milk following maternal mucosal infection are not well understood[9,10].

Live attenuated influenza vaccine (LAIV) is inhaled and thus directly stimulates the antigen recognizing cells in the respiratory tract mucosa[11], which together with the lactating mammary gland, the intestine and other mucosal sites constitute the integrated mucosal immune system[12,13]. Nasal secretory IgA is induced following LAIV administration in children and adults, and higher levels correlate with protection from culture-confirmed influenza[14-16]. Stimulation of one mucosal site may result in antibody production at other mucosa-associated lymphoid tissues[17]. Goldblum, Ahlstedt and colleagues[18] showed that ingestion of a certain E. coli strain during pregnancy resulted in high levels of inoculum-specific IgA levels and IgA producing cells in human milk, but not in saliva or blood. Vaccination against Streptococcus pneumoniae during pregnancy has been demonstrated to increase the levels of IgA in human milk, which may be associated with protection in the infant[19,20]. Limited data is available on milk IgA induced by the mucosal LAIV[21]. Intramuscular vaccination of mothers, including against influenza, during pregnancy induces pathogen-specific IgA in milk[22-25]. Human milk immune responses after vaccination beyond IgA have not yet been studied.

In this study, we simulated mucosal infection through LAIV vaccination of mothers to evaluate the subsequent immune responses. Intramuscular vaccination with IIV was used as a control. We further characterized the molecular gene signatures and pathways in human milk cells that may predict an immunogenic response.

Results

We enrolled 78 lactating women during the months of influenza vaccine availability from 2011 to 2016. Participant characteristics are shown in Table 1. The mothers were post-partum for a median of 7.5 days (interquartile range, 5–78 days) at the time of influenza vaccination. None had received a prenatal influenza vaccine during the influenza season of participation until study enrollment. Many (48 [61.5%]) reported never having received an influenza vaccine in a previous season in their lifetime. Figure 1a shows an overview of the enrollment and analyses performed in this study.

Table 1

Participant characteristics.

Characteristic	IIV 32 (41.0%)	LAIV 46 (59.0%)	P value
Race/Ethnicity, N (%)			NS
African American	0 (0)	1 (2.2)
Asian/Pacific Islander	5 (15.6)	2 (4.3)
Hispanic/Latino	24 (75.0)	37 (80.4)
White	2 (6.3)	5 (10.9)
Mixed	1 (3.1)	1 (2.2)
Postpartum day at vaccination, Median (IQR)	9.5 (5–78)	6.5 (5–50)	NS
Exclusive Breastfeeding, N (%)	17 (53.1)	18 (39.1)	NS
Influenza Season, N (%)			NS
2011–12	4 (12.5)	4 (8.7)
2012–13	0 (0)	2 (4.3)
2013–14	13 (40.6)	19 (41.3)
2014–15	10 (31.3)	14 (30.4)
2015–16	5 (15.6)	7 (15.2)
Ever received influenza vaccine in lifetime, N (%)			NS
Yes	14 (43.8)	14 (30.4)
No	18 (56.2)	30 (65.2)
Unknown	0 (0)	2 (4.4)

IQR interquartile range.

Fig. 1

Study design and immune responses induced by vaccination.

a Study overview: Total of 78 breastfeeding mothers were recruited and received doses of either LAIV (n = 46) or IIV (n = 32). Blood and milk samples were obtained for transcriptional profiles at days 0, 1, and 7; milk IgA and blood HAI immune responses were measured on days 0, 30, 90, and 180. b Human milk influenza-specific IgA increases by day 30 and persists above baseline (day 0) through at least 180 days in response to both LAIV and IIV. c Serum HAI increases by day 30 and persists through at least 180 days in serum after IIV but not LAIV. d Circos plot summarizing the antibody response induced by vaccination. The first and second outer tracks represent the vaccination and sample type, respectively. The middle track shows a barplot with max log2[fold change] of influenza-specific HAI or IgA for each sample of blood or milk, respectively. The three inner tracks show a heat map of log2[fold change] of influenza-specific HAI or IgA at days 30, 90, and 180 days. Inner ribbons connect mothers in which both milk and serologic responses are present.

Forty-six (59.0%) mothers received LAIV and 32 (41.0%) received IIV at enrollment. Milk influenza A(H1N1)pdm09-binding IgA responses were statistically similar between those stimulated by the mucosal vs. systemically-delivered vaccines at all time points post-vaccination. Both LAIV and IIV induced an influenza-specific IgA response in milk by day 30 that remained elevated as compared to pre-vaccination baseline titers through at least 180 days (Fig. 1b). Individuals with low A(H1N1)pdm09-specific IgA titers prior to vaccination had a higher increase in response to vaccination compared to those who had a high pre-existing influenza-IgA titers (Supplementary Fig. 1A, B). Serum A(H1N1)pdm09 HAI titers were higher for IIV recipients at all three post-vaccination time points and persisted through at least 180 days; no elevation in mean HAI titers were observed in LAIV recipients (Fig. 1b). Figure 1d shows a summary of the paired milk and serologic responses. Human milk IgA responses did not correlate with serum HAI responses in recipients of either vaccine (p > 0.05).

Human milk samples showed a significantly higher degree of molecular perturbation compared with blood samples in the same individual following either LAIV and IIV delivery (both p < 0.001, Fig. 2a). We used differential expression analysis to normalize sample variations to uncover genes consistently altered across all participants, which allows comparison of the vaccine types at different time points. Herein, differentially expressed genes (DEGs) were analyzed within the vaccination groups in days 1 and 7 compared with baseline for blood and human milk samples. We found 984 upregulated and 584 downregulated genes in human milk following live mucosal vaccination compared with only 8 upregulated and 8 downregulated by systemic vaccination on day 1. The opposite findings were seen in blood where IIV induced more differential expression compared with LAIV (Fig. 2b and Supplementary Fig. 2A–H). There are no statistically significant differences between timepoints within a sample type. An upregulation in genes related to inflammatory and antimicrobial responses in blood after IIV was comparable to blood transcriptional modules seen in previous studies (Supplementary Fig. 3)[26]. LAIV did not induce as many changes in immunological genes in blood transcriptional modules compared with IIV.

Fig. 2

Vaccines perturbation on transcriptome profiles for blood and milk.

a Scatter plot of sample molecular degree of perturbation scores. Scores were calculated for each sample summing the absolute log2 [fold change] >1 divided by the number of perturbed genes. Scores were compared between blood and human milk in both vaccines at day1 and day 7. Line in the scatter plot represents the median value. b Radar chart showing the number of differentially expressed genes for LAIV and IIV in blood and human milk. c Dotplot showing the gene set enriched score in day 1 and day7 for IIV and LAIV. Black borders indicate scores with P adjusted <0.25. d Volcano plot for meta analysis of Spearman correlation between human milk IgA fold-change and gene expression fold-change (at days 1 and 7) and NES from GSEA analysis.

In milk, DEGs upregulated after LAIV on day 1 comprised genes with roles in innate immunity including IL(Interleukin)-6, TLR (Toll-like receptor)-4, TLR6, CD14, EIF2A (eukaryotic translation initiation factor 2α), BST2 (Bone Marrow Stromal Cell Antigen 2), IFI-16 (Interferon Gamma Inducible Protein 16) and TRIM27 (tripartite motif-27). Upregulated DEGs with roles in bridging innate to early adaptive immunity included SPP1 (Secreted phosphoprotein 1), IL18, and IL12. Multiple upregulated genes were leading edge genes in the ATP synthesis pathways including Respiratory complex I (NDUFV1, NDUFA7, NDUFA2, NDUFB8, NDUFA5, NDUFA6, NDUFS3, NDUFA9, NDUFA8, NDUFA12, NDUFA10, NDUFS2, NDUFS4, NDUFA13, NDUFB5, NDUFB6, NDUFAF1, SDHAF2); complex II (SDHAF2); complex III (UQCRQ) and complex IV (COX6B1 COX8A, COX5B, COX7A2L) and a mitochondrial peptidase PMPCB, which cleaves the mitochondrial sequence off newly imported protein precursors.

The number of upregulated DEGs for LAIV in human milk dropped considerably at day 7 (from 984 to 27) in contrast with IIV, which increased in number (from 8 to 52). LAIV up-regulated genes include RORA and SPP1 which plays roles in regulation of the inflammatory response[27]. Both LAIV and IIV shared four up-regulated genes at day 7 (CREB3L1, FOLR1, GPX3, CDR1). CREB3L1 is activated in response to viral infection to inhibit proliferation of virus-infected cells, FOLR1 codes for a folate receptor involved in cell proliferation. A detailed list of all genes uncovered by the differential expression analysis is summarized in supplementary table 1.

Gene set enrichment analyses (GSEA) were performed at subject (Supplementary Fig. 4) and group (Fig. 2c and Supplementary table 1) levels using the sample-wise and the differential expression analysis fold-changes respectively. Individual variation and heterogeneity within each vaccination group was observed for blood and milk at all time points (Supplementary Fig. 4). Regulation of type I interferon production, toll-like receptor, and pattern recognition receptor signaling pathways were highly upregulated in milk on day 1 following LAIV but not after IIV vaccination. Neutrophil mediated immunity; T-cell receptor signaling pathway, and Fc-receptor signaling pathway had stronger positive normalized enrichment scores (NES) in milk at day 1 compared to day 7. The T-cell receptor and Fc receptor signaling pathways were not upregulated until day 7 after IIV. Mitochondrial respiratory related pathways were highly upregulated following LAIV on day 1, and then only mildly upregulated on day 7 in milk. The same pathways showed mild to no upregulation on day 1 but showed high upregulation on day 7 following IIV (Fig. 2c).

To investigate the association of enriched pathways with IgA responses, we performed correlation tests between DEGs and individual normalized enrichment score (NES) and IgA fold-change (Fig. 2d). Both innate and early adaptive immune pathways, including response to type I interferon, positive regulation of phagocytosis, and B cell activation were found to be positively correlated with IgA in response to LAIV vaccine on days 1 and 7 (Fig. 2d). The only pathway significantly upregulated with correlation to IIV was positive regulation of leukocyte chemotaxis, which was seen on day 7 (Supplementary Table 1).

Co-expression analysis of human milk expression profiles revealed the presence of three (M1-3) and four (M1-4) modules of co-expressed genes at day 1 and day 7, respectively (Fig. 3a). The M2 and M3 modules found at day 1 were enriched in immunological pathways including I-kappaB kinase/NF-kappaB signaling, and inflammatory response, while M1 was enriched with mitochondrial and protein processing pathways but include antigen processing and presentation and viral processing. On day 7, three modules enriched in immune pathways (M1, M3 and M4) were seen, while M2 (not shown) was enriched only in lipid biosynthetic processes. M1 is enriched mainly in neutrophil, mitochondrial and protein processing pathways. M3 and M4 were enriched mainly by cytokines and interferon I pathways, respectively. The complete list of genes and pathways enriched for each module is available in Supplementary Table 2.

Fig. 3

Co-expressed modules in human milk.

a Gene co-expression analysis of differentially expressed genes in LAIV and IIV. Functional over representation interaction plot of gene co-expression module using each module as gene list and GO biological process (2018) as gene set. b, c Network obtained from protein-protein interactions (InnateBD) for modules M2 (day 1) and M1 (day 7). Plot highlights gene nodes that act as potential hubs. Genes in bold are DEGs. Nodes were sized and colored proportionally to its degree.

The protein-protein interaction network (Fig. 3b) revealed highly connected genes that were strongly related to innate immunity including EIF2AK2, STAT1, DDX5, TLR2, and TLR4 on M2 day 1. Of those, TLR4 was a significantly upregulated DEG following LAIV on day 1 (Supplementary Table 1). Five DEGs functioned as hubs on day 7 including TLR4, NEDD8, COPB1, USO1 and CCT2 (Fig. 3c).

A deconvolution analysis into cell type proportions revealed no significant differences in cell composition in human milk or blood between the day 0, 1, and 7 timepoints or between vaccine type received. (Supplementary Table 1, FDR < 0.05). We evaluated the association between the top 100 cell type expressed genes from single cell data and DEGs. The upregulated DEGs in milk in response to LAIV on day 1 were significantly enriched with genes expressed by macrophages, monocytes and epithelial cells. In response to IIV on day 7, the upregulated DEGs were enriched with genes expressed by epithelial cells from matured secretory mammary epithelial cells (secMECs). No enrichment associated with a specific cell type was seen in response to IIV on day 1 (Supplementary Table 1).

Discussion

In this study, we combined gene expression profiling with antibody response quantification to identify key genes and pathways associated with reactogenicity and immunogenicity in human milk following mucosal or systemic influenza vaccination. Human milk is a highly immunologically dynamic, functional substance that actively adapts to antigenic stress[28]. We found significantly higher perturbation in human milk compared to blood following both types of vaccinations, suggesting a remarkable maternal investment in the defense of her young infant[29]. Both LAIV and IIV vaccines induced influenza-specific IgA that persisted for at least 6 months. Yet, the expressed genes and enriched pathways suggest that mucosal antigen introduction induces innate immunity pathways in human milk not seen after systemic vaccination.

After mucosal vaccination, TLR4 was upregulated in milk and also found to be a hub for gene co-expression. Hubs have a high degree of centrality such that up or downregulation of the gene has a significant impact on downstream genes and vaccine outcome[30]. TLR4 plays a role in immediate innate recognition of microorganisms by modulating TNF alpha secretion through macrophage-like cells[31]. Viruses including RSV and SARS-CoV have been shown to activate TLR4[32]. Influenza A has been shown to activate TLR4 through interaction with host damage-associated molecular patterns (DAMPs)[33,34]. TLR4 agonists have been studied as adjuvants to be given with inactivated vaccines because of their ability to link the innate with the adaptive immune response and thereby induce robust and persistent antibody responses seen with live vaccines[35]. Additional upregulated DEGs in milk after LAIV on day 1 had roles in innate immunity (BST2, IFI-16)[36,37], interference with viral replication (eIF2α, TRIM)[38-40], or bridging the innate and adaptive arms of the immune system (IL-18, SPP1, NEDD8)[41-43]. Upregulation of type 1 interferon and B cell activation pathways were associated with higher IgA responses in human milk.

In response to mucosal vaccination, we also found multiple upregulated genes in milk associated with the ATP synthesis pathways including the respiratory complex 1 which plays a central role in energy metabolism. Viruses have the ability to induce metabolic reprogramming in host cells similar to the Warburg effect in cancer cells to generate ATP rapidly[44,45]. Targeting of complex I is important for maintaining the mitochondrial membrane potential and ATP production. ATP can serve as a signaling molecule in the regulation of innate immune responses to viral infection[46]. Evaluation of the pathways using GSEA and co-expression analyses revealed that the host immune system appears to sense the influenza antigen as pattern recognition receptors and toll-like receptor pathways were both upregulated, likely leading to production of type 1 IFN. T cell signaling pathways are also upregulated. Mucosal vaccination triggered the mitochondrial and innate immune pathways earlier and stronger (day 1) compared to intramuscular vaccination which showed late and weak upregulation of these pathways on day 7.

Strong induction of innate pathways with LAIV but not IIV likely resulted from live virus vaccination as well as mucosal introduction. LAIV mimics natural infection by replicating within the nasopharynx epithelial cells[47]. Furthermore, mucosal stimulation induces B cells that express mucosal homing molecules that enhance trafficking of memory immune T and B cells to various remote effectors sites[48,49]. Our data suggests that there is a similar communication with innate immunity between mucosal organs such that nasal viral mucosa stimulation upregulates innate pathways distally in human milk. Epithelial cells are an important source of cytokines that shape the response of recruited innate cells such as neutrophils and NK cells that coordinate local and distal responses[50]. Our data suggests that mammary epithelial cells contribute to the upregulated DEGs associated with innate immune responses, in addition to the macrophages, neutrophils and lymphocytes thought to infiltrate milk from blood. Previous studies have found that mammary epithelial cells are the first milk cells to respond during maternal infection and recruit other immune effectors[51].

The study has several limitations. The analyses focused on influenza-specific responses measured against the 2009 H1N1 influenza strain present in the vaccine across all study seasons; responses to other influenza strains in the vaccines are possible. Only transcriptional signatures were assessed rather than protein and cytokine measurement or cellular characterization via flow cytometry. There was a time gap in the conduct of the study and the analysis; this allowed for deeper data interrogation techniques with the advancement of bioinformatic methods. We did not assess for immune responses in the human milk-fed infant.

In conclusion, mucosal vaccination results in a rapid upregulation of innate immune responses in human milk not seen with intramuscular vaccination, providing timely protection until the antigen-specific immunity develops. This may contribute to how human milk protects against mucosal infections in the human milk-fed infant. Nine intranasal vaccines are currently in clinical development to prevent infection and transmission of SARS-CoV-2 with promising data in animal models[52-54]. LAIV is one of the few mucosally delivered vaccines; our findings likely apply to other mucosal vaccines and infections. The human milk responses are not only important for infant protection but also provide a window into mucosal responses delivered elsewhere, extending findings previously described in blood[47].

Methods

Study participants

Healthy mothers were recruited at Children’s Hospital Los Angeles (CHLA) and Hollywood Presbyterian Medical Center (HPMC), Los Angeles, California during the influenza seasons between 2011 and 2016. Mothers of healthy infants born at ≥35 weeks’ gestation and under 6 months of age at enrollment who had not received influenza vaccine during the current influenza season were eligible for enrollment and randomization to receive LAIV or IIV post-partum. The institutional review boards at CHLA and HPMC approved the study; informed consent was obtained from all mothers.

At the initial study visit, study personnel interviewed mothers to collect clinical data including basic demographics, pregnancy history, and immunization history. We recorded mother and infant health and infant feeding characteristics at each visit using the Research Electronic Data Capture (REDCap Consortium, Vanderbilt, Tennessee, USA) database. We obtained maternal blood and milk on days 0 (prior to vaccination), 1, 7, 30, 90 and 180. Samples were stored at −80 °C until further processing.

Hemagglutination inhibition assay

Influenza HA-specific HAI titers were measured in serum. Reference sera and donor plasma were treated with receptor-destroying enzyme (RDE; Denka Seiken Co., Ltd., United Kingdom), and serial 2-fold dilutions of the pretreated and inactivated serum (1:40) were incubated with a standardized amount, or 4 agglutinating doses (4AD HA units), of the designated influenza virus for 1 h at room temperature to allow antigen-antibody binding. An equal volume of 0.5% chicken red blood cell suspension was then added to the mixture. HAI titers were determined after incubation for 30 min at room temperature as the reciprocal of the highest plasma dilution that completely inhibited hemagglutination. Reference sera used as controls included A/California/04/2009 (H1N1).

Human milk IgA quantitation

IgA was quantified in human milk at baseline and post-vaccination time points using the Human IgA ELISA development kit (MabTech 3860-1AD-6). Samples were tested in duplicate. Samples that showed more than twofold variation in final concentration within each replicate were retested. The influenza specific endpoint was determined using detergent lysed concentrated, purified whole virus preparation of A/California/04/2009 (H1N1) (CA04). Normalized CA04 titers refers to the CA04-specific endpoint titer divided by the total IgA.

RNA extraction and Transcriptomic analysis

We compared gene expression profiles in blood and human milk at days 1 and 7 after vaccination to baseline profiles prior to vaccination[55,56]. RNA was extracted from whole blood and milk samples using Maxwell® Simply RNA purification kit (Promega). RNA integrity was assessed using Agilent 2100 Bioanalyzer (Agilent, Palo Alto, CA). The NuGEN Ovation system was used to amplify the RNA from all samples. The mRNA transcripts were profiled using Affymetrix Human Exon 1.0 ST Array (Affymetrix, Santa Clara, CA) at the CHLA Genomics Core Laboratory. Signal intensities for the probe selection regions (PSR) were quantile normalized by robust multichip averaging. The expression signals from the 1.4 million PSRs were summarized into 17,880 transcripts that correspond to known or predicted genes. The gene symbol for each probe was obtained from the annotation Affy_huex_1_0_st_v2. Probes that matched the same gene symbol were collapsed, and the one with the lowest P value was selected. Gene expression fold-change was calculated for each time point in a subject-wise manner and a sample perturbation score was obtained summing the gene expression (absolute values) of all perturbed genes (log2 fold change >1 or <−1) and values were normalized by the number of perturbed genes per sample. Differential expression analysis for vaccine groups between day 1 and day 7 from baseline was conducted using the Limma R package v3.40.6[57].

Publicly available data used for LAIV and IIV Blood comparison

PBMCs or whole blood microarray data from young adults were collected from GEO and processed as described by de Lima, Cardozo and colleagues[58]. Studies data are available to download at GEO (https://www.ncbi.nlm.nih.gov/geo/) and ArrayExpress (https://www.ebi.ac.uk/arrayexpress/).

Gene over-representation and enrichment analysis

For the gene set enrichment analysis (GSEA), the log2[fold-change] gene expression was calculated for each time point in a subject-wise manner for blood and human milk samples. The fold-change for each participant was used as a pre-ranked gene list. Enrichment analysis was performed using fgsea Bioconductor package BTM[59,60] and GO Biological Process 2018[61] were used as gene sets for blood and human milk expression respectively. Redundant GO pathways collapsed using coverage of the gene set of 50%[62].

Analysis of gene co-expression modules and networks

To analyze gene co-expression modules, the log2[fold-change] gene expression in human milk was calculated for each time point in a subject-wise manner and analyzed using the Co-Expression Modules identification Tool (CEMiTool) Bioconductor package[63]. All parameters were left as default except for apply_vst = F and force_beta = T (for day 7). Module networks were created using InnateDB curated Protein-Protein Interactions[64].

Correlation between gene/pathway expression and IgA

To identify signatures associated with immunogenicity, individual pairwise log2-transformed gene expression data (days 1 and 7 post-vaccination) and normalized enrichment score (NES) from GO Biological Process were correlated with influenza-specific IgA enzyme-linked immunosorbent assay antibody titers in human milk (determined at days 30, 90 and 180 post-vaccination) using Spearman correlation. Effect sizes of gene expression and antibody titer correlation analyses obtained in single cohorts were subsequently integrated using the random-effects models for meta-analyses, as implemented in the metafor R package[65]. For the correlation-based analyses, sampling variances of Spearman Correlation Coefficient (SCCs) from each cohort were estimated using the escalc function. To estimate the summarized effect sizes, SCCs and variances were submitted to the rma function using Hunter and Schmidt (HS). Studies were weighted by the inverse of variance.

Cell composition profiling

The deconvolution of bulk gene expression data into cell type proportions was conducted using CIBERSORTx[66]. Bulk whole blood microarray data was deconvoluted using the default LM22 cell type signature matrix, quantile normalization and 500 permutations for significance analysis. Human milk bulk data was deconvoluted using the reannoted single cell collected from GEO (accession number GSE153889)[67], S-Mode batch normalization and 500 permutations for significance analysis. Single-cell analysis was performed using Seurat R package v4.1.0[68]. Genes expressed in less than 3 cells and cells with less than 200 detected genes were removed. Gene expression for each cell was normalized based on total expression and log transformed. Highly variable genes were identified by variance stabilizing transformation (vst) method, and principal component analysis was performed in the Seurat R package. SingleR v1.8.1[69] was employed to annotate cell clusters. The package utilizes pre-annotated reference samples with known annotations to assign those annotations based on similarities in their expression profiles. We used Human Primary Cell Atlas[70] as reference to annotate cell clusters.

Statistical analysis

Data were analyzed using SPSS version 25 (IBM Corp), GraphPad Prism version 7.0 and R statistical program v3.5.3 (2019). Statistical analysis was performed using analysis of variance (ANOVA) with Bonferroni’s post-test, with p < 0.05 deemed as statistically significant unless otherwise stated.

Supplementary Figure

Supplementary Table 1

Supplementary Table 2