Construction of a Transcriptome-Driven Network at the Early Stage of Infection with Influenza A H1N1 in Human Lung Alveolar Epithelial Cells.

We aimed to understand the molecular changes in host cells that accompany infection by the seasonal influenza A H1N1 virus because the initial response rapidly changes owing to the fact that the virus has a robust initial propagation phase. Human epithelial alveolar A549 cells were infected and total RNA was extracted at 30 min, 1 h, 2 h, 4 h, 8 h, 24 h, and 48 h post infection (h.p.i.). The differentially expressed host genes were clustered into two distinct sets of genes as the infection progressed over time. The patterns of expression were significantly different at the early stages of infection. One of the responses showed roles similar to those associated with the enrichment gene sets to known 'gp120 pathway in HIV.' This gene set contains genes known to play roles in preventing the progress of apoptosis, which infected cells undergo as a response to viral infection. The other gene set showed enrichment of 'Drug Metabolism Enzymes (DMEs).' The identification of two distinct gene sets indicates that the virus regulates the cell's mechanisms to create a favorable environment for its stable replication and protection of gene metabolites within 8 h.

INTRODUCTION

Influenza A H1N1 Virus (INV) is a RNA virus belonging to the family Orthomyxoviridae. Commonly referred to as the flu virus, H1N1 is specific to birds and mammals. The latest worldwide flu pandemic occurred in April 2009 and accounted for an estimated 284,500 global deaths (Dawood ).

In general, a variety of host responses during viral infection have been identified, including activation of numerous cell death and survival pathways. These pathways include: 1) programmed cell death I (apoptosis), 2) programmed cell death II (autophagy), and 3) endoplasmic reticulum stress with subsequent unfolded protein response (UPR). There has been extensive research on the regulatory roles of these pathways during the influenza virus life cycle (Yeganeh ). Pathogen sensory pathways (e.g., RIG-I) displayed long-lasting associations with cytokine/chemokine signaling through day 8 (Dimitrakopoulou ). Transcriptome research for H1N1 infection revealed that host cell innate immunity is induced at 3 to 5 days post infection (Zou ). Another report described that the influenza induced an early immune response 3 to 7 days after infection (Park ). However, influenza symptoms typically begin two days after exposure to the virus (Nair ). In order to understand host cell pathway hijacking in H1N1 infection within two days, transcriptome changes during the first 48 h post infection could provide important insights into the mechanisms underlying the pathogenicity of influenza A viruses. The host cell transcriptome profile was rapidly changed in the early stages of infection (within 8 h of infection). Several mediators of apoptosis were activated in this stage. After 8 h, the most significant change was found to be that in expression of genes encoding drug-metabolizing enzymes.

MATERIALS AND METHODS

Host cell preparation

A549 cells are adenocarcinoma human alveolar basal epithelial cells. The A549 cell line is widely used as an in vitro model as a type II pulmonary epithelial cell model for studies of drug metabolism and as a transfection host. A549 cells were maintained in Dulbecco’s modified Eagle’s medium (DMEM). All media contained 10% fetal bovine serum and antibiotics (penicillin streptomycin).

Viral infection and total RNA extraction

A549 cells were infected with influenza A virus at a multiplicity of infection (MOI) of 3 for 30 min at 37°C. The virus was propagated at 37°C in 11-day-old chicken embryos. After being washed with phosphate-buffered saline (PBS), cells were infected with influenza virus (IV) at an MOI of 3 for 30 min at 37°C. Infection media consisted of 0.1% glucose, 0.05% vitamin solution, and 0.5 mg/ml L-(tosylamido-2-phenyl) ethyl chloromethyl ketone (TPCK)-treated trypsin in DMEM. Infected cells were then washed and incubated in medium without serum for varying lengths of time. All infected cells were inoculated at less than 37 CT (cycle of threshold) for RT-PCR validation experiments. Total RNA was extracted from A549 cells with the RNeasy Mini Kit (74104, QIAGEN, MD, USA).

Sample QC & RNA quantification

RNA integrity was assessed using an Agilent Technologies 2100 Bioanalyzer with an RNA Integrity Number (RIN) value greater than or equal to 7. Cell-to-cell variability of influenza viral infection was quantified by subjecting the cell supernatant to a plaque assay to determine the load of infectious viral particles and then performing real-time reverse transcription quantitative PCR (RT–qPCR) on cell lysates to quantify the level intracellular viral RNA (vRNA) of individual genome segments. The size of PCR-enriched fragments was verified by assessing the template size distributions on an Agilent Technologies 2100 Bioanalyzer using a DNA 1000 chip (Agilent, CA, USA).

