Sex disparities in influenza: A multiscale network analysis.

Sex differences in the pathogenesis of infectious diseases because of differential immune responses between females and males have been well-documented for multiple pathogens. However, the molecular mechanism underlying the observed sex differences in influenza virus infection remains poorly understood. In this study, we used a network-based approach to characterize the blood transcriptome collected over the course of infection with influenza A virus from female and male ferrets to dissect sex-biased gene expression. We identified significant differences in the temporal dynamics and regulation of immune responses between females and males. Our results elucidate sex-differentiated pathways involved in the unfolded protein response (UPR), lipid metabolism, and inflammatory responses, including a female-biased IRE1/XBP1 activation and male-biased crosstalk between metabolic reprogramming and IL-1 and AP-1 pathways. Overall, our study provides molecular insights into sex differences in transcriptional regulation of immune responses and contributes to a better understanding of sex biases in influenza pathogenesis.

Introduction

Sex-related differences shaped by multidimensional biological characteristics that define females and males exert considerable influence on the pathogenesis of various human diseases (Klein and Flanagan, 2016), including autoimmune diseases (Jacobson et al., 1997), cancers (Cook et al., 2009, 2011; Kim et al., 2018), and infectious diseases caused by diverse pathogens (vom Steeg and Klein, 2016; Fischer et al., 2015; Sawyer, 2012). Each infectious disease exhibits a distinct pattern of sex bias in the prevalence, intensity, and outcome of infections (vom Steeg and Klein, 2016; Giefing-Krôll et al., 2015), as well as in the responses to antiviral drugs and vaccines (Klein, 2012; Morgan and Klein, 2019). For example, females in the human population as well as in other mammals have a higher fatality rate following exposure to influenza A viruses (IAV) (World Health Organization, 2010) and more robust antibody responses and adverse reactions after vaccination than males (Lorenzo et al., 2011; Fink et al., 2018; Potluri et al., 2019), while accumulating evidence suggests a male bias in COVID-19 mortality (Scully et al., 2020; Takahashi et al., 2020). Thus an in-depth understanding of sex differences in the pathogenesis of diseases and consideration in the rational design of prophylactic and therapeutic strategies represent an important step toward segmented medicine.

Sex differences in organisms are as old as the development of sex itself. Sexual reproduction has been reported in early eukaryotes, such as metazoa, fungi, amoebozoa, euglenozoa, or haptophyte (Speijer et al., 2015). According to Rhen, four mechanisms contribute to sexual differentiation: (i) convergent evolution of genetic differences between sexes, (ii) sex-limited or differential expression of autosomal loci, (iii) transgenerational epigenetic effects, and (iv) phenotypic plasticity for sexual traits (i.e., environmental influences on sexual development) (Rhen, 2007). The observed sex differences in disease pathogenesis have been primarily attributed to the differences in the innate and adaptive immune responses between female and male mammals. In both arms of immunity, the sexes differ in multiple aspects, including the detection of pathogen nucleic acids by pattern recognition receptors (PRRs) such as the toll-like receptors (TLRs), the number and activity of immune cells, and the production of cytokines and chemokines (reviewed in (Klein and Flanagan, 2016; Galligan and Fish, 2015)). These intricate differences in immune functions between sexes have a powerful impact on infectious disease pathogenesis. For example, an augmented response to pathogens in females allows better control and clearance of pathogens while promoting increased immunopathology (vom Steeg and Klein, 2016; Klein, 2012; Fischer et al., 2015). During influenza infections, female mice exhibited a more robust induction of pro-inflammatory cytokines and chemokines in their lungs, including TNF-, IFN-, IL-6, and CCL2, accompanied by greater weight loss, hypothermia, and mortality than male mice (Robinson et al., 2011).

The etiology underlying sex differences in immunity involves intrinsic and extrinsic factors that exert a combinatorial effect on the immune system’s functioning. Although the microbiome (Markle et al., 2013; Yurkovetskiy et al., 2014, Vom Steeg and Klein, 2017; Vemuri et al., 2019) and nutritional status (Khulan et al., 2012; Tobi et al., 2009; Sinha et al., 2003; Kawai et al., 2010; Osrin et al., 2005; Christoforidou et al., 2019) have been implicated in modulating immune responses, sex hormones and genetic mediators are the most widely appreciated factors shaping differential immunity between females and males. In addition to the profound effects of sex hormones that have been extensively demonstrated, genetic differences attributed to immune-related genes and microRNAs (miRNAs) that are located on the sex chromosomes also play an important role in determining the distinct immune responses between sexes, especially in prepubertal children, postmenopausal females, and age-matched males (reviewed in (Galligan and Fish, 2015; Klein and Flanagan, 2016; Fischer et al., 2015; Fish, 2008; Klein, 2000; vom Steeg and Klein, 2016; Schurz et al., 2019; Bianchi et al., 2012)). However, the pathways and cellular responses that mediate the differences in response to influenza infection have not been well elucidated. Moreover, the confounding effects of hormone and genetic factors in intact animals and human populations impose a challenge in resolving the mechanism of sex differences in influenza pathogenesis.

To address the molecular mechanisms underlying sex differences in influenza infection responses, we assessed the transcriptomic dynamics of blood cell responses in neutered female and male ferrets throughout the infection. We detected a temporal shift in mounting immune responses against viral infection between females and males. Using a multiscale network-based approach, we further identified pathways and cellular processes commonly induced in both sexes and those uniquely regulated in females or males. Our data revealed gene regulatory pathways that were differentially regulated between sexes and likely involved in marked sexual differences in influenza virus-induced pathogenesis, providing molecular and functional insights for functional evaluation and development of therapeutic strategies.

