Gene expression was profiled using species-specific commercial oligonucleotide arrays. Each array contained a different set of transcripts, with greater redundancy for macaque and swine compared to mouse. There are a total of 43,379 probes represented on the Mouse Whole Genome Gene Expression Microarray, a total of 20,217 probes represented on the Rhesus Macaque Gene Expression Microarray, and a total of 47,813 probes represented on the Porcine Gene Expression Microarray V1. Due to differences in gene annotation across species, differential gene expression analysis was performed for the transcripts associated with 4118 unique genes common to all three arrays using Ensembl gene identifiers associated with each array probe. This strategy accounted for differences in the number of probes with annotated transcripts and reduced the likelihood of falsely identifying differences in gene expression due to gaps in annotation or gene representation. Within this set of genes, there was a large representation of genes associated with Cell Death, Cancer, and Cellular Growth and Proliferation functional categories, as well as genes associated with Glucocorticoid Receptor Signaling, IL-12 Signaling and Production in Macrophages, and Acute Phase Response Signaling canonical pathways.

We investigated the functional classes of underrepresented gene sets present on each specie-specific array, as the focus on unique genes common to all three arrays may introduce bias in our analyses. There was enrichment of Molecular Mechanisms of Cancer, Axonal Guidance Signaling and G-Protein Coupled Receptor Signaling canonical pathways, and Gene Expression functional annotations related to transcription, organismal death and abnormal morphology of cells. Data integration and interpretation with a cross-species transcriptomic analysis brings its own challenges and although complete physical maps have been developed for mouse, macaque and swine genomes to support genome sequencing and comparative genomics, functional annotation to date is mostly based on human, mouse, and rat literature. As annotation improves, particularly for less characterized species such as swine, we will likely be able to more fully understand host responses to influenza virus using microarray and next-generation sequencing technologies.

Of the 4118 unique genes common to all three species arrays, we identified a total of 696 differentially expressed (DE) genes in mice, 771 DE genes in macaque, and 611 DE genes in swine that significantly changed in response to CA04 virus on at least one day (Student’s t-test P < 0.05, average fold-change ≥ 2). The Venn diagram shown in Figure 1 illustrates the overlap of CA04 virus-induced DE genes among mouse, macaque and swine infection models, and only 53 genes were differentially expressed in all three species. Functional analysis of each specie DE gene set revealed significant enrichment of genes associated with Acute Phase Response Signaling, LXR/RXR Activation, and Atherosclerosis Signaling canonical pathways in infected mice, VDR/RXR Activation, Hematopoiesis from Pluripotent Stem Cells, and Communication between Innate and Adaptive Immune Cells canonical pathways in infected macaques, and LXR/RXR Activation, Nitrogen Metabolism, and Antigen Presentation Pathway canonical pathways in infected swine (Table 1). The genes associated with each canonical pathway shown in Table 1 are reported in Additional file 2.

Inflammatory Response was significantly enriched in all three CA04 infection models, with predicted increased activation in macaques and swine (corrected z-scores of 3.437 and 4.02, respectively) (Additional file 3: Table S2). A similar number of DE inflammatory response genes characterized CA04 infection in mice, macaques and swine (175 DE genes in mice; 166 DE genes in macaques; and 160 DE genes in swine), though the kinetics of the host response was distinct in each species (Figure 2). In general, mice showed an increase in gene expression as infection progressed, whereas the majority of these genes were consistently upregulated or downregulated on days 1 and 6 p.i. in CA04 virus-infected macaques. In contrast to mice and macaques, swine exhibited enhanced gene expression on days 3 and 5 p.i. and the host response tapered as the infection resolved. The dissimilarity in the kinetics of gene expression across species may be in part due to differences in virus replication and timing of acute phase responses. Therefore, a more complete kinetic time course, including a range of inoculation doses, would be necessary to fully appreciate the impact of differences in acute phase responses to pH1N1 virus infection.

