The responses of lungs and adjacent lymph nodes in responding to Yersinia pestis infection: A transcriptomic study using a non-human primate model

Initiation of treatment during the pre-symptomatic phase of Yersinia pestis (Y. pestis) infection is particularly critical. The rapid proliferation of Y. pestis typically couples with the manifestation of common flu-like early symptoms that often misguides the medical intervention. Our study used African green monkeys (AGM) that did not exhibit clear clinical symptoms for nearly two days after intranasal challenge with Y. pestis and succumbed within a day after showing the first signs of clinical symptoms. The lung, and mediastinal and submandibular lymph nodes (LN) accumulated significant Y. pestis colonization immediately after the intranasal challenge. Hence, organ-specific molecular investigations are deemed to be the key to elucidating mechanisms of the initial host response. Our previous study focused on the whole blood of AGM, and we found early perturbations in the ubiquitin-microtubule-mediated host defense. Altered expression of the genes present in ubiquitin and microtubule networks indicated an early suppression of these networks in the submandibular lymph nodes. In concert, the upstream toll-like receptor signaling and downstream NFκB signaling were inhibited at the multi-omics level. The inflammatory response was suppressed in the lungs, submandibular lymph nodes and mediastinal lymph nodes. We posited a causal chain of molecular mechanisms that indicated Y. pestis was probably able to impair host-mediated proteolysis activities and evade autophagosome capture by dysregulating both ubiquitin and microtubule networks in submandibular lymph nodes. Targeting these networks in a submandibular LN-specific and time-resolved fashion could be essential for development of the next generation therapeutics for pneumonic plague.


Introduction
Time-sensitive diagnosis is the foremost challenge in the management of infection with Yersinia pestis (Y. pestis) [1,2], the causative agent of pneumonic plague [3]. Rapid

Aerosolized Y. pestis exposure to non-human primate (NHP) population
The experimental protocol is described elsewhere [4]. All animal experiments were approved by the Institutional Animal Care and Use Committee (IACUC) at the Walter Reed Army Institute of Research (WRAIR), Silver Spring, MD, and were performed in a facility accredited by the Association for the Assessment and Accreditation of Laboratory Animal Care International (AAALAC). The adult male African green monkey (Chlorocebus aethiops) model used was identified previously as reliable for the study of the pathogenesis of plague [5]. Briefly, the AGMs obtained from the primate colony of WRAIR were 4.8-7.0 kg body weight and determined to be negative for tuberculosis (TB), simian immunodeficiency virus (SIV) and simian retrovirus (SRV). The typical Animal Biosafety Level 3 (ABSL-3) housing conditions and pathogen exposure protocol were reported earlier [4].
A target dose of 100 ± 50 LD 50 aerosolized Y. pestis strain CO92 was given to animals that had fasted for the previous 6 hours and were anesthetized using 4 mg/kg Telazol (Fort Dodge Animal Health, Fort Dodge, IA). Our previous report described the methods involved in bacterial inoculation and aerosol delivery [4]. Individual AGMs received Y. pestis within the range from 0.33 x 10 6 to 3.55 x 10 6 CFU via the aerosol particles with an estimated size (mean mass aerosol diameter or MMAD ± geometric standard deviation) of 1.03 μm ± 1.46.
A group of three animals was anesthetized but not exposed to Y. pestis. This control cohort was euthanized 2 h after exposure to sham aerosol. At the 78 h post-exposure time point, we euthanized one Y. pestis-challenged moribund animal, which was randomly chosen to test the lethality of the given pathogenic load. Between the Y. pestis exposure and 78 h post-exposure, sub-cohorts (N = 3) were sequentially euthanized at 6 h, 9 h, 12 h, 18 h, 24 h, 32 h and 42 h p.i. [4]. The present study is focused on 9 h, 12 h, 24 h, 32 h and 42 h p.i. results. A number of organs, including submandibular LN, mediastinal LN and lungs, were collected from every animal including the controls (Fig 1).
The protocol associated with enumeration of bacterial colonies of blood was discussed earlier [4]. The tissue samples were aseptically collected and portions approximately 1 gm in weight were removed. These portions were placed into properly-labeled, disposable sterile sample bags, heat-sealed, weighed, and placed on ice. These samples were dissociated using a handheld tissue homogenizer (Tissue Tearor, Inc.) for 120 seconds, serially diluted in sterile saline, and cultured for quantitative cultures. Congo red agar plates were incubated at 28˚C for 72 h and the CFUs were enumerated based on the geometric mean of CFU/mg of the dilution factor and colony counts (Fig 2, S1 Table).
Protein was isolated from the phenol-ethanol supernatant saved during the RNA extraction protocol [23]. Briefly, 100% EtOH was added to precipitate the protein. Protein precipitate was washed with 0.3 M guanidine hydrochloride in 95% EtOH. The protein precipitate was dissolved in 1% SDS, and isolation and purification were completed using the Qiagen protocol.