To evaluate the genetic factors-mediated sex differences in the immune response to influenza virus infection, we conducted a study with neutered female and male ferrets. Ferrets are considered the gold standard in influenza research (Skarlupka and Ross, 2020; Belser et al., 2011). The ferret is currently the only small mammalian model available that is equally well suited to studying both influenza pathogenesis and transmission. Ferrets and humans show very similar clinical and pathological signs of influenza infection, depending on the host age, viral strain, environmental conditions, degree of secondary bacterial infection, and other variables. Clinical features in both organisms involve fever, nasal secretion, coughing, gastrointestinal complications, serum abnormalities, weight loss, lethargy, lymphopenia, hypercytokinemia, and transmission to susceptible contacts (Belser et al., 2011).

Concerning cellular immune response, cytotoxic lymphocytes are the primary cells important for attacking the immune system against infected cells. In the case of the humoral immune response, IgA antibodies seem to protect the upper tract in both humans and ferrets (Maher and DeStefano, 2004). Overall, the ferret model has played an essential part in mapping out the pathogenesis of the typical influenza infection, transmission, and reassortant variants of the virus, Reye’s syndrome, drug treatment, and resistance of the influenza disease to many antiviral drugs (Burlington et al., 1981; Herlocher et al., 2003).

Results

An integrated network approach for systematic characterization of the immune response to influenza infection in females and males

We examined the transcriptional response of neutered female and male ferrets in blood cells over the course of the infection (Figure 1A). Briefly, whole blood samples were collected from spayed adult female and castrated adult male ferrets infected with the influenza virus A/CA/07/2009 (H1N1pdm09) strain over a period of 1 to 8 (for males) or 1 to 14 (for females) days postinfection (dpi). Samples from infected animals and uninfected control samples, including those collected before the infection from infected animals or before and after infection from uninfected individuals, were used to serve as the baseline for the comparison of infection responses within each sex. The whole blood transcriptome was profiled to assess the impact of genetic sex on immune response regulation without the additive influence of sex hormones. Given the time frame of sample collection, this study primarily focused on identifying sex differences in innate immunity and the early phase of adaptive immunity (Rowe et al., 2010).

Figure 1

Overview of study design and analysis scheme

(A) The sample collection strategy of adult female and male ferrets infected with influenza A/CA/07/09 (H1N1pdm09) virus. In the serial sacrifice study of female ferrets (n = 45), blood samples were collected at the uninfected baseline (n = 26) and days 1 (n = 10), 3 (n = 12), 5 (n = 8), 8 (n = 7), and 14 (n = 8) postinfection from animals sacrificed at a given day. For the standard study of female ferrets (n = 9), samples were collected at baseline and from naive (n = 4) and infected (n = 5) animals at days 3, 5, and 8 postinfection. Similarly, for the standard study of male ferrets (n = 34), samples were collected at baseline (n = 33) and at days 1, 3, 5, and 8 postinfection (n = 30).

(B) Framework for systematic characterization of sex difference in immune response to influenza infection with whole blood transcriptome. Transcriptomic data generated from infected and control samples of each sex were used for differential expression analysis and Multiscale Embedded gene co-expression network analysis (MEGENA). Disease-relevant modules of co-expressed genes were then identified by the enrichment of differentially expressed genes (DEGs) and the association with physiological traits. Key regulators in those selected modules were predicted using the adopted Fisher’s inverse Chi-square approach in MEGENA (see STAR Methods). Pairwise comparisons were performed to further determine similarities in module composition, module key regulators, and intra-module connectivity (see STAR Methods).

To systematically characterize sex differences in the immune response to influenza virus infection, we used a co-expression network-based approach that integrated transcriptional patterns with physiological traits and enabled the identification of sex-specific gene expression signatures (Figure 1B). Specifically, differential expression analysis was carried out within each sex by comparing the infected samples with the uninfected controls at each time point after infection. A Multiscale Embedded Gene co-Expression Network Analysis (MEGENA) was then performed (Song and Zhang, 2015) on the transcriptomic data for each sex separately. We further prioritized the co-expressed gene modules in each sex identified by MEGENA based on their enrichment for the respective sex-specific differentially expressed gene (DEG) signatures that captured molecular responses to influenza infection. We characterized these modules with the Molecular Signature Database (MSigDB) and blood cell-type-specific marker gene signatures. Disease-relevant modules were then correlated with pathological traits, including weight loss and viral load, and selected for comparison of response signatures and co-expression patterns between sexes. Key regulators determined by multiscale hub analysis (MHA, see STAR Methods), as well as module compositions and within-module correlation patterns in selected modules, were further investigated to elucidate shared or sex-unique processes and regulatory patterns.

The pathophysiological response during infection is significantly different between the sexes

We were first interested in the overall physiological responses, such as weight loss and viral load (Figure S1) and potential sex difference during influenza infection. We observed inter-ferret variability and identified significant differences in weight loss and viral load, using a 95% CI. Concerning the weight loss of infected ferrets, the sickest female maintained 77% of her original body mass after eight days (Figure S1A). In comparison, the healthiest infected females kept their body mass at 97%. For males, we observed a slightly narrower range of weight loss. The sickest male lost weight by 22% after eight days, whereas the healthiest infected male maintained his weight at 94%. Weight loss at day 8 postinfection is significantly different between females (mean 87.8%) and males (mean 85.4%; p = 0.011; Figure S1A). With respect to viral load, we have observed the largest variation at 1 dpi. At this time point, females experienced a viral load with a range between 3.05 × 106 and 3.10 × 108 PFU/mL (mean 9.70 × 107 PFU/mL) and males with a range between 150 and 1.07 × 108 PFU/mL (mean 4.11 × 106 PFU/mL). The difference in viral load at 1 dpi between females and males is significant (p = 0.0011; Figure S1B). Overall, although females experienced a significantly larger viral load at 1 dpi, their weight loss during infection is less severe than the males’.