Among the inflammatory response genes there was representation of Acute Phase Response Signaling, LXR/RXR Activation, and Dendritic Cell Maturation canonical pathways in infected mice, Role of NFAT in Regulation of the Immune Response, Role of Macrophages, Fibroblasts and Endothelial Cells in Rheumatoid Arthritis, and G-Protein Coupled Receptor Signaling canonical pathways in infected macaques, and enrichment of inflammatory response genes associated with Glucocorticoid Receptor Signaling, Dendritic Cell Maturation, and LXR/RXR Activation canonical pathways in infected swine. When we examined the 53 genes common to mice, macaques and swine, we noted the majority of genes were upregulated in response to infection in mice on days 3 and 5 p.i., in macaques on days 1 and 6 p.i., and in swine on days 3 and 5 p.i. (Figure 3). Inflammatory response genes (highlighted yellow) included interferon (IFN) signaling molecule, STAT1, IFN-regulated antigen presentation and immunoproteasome components, TAP1, PSMB8 and PSMB9, Toll-like receptor 4 coreceptor, CD14, and antiviral effectors, IFIH1 and IFIT2, previously shown to play essential roles in innate immune response to influenza infection. We also found several leukocyte-specific genes commonly differentially expressed during CA04 infection. Members of the immunoglobulin receptor superfamily upregulated across species included CD72, CD274 (also known as programmed cell death ligand 1), and CD180, a TLR homologue expressed on B cells, macrophages and dendritic cells. Macrophage-restricted receptor, SIGLEC1, neutrophil factors, PLUNC and NCF4, and T lymphocyte coreceptor, CD8A, were also differentially expressed during CA04 infection in each species (Figure 3, Additional file 4: Table S3).

Several models show an association between enhanced immune cell infiltrate and severe lung immunopathology. Excessive macrophage and neutrophils are observed in the lung of mice following H5N1 and 1918 infection [21], and in pregnant animals infected with pH1N1 virus [22]. In mice, macaques, and swine infected with CA04 virus, leukocyte and lymphocyte responses were evident in all three species based on immune cell-specific gene expression changes detected in the lung, though the gene expression patterns varied across species. For example, there was greater upregulation of macrophage factor, SIGLEC1, observed in mice (days 3 and 5 p.i.) and macaques (days 1 and 6 p.i.) compared to swine on day 3 p.i. Enhanced expression of leukocyte associated immunoglobulin-like receptor 1, LAIR1, was observed in macaques and swine infected with CA04 virus, as compared to infected mice. In contrast to mice and macaques, swine exhibited strong expression of neutrophil factor, PLUNC, which we noted was downregulated during infection in the other species (Figure 3). These results reflect a rapid inflammatory shift in the lungs of mice, macaques, and swine during CA04 infection that involves different immune cell responses, but unlike H5N1 and 1918 viruses, these immune cell responses do not cause immunopathology. A cross-species transcriptomic comparison of the host response to H5N1 and 1918 viruses would be necessary to further explore the potential role of specific immune cell responses to influenza pathogenesis.

We further sought to evaluate CA04 virus-induced host responses that were unique to each species by investigating non-overlapping gene sets shown in Figure 1. The most significant canonical pathways represented in each of these three gene sets are summarized in Additional file 5: Figure S2. The 207 DE gene set unique to swine was enriched for genes associated with Role of JAK1 and JAK3 in γc Cytokine Signaling canonical pathway that included JAK3 and STAT5A genes. Differential regulation of JAK3/STAT5 signaling in CA04 virus-infected swine may impact the homeostasis and activation of peripheral T lymphocytes. For example, STAT5A/B defects results in enhanced apoptosis of T lymphocytes in mice [23]. Differential expression of JAK3 and STAT5A genes in swine during CA04 infection may also suggest a balance between inflammatory host defenses and glucocorticoid receptor (GR)-mediated cellular growth and survival. Closer inspection of each species inflammatory response DE gene set showed noteworthy representation of GR Signaling related genes, such as A2M, FGG, HSPA5, IL4, and MAPK13 genes in mice, FKBP4, JAK1, NCOA2 (also known as GRIP-1), and NRAS genes in macaques, and HSPA4 (also known as HSP70) in addition to JAK3, and STAT5A genes in swine. Alpha-2-macroglobulin (A2M) identified by mass spectrometry in human saliva was found to exhibit antiviral activity against pH1N1 virus [24].

