Transcriptome Analysis Reveals Unfolded Protein Response Was Induced During the Early Stage of Burkholderia pseudomallei Infection in A549 Cells

Burkholderia pseudomallei is a zoonotic pathogen that usually affects patients' lungs and causes serious melioidosis. The interaction of B. pseudomallei with its hosts is complex, and cellular response to B. pseudomallei infection in humans still remains to be elucidated. In this study, transcriptomic profiling of B. pseudomallei-infected human lung epithelial A549 cells was performed to characterize the cellular response dynamics during the early infection (EI) stage. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed by using the online databases DAVID 6.8 and KOBAS 3.0. Real-time quantitative PCR and western blot were used for validation experiments. Compared with the negative control group (NC), a set of 36 common genes varied over time with a cut-off level of 1.5-fold change, and a P-value < 0.05 was identified. Bioinformatics analysis indicated that the PERK-mediated unfolded protein response (UPR) was enriched as the most noteworthy biological process category, which was enriched as a branch of UPR in the signaling pathway of protein processing in the endoplasmic reticulum. Other categories, such as inflammatory responses, cell migration, and apoptosis, were also focused. The molecular chaperone Bip (GRP78), PERK, and PERK sensor-dependent phosphorylation of eIF2α (p-eIF2α) and ATF4 were verified to be increasing over time during the EI stage, suggesting that B. pseudomallei infection activated the PERK-mediated UPR in A549 cells. Collectively, these results provide important initial insights into the intimate interaction between B. pseudomallei and lung epithelial cells, which can be further explored toward the elucidation of the cellular mechanisms of B. pseudomallei infections in humans.


