Sex-specific differences of humoral immunity and transcriptome diversification in older adults vaccinated with inactivated quadrivalent influenza vaccines

Different humoral immune responses to QIVs in old adults

To explore the influences of sex differences on humoral immune response to QIVs, HAI assay was performed to compare humoral immune responses between old females and males. The HAI antibody titers were detected in 58 participants (43.1% female) at 0 and 28 days after inoculation (Figure 1). While male and female subjects had similar healthy state and age (Table 1), we found significant difference (p=0.0132) of HAI antibody titer between females and males at 28 days post-vaccination against vaccine H1N1 strain of QIVs formulation (Figure 1A).

Females demonstrated higher median antibody titer after inoculation with QIVs against homologous vaccine H1N1 and H3N2 strains, with significant statistical difference (p=0.0132) against H1N1 strain at 28 days post-vaccination. The humoral immunity against H1N1 and H3N2 vaccine strains was equivalent between females and males before vaccination. Nevertheless, males had stronger immune responses to B Yamagata lineage wild type virus and B Victoria lineage wild type virus strains compared to the females before vaccination (p=0.0069 and p=0.0015, respectively) (Figure 1C, 1D), consistent with the investigation that males could sustain higher antibody titer after influenza infection or immunization [18]. No statistical difference of immune response was observed in the different sex groups to the B_YAM and B_VIC wild type strains after inoculation with the QIVs (p>0.05).

Differential Expression Genes (DEGs) in old adults before and after vaccination

To investigate sex-based differential transcriptome related to increased adverse reaction in females and lower immune response in males after immunization with QIVs, we analyzed whole blood transcriptome of sixteen older adults using next-generation sequencing. DEGs between two sex groups were distinguished by hierarchical clustering method (Figure 2).

At the Day 0 pre-vaccination, 34 upregulated genes and 17 downregulated genes were identified in female compared to male group. At the Day 3, 20 upregulated genes and 14 downregulated genes were identified in female compared to male group. We found 41 and 25 DEGs in the female group, compared to the male group at Day 28 and Day 180, respectively (Supplementary Figure 1). Among thirteen DEGs identified at four time points, only one gene FAXDC2 was on the autosome 5, contributing to the biological processes including oxidation-reduction process, lipid biosynthetic process and sterol biosynthetic process (Supplementary Figure 2).

Multiple DEGs were identified at each time point pre- and post-vaccination (listed in Supplementary Figure 1). More genes were significantly upregulated on Day 0, Day 3 and Day 28 in female group compared to the male counterpart, which may result in variant humoral immune responses in sex groups (Figure 2).

Based on the threshold standard, 80 DEGs were filtered for subsequent analysis. Sex-specific gene transcription was minor in the analysis, with 75% of 16 gene detected only in males with FPKM mean values of transcription < 2, consistent with early reports that significantly fewer genes and pathways were affected in males after influenza infection or vaccination [4, 7, 15].

The schematic diagram of DEGs with sex-bias was shown in Figure 3. The female-bias GYPB gene was highly expressed at four time points, which was activated after immunization and reached the highest expression level at Day 3 (Table 2, females vs. males (Log₂ |FoldChange|) =2.14). The DEGs with sex-bias, dynamically expressed at four time points, which might mirror the development of female-biased higher HAI titers to influenza A/H1N1 vaccine virus (Table 2).

Based on gene clustering and differential analysis, 19, 14, 17 and 12 genes showed differential expression over ten-fold changes between two sex groups at Day 0 pre-vaccination, Day 3, Day 28, as well as Day 180 post-vaccination, respectively. The majority of those genes were upregulated in male group, which are distributed on the Y chromosome (Table 2). The gene transcription intensity detected as sexual dimorphism varied at each time point, with moderate fold-changes detected on the autosomes, and wide variant fold-changes detected on X chromosome.