Since GR signaling appears to play a role during CA04 infection in mice, macaques, and swine, we used a transcription factor (TF) analysis strategy that combined transcriptomic and DNA sequence information to investigate GR and other potential regulators controlling host responses during CA04 infection (described in Methods). As summarized in Table 2, we found predicted inhibition of hepatocyte nuclear factor 1alpha, HNF1A, in infected mice on all three days, consistent with our previous report [10], and predicted activation of pro-inflammatory regulators, STAT1, NFκB, and IRF7 in at least two species. Examination of downstream gene targets of GR determined using IPA Upstream Regulator Analysis showed differences in expression across species with more pronounced transcriptional changes observed in mice on days 3 and 5 p.i. and on the early time points in macaques and swine (Figure 4). Through a DNA sequence-based approach, we assessed the representation of DNA-binding motifs in the promoter region of DE genes targeted by HNF1A and STAT1 using DNA sequences from each species genome and regulator DNA-binding preferences, modeled as affinity Position Weight Matrices (PWM) (Additional file 6: Figure S3).

We have identified enriched TFs and targets from our dataset; however, the analyses reported here coordinately investigate the expression and sequence-specific TF DNA binding sites of the DE genes, and it does not investigate the combinatorial effects of multiple factors binding to the genomic regions or molecular determinants of transcriptional responses such as receptor ligation or histone acetylation that can collectively contribute to distinct transcriptional profiles [25]. In addition, restriction of the search space to -450 to +50 nucleotides relative to the transcription start site (TSS) can preclude the identification of more distal TFs effecting expression, which may explain why GR is not identified as a predicted regulator (Table 2), although many target genes are differentially expressed in the data set (Figure 4). Future investigations of GR-mediated host responses to CA04 virus would need to take into account these considerations as well as the physical interactions between GR and STAT5, for example, and chromatin modifications known to occur [26] to provide a more comprehensive view of transcriptional regulation during CA04 infection across different species.

Vitamin D receptor (VDR) and liver X receptor (LXR) are ligand-activated transcription factors of the nuclear receptor superfamily that heterodimerize with retinoid X receptor (RXR) to modulate an array of immune and metabolic programs (reviewed in [27,28]). Perturbation of RXR-mediated signaling pathways have shown unique roles for VDR in macrophage responses to Mycobacterium tuberculosis[29], and LXR in macrophage responses to Listeria monocytogenes[30]. In mice and swine infected with CA04 virus, we found LXR/RXR Activation was among the most significant canonical pathways differentially regulated during acute infection, while VDR/RXR Activation was the most significant canonical pathway differentially regulated in CA04 virus-infected macaques (Table 1). As shown in Figure 5, gene expression profiles of LXR/RXR and VDR/RXR associated pathway molecules, as well as genes enriched for Communication between Innate and Adaptive Immune Cells and Acute Phase Response Signaling canonical pathways strongly emphasize the distinct nature of the host response to CA04 virus in these three animal models.

Macrophages have been shown to increase in mouse lungs following influenza infection and these cells are specific targets of viral infection [21,22]. In response to CA04 virus, we observed a number of genes related to macrophages and with roles in lipid metabolism, such as cholesterol efflux, when examining LXR/RXR Activation in mice and swine. In mice, there was marked decreased expression of apolipoprotein genes, including APOE and APOA1, and moderate downregulated expression of macrophage-specific genes involved in recognition of oxidized lipids, such as macrophage scavenger receptor genes MSR1 (SR-A1) and CD36 (Figure 6, top panel; Additional file 7: Table S4). The ABCA1 gene, known for its role in cholesterol transport, and the LPL gene, involved in lipid catabolism, were also downregulated in mice. In contrast to mice, swine showed strong upregulation of these genes associated with cholesterol efflux on day 3 p.i. Previous studies have shown a decrease in ABCA1 mRNA expression in LXR-activated primary macrophages via an IRF3-dependent mechanism following influenza infection [31]. Here, we observe that ABCA1 and other genes associated with cholesterol efflux are significantly downregulated in response to CA04 infection in mice, but upregulated in CA04 virus-infected swine lung, suggesting that mice may be more efficient at suppressing cholesterol efflux, and in turn, possibly preventing the formation of lipid-laden macrophages (also known as foam cells). Foam cells are the product of inflammatory responses ranging from atherosclerosis (reviewed in [32]) to infection by Chlamydia pneumonia[33] and Toxoplasma gondii[34]. While the link between influenza-induced innate immunity and cholesterol metabolism is poorly defined, our findings suggest a role for LXR/RXR signaling during CA04 infection and indicate that this cellular pathway is impacted differently in two distinct hosts.