INTRODUCTION
Melioidosis is a potentially fatal infectious disease in tropical and subtropical countries worldwide caused by the Gram-negative bacteria Burkholderia pseudomallei (Wiersinga et al., 2012). Humans and animals usually acquire this disease through broken skin, inhalation, or ingestion of the Tier-1 (top tier) select agent B. pseudomallei (Wiersinga et al., 2006). It is estimated that there are ∼165,000 human melioidosis cases and ∼89,000 deaths (54%) worldwide per year (Limmathurotsakul et al., 2016). Due to the globalization of the world's tourism and trade, the epidemic areas are expanding (Wiersinga et al., 2012(Wiersinga et al., , 2018Limmathurotsakul et al., 2016). In addition, there is no licensed vaccine for the prevention of disease (Titball et al., 2017). It is particularly important and prominent to make a profound study on the pathogenic mechanism of B. pseudomallei.
B. pseudomallei is a facultative intracellular bacterium, and it can adhere to and invade a number of mammalian cell types (Jones et al., 1996;Williams et al., 2014) and persist in vivo for many years (Mays and Ricketts, 1975). During the early infection (EI) stage, B. pseudomallei can escape from the endocytic vacuole into the cytosol, next modulate different cellular responses, and evade intracellular killing in infected cells (Willcocks et al., 2016). Previous studies have shown that B. pseudomallei modulated host iron homeostasis to facilitate iron availability and intracellular survival (Schmidt et al., 2018). It could also evade LC3-associated phagocytosis (Gong et al., 2011) and induce macrophages pyroptosis (Bast et al., 2014). Our previous research demonstrated that B. pseudomallei could interfere with the progress of Rab32 GTPase-mediated B. pseudomallei-containing phagosomes maturation and escape into the cytosol (Hu et al., 2019). B. pseudomallei could evade autophagy by regulating ATG10 in A549 cells (Li et al., 2015). Moreover, B. pseudomallei acquired the evolutionary ability to subvert autophagy by hijacking host lipid metabolism for intracellular survival (Tang et al., 2020). However, the interaction of B. pseudomallei with its hosts is more complex, and cellular response to B. pseudomallei infection in humans still remains incomplete.
Transcriptome profiling is an effective strategy for understanding the molecular events in host-pathogen interactions in recent years (Tuanyok et al., 2006;Chin et al., 2010). Previous transcriptome studies about B. pseudomallei have shown that it could up-or down-regulate numerous genes to adapt rapidly to the intracellular environment in human macrophage-like U937 cells during the EI stage (Chieng et al., 2012). B. pseudomallei could up-regulate two component signal transduction systems and a denitrification enzyme pathway for biofilm production and virulence (Chin et al., 2015;Wong et al., 2015). The small colony variant (SCV) could up-regulate many virulence and survival factors pre-and post-exposed to A549 Abbreviations: MNGCs, multinucleated giant cells; EI, early infection stage; LI, late infection stage; RT-qPCR, real-time quantitative PCR; IL-1A, interleukin-1A; IL-6, interleukin-6; TNF-α, tumor necrosis factor-alpha; p-eIF2α, phosphorylation of the alpha subunit of eukaryotic initiation factor 2; ATF4, cyclic AMP-dependent transcription factor 4; BPC006, B. pseudomallei C006 strain; DEG, differentially expressed gene; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; WHO, World Health Organization. cells (Al-Maleki et al., 2020). Previous studies are limited by a narrow dynamic expression range, and little is known about the dynamic changes of host response to B. pseudomallei infection.
Human lung epithelial cells are particularly susceptible following exposure by inhalation (Zueter et al., 2016). Multinucleated giant cell (MNGC) formation is an important characteristic feature of the intracellular life cycle of B. pseudomallei during the late dissemination stage (Whiteley et al., 2017). The persistent infection of B. pseudomallei promotes the fusion of an infected mononuclear cell with one or more neighboring cells (Kespichayawattana et al., 2000;McNally and Anderson, 2011). Identifying the differentially expressed genes (DEGs) of host cells and digging their inter-relationships, such as enriched biological processes and pathways, will help to provide insight into the combined influence of B. pseudomallei on body function. Therefore, we divided the EI stage according to no more than 50% MNGC formation and dynamically investigated the global transcriptional response of human lung epithelial cell A549 infected with B. pseudomallei using microarray analysis during the EI stage to get a better understanding of genes association. Analysis of the DEGs revealed that the response genes of A549 cells to B. pseudomallei infection were mainly involved in unfolded protein response (UPR). UPR affects host's proteins synthesis and maturation, and the finding may provide a new direction for future work to reveal the pathogenic mechanism of B. pseudomallei in manipulating UPR signaling.

Cell Lines and Bacterial Strains
Human lung epithelial cell line A549 (ATCC, CCL_185) was grown in Dulbecco's Modified Eagle Medium (DMEM, Gibco, 11965-092) containing 10% fetal bovine serum (Gibco, 10099-141) without the addition of antibiotics at 37 • C with an atmosphere containing 5% CO 2 . B. pseudomallei C006 (BPC006), as a representative of clinical virulent strain, which was the first sequencing strain in China, was used for bacterial infection in this study (Fang et al., 2012). The genome sequence of BPC006 was very close to the international standard strain K96243, and the mapped reads were ∼94% (Supplementary Table 1). BPC006 was cultured in Luria-Bertani (LB) broth at 37 • C with shaking at 200 rpm, and the overnight culture was washed twice in phosphate-buffered saline (PBS) and adjusted to an appropriate concentration by measurement of the optical density at 600 nm for infection. Live B. pseudomallei was handled under standard laboratory conditions (biosafety containment level 3).