We further investigated the cell content of immune cells during infection. As we did not have measured immune cell data (e.g., from flow cytometry measurements), we utilized publicly available single-cell data. Owing to the limited availability of ferret single-cell information, in particular from ferrets infected with the influenza virus, we settled for a recently published dataset by Lee et al. from bronchoalveolar lavage fluid (BALF) of COVID-19 infected ferrets (Lee et al., 2021) (GEO: GSE171828). Although this was suboptimal, we hypothesized that influenza and COVID-19 had similar immune cell-type-specific transcriptional signatures in BALF and blood during the first few days after infection. However, a more appropriate, single-cell-based approach to this infection scenario is required for detailed analysis. With these limitations in mind, we deconvoluted the ferret bulk dataset with the Lee et al. single-cell data using MuSiC (Wang et al., 2019). The originally assigned cell types by Lee et al. included dendritic cells, macrophages, granulocytes, mast cells, natural killer (NK) cells, γδ-T cells, CD8+ T cells, CD4+ T cells, proliferating T cells, B cells, plasma cells, and epithelial cells.

Overall, red blood cells are overly abundant in both female and male ferrets (Figures S2A and S2C). But other cell types such as granulocytes and M1 macrophages are similarly abundant. We were further interested in the relationship between cell types and physiological traits. Most pronounced in both sexes are the negative correlations between weight and cell types, such as granulocytes, M1 macrophages, mast cells, plasma cells, and proliferating T cells (Figures S2B and S2D). Granulocytes and macrophages are the first line of immune cell response. So are mast cells that identify viral threats through several different receptors, such as TLR3 or Rig-I (DDX58). Interestingly, these three cell types are strongly positively correlated with viral load in females but not as much in males. On the other hand, immune cells of the adaptive immune system, such as proliferating T cells or plasma cells, have strong negative correlation with weight in males and weaker correlation in females. Generally, females mount faster and more robust innate and adaptive immune responses than males. However, our temporal window only extends to 8 days (14 days for females), which does not cover the full extent of adaptive immune response. Thus, we identified a strong correlation between innate immune cells and viral load in females. On the other hand, we observed a quenched immune response by immune cells of the adaptive immune system. However, this restricted response of the adaptive immune system could also be attributed to the limited temporal window of sample selection, as a more robust adaptive immune response may occur at a later time.

We were further interested in the detailed temporal response of immune cells (Figures S2E–S2H). Most noteworthy is the female response of dendritic cells compared to males. At 3 dpi, the relative abundance of dendritic cells in females is almost three times higher than in males (1.94 vs. 0.66%, resp., p = 1.11 × 10−4). In contrast, dendritic cells in males remain almost constant during the entire course of infection. Berghöfer et al. discovered in healthy subjects that TLR7 ligands induce a higher IFN-α production in females (Berghöher et al., 2006). However, the authors did not differentiate cell abundance between the sexes. Thus, other processes induced by the viral infection may cause the increased production of dendritic cells in females compared to males. Although taken with caution as cell type abundance is relative and inferred by transcription-based single-cell deconvolution, overall, a higher relative abundance of immune cells was observed in females compared to males.

Differential expression analysis reveals global temporal differences in immune response dynamics between sexes

Following the recognition of viruses, various immune cells equipped with versatile strategies to mount a potent antiviral response are mobilized to establish physical and chemical barriers against viral infection. To determine the global dynamics of blood cell response in females and males over time, we first carried out differential expression analysis of the whole blood transcriptome data collected over eight- or 14-day periods after influenza infection in comparison with sex-matched uninfected controls. Among different sets of DEGs identified over time in females and males (Tables S1 and S2), the most significant overlaps between any two given sets were observed at 3 and 5 dpi within or between the sexes, followed by those detected between 3 to 5 dpi in each sex and 1 dpi in females or eight dpi in males (Figures 2A and S3). It is noteworthy that DEGs at 14 dpi in females did not share significant overlaps with any other sets (Figure 2A), which suggests a largely subsided initial immune response by 14 dpi. This finding is consistent with the reported kinetics of an early innate phase in biphasic immune responses during influenza infections in ferrets (Rowe et al., 2010). To further understand the biological processes associated with DEGs during viral infection, we explored the enrichment of MSigDB gene sets in the upregulated and downregulated DEG signatures at each time point over 8 days in each sex. As expected, immune response related pathways are enriched in the signatures in both sexes and they included the innate immune response and cytokine signaling — particularly type I () and II () interferon (IFN) signaling — hemostasis, the urokinase-type plasminogen activator (uPA) and its receptor (uPAR)-mediated signaling, DNA replication, and the cell cycle (Figure 2B). These processes exhibited differential temporal kinetics of activation, with IFN signaling being initiated immediately upon infection. However, cell cycle-related processes were revealed later in the infection (Figure 2B). We also detected enrichment of downregulated genes during infection in T cell functioning and IL-17 pathways in both sexes, and several generic cellular processes were related to transcription and translation in males (Figure 2B, Tables S3 and S4).

Figure 2

Global dynamics of blood cell response following infection in both sexes

(A) Pairwise comparison of differentially expressed genes (DEGs) at each time point in both sexes. DEGs were identified by comparing the infected samples with the respective controls for each sex. The solid bars are colored by -log10(p value) generated with the super exact test and the numbers at the top show the observed intersects of DEGs between a pair of DEG sets. The semitransparent blue bars overlaying the solid bars indicate the expected intersects. The green dots below denote two sets of DEGs used in the comparisons. The numbers to the right of the dot matrix indicate the number of DEGs in a given set. Detailed comparisons of upregulated and dowregulated DEGs between sexes over time can be found in Figure S3. DEG lists can be found in Tables S1 and S2.

(B) Molecular Signatures Database (MSigDB) gene sets enriched in the DEG signatures at each time point in both sexes. The size of dots indicates adjusted p values from the enrichment tests, and the color denotes the sex. The top five gene sets enriched in each DEG signature were shown and highlighted with strokes around the dots. Detailed information can be found in Tables S3 and S4.