Next generation sequencing

The sequencing library was prepared by random fragmentation of the DNA or cDNA sample, followed by 5’ and 3’ adapter ligation. The Illumina Hiseq 4000 generated raw images utilizing HiSeq Control Software (HCS, Illumina, CA, USA) v3.3 for system control and base calling through the integrated primary analysis software Real Time Analysis (RTA, Illumina) v2.5.2. The binary BCL (base calls) value was converted into FASTQ utilizing Illumina package bcl2fastq (V2.16.0.10, Illumina).

Sequence alignment

The alignment software application “New Tuxedo Protocol” was used for transcript assembly (Pertea ). In the alignment step, HISAT2 was used to generate the alignment data (Kim ), which is a fast alignment program for mapping next-generation sequencing to read segments of the human genome (Siren ).

Transcript assembly

In the transcript assembly, StringTie (Pertea ) is known a fast and highly efficient assembler of RNA-Seq alignments into potential transcripts. It uses a novel network flow algorithm and an optional de novo assembly step to assemble and quantitate full-length transcripts representing multiple splice variants for each gene locus. StringTie estimates gene-level measurements by appropriately combining fragments per kilobase of exon per million reads (FPKMs) from the transcripts of each gene.

Time-series expression analysis

Ballgown (Frazee ) is a bioconductor-based suite that hosts a comprehensive set of bioconductor communities. The “stattest” function of Ballgown is available with a gene feature that prevents data skewing in the time series and was used for statistical tests of time-series expression. Gene-level measurements were performed by appropriately combining FPKMs from the transcripts comprising the gene (Li and Dewey, 2011) and it proposed a directed graph model combined with an expectation-maximization algorithm to estimate abundance. Ballgown generates a reference index by preparing reference transcript data and calculated RSEM values by inputting RNA-Seq data. A natural splines model with default settings was used to fit time profiles and six degrees of freedom (df) were used to fit the splines. “Timecourse=TRUE” parameters were assigned for time profiling using the natural splines function (Supplementary Table 1).

Functional annotation

DAVID bioinformatics resources consist of an integrated biological knowledgebase and analytic tools aimed at systematically extracting biological meaning from large gene/protein lists (Huang ).

Gene set variation analysis (GSVA)

GSVA is a gene set enrichment method which estimates variation of pathway activity over a sample population in an unsupervised manner. Program operation was employed using basic parameters in ‘RNA-seq’ mode and ‘kernel=TRUE’ using ‘gsva’ function (release version 3.5). The known-gene sets database was updated with the ‘Hallmark class’ and ‘c7 class’ by extending the default class (Hänzelmann ).

Weighted correlation network analysis (WGCNA)

Time-series differentially expressed genes (DEGs) were categorized by WGCNA. Categorized data were divided into 300–1000 gene categories with respect to their biological annotation-interpretation (Supplemetal Table 2). Constructing a weighted gene network entails calculation of the soft thresholding power β, to which co-expression similarity is increased to calculate adjacency. We did this with a power of 16, which is the lowest power for which the scale-free topology fit index curve flattens upon reaching a high value (in this case, approximately 0.95), a relatively large minimum module size of 30, and a moderate sensitivity (Supplementary Table 3).

Biological metabolic network in-silico validation