DEGs with sex-bias over two-fold change and padj <0.05 were analyzed by various databases including Kyoto Encyclopedia of Genes and Genomes (KEGG), DisGeNET, Reactome, Disease Ontology (DO), and Gene Ontology (GO) to investigate biological pathways involved in disease susceptibility and physiological distribution. No significant gene enrichment to the functional pathway was identified by mapping to KEGG database. Analysis based on DisGeNET showed that DEGs upregulated in females were primarily associated with anemia related diseases such as Hemoglobinopathies, Reticulocytosis, Increased bilirubin level, Hereditary spherocytosis and Anemia Hemolytic. DEGs upregulated in males were associated with reduced fertility such as azoospermia and Non-obstructive azoospermia. REACTOME analysis showed two genes KDM5D and UTY on the Y chromosome related to HDMs demethylate histones had lower expression in male group from Day 0 to Day 180 post-vaccination. DEGs upregulated in females included TRIM58, TRIM10, KLF1, ALAS2 and EPB42, which are mainly associated with erythrocyte differentiation, erythrocyte homeostasis, myeloid cell homeostasis and homeostasis of number of cells, while the genes highly expressed in males as KDM5D and UTY were enriched to histone lysine demethylation and protein dealkylation pathways.

Analyzing the correlation of these eighty DEGs collected by sex pairing comparison from four time points with female antibody against influenza A/H1N1 vaccine virus, to identify the hub genes contributing to the female-bias humoral immunity. No DEGs was identified significant statistical correlation with female-biased HAI titers to influenza A/H1N1 vaccine virus.

The robustness of RNA-Seq data was validated by the quantitative real-time PCR (qRT-PCR), detecting the relative expression level of interested genes. The RNA-Seq result validation was conducted referring to previous study, and the gene expression patterns of CD274 and SOCS1 genes were shown to be consistent with the transcriptome results higher expression in female group (Supplementary Figure 3).

DEGs associated with innate immune responses to QIVs vaccination in old adults

To further exploit the differences of DEGs patterns between older sex groups, we compared Day 3 transcripts with the baseline Day 0 in females and males, respectively. There were 16 activated and 24 inhibited genes expressed in females in the early stage after immunization, while no genes with significantly different expression were identified in males (Table 3).

GO enrichment analysis showed that compared with the males, upregulated genes in females were involved in innate immune process, such as type I interferon signaling pathway (IFI35/USP18/IFITM3/IFI6, GO:0060337), complement activation of classical pathway (SERPING1/C1QA/C2, GO:0006958), but down-regulated genes were not enriched in any biological function. By KEGG analysis, the activated genes (SERPING1/C1QA/C2) involved in biological process of complement and coagulation cascades (hsa04610, padj =0.0001), which may contribute to early vaccination responses such as muscle contraction, chemotaxis, phagocyte recruitment and inflammation (Supplementary Figure 4; Complement and Coagulation Cascades, hsa04610).

Hub genes IFITM3 (cor= 0.563, p= 0.023) and C1QA (cor= 0.665, p= 0.005) were positive correlated with female-biased HAI titers to influenza A/H1N1 vaccine virus, would play critical role in the innate immune responses to QIVs.

DEGs associated with adaptive immune responses to QIVs vaccination in old adults

DEGs associated with adaptive immunity to QIVs vaccination were mostly downregulated both in females and males, at Day 28 post-vaccination compared to the baseline Day 0. Total 208 and 101 DEGs were identified in females and males at Day 28, respectively (Figure 4), with 27 downregulated DEGs shared by females and males. Only 6 genes in females including IGKV1-12, AP001000.1, TIPARP-AS1, HSPD1P1, IGHV2-26 and AC008429.2 and 4 genes in males including TRIM36, BIRC5, C2 and C1QB showed higher expression compared to baseline Day 0. The highly expressed genes (BIRC5/C2/C1QB) involved in the regulation of protein processing (GO:0070613, padj=3.66E-05) and complement activation of classical pathway (GO:0006958, padj= 0.0004), and enriched in Complement and coagulation cascades (C2/C1QB, hsa04610, padj=0.017) KEGG pathway term were confirmed in males.