To further investigate the temporal dynamics of infection-responsive processes in both sexes, we examined the fold change (FC) of DEGs associated with the enriched MSigDB gene sets over time. We found more robust activation of type I and II IFN signaling pathways in females than males, as many genes were immediately and more strongly activated at 1 dpi in females (Figures 3A and 3B). A temporal shift in the expression kinetics between sexes was also evident in platelet activation, signaling, and aggregation, which were important in maintaining hemostasis and modulating inflammation (Le et al., 2015). Many genes involved in platelet activity were immediately activated at 1 dpi and substantially subsided by 8 dpi in females, in contrast to an attenuated activation at 1 dpi and continuous expression beyond 8 dpi in males (Figure 3C), further suggesting a more rapid response in females than males. Interestingly, we also detected a DEG, GLIPR1L2 (GLI pathogenesis-related one like 2), with extreme temporal shift as manifested by inverse alteration kinetics between the sexes (r = −0.9993, p = 0.0238). GLIPR1L2 was downregulated at 1 to 3 dpi and returned to baseline expression by 5 dpi in the females but was continuously downregulated from 3 dpi onward in males (Figure 3D). Although the function of GLIPR1L2 in viral infections remains unclear, it has been shown to reside in the same module with GLIPR1 (GLI pathogenesis-related protein 1) and GLIPR1L1 and can be targeted by the tumor suppressor p53 during cancer development and progression (Ren et al., 2006). Taken together, these results indicate temporal differences in the immune response dynamics between the sexes, marked by a prompt response in females compared with a lagged response in males.

Figure 3

Temporal differences in dynamics of response to infection between the sexes

(A–D) log2(fold change (FC)) of DEGs involved in (A) interferon (IFN) signaling, (B) IFN signaling, and (C) platelet activation, signaling, and aggregation over the course of the infection in both sexes. DEGs that were considered significant (i.e., with a log2(FC) ≥ log2(1.5) or ≤ -log2(1.5)) at one or more time points in females or males are shown. (D) log2(FC) of GLIPR1L2 over time. The colors denote different comparisons within each sex.

Co-expression network analysis identifies shared and sex-unique expression patterns in the immune response

Although differential expression analysis enabled a global overview of the immune response in both sexes over the course of the infection, it did not provide the resolution for gene sets with sex-unique expression patterns. To dissect the sex-biased factors, we next sought to employ a network-based approach to identify consistent and unique gene expression patterns between sexes with increased granularity. We performed MEGENA on the transcriptomic data corresponding to each sex (Tables S5 and S6) and prioritized the modules of co-expressed genes in each sex based on their enrichment of DEGs. To assess the direct link between modules and pathogenesis, we also examined correlation between modules and physiological features, including viral titers in the nasal washes and weight change (Figure S1). Overall, we identified modules enriched for upregulated or downregulated DEGs, or both, at one or more time points in each sex (Figures 4A, 4B, S4, Tables S7, and S8). Those modules were enriched for a diverse array of gene set signatures involved in distinct but related biological processes, such as IFN signaling, inflammation, immune cell migration, hemostasis, and apoptosis induction and clearance exerted by various blood cells (Figure 5), revealing more specific processes than those associated with DEGs globally. Some modules were also positively associated either with viral load in one sex or both or with weight loss in males (Figures 5 and S5). For example, viral load-related modules in both females and males were consistently enriched in the IFN signaling and primarily type I IFN signaling pathways. In contrast, viral load-related modules in females alone were associated with specific inflammation-related pathways, lysosome-related processes, the unfolded protein response (UPR), and hemostasis. Those only in males were associated with apoptosis-related pathways (Figures 5 and S5). Although we did not detect any weight loss-related modules in females, the weight loss-related modules in males were associated with generic processes and pathways such as transcription, translation4-1BB-dependent immune response and porphyrin metabolism (Figures 5 and S5), suggesting a complex relationship between weight loss and underlying molecular processes. Moreover, several modules in each sex were enriched for blood cell-type-specific signatures: (i) the female modules M35, M196, and M524 and the male module M127 were enriched for the T cell signature which was associated with either the translocation of ZAP-70 or the immunological synapse, (ii) the male module M15 was enriched for the generation of second messenger molecules pathway, and (iii) the male modules M17, M134, M356, and M523 were enriched for the CD19+ B cell signature which was related to primary immunodeficiency (Tables S7 and S8). These results elucidate a comprehensive landscape of complex immune responses elicited by diverse types of immune cells. Although both sexes have a common arsenal of strategies as immune system response, females and males further employ different pathways and processes contributing to the fine granularity of sex disparity in such a host-pathogen system.

Figure 4

Co-expressed gene modules enriched in infection-responsive genes in both sexes

(A and B) Sunburst plots showing the hierarchy of all the modules and those enriched in differentially expressed genes (DEGs) in the (A) female and (B) male co-expression networks, respectively. Modules are colored by the enrichment in upregulated (red) or downregulated (blue) DEGs or both (yellow). Modules enriched in the MSigDB gene sets observed in both sexes are highlighted. Detailed information about each module’s composition can be found in Tables S5 and S6, and information about modules enriched in DEGs at each time point in Figure S4 and Tables S7 and S8.

Figure 5

Co-expressed gene module functions and temporal responses

Bubble plot showing the MSigDB gene sets associated with modules enriched in DEGs in both sexes. The size of the bubbles indicates the number of modules enriched in an MSigDB gene set, and the color denotes the sex. The modules and their hierarchy were shown on the right side of the plot and colored by sex. The asterisks and squares next to the MSigDB gene sets indicate the presence of correlation between at least one module among those associated with the corresponding MSigDB gene set and physiological traits (i.e., asterisks for viral load and squares for weight loss). Only the correlations with p value ≤ 0.05 and coefficient ≥ 0.5 were included. Detailed information can be found in Figure S5. Detailed time-resolved development of the IFN signaling modules and the AHSP pathway modules for females and males can be found in Figures S6 and S7, and in Figures S8 and S9, respectively. Detailed information about key regulators in modules can be found in Figure S10.