In human macrophages and respiratory epithelial cells, vitamin D plays a role in innate immune responses by increasing TLR coreceptor, CD14, and stimulating expression of antimicrobial peptides, such as cathelicidin during respiratory infection [35,36]. Vitamin D has also been investigated in epidemiological studies of reduced risk of influenza infection in different human patient cohorts [37]. In contrast to mice and swine, IL2 gene expression was more strongly induced in infected macaques on both days. This may indicate macaques mount enhanced Th1 cell responses during CA04 virus infection compared to mice and swine (Figure 6, lower panel). As shown in Figure 6, we also found increased, expression of VDR/RXR Activation associated genes such as CYP24A1, encoding a hydroxylase involved in vitamin D catabolism, RXRG, and HR, a transcriptional corepressor of vitamin D receptor. Intriguingly, antimicrobial peptide, CAMP (also known as LL-37), was found to be downregulated in mouse and macaques infected with CA04 virus, and the decreased expression of CAMP in CA04 virus-infected macaques may be explained in part by the increased expression of APOA1 (Figure 6, top panel; Additional file 7: Table S4), which has been shown to bind and inhibit CAMP [38]. While there are several genes associated with VDR/RXR Activation that are differentially expressed in macaques infected with CA04 virus, vitamin D control of innate immune responses does not appear to play a major role during infection.

2009 pandemic H1N1 influenza virus A/California/04/2009 (CA04) was isolated from a nasal swab of a <18 y.o. boy from San Diego, California [43].

Female BALB/c mice (Mus musculus), 6 to 8-week-old, were intranasally inoculated with 10⁶ plaque-forming units of CA04 virus in 50 μl (n = 9) or inoculated with 50 μl of phosphate-buffered saline (control; n = 8). Whole lungs were harvested from infected animals at days 1, 3 and 5 post-inoculation (n = 3 per time point) and from time-matched control animals (n = 3 on days 1 and 3 and n = 2 on day 5) for extraction of total RNA as previously described [10]. Crossbred pigs (Sus Scrofa), 4-week-old, were inoculated intratracheally with either 10⁶ TCID₅₀/pig egg-derived CA04 virus (n = 15) or mock inoculated with non-infectious cell culture supernatant (control; n = 15) as described elsewhere [44]. Animals were euthanized on 3, 5, and 7 dpi (n = 5 per time point). Cynomolgus macaques (Macaca fascicularis), 4 to 15 y.o., weighing 3.0-8.7 kg, were infected with CA04 virus (n = 6) under anesthesia through a combination of intratracheal (4 ml), intranasal (0.5 ml per nostril), conjunctival (0.5 ml per eyelid) and oral (1 ml) routes with a suspension containing 1×10⁶ TCID₅₀/ml (total infectious dose was 7×10⁶ TCID₅₀) as described elsewhere [7]. Animals were euthanized on 1 and 6 dpi (n = 2 per time point).

Mouse infection experiments were completed at the CDC under the guidance of the CDC’s Institutional Animal Care and Use Committee in an Association for Assessment and Accreditation of Laboratory Animal Care International-(AAALAC)-accredited animal facility. Swine infection experiments were completed at the Central States Research Center (CSRC), Inc BSL-3 facility (Oakland, NE) in compliance with the CSRC’s Institutional Animal Care and Use Committee studies. Macaque infection experiments were approved by the RML Institutional Animal Care and Use Committee (IACUC), and performed following the guidelines of AAALAC by certified staff in an AALAC approved facility.

Real-time PCR was performed using a Custom TaqMan Gene Expression Assay (Applied Biosystems) designed for CA04 HA sequence (forward primer: AGCTCAGTGTCATCATTTGAAAGGT; reverse primer: GGACATGCTGCCGTTACAC; reporter: TTGGGCCATGAACTTG). cDNAs were generated using a QuantiTect reverse transcription kit (Qiagen). Samples from individual animals were run in quadruplicate. rRNA (18S) was used to normalize quantification of the target and quantification of normalized target was performed using the 2^-ΔΔCt calculation [45]. HA expression was quantified relative to expression of an endogenous control for each specie sample set that did not change with infection and expression in an uninfected lung sample. The following TaqMan Gene Expression Assays (Applied Biosystems) were used: Mfap1a (Assay ID Mm00849648_gH) served as the endogenous control for mouse samples, B2M (Assay ID Rh02847368_m1) served as the endogenous control for macaque samples, and RPS6 (Assay ID Ss03374061_g1) served as the endogenous control for swine samples. Average log10RQ expression is shown for each species at each time point ± standard deviation.