Elsevier’s Pathway Studio Mammalian, ChemEffect, DiseaseFX analysis: The set of biological networks in this study was created using Pathway Studio 11.0 (Ariadne) software (https://mammalcedfx.pathwaystudio.com/). Elsevier’s Pathway Studio® enables users to explore molecular interactions and cause-and-effect relationships associated with biological processes by integrating a vast knowledge base of biological relationships with analytical and visualization tools.

Ingenuity Pathway Analysis (IPA) metabolomics analysis: DEGs were analyzed using QIAGEN’s IPA software. The canonical pathways and functional processes with the most significant biological importance were identified using the list of DEGs identified with RNA-seq and the Ingenuity Pathways Knowledge Base. Pathway enrichment p-values (Fisher’s exact test) and activation z-scores were calculated using IPA.

RESULTS

The biological functions are divided into two groups approximately 8 h after infection

In hierarchical clustering, all samples were divided two groups (Supplementary Fig. 1). Early infection and late infection state samples were clustered into separate groups (Supplementary Fig. 2). Known biological function pathways were calculated using GSVA and all samples were divided into two groups by pathway activity score (Fig. 1). GSVA provides increased power to detect subtle pathway activity changes with the ‘Hallmark class’ and ‘c7 class’: ‘Apoptosis (p<0.002)′, ‘G2M checkpoint (p<0.002)′, ‘Interferon alpha response (p<0.002)′ and other biological functions were highly enriched in samples collected from cells infected for less than 8 h. ‘Arachidonic acid metabolism (p<0.001)′, ‘glutathione metabolism (p<0.001)′ and others were highly enriched in samples collected from cells infected for more than 8 h. The clustering result for whole transcriptome profile and pathway activity showed the same pattern, indicating that 8 h is a critical point for functional changes in host cells.

Fig. 1.

Unsupervised hierarchical clustering of gene set variation analysis. Gene set variation analysis (GSVA) is a non-parametric, unsupervised method for estimating variation of gene set enrichment through the expression data of infected cell. GSVA transformed each line gene-based expression matrix of gene by sample to matrix of a gene-set by sample. The matrix shows clear separation between the groups.

Changes in gene expression were observed in the two groups at 8 h

To identify host cell pathway hijacking, DEGs were selected by statistical testing of time-series expression data (Fig. 2). Gene expression changes were divided into two major groups. Fig. 2 shows 1,031 genes of the time-series DEGs, which were filtered by a threshold Q-value<0.05 (Supplementary Table 1). To find gene expression patterns in a functional module, WGCNA was used to predict groups of highly correlated genes. The results of the WGCNA analysis package analysis divided the DEGs into two groups. The biological functions are divided into two groups at approximately 8 h after infection (Supplementary Table 2). The blue module (blue line) shows increases expression after 4 h and decreases after 8 h, while the turquoise module (turquoise line) decreases at 4 h and then increases after 24 h.

Fig. 2.

Visualization of differential expression gene sets (DEGs) in the two groups. (A) early stage expression from 30 min to 8 h and (B) the late stage expression pattern, which progressively increased after 8 h of infection. (C) A network graph showing correlated genes. The graph shows clear separation between the groups.

Gene ontology (GO) enrichment analysis (Dennis ) results of these two modules are displayed in Table 1. The envelop surface glycoprotein gp120 (p<0.007) was enriched in group 1. The neuroactive ligand-receptor interaction (p<1.97E−04) and metabolism of xenobiotics by cytochrome P450 (p<2.60E−04) were enriched in group 2.

Table 1.

Time series DEGs: GO annotation and pathway menrichment analysis for all timepoints and early- and late-specific timepoints

Gene groups	Category		Term	Count	%	p-value	Fold enrichment	Bonferroni	Benjamini	FDR
Early stage changed gene sets (<8 hour)	HIV_Interaction	env	env:Envelope surface glycoprotein gp120	20	5.68	0.007	1.65	0.05	0.05	3.54
	Molecular function	GO:0015267	Channel activity	28	5.53	2.24E-05	2.49	0.01	0.01	0.03
		GO:0022803	Passive transmembrane transporter activity	28	5.53	2.34E-05	2.49	0.01	0.01	0.03
		GO:0030955	Potassium ion binding	14	2.77	4.50E-05	4.01	0.03	0.01	0.07
		GO:0022838	Substrate specific channel activity	26	5.14	9.05E-05	2.40	0.05	0.01	0.13
	KEGG Pathway	hsa04080	Neuroactive ligand-receptor interaction	19	3.75	1.97E-04	2.68	0.02	0.02	0.23
		hsa00980	Metabolism of xenobiotics by cytochrome 450	9	1.78	2.06E-04	5.41	0.02	0.01	0.24
		hsa00982	Drug metabolism	9	1.78	2.60E-04	5.24	0.03	0.01	0.30

Category refers to the GO functional category; Count refers to the number of enrichment DEGs.

Biological features of host cell before 8 h after infection: apoptosis

The protein gp120 is essential for viral entry into cells as it plays a vital role in attachment to specific cell surface receptors (Gram and Hansen, 1998). Viruses, such as HIV, with gp120 in their genome induce cellular apoptosis (Kapasi ). Apoptosis is important in an infection response and is activated at the early infection stage (Fig. 1). Supplementary Fig. 3 shows that apoptosis marker genes were expressed in early infection stage (Subramanian ). The viral gene gp120 is known to be associated with BCL-2 expression (Gaubin ). BCL-2 and other apoptosis-related genes are down-regulated in the early stages of influenza A infection. At timepoints before 8 h post-infection, expression of apoptosis marker genes, including BCL2, was decreased (Supplementary Fig. 4A).

Specific biological features after 8 h: DME-centric pathways

DMEs are the set of metabolic pathways that modify the chemical structure of xenobiotics, which are compound of foreign materials such as replicated viral RNA and other metabolites (Mackenzie ). Suppression of DMEs is necessary for viruses to maintain their presence in a cell and replicate themselves. In general, ‘endoplasmic reticulum (ER) stress response’ and ‘cytochrome P450’ are typical examples of DMEs. Both were found to be DEGs in the experimental results (Supplementary Fig. 4B, 4C).

In the Supplementary Fig. 4B, the endoplasmic reticulum (ER) stress response, also known as the unfolded protein response (UPR), is a primitive, evolutionary conserved molecular signaling cascade which has been implicated in multiple biological processes, including innate immunity and the pathogenesis of certain viral infections. Influenza A virus induces ER stress in a pathway-specific manner (Hassan ). In the Supplementary Fig. 4C, various cytochrome P450 (p< 2.06E−04) genes were differentially expressed 8 h post infection. The major role of cytochrome P450s is in steroid and bile acid syntheses (Nelson, 2009; Lorbek ). They are essential for the solubilization and absorption of lipids and fat-soluble vitamins in the intestine and they contribute to the major DME pathways by which cholesterol is excreted from the body (Lu ).

DISCUSSION

The timepoint of 8 h post infection is an important inflection point for gene expression and biological functions

Kim reported that relatively few host cell genes were differentially expressed at 8 h post infection. In this study, our data are in agreement with the report by Kim . This study indicates that 8 h post infection is an inflection point for host cell function and viral transition from the early to late infection stage.

In the early infection stage ‘Apoptosis’, ‘G2M checkpoint’, ‘Interferon alpha response’, and other biological functions are enriched. One study reported that apoptosis of virus-infected cells is one important host strategy used to limit viral infection (Scott and Norris, 2008). Cell cycle arrest may inhibit early cell death of infected cells, allow the cells to evade immune defenses, or help promote virus assembly (Bagga and Bouchard, 2014). Innate cytokine responses, such as alpha interferon (IFN-alpha) have roles in determining the rate of virus replication in the initial stages of infection (Price ). These results indicate that influenza virus naturally establishes an environment for self-reproduction in early infection stage and hijacks host cell pathways, termed ‘cytopathic viral effects’ (Heaton, 2017) (Fig. 1).

In the late infection stage, ‘Arachidonic acid metabolism’ and ‘glutathione metabolism’ were found to be enriched. Lu reported results from metabolic profiling and biochemical detection which indicated significant metabolic changes in the arachidonic acid metabolic pathway. Glutathione is partially depleted in influenza A virus infection (Alsuwaidi ). These result suggest that late expressed genes create circumstances necessary to hijack host cell pathways by inhibiting the inflammatory response and the oxidation required for metabolic breakdown of xenobiotic materials by living organisms, termed the ‘Drug metabolism pathway’ (Commandeur ).

However, driver modules in host cell function hijacking was still unclear. To identify driver function in the early infection stage, WGCNA used co-expression to predict functional modules. We predicted an early infection module (blue module) and found enriched ‘GP120 pathway in HIV’ (Fig. 3). The MAVS gene (IPS-1, VISA or Cardif) is critical for host defenses to viral infection by inducing type-1 interferons (IFN-I), though its role in virus-induced apoptotic responses has not been elucidated (Lei ). Like gp120 in HIV, the PB1-F2 gene has been shown to induce BCL2 by changing the level of MAVS gene expression (Varga ). Unlike HIV, influenza A requires a mediator such as gp120. The PB1-F2 viral gene in the influenza genome is a candidate to be a mediator of gene sets regulated by gp120 (Fig. 3). The PB1-F2 viral gene is known to induce cell death (Chen ) and apoptosis (Nicholson ), similar to gp120.

Fig. 3.

The reconstruction signaling and regulatory network of gp120 at apoptosis pathway at the early stage. Schema of the gp120 protein regulated pathway in influenza-infected cells. The genes highlighted in blue are the genes related to “glycoprotein gp120” from the enrichment gene analysis. The remaining genes were identified through network analysis (PATHWAY studio).

In the WGCNA results for late stages of infection, we predicted a late module (turquoise module) and found the ‘Drug metabolism enzyme pathway (DMEs)’ to be enriched (Supplementary Fig. 5). DMEs were additionally found to be enriched in the late stage of infection by DAVID: ‘metabolism of xenobiotics by cytochrome 450 (p<0.004)′ and ‘drug metabolism (p<0.0003)′ (Supplementary Table 3). P450 converts exotic materials to chemically reactive metabolites, which, if not detoxified, may lead to various forms of hepatic and extrahepatic toxicity, including cellular necrosis (Park ). In Supplementary Table 4, we show that CYP2C9 (P-value=7.46E−08), CYP2C18 (p-value=9.41E-06), CYP2F1 (p-value=0.0012), and CYP3A7 (p-value=9.95E−05) are involved in drug and steroid metabolism (Nelson, 2009).

The hub gene of the transcriptome-driven reconstruction network is a candidate drug targeting for influenza A

A therapeutic strategy that might escape viral resistance is to target host cellular mechanisms involved in viral replication and pathogenesis. The DGIdb integrates drug-gene interactions from 15 different sources (Wagner ) (Table 2). Reconstruction of biological networks has been vital to understanding complex biological activities (Yu ). In the gp120-regulated pathway of the early stages of infection, two genes (BCL-2, TERT) that surpassed the score thresholds (Sum score of DGIdb >2) were found. BCL2 was enriched in Paclitaxel, Oblimersen, Docetaxel, and Obatoclax treatments, while the TERT gene was only enriched in GV1001-treated cells.

Table 2.

Drug-gene interaction database for candidate genes

Gene	Drug	CFB	CFCT	DB	GPI	PG	TALC	TE	TT	TC	#of DB	# of PMID	Sum score (DB+PMID)
BCL2	PACLITAXEL	-	TRUE	TRUE	-	-	-	TRUE	-	TRUE	4	5	9
BCL2	DOCETAXEL	-	-	TRUE	-	-	-	-	-	TRUE	2	6	8
BCL2	RASAGILINE	-	-	TRUE	-	-	-	-	-	-	1	5	6
BCL2	OBLIMERSEN	-	-	-	-	-	TRUE	-	TRUE	TRUE	3	0	3
FPR1	NEDOCROMIL	-	-	TRUE	-	-	-	-	-	-	1	3	4
NR1H4	CHENODEOXYCHOLIC ACID	-	-	TRUE	TRUE	-	-	-	-	-	2	2	4
NR1H4	FEXARAMINE	-	-	TRUE	TRUE	-	-	-	-	-	2	2	4
TERT	GV1001	-	-	TRUE	-	-	-	-	-	TRUE	2	1	3
TNF	THALIDOMIDE	-	-	TRUE	-	-	-	TRUE	TRUE	TRUE	4	7	11
TNF	AMRINONE	-	-	TRUE	-	-	-	-	-	-	1	5	6
TNF	CHLOROQUINE	-	-	TRUE	-	-	-	-	-	-	1	5	6
TNF	CLENBUTEROL	-	-	TRUE	-	-	-	-	-	-	1	5	6
TNF	PSEUDOEPHEDRINE	-	-	TRUE	-	-	-	-	-	-	1	2	3

CFB: Clearity Foundation Biomarkers, CFCT: Clearity Foundation Clinical Trial, DB: Drugbank, GPI: Guide to Pharmacology Interactions Version: 04-March-2015, PG: PharmGKB - The Pharmacogenomics Knowledgebase, TALC: Targeted Agents in Lung Cancer, TE: Trends in the exploitation of novel drug targets, TT: Therapeutic Target Database, Version 4.3.02, TC: The Druggable Genome: Evaluation of Drug Targets in Clinical Trials Suggests Major Shifts in Molecular Class and Indication.

Brown et al. reported that Docetaxel should be further studied for its use in influenza vaccine production (Chen ; Brown ). Obatoclax is a novel inhibitor of endosomal acidification which prevents viral fusion and could be pursued as a potential broad-spectrum antiviral candidate (Varghese ). The TERT gene decreased rapidly at 8 h post infection (ENSG00000164362, p-value=3.34E−05). TERT is the major catalytic subunit of telomerase. TERT has additionally reported to be associated with apoptosis (Fu ).

The final goal of studying host-microbe interactions is to inform in silico models in order to identify crucial pathways, hubs, and bottlenecks, which will lead to development of new targets and strategies for prevention, diagnosis, risk assessment, and treatment of severe influenza infection in humans (Kollmus ). This study adopted functional transcriptome analysis with transcriptome sequencing and identification of genes after focusing on key time-series processes, providing several gene set enrichment analysis methods to identify genes related to seasonal influenza strains. In the future, a series of functional validation experiments need to be performed to compare more virulent influenza strains.