We also investigated the temporal response of the co-expression networks. For this purpose, we constructed separate co-expression networks for each time point. We assessed the similarity (or difference) of modules between time points by using the Jaccard index, focusing on the gene content by disregarding the network connectivity. This allowed more robust temporal trajectories between modules with significant overlap measured by Fisher’s exact test (FET). For our comparative analysis, we focused on overall common responses, such as the upregulated type I IFN signaling and the downregulated AHSP pathway. Figures S6 and S7 show the time-resolved module development of type I IFN signaling in females (relating to modules M104 and M400) and males (relating to modules M3, M26, M165, and M478), respectively. Figures S8 and S9 depict the AHSP pathway in females (relating to modules M414) and males (relating to modules M6, M60), respectively.

We observed a similar temporal shift between the sexes in type I IFN signaling module conservation as we observed with the gene expression of immune genes and corresponding functions. The female type I IFN response is established and conserved early between 1 dpi and 3 dpi (Figures S6A and S6C), whereas males muster their IFN defense much later (Figure S7A). Their IFN response pathway is organized at 3 dpi and 5 dpi (Figure S7C), with an established network at 5 dpi continuing to eight dpi (Figure S7D). However, the male IFN signaling network displayed an early organized phase between 0 dpi and 1 dpi but was disorganized at 3 dpi (Figure S7A). On the other hand, interferon-induced helicase C domain 1 (IFIH1) is central in IFN signaling for both females and males. The overall downregulated AHSP pathway shows similar behavior (Figures S8 and S9). We observed the strongest conservation in module membership and connectivity between 5 dpi and 8 dpi in males (Figures S8A and S8D). On the other hand, females have a highly conserved temporal response (especially the interaction between hemoglobin subunit beta (HBB) and the novel gene ENSMPUG00000010158) across all time points (Figure S9A).

To identify the molecular factors mediating sex differences in the immune response to infection, we performed an in-depth comparison of the structural and functional relationships between modules across the female and male networks. Given that some processes and pathways were consistently detected over the 8-day period in the female and male networks, we first asked if the modules enriched in the same gene set signatures in both sexes shared similar members and key regulators. These shared signatures depicting temporal dynamics of the immune response included (i) IFN signaling (primarily type I IFN), lysosome-related processes, hemostasis, mitotic cell cycle, and DNA replication, and graft versus host disease signature, which exhibited gradual activations over time, and (ii) primary immunodeficiency signature, AHSP pathway, and translocation of ZAP-70 to the immunological synapse, which were generally suppressed in both sexes (Figure 4), similar to those observed in DEGs. We conducted pairwise comparisons between the modules in the female networks and those in the male one by using FET to determine the similarity of modules and their key regulators (Table S9). Not surprisingly, the modules enriched for the same signatures in both sexes had a significant overlap (Table S9). However, the key regulators in modules with matched signatures shared limited overlaps (Figure S10 and Table S9). These results suggest similar transcriptional programs in both sexes upon infection are coordinated by different regulators in females and males.

To further pinpoint the sex-specific gene co-expression patterns that implicated distinct functional signatures, we next focused on the modules enriched for sex-unique gene signatures to determine the uniqueness of those modules to each sex. We examined module-module overlap by using FET and performed differential gene correlation analysis (DGCA) (McKenzie et al., 2016) of the modules enriched for the signatures uniquely observed in one sex. DGCA assesses if those genes possessed a different regulatory relationship in their counterpart module initially determined by using FET in the other sex (Table S10). We identified one female module, M139, which was unique to females determined by FET (Table S9) and was associated with the activation of chaperone genes by XBP1(s) (Figure 5 and Table S7). DGCA identified two male modules, the AP-1 (activator protein-1)-associated M152 and the IL-1-associated M569, which were specific to males (Tables S8 and S10; Figure 5). XBP1 (X-box binding protein 1) was the key regulator of the female-specific module M139 (Figure 6A). The majority of the genes within the female module, M139, exhibited a rapid and strong induction immediately following infection at one dpi in females (Figure S11A) compared with a milder induction in males (Figure 6B), suggesting a prompt and potent activation of the IRE1/XBP1 branch of the UPR (So, 2018; Hetz, 2012) during infections, specifically in females. The male-specific module M152 was enriched for later-activated DEGs, including two key regulators HPGD (15-Hydroxyprostaglandin dehydrogenase) and MS4A2 (the chain of the high-affinity IgE receptor (FcRI)), which were downregulated at three dpi and upregulated at 5 to 8 dpi in both sexes (Figures 6C, 6D, and S11B) and likely linked to FcRI-mediated and prostaglandin E2 (PGE2)-mediated inflammatory responses (Ugajin et al., 2011; Inage et al., 2014; Mkaddem et al., 2019; Rastogi et al., 2020; Tai et al., 2006). The other male-specific module, M569, enriched for the upregulated DEGs at 3 to 8 dpi in males (Figure S11C) contained two key network regulators, ACSL3 (acyl-CoA synthetase long-chain family member 3) and IL18RAP (IL-18 receptor accessory protein) (Figures 6E and 6F), implying a link between lipid metabolism and inflammation (Yan et al., 2015; Roelands et al., 2019; Remmerie and Scott, 2018; Batista-Gonzalez et al., 2020; Kobori et al., 2018; Yasuda et al., 2019). Taken together, these results revealed sex-dependent regulatory relationships and the extent of the modulation of genes involved in IRE1-XBP1, AP-1, and IL-1 pathways upon infection, suggesting sex differences in specific aspects of the UPR, lipid metabolism, and inflammation independent of sex hormones.

Figure 6

Networks of sex-unique modules and the temporal expression pattern of their key regulators and members

(A, C, and E) Networks of co-expressed genes in (A) a female-specific module M139 and (C and E) two male-specific modules M152 and M569. Node shapes denote module membership (diamonds for key regulators and circles for regular members), and their sizes are proportional to the node strength. The solid purple lines in the center of nodes show the log2 (fold change (FC)) of a given DEG, and the dashed gray lines indicate zero.

(B, D, and F) log2(FC) of (b) DEGs in the female-specific module M139, or that of the key regulators in the male-specific modules (d) M152 and (f) M569. The colors indicate the comparison within each sex (i.e., pink for females and blue for males). Detailed information about DEGs in each module over time can be found in Figure S11.