Total RNA isolated from lung tissue from individual animals on each day of euthanasia was used for oligonucleotide array experiments. For swine and mice, RNA isolated from mock-infected animals at each time point served as an uninfected reference. For cynomolgus macaques, a pool of RNA from the lungs of eight uninfected animals matched for age and sex was used as the uninfected reference. NanoDrop ND-1000 and Agilent 2100 Bioanalyzer instrumentation was used to determine the concentration and quality of all RNA samples. A total of 1500 ng Cy3-labeled probe was used for microarray slide hybridizations, thereby normalizing for the input RNA amount. Mouse, macaque and swine samples were measured with 4×44K commercial arrays from Agilent Technologies designed for each species, Mouse Whole Genome Gene Expression Microarray (G4122F), Rhesus Macaque Gene Expression Microarray (G2519F; Design ID: V2: 026806), and Porcine Gene Expression Microarray V1 (G2519F; Design ID: V1: 020109).

The background corrected data from the Feature Extraction output were normalized across replicates within each species using central tendency normalization (75% percentile, target 1000) within Genedata (Analyst 7.0). All primary microarray data have been deposited in NCBI's Gene Expression Omnibus (GEO) under GEO Series (GSE) accession number GSE40092. The primary microarray data are also available at the University of Washington’s Public Microarray Data Download site ( http://expression.microslu.washington.edu).

For mice and swine, Student’s t-test was performed on background corrected, normalized log-intensity data comparing CA04 virus-infected lung gene expression to time- and species-matched mock-infected lung gene expression at each time point (Unadjusted P-value < 0.05). For cynomologus macaques, Student’s t-test was performed on background corrected, normalized log-intensity data comparing CA04 virus-infected lung gene expression to an uninfected lung reference pool (n = 8) matched for age and sex at each time point (Unadjusted P-value < 0.05). Differentially expressed genes were then filtered to include only genes that changed at least two-fold compared to mock on at least one day within each species. Statistical comparison using Student’s t-test with the Benjamini-Hochberg multiple testing correction [46] (Adjusted P-value < 0.05) and no fold-change parameter resulted in very few macaque genes that were statistically significant, though there was considerable overlap in the DE gene sets identified from the two tests for each species.

Functional analysis of DE genes was performed using Ingenuity Pathways Analysis (Ingenuity Systems), which analyzes the experimental dataset in the context of known biological response and regulatory networks in the Ingenuity Pathways Knowledge Base (IPKB). Ensembl human gene annotations were used for functional analyses of macaque and swine gene sets (Additional file 8: Table S5). The right-tailed Fisher’s Exact test was used to determine the statistical significance of each biological function assigned to the gene expression data, and the Benjamini-Hochberg (B-H) Multiple Testing Correction was applied to p-values to reduce the likelihood that statistical associations were due to random chance.

Upstream Regulator Analysis in IPA incorporates expression of downstream target genes from the experimental dataset and compiled knowledge of reported relationships between regulators and their known target genes within IPKB. This analytical tool was used to predict upstream regulators and infer their activation state by calculating a z-score that determines whether gene expression changes for known targets of each regulator (z > 2, regulator predicted as “activated” and z < -2, regulator predicted as “inhibited”. Transcription factor binding motif enrichment was performed with PSCAN [47] using Position Weight Matrices (PWM) for human and mouse species obtained from JASPAR CORE database [48] and the promoter of each target gene defined from -450 to +50 nucleotides relative to the TSS. PSCAN computes a z-test P-value for each regulator, which is an assessment of whether there is significant representation (P < 0.001) of the regulator DNA-binding motif in promoters of the queried genes.

Transcription factor (TF) DNA-binding analysis was performed using STAT1 (human MA0137.2) and HNF1A (mouse MA153.1) PWMs obtained from the JASPAR CORE database [48] and STAT1 (mouse M00224) and HNF1A (human M00206) PWMs obtained from the TRANSFAC database 7.0 [49]. Promoters sequences of genes analyzed were defined from -450 to +50 nucleotides relative to the TSS and were retrieved from the Ensembl database (release 66). Macaca mulatta sequences were used instead of Macaca fascicularis sequences for the TF DNA-binding analysis. P-values associated for each target gene were calculated using the TFM-Pvalue analytical tool described in [50] and distribution of nucleotides amongst the promoter sequences have been taken into account for background correction. Significant target genes have been identified as having a P-value < 10^-4.