cDNA microarray hybridization and post-processing
Agilent Technologies (Agilent Technologies, Inc., CA) was the source of the Rhesus Monkey Oligo Microarray slides containing~40,000 probes. BioChain's monkey universal reference RNA (BioChain Institute, Inc., Newark, CA) was labeled with Cy3 dye, and transcribed mRNA (from either control or Y. pestis-exposed samples) was labeled with Cy5 dye using Agilent's Quick Amp Labeling Kit (Agilent Technologies, Inc., CA). The cDNA array was cohybridized using the Cy3-and Cy5-labeled RNAs, the scanned images were visualized and normalized using Feature Extraction (Agilent Technologies, Inc., CA), and the ensuing data were analyzed using GeneSpring software (Agilent Technologies, Inc., CA).

Protein ELISA
Detection and quantitation of total proteins extracted from all three organs were determined using the Thermo Scientific (Rockford, IL) Pierce BCA Protein Microplate Assay Kit. The microplate reader was set at 562 nm absorbance detection wavelength to determine the quantities of proteins present in the tissue samples. MyBioSource (San Diego, CA) monkey ubiquitin, and microtubule-associated protein Tau (MAP Tau) and NFκB ELISA kits were used. The vendor-recommended protocol and analytic  procedure were used to compare the result with a standard curve to determine the quantitative measure of the expressed proteins.

Transcriptomic assay validation
The detail of the QuantiGene-Plex protocol was described elsewhere [24]. From Thermo Fisher Scientific (Frederick, MD), we purchased the customized probes of the following NHP genes: ADNP, CYLD and TRIM9 in lungs, and USP12, USP10, SKP2 and PSMA4 for Genes clustered here showed alteration in at least one of the two early time points, 9 h and 12 h p.i. Ubiquitin-associated genes located in submandibular LN are shown in the left hierarchical tree. These genes linked to ubiquitination weren't expressed in the lungs or mediastinal LN. All of these genes are listed in Table 3, which is in the same order as that in the figure. From the pool of genes linked to microtubule network, only those were mined that ranked in the primary and secondary tiers computed by GIH algorithm. Genes showing early (9 h or 12 h p.i.) alteration in expression are clustered in the right column. All of these genes are listed in Table 4 is in the same order as that in the figure. A scale depicting the range of color corresponding to the gene expression levels is shown at right. submandibular LN. Beads linked to the probes were hybridized with total RNA selected for individual time points of corresponding tissue types. The hybridization plate was incubated for 22 hours at 54˚C, while shaking at 600 RPM in a VorTemp 56 shaking incubator. Posthybridization, the samples were captured in magnetic separation plates, washed, signal-amplified and scanned using the BioRad BioPlex 100 instrument. The BioPlex instrument settings were sample size 100 ul, timeout 60 seconds, and Bead Events/Bead region 100.
Fluorescent readings from blank wells were subtracted from fluorescent values for each mRNA of interest. These values were then normalized against the geometric mean expression of two control genes for each sample: GAPDH and β-Actin (ACTB). The array data were compared with QuantiGene Plex data in Table 1.
In addition, the longitudinal expression dynamics of a number of genes in blood were validated using the qPCR assay we published earlier [4]. The expression of genes linked to the functions relevant to this manuscript was reported, including INFγ, XIAP, ELP2, UBE2D1, ADAMTS12, NCR1, SOCS1, NEDD4, ALDH1L1 and IL6.