Discussion

A comprehensive understanding of the molecular mechanism underlying sex differences in response to viral infections has been either largely limited to sex chromosome-linked genes or hindered by the presence of confounding factors, especially in human populations. In this study, we systematically characterized sex differences in the transcriptional landscape of blood cell responses to influenza virus infection in neutered ferrets, thus precluding the impact of other factors. We detected temporal and regulatory differences in the processes commonly induced in both sexes as evidenced by a more rapid response in females than males and different key regulators orchestrating those shared responses. With a network-based approach, we further identified sex-specific modules of co-expressed genes that exhibited unique regulatory relationships in each sex. Those sex-unique modules revealed a rapid UPR activation via the IRE1/XBP1 pathway upon infection in a female-biased manner and a male-specific regulation of inflammatory responses mediated by a cascade involving FcRI, AP-1, and PGE2, and the IL-1 family signaling associated with lipid metabolism. These data provided a comprehensive picture of genetic factors-mediated molecular differences in the immune response to influenza infection between sexes.

Sex-differentiated gene expression in normal and disease states of human tissues has been identified in many studies (Kassam et al., 2019; Trabzuni et al., 2013; Lopes-Ramos et al., 2020; Bongen et al., 2019; Jansen et al., 2014). However, the interpretation of underlying molecular mechanisms of infectious disease pathogenesis, especially genetic factors-mediated baseline differences between the sexes, remains challenging because of the inevitable presence of confounding factors such as sex hormones, age, microbiota composition, and nutritional status in human studies. With the neutered adult ferrets, our study provided a unique opportunity to interrogate sex-biased transcriptional response to influenza virus infection driven by genetic mediators. From the blood transcriptomic data, we detected a temporal difference in the dynamics of immune responses between sexes, manifested by a more rapid activation in females than males upon infection. The temporal difference between the sexes was observed at the gene co-expression network level at different time points. The immediate response in females is likely attributed to germline-encoded bias in innate immune sensing, as TLR7 that is located on the X chromosome and encodes an innate PRR-recognizing single-stranded RNA can not only have a higher expression level in females than males because of X-inactivation escape (Pisitkun et al., 2006) but also induce more robust IFN- production in peripheral blood mononuclear cells (PBMCs) isolated from females than those from males in vitro (Berghöher et al., 2006). The sex disparity in TLR7 expression may also relate to the observed difference in the response of dendritic cells. After cell type deconvolution using a public single-cell dataset, we identified a significant increase in relative abundance of dendritic cells in females at 3 dpi, whereas the dendritic cell abundance in males remained constant.

We further found a limited overlap in key network regulators modulating commonly altered processes upon infection between the sexes, suggesting differential regulation of various aspects of immune responses in females and males. This observation is consistent with a previous report of sex differences in regulatory targeting patterns identified across multiple healthy tissues in humans (Lopes-Ramos et al., 2020). We further identified a peculiar gene, GLIPR1L2, with an extreme temporal shift as manifested by inverse alteration kinetics between the sexes. Although its potential role during infection remains obscure, GLIPR1L2 is a target of tumor suppressor p53 during cancer development and progression (Ren et al., 2006).

Using a network-based approach, we identified sex-specific modules that represented distinct pathways differentially impacted by sex and exhibited unique regulatory structures and co-expression patterns. These modules revealed several sex-biased processes involved in UPR, lipid metabolism, and inflammation. We found a robust early activation of IRE1/XBP1-dependent UPR in the female ferrets. Activation of IRE1 can lead to the splicing of XBP1 mRNA (XBP1u) and subsequent expression of spliced XBP1 (XBP1s) that acts as an active and stable transcription factor for UPR-related genes (So, 2018; Hetz, 2012). Extensive studies on the role of UPR in immune cells demonstrated a fundamental impact of the IRE1/XBP1 pathway on promoting innate pro-inflammatory responses following TLR signaling in macrophages (Martinon et al., 2010) and on modulating differentiation of several cell types (reviewed in (So, 2018)), including stimulating the differentiation of B cells into antibody-secreting plasma cells (Iwakoshi et al., 2003). For example, the UPR is known to play a role in influenza virus infection (Mehrbod et al., 2019). Although predominantly observed in the lung, influenza virus infection induces ER-stress and activates the IRE1 pathway with subsequent XBP1 splicing (Hassan et al., 2012). Interestingly, emerging evidence from glycomic profiling of lung tissues collected from influenza-infected adult female ferrets indicates a link between the UPR, glycosylation, augmented immune response, and severe infection outcome; the activation of the IRE1/XBP1 pathway induces high mannose immediately upon infection, which can be recognized by the innate immune activator mannan-binding lectin (MBL2) and is associated with alveolar severity (Heindel et al., 2020). Moreover, a higher level of XBP1s expression in rat livers from prepubertal females compared with age-matched males has been reported, suggesting genetically determined female bias in XBP1s expression (Rossetti et al., 2019). Taken together, rapid activation of the IRE1/XBP1 axis of UPR observed specifically in females, likely because of genetically encoded bias in XBP1s expression, may trigger potent immune and inflammatory responses and eventually could have contributed to a more severe disease seen in females than in males. Future studies will need to determine the functional consequences of XBP1s activation on a specific immune cell type during influenza infection and evaluate the impact of sex hormones on the pathways outlined.