Immunofluorescence and MNGC Formation
A549 cells were seeded in 12-well plates with 1 × 10 5 cells per well and grown on glass coverslips (NEXT, China) overnight. Subsequently, the cells were infected with B. pseudomallei at a multiplicity of infection (MOI) of 10 and incubated for 6 or 12 h. The cells were then fixed with 4% paraformaldehyde (PFA, Electron Microscopy Sciences) for 10 min and permeabilized with 0.05% Triton X-100 (Sigma, T8787). B. pseudomallei-infected fixed cells were incubated with a rabbit polyclonal antibody of B. pseudomallei (1:500) from our previous research (Hu et al., 2019) and subsequently with the Alexa Fluor 568 anti-rabbit secondary antibody (Abcam, UK) (1:2,000). The cells were stained with Actin-Tracker Green (Biotechnology, Shanghai, China). Nuclei were counterstained with DAPI (Life Technologies, 300 nM). Cover glasses were mounted on glass slides using a fluorescence mounting medium (DakoCytomation). Images were acquired on a fluorescence microscope (Zeiss, Germany) using a 40× (60× for B. pseudomallei) lens objective.
MNGCs were defined as cells containing at least three nuclei. Counting statistics of MNGCs was conducted according to the method outlined in previous reports (Whiteley et al., 2017). In brief, MNGC formation efficiency (as a percentage) was determined with a 20× objective using the following formula: (N within multinucleated giant cells / total N) × 100, where N is the number of nuclei. A minimum of 2,300 nuclei were counted for each condition. Three independent replicates were performed at each time point in the test.

Intracellular Survival of Bacteria
Intracellular survival of B. pseudomallei in A549 was estimated as previously described with some modifications (Li et al., 2015). In the initial 2 h infection with B. pseudomallei, the bacteria and cells were co-incubated for 1 h, and then the cells were overlaid with DMEM containing 250 µg/ml kanamycin for 1 h to kill extracellular bacteria. Subsequently, the cells were washed with PBS and maintained with 20 µg/ml of kanamycin in DMEM for the next 1, 4, 8, 12, 16, and 20 h. After incubation, infected A549 cells were washed three times with PBS and lysed with 1 ml of 0.1% Triton X-100 (Sigma, T8787). Cell lysates were collected and serially diluted 10-fold in PBS, and aliquots of 10 µl were plated onto LB agar to assess viable bacterial counts. We performed the experiment three replicates independently.

GeneChip TM Microarray Assay
The A549 RNA extracted from 6 to 12 h post-infection (hpi) was isolated by means of the RNeasy Micro Kit (Qiagen) according to the manufacturer's protocol for transcriptome analyses using GeneChip Human Gene 2.0 ST Arrays (Affymetrix). The control wells at 0 hpi without B. pseudomallei were used as reference group for normalization. Sample preparation for microarray hybridization was carried out as described in the Affymetrix GeneChip R Whole Transcript (WT) Sense Target Labeling Assay manual (Affymetrix, Inc., Santa Clara, CA, USA). In brief, 300 ng of total RNA was used to generate double-stranded cDNA. First, cRNA was synthesized (WT cDNA Synthesis and Amplification Kit, Affymetrix), purified, and reverse transcribed into singlestranded DNA. Purified ssDNA was then fragmented and labeled with biotin (WT Terminal Labeling Kit, Affymetrix). Finally, 2.3 µg DNA was hybridized to GeneChip Human Gene 2.0 ST Arrays (Affymetrix) for 16 h at 45 • C in a rotating chamber. Hybridized arrays were washed and stained in the Affymetrix Washing Station FS450 using Hyb, Wash & Stain Kit (Affymetrix), and the fluorescent signals were measured in the Affymetrix GeneChip R Scanner 3000-7G. Sample processing was performed at the Affymetrix Service Provider and Core Facility, "KFB-Center of Excellence for Fluorescent Bioanalytics" (Regensburg, Germany; http://www.kfb-regensburg.de). For microarray data analysis, the RMA algorithm in the Affymetrix GeneChip Expression Console software 4.0.1 was used to summarize probe signals. Then, they were exported to Microsoft Excel, and average signal values, comparison fold changes, and significance P-values were calculated. Probe sets with a fold change above 1.5-fold and a Student's t-test P-value lower than 0.05 were considered as significantly regulated according to the criteria in a previous report (Manalo et al., 2005).