In females, G-protein coupled peptide receptor activity (GO: 0008528) and peptide receptor activity (GO:0001653) associated genes (LTB4R2/RXFP4/MC1R/GPR75) were downregulated at Day 28 compared to the baseline Day 0. Notch receptor 3 (NOTCH3) which contributes to the activation of Th1 cellular immunity, and linker for activation of T cells (LAT) and Fos proto-oncogene (FOS) which directly activate Th2 humoral immunity were significantly repressed in females at Day 28. In males, two genes AOC3 and AOC2 were downregulated, associated with molecular function involved in oxidoreductase activity (acting on the CH-NH2 group of donors, oxygen as acceptor, GO:0016641), which impacted biological process of cellular response to xenobiotic stimulus (GO:0071466).

Hub gene IGHV2-26 (cor= 0.524, p= 0.037) was positive correlated with female-biased HAI titers to influenza A/H1N1 vaccine virus, would play critical role in the adaptive immune responses to QIVs.

Variability of DEGs and biological function pathways associated with decline of humoral antibody titers to QIVs vaccination

Previous studies reported that sustained capacity of protective antibody titers specific to influenza vaccine strains was poor in older adults [19, 20]. In females, GO enrichment analysis for DEGs at Day 180 compared to baseline Day 0 showed that primary DEGs were involved in the process associated with endoplasmic reticulum activity, including protein targeting to ER (GO:0045047), rRNA processing (GO:0006364), polysomal ribosome (GO:0042788) and cytoplasmic side of endoplasmic reticulum membrane (GO:0098554) (Table 4). Differential analysis at Day 180 compared with Day 28 in females showed that DEGs were enriched to multiple terms related to regulation of leukocyte activity including the biological process of positive regulation of myeloid leukocyte differentiation (GO:0002763) and Fc-epsilon receptor signaling pathway (GO:0038095) during the decline stage of specific antibody (Table 4). On the other hand, DEGs in males at Day 180 compared to baseline Day 0 were enriched for GO terms including protein targeting to membrane (GO:0006612), oxidative phosphorylation (GO:0006119), cellular protein complex disassembly (GO:0043624) and protein degradation (Table 4).

KEGG pathway enrichment analysis in males at Day 180 compared to baseline Day 0 showed that DEGs were involved in regulation of cell cycle, cell differentiation and signal transduction pathways including oxidative phosphorylation, thermogenesis, proteasome and pyrimidine metabolism pathways, but no significant enriched immune response pathway was identified in the females.

Ethics statement

Before the activity of this clinical study was employed, the written informed consent was obtained from each subject. The study was conducted in accordance with the Declaration of Helsinki, the clinical protocols were approved by the Medical Ethics Committee of Guangdong Centers for Disease Control and Prevention, and conducted according to the local institutional ethics committee guidelines. The trial was registered with China Drug Trials.org.cn (Registration Numbers: CTR20190846; subjects over 60 yrs). This study population and laboratory detection protocols were described previously [36–39].

Clinical study organization

The initial study cohort recruited 60 physically healthy older adults of Han Chinese (median age 67 years and age range 60-80 years, 43.3% female) who were inoculated by intramuscular injection at nondominant arm with the 2018-2019 seasonal QIVs (Wuhan Institute of Biological Products Co., Ltd., lot SH201805649) or control vaccine (Hualan Biological Engineering Co., Ltd., lot 201809B033; containing the A/Michigan/45/2015 NYMC X-275A [H1N1; an A/Michigan/45/2015(H1N1)pdm09 like virus], A/Singapore/INFIMH-16-0019/2016 IVR-186 [H3N2; an A/Singapore/INFIMH-16-0019/2016(H3N2) like virus], B/Phuket/3073/2013 wild type virus [B Yamagata lineage], and B/Colorado/06/2017 wild type virus [B Victoria lineage]). Whole blood samples (approximate 4 ml each) of the subjects were promptly collected at Day 0 (pre-vaccination), 3-, 28- and 180-days post-vaccination with EDTA-K2 contained tubes. Hemagglutination-inhibiting (HAI) antibody titers were detected from 58 participants (43.1% female) at 0 and 28 days post-inoculation, and physical examination records of all subjects were collected. Subsequently, 16 of 58 subjects were selected for transcriptome analysis of whole blood samples, and all the data were applied to statistical correlation analysis. The characteristics of subjects enrolled was listed in Table 1.