Statistical and functional analysis
For all transcriptional and protein analyses, the untreated controls were used as the baseline (or threshold). Transcripts altered at individual time points (moderated t-test analysis values of p < 0.01 and fold change >|2|) were mined from each organ type. Subsequent network analysis using the significantly altered genes was conducted using the Ingenuity Pathways Analysis platform (IPA, QIAGEN, Inc., CA). Networks meeting the hypergeometric threshold (p < 0.05) were considered for subsequent analysis. Functional analysis identified those genes that were associated with various networks of interest. The Molecule Activity Predictor (MAP) toolkit of IPA was used to predict the temporal regulation patterns of networks of interest. The MAP algorithm suggests the state of activation or inhibition of the network of interest by considering the observed expression changes of the molecules associated with this network.
Since a large gene set was identified linked to the microtubule network, it became essential to sort the genes based on their involvement in this biological process. As a result, we were able down-select to obtain those genes which were primarily involved with Y. pestis-mediated manipulation of the microtubule network. We used the Gene Interaction Hierarchy (GIH) algorithm [25] to sort the genes based on their involvement in the process by which Y. pestis manipulates microtubule networks. The genes were sorted based on the number of interactions they made with their neighbors within the microtubule networks domain. Candidates  50] ranked within the top 90 percentile were placed into the primary tier, and candidates ranked between 75 and 90 percentile were placed into the secondary tier. The remaining genes that made at least one interaction were pooled within the peripheral group. Some genes did not interact with their neighbors, and they were clustered as the orphan group. Hierarchical clustering of each subgroup was carried out using a Euclidian algorithm. Protein expression was calculated for every time point using the untreated control samples as the baseline. Fold change cutoff was set at ±1.5. GraphPad Prism was used for visualization and calculation of the statistical significance using Welch's t-test ( # p<0.1; � p< 0.05; �� p<0.01).

Results
AGMs were exposed to aerosolized Y. pestis and sequentially euthanized at 6 h, 9 h, 12 h, 18 h, 24 h, 32 h and 42 h p.i. Three untreated animals were used as the baseline cohort. Two animals marked for lethality testing died at 72 h-78 h p.i. Post euthanasia, we extracted 13 tissues that included blood, lungs, submandibular LN and mediastinal LN from individual NHPs. [4].
Early Y. pestis colonization detected in lungs, submandibular LN and mediastinal LN S1 Table documented the number of Y. pestis colonization after the pathogenic exposure. Our previous publication presented the enumeration of bacterial colonies of all the tissue types beyond blood, lungs, submandibular LN and mediastinal LN; in addition, we justified the selection of the threshold [log 10 (1.5 CFU/mg)] and reasons for reporting the sub-threshold number [4].
At 9 h p.i., two animals out of three showed bacterial colonization in the submandibular LN with an average load ten times higher than the threshold. Subsequent time points showed some loss of bacterial load in submandibular LN. A very consistent above-threshold bacterial colonization was observed in the lungs. In mediastinal LN, Y. pestis colonization that met the threshold was detected at relatively late time point (between 12-24 h p.i). The number of AGMs detected with positive infection in mediastinal LN were found to increase at subsequent time points. As reported earlier, blood showed further delayed detection, with Y. pestis colonization first observed to meet the threshold value at around 32 h p.i.