Functional Annotation and Enrichment Analysis of DEGs
Gene Ontology (GO) is a functional annotation of the DEGs. The mRNAs exhibiting differential expression were entered into the integrated discovery online server (http://david.abcc.ncifcrf. gov) to be annotated and visualized. The analyses included classifications of cell constituents and molecular function, as well as biological processes, with a confidence level of 95%. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis using the KEGG Orthology-Based Annotation System (KOBAS; http://kobas.cbi.pku.edu.cn/home.do) was also conducted, and the signaling pathways were verified using the Fisher's test and an false discovery rate (FDR) of ≤0.05. Significant DEGs were used to populate human pathways in KOBAS, with statistical significance determined by Fisher's exact test or hypergeometric test.

Validation of Microarray Assay
To validate the data generated from GO analysis and pathway analysis, 11 DEGs were selected from different functional categories for real-time quantitative PCR (RT-qPCR) analysis. Supplementary Table 2 shows the respective primer sequences of selected mRNA transcripts. β-Actin and GAPDH were used as reference genes for normalization. The total RNA used for RT-PCR was extracted by TRIzol (Invitrogen Life Technologies) according to the manufacturer's recommendation. Next, complementary DNA (cDNA) was retro-transcribed from 1 µg of total RNA using PrimeScript TM RT reagent kit with gDNA Eraser (Takara, Dalian, China). CFX96 Touch Real-Time PCR detection system (Bio-Rad, USA) was used for PCR detection. The reaction system (total volume: 25 µl) consisted of the following: 12.5 µl SYBR Premix Ex Taq II (Takara, Dalian, China), 1 µl 10 mmol/L upstream primer and 1 µl 10 mmol/L downstream primer, 2 µl of ∼50 ng/µl cDNA template, and supplemented with ddH 2 O to a final reaction volume of 25 µl. The PCR program was set as 95 • C for 2 min and then 40 cycles of 95 • C for 5 s and 60 • C for 30 s. To confirm that only one PCR product was amplified, the PCR product was subjected to dissociation curve analysis at the end of each PCR reaction. Each assay was conducted in three replicates. The method of 2 − CT was used to calculate the relative expression levels of mRNAs in cells after B. pseudomallei infection, which were expressed as the relative fold change in the expression level in infected cells divided by that in control cells.
Membranes were subsequently incubated with HRP-linked antirabbit IgG secondary antibody (CST 7074S, 1:5,000) at the room temperature of 25 • C for 2 h. The bound antibodies were monitored with chemiluminescence using an electrogenerated chemiluminescence detection system (ChemiDoc XRS System, Bio-Rad, USA), and the density of each protein band was quantified using Image Lab software 6.0 (Bio-Rad Laboratories). GAPDH (CST 8884, 1:1,000) and β-actin (CST 3700, 1:5,000) were used as loading controls, and three independent replicates were performed.

Statistical Analysis
Statistical analysis of gene expression was performed by means of Student's t-test. SPSS software for Windows (IBM, V.20.0) was used. Statistical significance was stated in case of P-values being lower than 0.05.

MNGC Formation
To observe the cell morphology, cytoskeleton protein actin was stained with Actin-Tracker Green, and confocal microscope was employed. A time series assessment of A549 phenotype challenged with B. pseudomallei was shown in Figure 1A. Three independent replicates were performed for each time point in the test. There were no observable microscopic reaction lesions at 0 hpi. By 6 hpi, MNGCs had appeared with a percentage of 5.67 ± 1.15, and most of the infected cells did not fuse at this period, maintaining normal cell morphology with clear cell boundary and intact outline. However, at 12 hpi, the MNGCs had developed much larger, and the percentage increased up to 53.33 ± 4.935  Figure 1B). Actually, the infection of B. pseudomallei was a continuous process. The early and late stages of B. pseudomallei infection were not strictly defined. We divide the early stage of infection based on the formation of MNGCs no more than 50% for the convenience of research. Therefore, the results indicated that the early stage of infection could be determined before 12 hpi.