HAI assay

HAI assay was performed in accordance with standard WHO procedures [40]. A 1% suspension of chicken erythrocytes and 4 HA units/25 μl of corresponding influenza virus antigens were employed in the detection of functional antibody titers. Serum samples were treated with receptor destroying enzyme and tested in duplicate in serial 2-fold dilutions initiating from 1:10. Each plate contained both negative and positive serum controls. The HI titer was defined as the inverse of the highest serum dilution to inhibit hemagglutination. For the convenience of statistical analysis, HI antibody titers below 10 were treated as 5.

RNA isolation and mRNA sequencing

Total RNA was extracted from whole blood samples using TRIzol (Bio Basic Inc.) according to the manufacturer’s instructions, and analyzed using RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA) and Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). mRNA sequencing libraries were built using 3 μg RNA and NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer’s recommendations. Poly-T oligo-attached magnetic beads were employed to isolate poly-A RNA. First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. After adenylation of 3’ends of DNA fragments, NEBNext Adaptor with hairpin loop structure were ligated to prepare for hybridization. At last, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBotCluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform (Illumina NovaSeq 6000) and 125 bp/150 bp paired-end reads were generated.

Hisat2 v2.0.5 was used for building index of the GRCh38 (hg19) reference genome, as well as the alignment of paired-end clean reads to the reference genome. The Conditional Quantile Normalization was applied to gene counts to adjust for GC content and gene length. Genes numbers mapped to each gene were counted by FeatureCounts v1.5.0-p3. Subsequent analyses were based on the normalized values.

Differential expression analysis

(For DESeq2 with biological replicates) Differential expression analysis of two conditions/groups (two biological replicates per condition) was performed using the DESeq2 R package (1.16.1). The P-values were adjusted using the Benjaminiand Hochberg’s approach for controlling false rate. Genes with an adjusted P-value <0.05 found by DESeq2 were assigned as differentially expressed.

GO and KEGG enrichment analysis of differentially expressed genes

Gene Ontology (GO) enrichment analysis of differentially expressed genes was implemented by the Cluster Profiler R package, in which gene length bias was corrected. GO terms with corrected P-value < 0.05 were considered significantly enriched by differential expressed genes. Cluster Profiler R package was used to test the statistical enrichment of differential expression genes in KEGG pathways (http://www.genome.jp/kegg/).

Protein-Protein Interactions (PPI) analysis of differentially expressed genes

PPI analysis of differentially expressed genes was based on the STRING database to predict Protein-Protein Interactions. Cytoscape version 3.7.2 software was used for gene network visualization.

Data availability statement

RNA-Seq Expression data in our study were deposited at the National Center for Biotechnology Information Gene Expression Omnibus (GEO) public repository (accession number: GSE151558), to review GEO accession GSE151558: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE151558.

Statistical analysis

Statistical analyses were conducted using IBM SPSS Statistics 26.0.0.0 (SPSS Software). Spearman’s correlation analysis was performed, and p values < 0.05 were considered significant. Independent-sample t test was employed to compare two groups, with a 95% confidence interval. A least significant difference (LSD) of the post-hoc test of one-way analysis of variance (ANOVA) was used to compare multiple groups at a confidence interval of 95%.