Our integrative network analysis also captured two male-biased signaling pathways related to inflammatory response and lipid metabolism. Genes involved in the IL-1 pathway, including IL1RAP, JUN, IRAK3, and IRAK4, as well as two key network regulators (ACSL3 and IL18RAP), exhibited stronger inductions from 3 dpi onwards in males. ACSL3, primarily expressed in the ER and cytosolic lipid droplets, is responsible for fatty acid (FA) activation by converting free FA to FA acyl-CoA esters, an initial step contributing to lipid biosynthesis and FA oxidation (Roelands et al., 2019; Yan et al., 2015; Tang et al., 2018). Interestingly, lipid metabolism has been shown to dramatically alter the phenotype and function of monocyte/macrophages and modulate inflammatory responses, as increased lipid synthesis fueled by glucose metabolism has been linked to inflammatory macrophages, and FA oxidation is associated with alternatively activated macrophages (Remmerie and Scott, 2018; Batista-Gonzalez et al., 2020; Yan and Horng, 2020). However, the impact of ACSL3 induction and lipid metabolism on inflammatory responses in the blood cells and the crosstalk between the metabolic reprogramming and IL-1 family member IL-18-mediated signaling require further investigation. In addition to genes involved in the IL-1 pathway, we found another module enriched for the AP-1 pathway with a male-specific regulatory pattern. However, the extent of change in expression of its key regulators upon infection is similar between sexes. The AP-1 pathway composed of members in the Jun, Fos, Maf, and ATF subfamilies (Eferl and Wagner, 2003; Shaulian and Karin, 2002) plays an important role in regulating immune and inflammatory responses such as T cell activation and cytokine production (Atsaves et al., 2019; Zenz et al., 2008; Foletta et al., 1998; Schonthaler et al., 2011). Intriguingly, two key regulators in this AP-1 pathway-related module — MS4A2 acting to amplify FcRI expression and signaling during allergic responses (Kraft et al., 2004) and HPGD responsible for inactivating potent lipid mediator PGE2 (Tai et al., 2006) — have been implicated in susceptibility to anaphylaxis by modulating PGE2 homeostasis in the mast cells (MCs) (Rastogi et al., 2020). PGE2 stabilization by inhibiting HPGD-catalyzed degradation can attenuate FcRI-mediated MC degranulation, preventing MC hyperresponsiveness and reducing susceptibility to anaphylaxis (Rastogi et al., 2020). Thus, these observations likely indicate genetically encoded male bias in the expression and regulation pattern of genes involved in specific aspects of lipid metabolism and inflammatory response in the blood cells, which requires further investigation.

Taken together, our data revealed an emerging picture in which genetic factors-mediated sex differences could manifest in the temporal dynamics and regulatory relationship of commonly induced immune and inflammatory responses, the activation of the IRE1/XBP1 pathway during the UPR, and the modulation of lipid metabolism and inflammatory response involving the IL-1 and AP-1 pathways. Although this proposed paradigm awaits further validation, especially the assessment of its generalizability in the presence of sex hormones, it provides insights into the molecular mechanism of sex differences in the cellular function and activity of blood cells, which can also be relevant to the pathogenesis of other infectious and inflammatory diseases. We noticed several limitations in this study, including the lack of granularity in differentiating cell type-specific responses and a restricted time course that does not allow a thorough investigation of adaptive immune responses, particularly humoral-mediated immunity. However, our study represents an important step toward a comprehensive understanding of sex differences in influenza virus-induced pathogenesis, which holds promise in developing precision medicine with improved efficacy for both sexes.

Limitation of the study

This study uses animal blood as a typical and convenient resource to analyze the host response against influenza virus. However, as a respiratory disease, the influenza virus primarily infects the respiratory system of the host. But the extraction of samples, in particular from the lung, usually requires sacrifice of the laboratory animal. A second limitation constitutes the fact that female and male ferrets were neutered because of animal welfare considerations. In particular, male non-neutered ferrets are quite aggressive against their fellow ferrets. Thus, sex hormone levels as contributing factors are lower in the studied neutered ferrets compared to the levels in non-neutered animals. However, with the neutered adult ferrets, our study provided a unique opportunity to interrogate sex-biased transcriptional response to influenza virus infection driven by genetic mediators.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Christian Forst (christian.forst@mssm.edu).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

Animal experiments

Female and male Fitch ferrets (Mustela putorius furo) were obtained from Triple F Farms (Sayre, PA, USA) and were seronegative to a circulating influenza A/California/07/2009 (A/CA/07/09; H1N1pdm09) virus. As previously described (Bissel et al., 2019), adult ferrets, 6 to 12 months of age, were pair-housed in stainless steel cages (Shor-Line, Kansas City, KS) that contained Sani-Chips laboratory animal bedding (P.J. Murphy Forest Products, Montville, NJ) and provided with food and fresh water ad libitum. To achieve a mean of 20% weight loss no earlier than 8 dpi, adult ferrets were infected intranasally with the H1N1pdm09 virus at a dose of 106 plaque-forming units (PFU). The animals were monitored daily for the severity of clinical disease using weight loss. Disease symptoms, including elevated temperature, low activity level, sneezing, and nasal discharge, were noted if present (data not shown). Any animal losing > 20% weight loss was humanely euthanized. Blood was collected from anesthetized ferrets via the anterior vena cava after infection. The University of Georgia Institutional Animal Care and Use Committee approved all experiments, which were conducted in accordance with the NIH’s Guide for the Care and Use of Laboratory Animals (Bissel et al., 2019), The Animal Welfare Act, and the Biosafety in Microbiological and Biomedical Laboratories guide of the Centers for Disease Control and Prevention and the NIH.

Ethics statement

All research studies involving the use of animals were reviewed and approved by the Institutional Animal Care and Use Committees of the University of Georgia. Studies were carried out in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals.

Method details

Viral plaque assay

Plaque assays were performed to determine the viral burden in snap-frozen nasal washes, as previously described in (Bissel et al., 2019). Supernatants were gently thawed on ice, forced through a cell strainer (70 μm) and syringe plunger in phosphate-buffered saline, then spun down (2500 rpm, 5 minutes, 4°C) to collect the supernatant. Supernatants were diluted in Iscove’s modified Dulbecco’s minimum essential medium. Madin-Darby canine kidney (MDCK) cells were plated (5 × 105) in each well of a six-well plate. Samples were diluted in Iscove’s modified Dulbecco’s minimum essential medium (final dilution factors of 100 to 10−6) and overlaid onto the cells in 100 μL of Dulbecco’s modified Eagle’s medium supplemented with penicillin-streptomycin followed by 1 hour incubation. Samples were removed, cells were washed twice, and medium was replaced with 2 ml of L15 medium plus 0.8% agarose (Cambrex, East Rutherford, NJ, USA) and incubated for 72 hours at 37°C with 5% carbon dioxide. Agarose was removed and discarded. The cells were fixed with 10% buffered formalin and then stained with 1% crystal violet for 15 minutes. After thorough washing in distilled water to remove excess crystal violet, the plates were dried, the number of plaques was counted, and the number of PFU per milliliter was calculated.