B. pseudomallei Intracellular Growth
Taking the infected time as the abscissa axis and the number of viable bacteria as the ordinate axis, growth kinetics curve was drawn in Figure 1C. The ability of B. pseudomallei to survive and replicate intracellularly demonstrated a slow increase from 1 to 12 hpi, with a growth rate of about 9 × 10 4 cfu/h. However, during 12-20 hpi, the growth of intracellular B. pseudomallei was rapid, with a high rate of about 4 × 10 5 cfu/h. Combined with intracellular B. pseudomallei in Figure 1A, the intracellular number at 6 hpi was much lower than that at 12 hpi. It should be mentioned that the morphology of the majority of infected cells at 6 hpi still remained normal. On the basis of all the above results, it could be deduced that the time of 6 hpi was similar to the lag phase of intracellular B. pseudomallei, that the bacteria were in the internalization stage, and that the time of 12 hpi was similar to the intracellular multiplication stage. We further choose the two time points for subsequent dynamic transcriptome research to appropriately reveal transcriptomic profiles of A549 with B. pseudomallei infection during the EI stage.  were used as the reference group for normalization. Applying a cut-off level of 1.5-fold change and a P-value < 0.05 on the mRNA level, we found that the cellular transcriptome profiles of duplicates were highly similar but diverged between time points, and that the number of DEGs at 6 hpi was greater than that at 12 hpi (Figure 2A). Figure 2B shows that the venn diagrams depict the overlap of DEGs from the time points of 6 and 12 hpi during the early stage of infection. A set of 306 genes at 6 hpi (Supplementary Table 3) and 103 genes at 12 hpi (Supplementary Table 2) compared with the uninfected groups (0 hpi) were observed. Thirty-six genes were involved in all two categories that indicated the 36 DEGs varied over time during the EI stage. The hierarchical clustering was performed for 306 genes (93 genes up-regulated and 213 genes downregulated) at 6 hpi and 103 genes (45 genes up-regulated and 58 genes down-regulated) at 12 hpi, using the Pearson correlation as a similarity matrix. The heat-maps were generated to visualize the quantitative differences in the expression levels of DEGs at 6 and 12 hpi compared with 0 hpi. Clustered samples in the columns and the highly up-regulated (red) and down-regulated (blue) genes in the rows were shown in Figures 2C,D, respectively.