Global transcriptional analysis of the organs of interest
Lungs, submandibular LN and mediastinal LN were selected for protein and transcriptomic analyses. Transcriptional profiles of blood samples were presented earlier, and based on the published results, we analyzed two early (9 h and 12 h), one middle (24 h) and two late time points (32 h and 42 h). In the subsequent text, we refer to 9 h and 12 h p.i. time points together as the early phase of infection.
In principal component analyses (PCA), the unexposed animals showed a clear separation from the Y. pestis-exposed animals in these three different organs (S1

Inhibition of ubiquitin indicated in submandibular LN
Functional analysis of the gene expression data implied that the regulation of the ubiquitin network was significantly altered in submandibular LN (Table 2). S2 Table lists the gene profile enriched for the ubiquitin network. Overall, we found 42 ubiquitination-associated genes with altered expression in the submandibular LN, and the majority (62%) of encoded proteins localized to the cytoplasm. S3 Fig depicts the entire ubiquitin network enriched by the differentially regulated genes; the genes were further segregated based on the locations of the encoded proteins. There were 29 genes linked to the ubiquitin network that showed significant changes in regulation during the early time points (Table 3). Since the objective of the present study is to understand the early pathogenesis, we depicted the hierarchical transcriptomic regulation of two early time points (Fig 3).

Cross-organ gene expression profiles predicted early inhibition of microtubule activities
All three major organs at their Y. pestis ports of entry showed signs of bacterial colonization and concurrent transcriptional changes, suggesting a shift in microtubule activities (Table 2). S3 Table lists the corresponding genomic profiles across these three organs, for which 278 genes were significantly altered and associated with 821 interactions with their neighbors. GIH algorithm [25] -based gene clustering was carried out to down-select those genes which were most critically involved with Y. pestis-mediated manipulation of microtubule functions. The GIH algorithm found 20 genes ranked in the primary tier that were involved in more than 37% of total interactions. There were 44 genes ranked in a secondary tier. One hundred twenty-nine peripheral genes were involved in at least one interaction with their neighbors. In addition, there were 87 genes with no interactions, defined as orphan genes (S3 Table). S4A and S4B Fig show the distributions of the primary and secondary ranked genes across the three organs, respectively. More than 50% of these genes were located in submandibular LN. MAP analysis showed early inhibition of the microtubule network in submandibular and mediastinal LN (Table 2). Table 4 lists all the genes related to microtubule functions that ranked in the primary and secondary tiers of GIH and showed early significant alteration in at least one organ of interest. Changes in gene expression of the hierarchical clustered genes listed in Table 4 are shown in Fig 3. Furthermore, we carried out network analysis seeding of only those genes which are documented in Table 4, and the list of networks is shown in S5 Fig. S4 Table lists organ-specific genes associated with apoptosis and inflammatory response. The MAP algorithm was used with the IPA platform to compute the regulation of these networks  Table 2. Interestingly, an early onset of apoptosis and early suppression of inflammatory response were predicted in all three organs. There were 473 apoptotic genes altered in the submandibular LN across all the time points. Among these genes, 204 (42 up-and 162 down-regulated) and 184 (108 up- Roles of submandibular LN in host-Y. pestis interaction and 76 down-regulated) were regulated at 9 h and 12 h p.i., respectively. Similarly, 160 inflammatory responding genes were altered in submandibular LN across all of the time points. Among them, 81 (12 up-and 69 down-regulated) and 69 (19 up-and 50 down-regulated) genes were regulated at 9 h and 12 h p.i., respectively. We also investigated the trend of the toll-like receptor (TLR) signaling networks across the time course (S4 Table). The most prominent responses were in the submandibular LN that included 14 genes linked to the TLR signaling networks. Among them, 9 (4 up-and 5 downregulated) genes were regulated at 9 h p.i.  Roles of submandibular LN in host-Y. pestis interaction MAP Tau protein abundance in the submandibular LN never met the threshold level. Decreased MAP Tau expression in the mediastinal LN crossed the threshold of detection only at the early time points. In the lungs, MAP Tau protein abundance became increasingly suppressed with time, crossing the detection threshold at 32 h p.i.

Ubiquitin, MAP Tau and NFκB protein loads were reduced during the early phase of infection
The submandibular LN had a low abundance of NFκB protein. The level of reduction at 9 h p.i. was nearly significant. Likewise, the load of NFκB protein remained at the sub-threshold level in mediastinal LN. In the lungs, the load of NFκB protein remained less than baseline level during the early phase of infection (from 9 h to 24 h p.i).

Discussion
A recent assessment of the challenges faced by the plague vaccine development program pointed out the scarcity of NHP data in the field of understanding plague pathogenesis [26]. Mizel et al., explored a multi-organ study using a NHP model where the assays were limited to an antigen-specific humoral immune response [27]. To our knowledge, our previous study [4] was the first to interrogate the global transcriptomic profile of NHP blood to understand the early pathogenesis of pneumonic plague. While the majority of the published literature used rodents to understand the temporal progression of Y. pestis [28-30], we deemed it a logical step to investigate NHPs since they are a species that is phylogenetically and physiologically near to humans. However, it is important to note that the study of Yersinia pathogenesis using rodent models has particular value, since these species are the natural hosts of this pathogen. Furthermore, the present study derived incentive for its focus from many rodent studies that suggested the critical roles of lungs and adjacent lymph nodes in the host response to Yersinia [31, 32], including a recent study that reported the particular affinity of Y. pestis in colonizing lymph systems post-subcutaneous infection, as detected using a whole body bioluminescent tracing protocol [33].
Like humans, AGMs are good at hiding their physical discomforts until very late [4,5]. No striking clinical symptoms were noticed until one day prior to their deaths. The cornerstone of our findings was that many immunological and host-response networks were perturbed in tissues much earlier than the bacterial colonization was observed in those tissues. Significant colonization of Y. pestis in blood was noted at 32 h p.i. although sub-threshold bacterial colonization were observed starting from 12 h p.i. (S1 Table). In contrast, lungs (Fig 2A, S1  Table) and submandibular LN (Fig 2B, S1 Table) manifested significant colonization of Y. pestis by 9 h p.i. These organs are located near the respiratory tract. Hence, early colonization of Y. pestis in these organs is rather expected since the NHPs were intranasally challenged with Y. pestis. Comprehension of those molecular events that take place in these organs is likely to elucidate the early host-pathogen interactions that ultimately failed to contain Y. pestis. Notably, mediastinal LN are in close proximity of the lungs and showed colonization between 12-24 h p.i. (Fig 2C, S1 Table), whereas the transcriptomic changes in mediastinal LN were noted starting at 9 h p.i. Hence, both blood and mediastinal LN underwent molecular perturbations before bacterial colonization was observed, which showcases the advantage of determining molecular markers as the early sign of infection.
Despite the anatomic differences between the lymph nodes and lungs, the trends of their transcriptomic responses to Y. pestis infection were very similar as suggested by their PCA plots (S1 Fig). Among all of these cases, the largest transcriptomic shift from the pre-infection state happened at the earliest time point of 9 h p.i. In comparison, the molecular shifts in the subsequent time points appeared negligible. PCA plots displayed a large displacement and essentially conglomeration of the post-infection time points away from the pre-infection time point for transcriptomic regulation. Therefore, one can posit that the impact of Y. pestis assault was immediate, and this early impact potentially caused the maximum damage in these organs.
Limitations of this study included using the Rhesus Monkey Oligo array to probe the AGM transcriptomic profile. The small sample size per time point also compromised the statistical power of the study. To mitigate this risk, a longitudinal study was designed which was held together by transcriptomics-proteomics analyses. The untargeted transcriptomic investigation was streamlined to focus on the networks of interest, and stringent statistical thresholds were in place to deliver a robust molecular interpretation. Time-specific NHP controls at every time point were replaced with a more cost-efficient solution: a baseline control was euthanized at the beginning of the study. The high throughput transcriptomic results were validated by targeted proteomics and transcriptomic data. The strict regulations imposed on BSL3 labs prohibited us from isolating cells critical for adaptive and innate immunity. Hence, the present study lacks cell-specific data and histopathology data.

Early highjack of the ubiquitination process in the submandibular LN potentially compromises proteolytic activity
Ubiquitin plays a versatile role in cell biology, and facilitating the host cellular defense is one of its key roles [8,9]. Our previous study showed early activation of the ubiquitin network in the blood transcriptomic profile [4]. The ubiquitin network became suppressed at the transcriptional level in the submandibular LN by 9 h p.i. In concert, a low abundance of ubiquitin was observed in submandibular LN. The consequences could be far reaching since the efficacy of the host immune mechanism essentially depends on the proteasome-mediated disposal of ubiquitinated proteins [34]. Ubiquitins play critical roles during infection by supporting the antigen-presenting cells involved in host immunity and in activation of the NFκB family of transcription factors. Ubiquitins also participate in the degradation of whole organelles and large protein aggregates via autophagy [34], which is discussed in a later section in the context of microtubule function.
Existing literature has identified YopJ, a virulence factor of Yersinia, as an efficient deubiquitinating agent [35] enabling deactivation of NFκB-mediated host defense mechanisms. Our study found a set of proteasome-encoding genes that included PSMA1, PSMD10, PSMD14 and PSMA4, that were all down-regulated during the early phase of infection. Our functional analyses predicted that the early suppression of the toll-like receptor network and the IκB kinase/NFκB cascade network (Fig 4) potentially triggered an early deactivation of the NFκB signal. In support of this hypothesis, we found a nearly significant reduced load of NFκB protein in submandibular LN at 9 h p.i. Impaired proteocatalytic activities were further showcased by a group of down-regulated ubiquitin-linked genes, USP1, USP3, USP8, USP10, USP10, USP12 and USP44, that encode proteins associated with peptidase. Hence, there is an early indication of dysregulated ubiquitin-proteasome activity caused by Y. pestis infection.
Notably, the ubiquitin network was not significantly perturbed in mediastinal LN and lungs, and this essentially highlights the complex inter-organismal relationship [36]. Nevertheless, the trend was very clear. In submandibular LN, the ubiquitin networks remained inhibited during the early and late phases of infection (Table 2). By integrating this information with other networks closely connected to ubiquitin, such as microtubule activities, inflammation, and apoptosis, we found a more compelling picture of the pathogenesis of plague.

Early dysregulation of the microtubule network potentially limited the protein sequestration process
Typically, pathogens trigger the microtubular activities in order to enhance the local spread of infections [16]. Alternatively, destabilization of host microtubules results in impairment of the autophagosome, permitting an easy access to the cytoplasm, and ultimately, evasion of the host defense. Our study showed early down-regulation of genes linked to microtubule networks in submandibular LN, mediastinal LN and lungs. The HSP90AA1 gene encoding the chaperone protein for MAP Tau was down-regulated at 9 h p.i. In concert, the load of MAP Tau protein, a critical constituent of stable microtubules [37], was reduced during the early phase of infection in all three organs.
In the event of proteosomal impairment, host cells are essentially compelled to use alternative routes to dispose of pathogens. Hence, the pathogen or its derivative is trafficked to aggresomes via retrograde transportation undertaken by microtubules. Aggresomes are controlled by GTPase-mediated vesicular trafficking and are eventually sequestered by autophagy [38]. Thus in the event of compromised microtubular activities, autophagy becomes dysfunctional and essentially facilitates the proliferation of pathogens.
The scientific literature suggested a somewhat contrasting picture. Y. pestis was found to block autophagosomes in vitro by altering the functions of AKT, AMPK and p53 signaling [39]. An in vivo study observed a slightly different mechanism, where Y. pestis evaded the host defense by entering mouse macrophages and continuing to proliferate within the autophagosomes [40]. The NHP model revealed a number of genes associated with autophagy, such as HSP90AA1, HSPA13, MYO6 and HMGB1 that were down-regulated during the 9 h-12 h p.i. period. In contrast, LMTOR3, the gene encoding a protein responsible for mTOR signal activation, was down-regulated during the early phase of infection, indicating positive regulatory effects of autophagy. As time went by, an increasing number of genes that potentially regulate autophagy were found down-regulated. In particular, a large set of genes associated with AMPK-mediated autophagy were down-regulated, including PK1A, PRKAA1, PRKACB, CREM (all at 24 h p.i); CREB3L3, PK1A, PRKAA1 (all at 32 h p.i.); and CREB3L3, CREB, PK1A, PRKAA1, PRKACB (all at 42 h p.i.). In addition, a number of genes associated with GTPase activities, which included GIMAP2, RACGAP1, AGBL5 and RAP1GDS1, were downregulated during the early phase of infection. Overall, a temporal trend was established, highlighting a possible mechanism adopted by Y. pestis to evade autophagy.
The functional analysis of the genes critically associated with microtubule activities (i.e., the genes enriching the two top tiers classified by the GIH algorithm) found that these genes coenriched the ubiquitin network, TLR signal, and apoptosis (S5 Fig). This result is supported by the emerging evidence that interlinks microtubule-protein networks and cell death [41]. Our results possibly indicated the existence of a parallel apoptotic mechanism induced by Y. pestis that further weakened the host defense.

Immunosuppression synchronized with the early onset of apoptosis results in the comprehensive dysregulation of the host defense
Our previous study indicated that Y. pestis triggered the onset of apoptosis in blood [4]. S5 Fig highlights the possible cross-talk between apoptosis with the ubiquitin-microtubule network. Interrogation of the three organs also revealed that apoptosis may have been triggered during the early phase of infection (Fig 4). The early trend of apoptotic risk was manifested during 9 h-12 h p.i. as nearly 60% of all of the apoptotic genes were altered during this early phase of Y. pestis infection. The early signal of apoptosis in submandibular LN was synchronized with putative inhibition of TLR signal (Fig 4). This observation is particularly interesting in light of a recent study that reported Yersinia species showed more virulent pathogenesis in TLR-deficient mice [42].
We found that suppressed TLR signaling and the early highjack of ubiquitin-microtubule networks were potentially synchronized with the suppression of the host inflammatory response. Although the genes encoding the pro-inflammatory cytokines such as IL-18 were up-regulated in submandibular LN, the overall inflammatory response was putatively suppressed during the early phase of infection. As the time went by, an increasing number genes encoding pro-inflammatory cytokines such as TGFA and IL-36 and anti-inflammatory cytokines such as IL-17 became down-regulated at the transcriptomic level.

Conclusion
Over time, the pathogenesis of Y. pestis involves multiple organs and affects many biological functions. To understand this process, we explored the conceivable relationships among the seemingly distinct networks. Aligned to the emerging hypothesis that complex disease biology is essentially driven by clusters of networks [18], we hereby showed the interconnectivity among a set of networks which was perturbed in a synchronized fashion across multiple organs during a pre-symptomatic time period. The present study added a time-resolved perspective to many of those molecular events, such as ubiquitin and microtubule functions that were reported previously in the context of Y. pestis pathogenesis [13,14,35,39,40]. In fact, our previous publication [4] upheld ubiquitin and microtubule functions as two major networks potentially manipulated by Y. pestis.
Submandibular LN emerged as a potential target organ for early intervention. Past records have shown high activity in submandibular LN in the event of aerosol challenge [43,44]. A particular time lag between submandibular LN and blood molecular profiles was found. For instance, TLR signal was activated in blood during the early phase of infection, and Y. pestis was able to gradually dysregulate this network. The ubiquitin network remained active across the time course in blood. However, both of these networks were dysregulated in the early phase of Y. pestis infection in submandibular LN. Thereby our study revealed how some of the key organs potentially cross-talk during the time course of pathogenesis. For an opportunistic pathogen like Y. pestis, which proliferates rapidly and is capable of turning the key host defense mechanisms to its favor, it is essential to target the appropriate organs and molecular networks to effectively confront and counter the pathogen. Our study identified submandibular LN as an early site for the host immune system to encounter pathogens, and we mined some potential pathways that were impaired by Y. pestis to create its permissive niche. A carefully designed in vitro study that is preferably based on synthetic organ clusters could be used to validate the present multi-organ based analysis. A bio bank obtained from well-regulated clinical efforts would be another valuable source for validating the results. This knowledge can potentially drive the strategies for next generation therapeutics.

Disclaimers
The views, opinions, and/or findings contained in this report are those of the authors and should not be construed as official Department of the Army position, policy, or decision, unless so designated by other official documentation.
Research was conducted in compliance with the Animal Welfare Act and other Federal statutes and regulations relating to animals and experiments involving animals and adheres to principles stated in the Guide for the Care and Use of Laboratory Animals (NRC 2011) in facilities that are fully accredited by the Association for the Assessment and Accreditation of Laboratory Animal Care, International.
Citations of commercial organizations or trade names in this report do not constitute an official Department of the Army endorsement or approval of the products or services of these organizations.