RNA-seq experiments

RNA was extracted from ferret blood using the Mouse RiboPure-Blood RNA Isolation Kit (Ambion). Both extraction methods followed the manufacturer’s protocols and incorporated a DNase treatment (QIAGEN) after passing the sample through the filter cartridge. Strand-specific total RNA-seq libraries from ribosomal RNA-depleted RNA were prepared using the TruSeq Stranded Total RNA Library Prep kit (Illumina) according to the manufacturer-supplied protocol. Libraries were sequenced 100 bp paired-end to a depth of approximately 40 million host genome reads on Illumina HiSeq 2500 instruments.

Quantification and statistical analysis

RNA-seq analysis

Paired reads were aligned to the ferret genome Ensembl version 1.0.80 using tophat (version v2.0.13). Following read alignment, featureCounts (v1.4.6) was used to quantify the mRNA expression levels based on the Ensembl gene model. mRNAs with more than 5 reads in at least 1 sample were considered expressed and retained for further analysis. Otherwise, they were removed. The mRNA read counts data was normalized using the trimmed mean of M-values normalization (TMM) method of the edgeR package (Robinson et al., 2010) to adjust for sequencing library size difference.

Identification of differentially expressed genes

Given the time series information for each different dataset, we used two distinct methods to identify differentially expressed genes. (i) The t-test was used to identify differentially expressed genes (DEGs) between case and control. (ii) We used a one-way ANOVA-like approach as implemented in the Limma-package (Smyth, 2004) to identify differentially expressed genes across the time series—termed significant response genes. Significance is defined based on a false discovery rate (FDR) of 5% or less. DEGs is used as a generic term for differentially expressed genes. Statistical significance of intersections between any two given DEG sets in both sexes was determined with SuperExactTest (Wang et al., 2015).

Gene co-expression network analysis

Multiscale Embedded gene co-expression network analysis (MEGENA) (Song and Zhang, 2015) was performed to identify host modules of highly co-expressed genes in influenza infection. The MEGENA workflow comprises four major steps: 1) Fast Planar Filtered Network construction (FPFNC), 2) Multiscale Clustering Analysis (MCA), 3) Multiscale Hub Analysis (MHA), 4) and Cluster-Trait Association Analysis (CTA). The total relevance of each module to influenza virus infection was calculated by using the Product of Rank method with the combined enrichment of the differentially expressed gene (DEG) signatures as implemented: , where, is the relevance of a consensus j to a signature i; and is defined as , where is the ranking order of the significance level of the overlap between module j and the signature.

Identification of enriched pathways and key regulators in the host transcriptome modules

To functionally annotate gene signatures and gene modules derived from the host transcriptome data, we performed an enrichment analysis of the established pathways and signatures—including the MSigDB and ARCHS4 tissues (Lachmann et al., 2018). The hub genes in each subnetwork were identified using the adopted Fisher’s inverse Chi-square approach in MEGENA; Bonferroni-corrected p-values smaller than 0.05 were set as the threshold to identify significant hubs.

Correlation analysis

The correlation between modules and physiological traits was performed using Spearman’s correlation. The correlation between the FCs of DEGs in females and males over time was estimated using Pearson’s correlation with the Benjamini-Hochberg correction.

Similarity comparison of modules and their key drivers

The pairwise similarity between a given female module and every male module or between their key drivers was estimated using Fisher’s Exact Test (FET) with Bonferroni correction. Modules in one sex without counterparts in the other sex (FET p ≤ 0.05) were considered sex-unique module compositions. Differential Gene Correlation Analysis (DGCA) was also performed with the “moduleDC” function implemented in the DGCA package (McKenzie et al., 2016) to further determine the unique regulatory relationship of the overlapped gene members between a pair of similar modules in females and males. Modules with significant differential correlation (p ≤ 0.05) with all of their counterparts were considered to exhibit sex-specific regulatory relationships.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Bacterial and virus strains
Influenza A/California/07/2009	International Reagent Resource	IRR: FR-201
Deposited data
Ferret RNAseq	https://www.ncbi.nlm.nih.gov/geo/	GEO: GSE168512
Data	https://www.synapse.org	https://doi.org/10.7303/syn26485897
Computer code	https://www.synapse.org	https://doi.org/10.7303/syn26485897
Experimental models: Cell lines
Madin-Darby canine kidney (MDCK) cell	International Reagent Resource	IRR: FR-242 (RRID:CVCL_0422)
Experimental models: Organisms/strains
Ferret: Ferret: -	Triple F Farms	Fitch ferrets
Software and algorithms
TopHat: v2.0.13	https://ccb.jhu.edu/software/tophat/index.shtml	RRID:SCR_013035
Ensembl Genomes: ferret: v1.0.80	https://www.ensembl.org	RRID:SCR_006773
featureCounts: v1.4.6	https://cran.r-project.org/	RRID:SCR_012919
LIMMA: v3.44.1	https://www.bioconductor.org	RRID:SCR_010943
SuperExactTest: v1.0.7	https://cran.r-project.org/	https://cran.r-project.org/package=SuperExactTest RRID:SCR_003005
MEGENA: v1.3.7	https://cran.r-project.org/	https://cran.r-project.org/package=MEGENA RRID:SCR_003005
ARCHS4	https://maayanlab.cloud/archs4/	RRID:SCR_015683
DGCA: v1.0.2	https://cran.r-project.org/	RRID:SCR_020964