Cellular Processes and Pathways Regulated by B. pseudomallei Infection
To explore cellular responses to B. pseudomallei infection and characterize DEG functions, we predicted enriched biological processes for up-regulated gene set and down-regulated gene set against the GO database, respectively, as a function of time following B. pseudomallei infection (Figure 3). No matter it was at 6 or 12 hpi, the most noteworthy enrichment biological process category was the PERK-mediated UPR with P-value = 2.48 × 10 −5 at 6 hpi and P-value = 2.75 × 10 −4 at 12 hpi in the upregulated gene set, including three common DEGs of EIF2S1, HSPA5, also known as GRP78 or BiP, and ATF4. The expression levels of these three DEGs at later time points were comparable to or higher than those at earlier time points (Table 1), suggesting that PERK-mediated UPR was induced and became stronger over time during the EI stage.
Certainly, many enrichment biological processes of DEGs in both up-and down-regulated gene sets pointed to inflammatory response, which included the processes of cellular response to tumor necrosis factor (TNF), cellular response to lipopolysaccharide, positive regulation of interleukin-8 and interleukin-6 production, and cellular response to interleukin-1.  The relative expression level of each mRNA transcript was represented as the n-fold change relative to NC (0 hpi). EIF2α and ATF4 were verified in PERK-mediated unfolded protein response; IL1A, IL6, TNFAIP3, NFκB2, and GKN2 were verified in inflammatory responses; ANXA9, CLDN7, and ICAM1 were verified in cell migration; and BCL2L15 were verified in apoptosis. All data are presented as mean ± SD (n = 3 in each time point). (L) The correlation of expression level between GeneChip TM microarray and RT-qPCR was high and significant (Pearson's r = 0.93; P < 0.0001).
including the signaling pathway of protein processing in the endoplasmic reticulum (ER) (Figure 4A). Differentially expressed proteins of BiP (also known as HSPA5/GRP78), eIF2α, ATF4, Hsp70 (including HSPA1A, HSPA1B, and HSPA2), HERPUD1, and EDEM2 were enriched in this signaling pathway, as shown in Figure 4B. What caught our attention was that the signaling pathway of protein processing in the ER was mainly induced by the branch of PERK-mediated UPR. The other top enrichment pathways were mainly involved in inflammatory response, including TNF signaling pathway and some viral-like pathogens infection signaling pathways, which were consistent with the results of GO enrichment analysis. In general, our results inferred that B. pseudomallei triggered the integration response of the host cell, but PERK-mediated UPR in the pathogenesis of B. pseudomallei should be paid more attention.

Validation of Microarray Data and UPR
To validate the expression results observed in microarray, we further randomly selected 11 genes for verification by means of RT-qPCR. Selected genes were, respectively, involved in PERK-mediated UPR (eIF2α and ATF4), inflammatory responses (IL1A, IL6, TNFAIP3, NFκB2, and GKN2), cell migration (ANXA9, CLDN7, and ICAM1), and apoptosis (BCL2L15) in Table 1. As summarized in Figure 5, our RT-qPCR experiments indicated that the mRNA expression level of the 11 selected genes was consistent with the expression of microarray The relative expression level of DEGs in every group. Cells without infection (0 hpi) was the control group. A t-test was performed to compare the means, and the asterisks represent the P-value as follows * <0.05, n.s = non-significant, and P-value > 0.05 is considered not statistically significant.
( Figures 5A-K), with the Pearson correlation r = 0.93 and Pvalue < 0.0001 (Figure 5L), which indicate that the microarray data were reliable. The PERK-mediated UPR is one branch of UPR in the signaling pathway of protein processing in the ER (Bettigole and Glimcher, 2015). It will be activated to maintain protein-folding homeostasis in the ER. Concomitantly, mRNA translation was transiently attenuated through the phosphorylation of eIF2α leading to ATF4 activation (Hotamisligil, 2010). To investigate B. pseudomallei-activated PERK arm, BiP, PERK, eIF2α/p-eIF2α, and ATF4 were selected for verification via western blotting analysis, and the result was shown in Figure 6. With B. pseudomallei infection prolonged in the EI stage, the expression level of the ER chaperone proteins BiP, PERK, and ATF4 was verified to be significantly increased. The relative expression of BiP, PERK, and ATF4 was 2.20, 3.82, and 1.68 at 6 hpi, respectively, and then it went up to 2.58, 7.60, and 3.38 at 12 hpi, respectively. There was no significant change about the expression level of eIF2α, but the phosphorylation of eIF2α, the key to the activation of PERK-mediated UPR (Pavio et al., 2003), was verified to be significantly increased with the relative expression from 1.40 to 1.74. Our present data suggested that B. pseudomallei infection induced PERK-mediated UPR in human lung epithelial cells.

DISCUSSION
The intracellular lifestyle of B. pseudomallei in host cells has been reported to impact the severity of melioidosis, ranging from acute fatal sepsis to chronic infection with or without clinical symptoms (Pilatz et al., 2006;Lazar Adler et al., 2009;Allwood et al., 2011). However, the pathogenic mechanism of B. pseudomallei is still not fully understood. Elucidating the cellular response to B. pseudomallei infection will aid in unraveling the pathogenic mechanisms of the pathogen. We previously reported that B. pseudomallei could modulate autophagy (Li et al., 2015) and Rab32-mediated phagocytic vesicle transport to evade intracellular killing in host cells (Hu et al., 2019).
Although previous studies have been carried out to determine the host transcriptional profiling of B. pseudomallei infection (Perumal Samy et al., 2015;Aschenbroich et al., 2019;Al-Maleki et al., 2020), limited by a narrow dynamic expression range, little is known about the dynamic changes of host response to B. pseudomallei infection. To further understand the cellular response dynamics caused by B. pseudomallei, we conducted a transcriptomic profiling of B. pseudomallei-infected A549 cells at 6 and 12 hpi during the EI stage. Our results generated a number of interesting observations. Upon B. pseudomallei infection, a set of 36 common genes that varied over time from 6 to 12 hpi were found. These DEGs are consistent with previous studies on the infection of B. pseudomallei, which involved inflammatory response (Wongprompitak et al., 2009;Perumal Samy et al., 2015), cell migration (Williams et al., 2014), apoptotic process (Hseu et al., 2014), and so on. The results indicated that these continuously regulated DEGs might be a sign of B. pseudomallei infection and played a key role in pathogenicity. Interestingly, the PERK-mediated UPR was enriched as the most noteworthy biological process category.
UPR is a cytoprotective response that prevents cytotoxicity effects caused by cellular accumulation of unfolded or misfolded proteins, which also invokes innate immune signaling in response to invading microorganisms (Celli and Tsolis, 2015). It is reported that various intracellular bacteria trigger the UPR to their advantage of replication during infection, such as Mycobacterium tuberculosis (Seimon et al., 2010), Brucella melitensis (Smith et al., 2013), and Listeria monocytogenes (Pillich et al., 2012). In our study, the expression of BiP, an ER chaperone protein, was increased as B. pseudomallei infection progresses, and B. pseudomallei infection increased PERK expression (Figure 6). This might be due to the accumulation of unfolded protein in the ER, resulting in cells needing more BiP to compensate for the dissociation of BiP from the transmembrane protein located in the ER, thereby helping the unfolding and assembly of protein structures (Bertolotti et al., 2000). Activated PERK could phosphorylate its target of the translation initiation factor eIF2α (Harding et al., 1999). Consistent with previous results, we confirmed that p-eIF2α increased with the time of B. pseudomallei infection. The p-eIF2α induced the activation of ATF4, which led to the transient weakening of gene translation involved in unfolded protein reaction (Choi and Song, 2019). The results showed that B. pseudomallei infection increased ATF4 expression (Figure 6). In conclusion, our present data suggested that B. pseudomallei induced PERK-mediated UPR in human lung epithelial cells. It was reported that UPR was exploited by energy and metabolically deficient pathogens to increase ATP and metabolites required for development, certain microbial agents, including viruses and bacteria, induce UPR (Celli and Tsolis, 2015), and it may be required for Brucella replication (Smith et al., 2013). Whether B. pseudomallei could hijack PERKmediated UPR for intracellular survival and replication remains to be delved deeper. Moreover, the expression levels of upregulated genes at later time points were comparable to or higher than those at earlier time points (Figure 6), suggesting that PERK-mediated UPR to B. pseudomallei infection became stronger over time. Although the UPR plays a major role in certain microbial infectivity, its role in B. pseudomallei pathogenesis is unknown. Other enrichment cellular processes, such as inflammatory responses, cell migration, and apoptosis, were enriched too. In general, our results inferred that B. pseudomallei triggered an integrated host cell response and the role of the UPR should be paid more attention to in the pathogenesis of B. pseudomallei.
Inflammatory response was also enriched according to microarray assay, as seen in Figures 3, 4. Microarray technology for transcriptome yields a large amount of data, and it is important to validate differential expression by independent methods, which conformed to the results of microarray assay. RT-qPCR analysis of DEGs included in inflammatory response, such as IL1A, IL6, TNFAIP3, NFκB2, and GKN2 (Figure 5), was therefore conducted. B. pseudomallei increased the mRNA expression of TNFAIP3 and IL-6 in human lung epithelial cell, which was consistent with the current research that proinflammatory mediators of TNF-α and IL-6 increased in the liver of mice following infection with B. pseudomallei and in serum samples of melioidosis patients (Ulett et al., 2000;Koosakulnirand et al., 2018). Activated UPR transcription factors can induce inflammatory response, which play an important role in the pathogenesis of inflammation (Bettigole and Glimcher, 2015). Previous studies showed that intracellular bacteria, such as M. tuberculosis, could induce p38 MAPK activation by TLR2 and TLR4, which in turn caused the accumulation of misfolded or unfolded TNF-α in the ER, and further induces the induction of macrophage ER stress (Oh et al., 2018;Choi and Song, 2019). Similarly, the MAPK signaling pathway was found in the process of B. pseudomallei infection in our manuscript (Figure 4A), and this might also be a possible way for B. pseudomallei to activate inflammatory response.
Additionally, our data yield that B. pseudomallei exposure could active the p-eIF2α/ATF4 pathway that was verified by western blotting in Figure 6. Accumulating evidences suggested that the p-eIF2α/ATF4 pathway mediated the induction of a gene expression program that referred to integrated stress response (ISR). ISR is a common adaptive pathway to restore cellular homeostasis, which could be triggered by integrating various types of stress signals, including ER stress, amino acid deprivation, virus infection, and oxidative stress (Harding et al., 1999;Zhang et al., 2002;Yerlikaya et al., 2008), while it is not clear whether this response benefits the host or the pathogen. However, considering the steady growth of intracellular B. pseudomallei at the early stage of infection in Figure 1C, it is likely that B. pseudomallei has evolved to modulate the ISR to its advantage during infection, but it still needs sufficient experimental evidence to reveal to what extent B. pseudomallei uses to exploit this intracellular niche to promote their replication. Additionally, the current evidence that bacteria may regulate UPR comes from the research on bacterial toxins. For instance, subtilase SubAB from Shiga-toxigenic Escherichia coli (STEC) can specifically degrade the ER chaperone BiP and further activated IRE1, ATF6, and PERK, the three known UPR signaling arms (Morinaga et al., 2008). B. pseudomallei has a strong virulence system (Wiersinga et al., 2006), whether UPR can be activated by it is not clear. Future work is likely to reveal to which novel mechanisms are employed by B. pseudomallei to manipulate UPR signaling.

CONCLUSIONS
In conclusion, this study has elucidated the dynamic DEG profiles following B. pseudomallei infection of A549 cells during the EI stage. Interestingly, a set of 36 common genes varied over time were found, and the PERK-mediated UPR was enriched as the most noteworthy biological process category. We have verified the DEGs with RT-qPCR and further confirmed that B. pseudomallei induced the expression of Bip (GRP78), PERK, p-eIF2α, and ATF4, the key proteins involved in the UPR. Overall, our analysis suggests that B. pseudomallei infection induced the PERK-mediated UPR in protein processing in the ER. Therefore, the data reported here are valuable resources for practical use to further characterize the mechanisms that are potentially involved in the intracellular living of B. pseudomallei.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
CR, CM, and YX are fully responsible for the study, including the study design, data collection, data analysis and writing. MZ, SY, WY and ZH help with experiment implementation and data collection. JY and XC provide experimental technical support. LD provides research funding management and reimbursement support. XM provides guidance and supporting contribution. QL and YL bear the fund, design, writing and critically reviewing the intellectual content of the article. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
This manuscript has been released as a pre-print at Research Square . We thank Gminix (Shanghai, China) for its assistance with the GeneChip TM microarray assay. Affymetrix microarray data are available as Excel files that we upload in the additional materials (Supplementary Table 6). NCBI accession numbers are not available due to a damaged CEL file.