The Transcriptomic and Bioinformatic Characterizations of Iron Acquisition and Heme Utilization in Avibacterium paragallinarum in Response to Iron-Starvation

Avibacterium paragallinarum is the pathogen of infectious coryza, which is a highly contagious respiratory disease of chickens that brings a potentially serious threat to poultry husbandry. Iron is an important nutrient for bacteria and can be obtained from surroundings such as siderophores and hemophores. To date, the mechanisms of iron acquisition and heme utilization as well as detailed regulation in A. paragallinarum have been poorly understood. In this study, we investigated the transcriptomic profiles in detail and the changes of transcriptomes induced by iron restriction in A. paragallinarum using RNA-seq. Compared with the iron-sufficiency control group, many more differentially expressed genes (DEGs) and cellular functions as well as signaling pathways were verified in the iron-restriction group. Among these DEGs, the majority of genes showed decreased expression and some were found to be uniquely present in the iron-restriction group. With an in-depth study of bioinformatic analyses, we demonstrated the crucial roles of the Hut protein and DUF domain-containing proteins, which were preferentially activated in bacteria following iron restriction and contributed to the iron acquisition and heme utilization. Consequently, RT-qPCR results further verified the iron-related DEGs and were consistent with the RNA-seq data. In addition, several novel sRNAs were present in A. paragallinarum and had potential regulatory roles in iron homeostasis, especially in the regulation of Fic protein to ensure stable expression. This is the first report of the molecular mechanism of iron acquisition and heme utilization in A. paragallinarum from the perspective of transcriptomic profiles. The study will contribute to a better understanding of the transcriptomic response of A. paragallinarum to iron starvation and also provide novel insight into the development of new antigens for potential vaccines against infectious coryza by focusing on these iron-related genes.


INTRODUCTION
Infectious coryza, caused by Gram-negative bacterium Avibacterium paragallinarum, is a severe respiratory disease of chickens that brings enormous economic losses for the poultry industry worldwide (Blackall, 1999). A. paragallinarum is classified into three serovars-A, B, and C-according to the Page schemes. In South Africa, it is considered as one of the most serious diseases, with C3 being the predominant serovar (Bragg et al., 1996;Boucher et al., 2014). In China, serovars A, C, and B were identified for the first time in 1987, 1995, and 2003, respectively (Zhang et al., 2003Sun et al., 2018;Xu et al., 2019). Several factors are responsible for the pathogenicity of A. paragallinarum. Among them, hemagglutinin (HA) has been regarded as the most key virulence factor because of its participation in tissue adhesion (Barnard et al., 2008). Vaccination remains the most effective preventive measure against A. paragallinarum. However, vaccination failures always occur with the emergence of new serovars or serovar variants (Blackall, 1999;Bragg, 2005). In addition, the advent of multidrug bacterial resistance is deemed one of the most alarming situations worldwide. Therefore, it is urgent for the development of novel broad-spectrum prophylactic agents against A. paragallinarum. To date, the molecular mechanisms of colonization, growth, and virulence of A. paragallinarum in chicken upper respiratory tract are still poorly understood.
Iron is an essential micronutrient for most bacteria, which plays essential roles in bacterial colonization, growth, and virulence. Unfortunately, the concentration of free iron in the host is not enough to support the growth of bacteria (del Rio et al., 2005). In order to achieve effective iron homeostasis, bacteria have evolved sophisticated iron-acquisition systems to scavenge iron from their surroundings (siderophores, hemophores, or host-molecule-binding proteins) to ensure adequate supplies (Parrow et al., 2013). The siderophore transport system, which is capable of acquiring free inorganic iron, consists of both siderophores and TonB-dependent transferrin receptors (TBDR) on the outer membrane of bacteria (Chu et al., 2010). In A. paragallinarum, Abascal et al. (2009) have identified three upexpressed outer membrane proteins acting as transferrin receptor and iron transport proteins in bacteria cultured in iron-restricted medium, which are regulated by the Fur protein (Abascal et al., 2009). They speculate that these proteins are associated with the colonization of bacteria in an iron-limited host environment. Besides transferrin and lactoferrin, enterobactin is also proved to be a crucial siderophore. Studies have confirmed that the growth ability of many bacteria has been enhanced in the presence of enterobactin, including not only enteric bacteria such as Escherichia coli and Campylobacter jejuni, but also respiratory bacteria such as Hemeophilus parainfluenzae and Bordetella pertussis (Williams et al., 1990;Thulasiraman et al., 1998;Anderson and Armstrong, 2004;Horiyama and Nishino, 2014). In terms of ferric enterobactin receptor, we have previously demonstrated that CfrA and CfrB contribute to iron acquisition in Campylobacter (Zeng et al., 2009;Xu et al., 2010). Until now, the mechanisms of iron acquisition by siderophores and the detailed regulation in A. paragallinarum are still on the preliminary stage and remain to be further explored.
Remarkably, the availability of free iron is strongly limited in vertebrate hosts, causing insufficient iron acquisition by siderophores, especially for pathogenic bacteria (Sekine et al., 2016). Thus, the heme utilization system has been well developed in bacteria to better obtain iron. Over the past decades, the utilization of heme is validated to be a common mechanism employed by pathogens. For example, heme is preferred as an iron source over ferrous iron for some Gram-positive bacteria such as Staphylococcus aureus and Streptococcus pyogenes (Rouault, 2004;Skaar et al., 2004;Zhu et al., 2008). For Gramnegative pathogens, proteins encoded within the heme transport operon have been studied such as Vibrio anguillarum (Sekine et al., 2016). To acquire heme from hemoproteins, bacteria have developed sophisticated heme acquisition machineries. Gramnegative pathogens can also utilize an outer TonB-dependent membrane receptor protein to transport heme into the periplasm, and then periplasmic heme-binding proteins transport heme to the cytoplasm (Contreras et al., 2014;Richard et al., 2019). In the cytoplasm, heme can be degraded by heme-degrading enzymes, which remove iron from the heme (Sekine et al., 2016). Nevertheless, the research on heme uptake and utilization by A. paragallinarum is rare.
Small regulatory RNAs (sRNAs) act as post-transcriptional regulators controlling bacterial adaptation to environmental changes (Koeppen et al., 2016;Gerrick et al., 2018). Their regulatory functions have become the center of attention in a variety of bacteria, such as E. coli, Salmonella, and Pseudomonas aeruginosa (Koeppen et al., 2016;Chareyre and Mandin, 2018;El Mouali et al., 2018). sRNAs play crucial roles in regulating outer membrane protein expression, iron ion balance, community induction, and bacterial pathogenicity. The regulatory mechanism involves binding to target bacterial mRNAs or directly interact with target proteins (Repoila and Darfeuille, 2009). For instance, sRNAs RyhB of enteric bacteria and PrrF of P. aeruginosa are highly upregulated during iron starvation, which are able to suppress mRNAs encoding ironcontaining proteins (Masse et al., 2005;Reinhart et al., 2015). Recently, MrsI is identified to be important for iron-sparing response in Mycobacterium tuberculosis (Gerrick et al., 2018). To date, sRNA studies in A. paragallinarum have been limited, let alone the interaction between any sRNA and its targets.
The identification and quantification of differentially expressed genes (DEGs) can be conducted at a transcriptomic level by using RNA sequencing (RNA-seq) technology. This transcriptomic approach has identified 13 upregulated genes encoding for proteins exposed to iron-restriction stress at the bacterial surface of Hemeophilus (Glässerella) parasuis (Alvarez-Estrada et al., 2018). In addition, RNA-seq can also be used to characterize sRNAs, including their secondary structure and predicted mRNA targets. Koeppen et al. (2016) have characterized differentially packaged sRNAs in outer membrane vesicle secreted by P. aeruginosa. Nevertheless, a comprehensive analysis of the transcriptomic changes in A. paragallinarum to iron-restriction stress has yet to be performed.
To get a deeper understanding of the biological pathways and main transcripts involved in iron-restriction stress as well as exploring the specific molecular mechanism behind iron acquisition and heme utilization in A. paragallinarum, we used RNA-seq for gaining information about gene and sRNA abundance in A. paragallinarum under iron-restriction condition. To the best of our knowledge, this is the first report regarding the transcriptomic and bioinformatic characterizations of iron acquisition and heme utilization in A. paragallinarum in response to iron starvation.

Bacterial Strain and Growth Conditions
The strain 3005 of the A. paragallinarum serogroup C was used in this study, which was isolated from China in 2018. As outlined in Supplementary Figure 1, tryptic soya broth (TSB) and tryptic soya broth agar (TSA) that were added with 10% chicken serum and 0.0025% nicotinamide adenine dinucleotide (NAD) were used for the propagation and maintenance of the strain. Then, A. paragallinarum was grown under in vitro culture conditions with iron restriction (R) or without iron restriction (C, as a control). For the ironrestriction condition, A. paragallinarum cells were inoculated at 1:100 in 200 ml minimal essential medium alpha (MEMα) (Invitrogen) with 0.0025% (w/v) NAD (Palyada et al., 2004). For the control condition, an extra 40 µM FeSO 4 was added in the cultures based on the related reference (Palyada et al., 2004). The cells were grown in a shaking incubator at 37 • C until reaching an optical density (OD) of 0.6-1.0 at 600 nm (exponential phase of growth), and samples were harvested by centrifugation, followed immediately by being resuspended in RNA later (Ambion, Carlsbad, CA, United States). Samples were then stored at −80 • C. Three replicates were made under each condition.

RNA Extraction and Sample Preparation
Total RNA for RNA sequencing was extracted using TRIzol R reagent (Invitrogen) according to the manufacturer's instructions, and the contaminating DNA was removed by treatment with DNase I (TaKaRa) and confirmed by negative PCR amplification results. Finally, RNA quality was determined by 2100 Bioanalyser (Agilent) and quantified using the ND-2000 (NanoDrop Technologies).

Library Preparation and Illumina Sequencing
Two micrograms of total RNA from each sample was used for RNA-seq transcriptome library preparation. After ribosomal RNA (rRNA) was depleted, mRNA was purified and fragmented. Subsequently, first-and second-strand cDNAs were synthesized and ligated with adapters. cDNA library preparation was performed sequentially with the TruSeq RNA Sample Prep Kit from Illumina (San Diego, CA, United States) following the manufacturer's protocol. After quantified by TBS380 (Picogreen), paired-end RNA-seq sequencing library was sequenced with the Illumina HiSeq × 10 (2 × 150 bp read length). The processing of original images to sequences, base-calling, and quality value calculations were performed using the Illumina GA Pipeline (version 1.6), in which 150 bp paired-end reads were obtained. A Perl program was written to select clean reads by removing low-quality sequences, reads with more than 5% of N bases (unknown bases), and reads containing adaptor sequences.

Mapping and Alignment of Sequence Reads
After the low-quality reads were trimmed and the adapters were removed, the clean reads were mapped and aligned onto the whole reference genome (A. paragallinarum strain AVPG2015, GenBank accession nr NC_011852.1) using Bowtie 2.0.

Differential Gene Analysis of RNA-Seq
The raw data were processed for the gene data analysis using the normalized transcripts per million mapped reads (TPM) values using RSEM, and differential gene expression was obtained by DESeq2. Based on the negative binomial distribution model, heterogenetic analysis among samples was obtained by controlling the fold change (FC) and p value statistical methods. The screening conditions of this study were | log2FC| ≥ 1 and p < 0.05. An analysis of variance (ANOVA), using the Benjamini-Hochberg multiple testing correction, was performed to identify genes significantly differentially expressed (p < 0.05). Significantly differential expressed genes with fold change ≥ 2 were subjected to further Gene Ontology (GO) and pathway analysis.

Bioinformatics and Statistical Analyses
Differentially expressed genes were classified for the categories using annotations of GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways with the DAVID online software 1 . The threshold of significantly enriched genes was determined according to the corrected p (Benjamini) (p ≤ 0.05). The default parameters (hypergeometric algorithm, Benjamini-Hochberg multiple testing correction) were performed to take overrepresentation analyses. GO terms (biological process, cellular component, and molecular function) and KEGG pathways terms with p-Values of ≤0.05 were considered to be significantly enriched with DEGs. The heatmaps of DEGs were achieved with the pheatmap packages from R.

sRNA Analysis
Small regulatory RNAs prediction was obtained by the current mainstream prokaryotic analysis software Rockhopper 2 , which was mainly based on the base sequencing coverage information. Besides, the RNAfold program in the Vienna RNA packages was used to predicate the secondary structure of sRNA based on the minimum free energy.

Data Deposition
The raw sequencing files from this study were deposited into the NCBI SRA repository and the SRA accession number was SRR13291129.

Validation of DEGs by Quantitative Real-Time PCR
Differentially expressed genes identified by the above-described method were validated using quantitative real-time PCR (qPCR). The gyrA was used as a reference control (Wen et al., 2016). The information on the primers is listed in Supplementary Table 1. The qPCR programs for VY92_RS06600, VY92_RS09125, VY92_RS09790, VY92_RS03730, and VY92_RS03735 were as follows: 1 cycle of 95 • C for 3 min; 35 cycles of 95 • C for 30 s, 60 • C for 30 s, and 72 • C for 40 s. The qPCR programs for VY92_RS01655, VY92_RS00335, VY92_RS05725, and gyrA were as follows: 1 cycle of 95 • C for 3 min; 35 cycles of 95 • C for 30 s, 56 • C for 30 s, and 72 • C for 40 s. Gene expression levels were calculated as previously described (Huo et al., 2020). Data analysis was conducted by using two-way analysis of variance with GraphPad Prism (ver. 5.0). p < 0.05 represents statistical significance. The results were shown as the means ± standard deviations of three independent experiments.

Quality Control of Sequencing Data and Mapping
As outlined in Supplementary Figure 1, the cells were grown and harvested at 5 h for RNA sequencing based on the measurement of OD values via bacteria growth assay (Supplementary Table 2). Then, the quality control of sequencing data was taken. As shown in Table 1, 13.70 to 16.68 M raw reads were generated for each sample, which the RawQ20 and Raw30 ratios were all above 90%. High-quality reads were obtained and used for further analyses after quality assessment using the FastQC tool and data preprocessing. After the quality control, clean reads with numbers averaging between 13.00 and 15.87 M were obtained per sample, in which the CleanQ20 and CleanQ30 ratios reached more than 95%. Besides, we also assessed the nucleotide distribution of both raw and clean data and found that sequences were all roughly distributed at about 150 bp (Supplementary Figure 2).
In order to quantify the gene expression in samples, sequence mapping and alignment analysis of high-quality clean reads was taken at first by comparing with the reference genome (A. paragallinarum strain AVPG2015, GenBank accession nr NC_011852.1). As shown in Table 2, more than 80% of reads from all experimental groups 14 R_1, R_2, and R_3 represent the three independent replicates in iron-restriction group, C_1, C_2, and C_3 represent the three independent replicates in control group. Q20 = bases of Q ≥ 20/all bases of sequencing. Q30 = bases of Q ≥ 30/all bases of sequencing. R_1, R_2, and R_3 represent the three independent replicates in iron-restriction group, C_1, C_2, and C_3 represent the three independent replicates in control group.
Frontiers in Microbiology | www.frontiersin.org could be aligned to the reference genome, and unique mapped reads ratio ranged from 46.32 to 76.1%. Besides, 81.69 to 87.63% of these clean reads were distributed in coding regions. The distribution and numbers of gene expression are shown in Figure 1A. In the density diagram of gene expression, high-and low-abundance genes could be shown and formed an obvious peak in each group. Furthermore, the TPM of genes was shown to be similar between these two groups, indicating the stability of the expression levels. Thus, the sequencing data conformed to the requirements of subsequently bioinformatics analysis.

Identification and Quantification of Differentially Expressed Genes
According to the quantification of gene expressions, we analyzed 2,566 transcripts, including 2,445 genes encoding proteins and 121 genes encoding sRNAs (Supplementary Data Sheet 1). By taking the Venn analysis, we found that there were 2,062 common genes in these two groups ( Figure 1B). In addition, 30 genes were uniquely expressed in A. paragallinarum under iron-restriction condition and 26 genes were uniquely expressed in A. paragallinarum under the control condition, respectively. Based on the criteria of | log2FC| ≥ 1 and p-adjust < 0.05, the numbers of genes that up-and downregulated their abundance were summarized. As shown in Figure 2 and Supplementary Data Sheet 1, 122 DEGs were identified significantly between the iron-restriction group and the control group. Compared with the iron-sufficiency control group, iron restriction could significantly downregulate 68 genes, including 65 genes encoding for proteins and 3 genes encoding for sRNAs, which include autotransporter outer membrane beta-barrel domain-containing protein and other proteins. Meanwhile, 54 genes were significantly upregulated in the iron-restriction condition than in the control, which consisted of 50 genes encoding for proteins and 4 genes encoding for sRNAs. Among these 54 differentially upregulated genes, 19 genes could directly evolve into sophisticated iron-acquisition systems, such as iron utilization receptors, ABC transporters, and TonB energy systems (Table 3). Additionally, several iron-regulated operons were identified for the gene clustering analysis, suggesting that a number of differentially up-expressed genes might be cotranscribed. There were several operons involved in ABC transporters: one that contained hemoglobin/transferrin/lactoferrin receptor (VY92_RS00325) and ABC transporter ATPbinding protein (VY92_RS00330 and VY92_RS00335), one that contained Fe(3+)-hydroxamate ABC transporter permease FhuB (VY92_RS05675) and iron-hydroxamate transporter ATP-binding subunit (VY92_RS05685), and another one that contained ABC transporter ATP-binding protein (VY92_RS05745) and ABC transporter permease (VY92_RS05750). An operon was shown to be involved in heme utilization systems that contained heme utilization protein HutZ (VY92_RS03730), heme utilization cytosolic carrier protein HutX (VY92_RS03735), and TonB-dependent receptor (VY92_RS03740). An operon was shown to be involved in iron transporter and iron permease that contained iron transporter (VY92_RS05765), iron permease (VY92_RS05770), and ironsulfur cluster carrier protein ApbC (VY92_RS05775). There was also an operon that was involved in TonB energy systems, containing energy transducer TonB (VY92_RS11145) and TonB system transport proteins ExbD and ExbB (VY92_RS11150 and VY92_RS11155). Interestingly, 98 genes encoding for   numerous domains of unknown function (DUF) domaincontaining proteins were discovered in A. paragallinarum. Among these genes, there were 57 (58%) upregulated and 17 (17%) downregulated DUF genes in the iron-restriction group. Furthermore, significant differences of gene expressions were shown in three upregulated DUF genes, which encoded DUF454, DUF4198, and DUF2318 domain-containing proteins, respectively. There was also one DUF gene that encodes DUF1345 domain-containing protein that was significantly downregulated in the iron-restriction group than in the control group.
In the case of 115 significantly DEGs encoding proteins, the hierarchical cluster analysis of these DEGs was carried out for understanding the relationships between genes across all 6 samples (Supplementary Figure 3). Above all, it suggests that iron restriction could result in significant changes of gene expressions in A. paragallinarum and seems to be likely to suppress expressions of the majority of genes in bacteria.

Functional Annotations and Bioinformatic Characterization of DEGs in A. paragallinarum Under Iron-Restriction Condition Using GO Enrichment Analysis
Based on the identification and quantification of DEGs, we further performed the functional annotations and bioinformatic characterization of DEGs in A. paragallinarum under iron restriction using GO enrichment analysis through the DAVID online analysis system. Firstly, we took the GO annotations analysis of these DEGs. As shown in Figure 3A and Supplementary Data Sheet 2, upregulated genes could be clustered into three categories, including six kinds of biological processes, five kinds of cellular components, and six kinds of molecular functions. Then, the GO enrichment analysis of these upregulated genes was carried out and the top 20 terms are summarized in Figure 3B and Supplementary Data Sheet 3. Twenty-two upregulated genes enriched in the 13 GO biological process terms, 22 upregulated genes enriched in the 5 GO cellular component terms, and 45 upregulated genes actively participated in the 27 GO molecular function terms were present in A. paragallinarum under the iron-restriction condition, respectively. Among them, three biological process terms that were involved in the response to macromolecule localization, protein localization, and establishment of protein localization were significantly enriched. Another three cellular component terms for protein binding were significantly enriched, which were mainly involved in the cell outer membrane, external encapsulating structure part, and membrane. Besides, biological process items for the transcription factor activity and sequencespecific DNA binding, nucleic acid binding transcription factor activity, and ATPase activity were significantly enriched.
In addition, the functional annotations and bioinformatic characterization of downregulated genes in A. paragallinarum under the iron-restriction condition were taken by GO enrichment analysis. As shown in Figure 3C and Supplementary Data Sheet 2, the downregulated genes could also be classified into three categories, which consisted of six kinds of biological processes, four kinds of cellular components, and five kinds of molecular functions. The GO enrichment analysis of these downregulated genes was also conducted and the top 20 terms are summarized in Figure 3D and Supplementary Data Sheet 4. The data showed that there were 391 downregulated genes enriched in the 75 GO biological process terms, 54 downregulated genes enriched in the 8 GO cellular component terms, and 185 downregulated genes actively participated in the 50 GO molecular function terms in A. paragallinarum under the iron-restriction condition, respectively. Additionally, 14 biological process items and 6 molecular function items were significantly enriched, which mainly participated in the oxidation-reduction process; electron transport chain; respiratory electron transport chain; generation of precursor metabolites and energy; "de novo" inosine-5'monophosphate (IMP) biosynthetic process; cellular respiration; IMP metabolic process; IMP biosynthetic process; sodium ion transport; cation transport; energy derivation by oxidation of organic compounds; nitrogen cycle metabolic process; single-organism metabolic process; metal ion homeostasis; oxidoreductase activity; acting on NAD(P)H, quinone, or similar compound as acceptor; heme binding; tetrapyrrole binding; acting on NAD(P)H; and electron carrier activity. However, no cellular component item was significantly enriched for these downregulated genes. Taken together, the identification of a larger number of DEGs between the iron-restriction group and the control group using DEGs GO enrichment analysis suggests that iron restriction can result in changes of a variety of biological processes and cellular component terms as well as molecular functions that are related to DEGs in A. paragallinarum.

Functional Annotations and Bioinformatic Characterization of DEGs in A. paragallinarum Under Iron-Restriction Condition Using KEGG Pathway Enrichment Analysis
In order to further detect the related pathways in which these DEGs in A. paragallinarum under the iron-restriction condition were involved in, the KEGG pathway analysis was performed. As shown in Figure 4A, upregulated genes could be classified into four KEGG categories, including four kinds of metabolism, one kind of genetic information processing, two kinds of environmental information processing, and one kind of cellular process. Then, KEGG enrichment analysis of these upregulated genes was carried out. As shown in Figure 4B and Supplementary Data Sheet 5, 17 upregulated genes were enriched in the 11 KEGG pathway items such as ABC transporters, pyrimidine metabolism, purine metabolism, quorum sensing, nicotinate and nicotinamide metabolism, and so on. However, no KEGG pathway item was significantly enriched.
By taking KEGG pathway analysis of downregulated genes in A. paragallinarum under iron-restriction condition, we found that the downregulated genes could also be classified into four KEGG categories, including six kinds of metabolism, one kind of genetic information processing, one kind of environmental information processing, and two kinds of human diseases ( Figure 4C). Furthermore, KEGG enrichment analysis showed that 75 downregulated genes were enriched in the 30 KEGG pathway items. Eight of them exhibited dramatic enrichment that were related to oxidative phosphorylation, citrate cycle (TCA cycle), nitrogen metabolism, two-component system, butanoate metabolism, methane metabolism, carbon fixation pathways in prokaryotes, and glyoxylate and dicarboxylate metabolism ( Figure 4D and Supplementary Data Sheet 6). Above all, the results of KEGG enrichment analysis demonstrate that iron restriction can result in changes of a series of biological metabolic pathways and signal transduction pathways that are associated with DEGs in A. paragallinarum.

Validation of DEGs in A. paragallinarum Under Iron-Restriction Condition Through qPCR Analysis
Based on the transcriptomic results, some crucial genes involved in the iron homeostasis showed significant difference between the iron-restriction group and the control group. In order to further validate the results, qPCR was carried out to detect the expressions of up-and downregulated genes in A. paragallinarum under iron-restriction condition. Here, three upregulated genes (VY92_RS03730, VY92_RS03735, and VY92_RS00335) and four downregulated genes (VY92_RS06600, VY92_RS09125, VY92_RS09790, and VY92_RS01655) were selected as representative DEGs for qPCR. Among them, VY92_RS03730, VY92_RS03735, and VY92_RS00335 were annotated as the heme utilization protein HutZ, heme utilization cystosolic carrier protein HutX, and ABC transporter ATPbinding protein, respectively. VY92_RS06600, VY92_RS09125, VY92_RS09790, and VY92_RS01655 were annotated as the H-type ferritin, cytochrome c nitrite reductase pentaheme subunit, succinate dehydrogenase/fumarate reductase ironsulfur subunit, and autotransporter outer membrane beta-barrel domain-containing protein, respectively. As shown in Figure 5A, the iron-restriction group showed significantly increased mRNA levels of all three upregulated genes than the control group. In addition, the expressions of four downregulated genes were lower in the iron-restriction group than in the control group, and a significant difference could be seen in VY92_RS06600 and VY92_RS09125 ( Figure 5B). Taken together, the results are in accordance with the above RNA-seq and suggest the distinct gene expressions of these iron-related DEGs in A. paragallinarum under iron restriction.

Identification and Bioinformatic Characterization of sRNAs in A. paragallinarum Under Iron-Restriction Condition
In previous results of the identification and quantification of genes in this study, 121 genes encoding sRNAs were verified. We further performed bioinformatic characterization of these sRNAs in A. paragallinarum under the iron-restriction condition.
Firstly, we summarized the length of these sRNAs. As shown in Figure 6 and Supplementary Data Sheet 7, the length of each sRNA was smaller than 500 nt, and most of them had the length of 101 to 150 nt. Subsequently, we determined if these sRNAs were likely to form stable secondary structures. Here, we selected four upregulated sRNAs (sRNA0022, sRNA0080, sRNA0084, and sRNA0103) and three downregulated sRNAs (sRNA0004, sRNA0036, and sRNA0099) in A. paragallinarum under the iron-restriction condition. As shown in Figure 7, these sRNAs could form desirable stem-loop secondary structures.
It is widely known that sRNA executes its regulatory function by pairing and binding to its target genes. We also identified potential interactions between these sRNAs and target mRNAs. As shown in Supplementary Data Sheet 8, a single sRNA could target various mRNAs, and the same target mRNA may also be regulated by different sRNAs. Among the four upregulated sRNAs, we found that they could target the genes mainly encoding the proteins such as DUF3277 domain-containing protein, Fic family protein, DUF1043 domain-containing protein, DUF1828 domain-containing protein, ABC transporter permease, ABC transporter ATPbinding protein, spermidine/putrescine ABC transporter permease, heme-binding protein, DNA helicase RecQ, DNAbinding protein, DNA repair protein RecO, cAMP-activated global transcriptional regulator CRP, transcriptional regulator MraZ, succinate dehydrogenase/fumarate reductase ironsulfur subunit, autotransporter domain-containing protein, MurR/RpiR family transcriptional regulator, met regulon transcriptional regulator MetJ, and so on. In terms of three downregulated sRNAs, sRNA0004 and sRNA0036 tended to target the genes which coded for the proteins such as DUF1367 domain-containing protein, Fic family protein, two-component system sensor histidine kinase QseC, and so on. Interestingly, we could see that genes encoding the same protein such as DUF1367 domain-containing protein could be targeted by both upregulated sRNAs and downregulated sRNAs. Unfortunately, no target gene was detected for sRNA0099 in this study. Considering that the gene encoding for the Fic protein showed no significance of transcriptomic changes of expression between the iron-restriction group and the control group as well as being targeted by the both up-and downregulated sRNAs, we also used qPCR to further detect the mRNA levels of the related gene. As shown in Supplementary Figure 4, the results were also consistent with the above RNA-seq and showed no difference of mRNA levels of VY92_RS05725 encoding for the Fic protein between the iron-restriction group and the control group, which suggested the stable expressions of Fic protein in A. paragallinarum under iron restriction.
Furthermore, the potential interaction between sRNAs and qPCR-validated DEGs was examined. Among those eight DEGs, VY92_RS09790 and VY92_RS00335 were shown to be targeted by various sRNAs, especially VY92_RS09790 (Table 4). VY92_RS09790, one of the crucial downregulated genes in A. paragallinarum under the iron-restriction condition, could be targeted by 13 sRNAs, including upregulated sRNA0022 and sRNA0103. Notably, upregulated VY92_RS00335 could also be targeted by sRNA0022, indicating the potential regulatory role FIGURE 5 | Validation of DEGs in A. paragallinarum under iron-restriction condition compared with the control condition by qPCR. (A,B) The mRNA levels of three upregulated genes and four downregulated genes in A. paragallinarum under iron-restriction condition compared with the control condition were determined by real-time PCR (N = 3), respectively. *p < 0.05, **p < 0.01, ***p < 0.001. of sRNA0022 in A. paragallinarum under the iron-restriction condition. Taken together, these results validate that sRNAs have a potential role in the regulation of iron homeostasis in A. paragallinarum.

DISCUSSION
Recently, transcriptomic and bioinformatic analyses, as a new frontier developed in life science, have been extensively used to study dynamic alterations in global transcriptomic profiles and explore the transcriptome response of bacteria to stimuli including nutritional stress (Ying et al., 2015). Therefore, we took RNA-seq to better understand gene expression profiles in A. paragallinarum under the iron-restriction condition. Here, we focused on the large quantities of bacterial genes regulated by iron restriction as well as possible functional activities. These analyses will be helpful for the in-depth realization of the molecular mechanisms involved in response to nutritional iron availability in A. paragallinarum to control cellular iron homeostasis when faced with iron starvation.
Bacteria are generally adaptable to changing circumstances such as various nutrient deficiencies or stresses. This adaption results in many physiological changes such as expression changes of global genes and induction of specific genes (Siegele and Kolter, 1992). Here, we also evaluated and compared the DEGs in A. paragallinarum under conditions with iron restriction and iron sufficiency as well as the related cellular functions and signaling pathways. Between these two conditions, we have revealed a number of distinct characteristics of global gene expression in bacteria. In comparison with iron sufficiency, iron restriction could lead to both upregulation and downregulation of genes, and it seems that the majority of genes showed decreased expression. Additionally, some genes were found to be uniquely present in the iron-restriction group. As previously reported, perhaps reduction of global gene expression and induction of unique genes are also helpful for A. paragallinarum to enter stationary phase by ceasing growth and retaining survivability for a long time (Miksch and Dobrowolski, 1995). Notably, several upregulated genes are demonstrated to be able to facilitate A. paragallinarum to effectively obtain iron for their growth (Abascal et al., 2009). In our study, when bacteria are in an ironrestriction environment, these upregulated gene products that we found can also be helpful for the release of iron from iron-binding proteins and compounds, promoting the acquisition and use of more iron sources by bacteria. The outer membrane receptors are generally induced by iron restriction and thought to enhance the rate of ferric-siderophore uptake, allowing bacteria to more effectively scavenge ferric-siderophores from their surroundings. Additionally, transport of ferric-siderophores through outer membrane receptors requires energy, and this energy is provided by the TonB energy system (Postle, 1993;Larsen et al., 1994;Higgs et al., 1998). Therefore, the increased expression of some proteins that are related to the TonB energy system in bacteria under the iron-restricted condition can facilitate the entry of ferric-siderophores across the outer membrane into the periplasm. TonB ExbB-ExbD proteins are found in many Gram-negative bacteria, mostly known to be found in the E. coli system (Braun et al., 1998;Palyada et al., 2004). In our study, the TonB system transport proteins ExbD and ExbB were found to be either iron stimulated and cotranscribed, and it is likely that the formation of the protein complex TonB ExbB-ExbD also existed. Due to that iron periplasm-to-cytosol transport requires ABC transporters on the cytosolic membrane, the transporter ATP-binding proteins are also upregulated to better contribute to the entry of iron into the cytoplasm as well as iron utilization for bacteria. Certainly, there are many complex regulatory mechanisms that were related to the use of iron, and we believe that our study will provide a good start for the iron regulation system in A. paragallinarum.
Furthermore, numerous cellular functions and signaling pathways were triggered or suppressed in A. paragallinarum under iron restriction in order to ensure iron homeostasis. Interestingly, several novel sRNAs were also elucidated in our study, while their expressions appeared to be different in A. paragallinarum under the iron-restriction condition compared with the iron-sufficiency condition. It suggests that these differences may reflect a survival strategy of A. paragallinarum in ensuring a rapid and accurate response to environmental iron changes. As far as we know, this is the first study of reporting comprehensive transcriptomic information in A. paragallinarum following iron restriction, with a bonus scene of the bioinformatic identification of new important sRNAs. In recent years, the iron-regulated proteins, especially the ferric siderophore outer membrane receptors, have been investigated as potential vaccine candidates (Alteri et al., 2009;Brumbaugh et al., 2013;Wang et al., 2019). Meaningfully, our study will also provide novel insight into the development of new antigens for potential vaccines against infectious coryza. The upregulated genes in A. paragallinarum under iron-restriction condition could be used for the preparation of gene-deleted bacterial strain and, thus, could be candidates for live attenuated vaccine and the genes themselves can also be regarded as candidates for subunit vaccines to control A. paragallinarum in poultry farms, just like other bacteria that have been previously reported (Zeng et al., 2009).

Hut Protein Is Required for Iron Acquisition and Heme Utilization in A. paragallinarum
It is well known that heme utilization systems have essential roles in bacterial iron acquisition and pathogenicity (Shi et al., 2019). A heme uptake system generally consists of the TonB-dependent outer membrane receptor protein, periplasmic binding protein, and inner membrane-associated ABC transporter . Despite extensive studies on the mechanism of heme transport from outside to the cytoplasm, proteins that participated in heme utilization remain unknown. Nowadays, Hut (heme utilization) protein has been identified in Vibrio cholerae according to genome sequence and bioinformatics-based predictions (Uchida et al., 2019). Uchida et al. have recently purified and characterized HutZ as a heme-degrading enzyme in V. cholerae, which releases iron into the cytoplasm (Uchida et al., 2012). As a unique heme storage protein, HutZ is critical for heme utilization and can cleave heme to biliverdin via verdoheme (Uchida et al., 2019). Besides, it is reported that HutZ is also required for biofilm formation and contributes to the pathogenicity of Edwardsiella piscicida (Shi et al., 2019). HutX, a cytoplasmic heme transport protein for HutZ, facilitates the transfer of heme to HutZ via a specific protein-protein interaction (Sekine et al., 2016). In A. paragallinarum, the function of the Hut protein family has not yet been elucidated. In this study, HutZ and HutX proteins were also identified in A. paragallinarum. RNA-seq analysis and qPCR successfully confirmed the preferential expressions of these proteins in A. paragallinarum under iron-restriction condition, indicating that heme uptake may be an important route for iron acquisition in A. paragallinarum. Herein, our results provide a clue that the heme transport system can play an essential role in iron acquisition as well as HutZ and HutX proteins are considered to be necessary for heme utilization in A. paragallinarum. Just like other bacteria mentioned above, this HutZ in A. paragallinarum can also be considered as a hemedegrading enzyme that releases iron into the cytoplasm and HutX is regarded as a cytoplasmic heme transport protein for HutZ, which facilitates the iron acquisition and heme utilization in A. paragallinarum. These heme acquisition-associated genes might be potential candidate genes for the development of subunit vaccines to control A. paragallinarum. Nevertheless, the regulatory mechanism of heme utilization by HutZ and HutX is worthy of further investigation.

Dominant Role of DUF Domain-Containing Proteins in A. paragallinarum Under Iron Restriction
Even when sequencing data is being increased, there are also a lot of genes or proteins that are uncharacterized and their domain functions are unknown (Nabi et al., 2020). These proteins are regarded as proteins having DUFs (Finn et al., 2016;El-Gebali et al., 2019). DUFs exist in a variety of species, such as around 1,500 in eukaryotes and 2,704 in bacteria (Goodacre et al., 2013). Evolutionary conservation suggests that many of the DUFs have vital functions. For example, DUF59, one of the most common DUFs, occurs in both bacteria and eukaryotes. Though the precise biochemical function of DUF59 remains unclear, previous studies have provided insights into its physiological roles. For instance, proteins containing DUF59 participate in many processes, which include maturation of iron-sulfur (FeS) proteins and intracellular iron homeostasis in S. aureus (Mashruwala et al., 2016). Recently, DUF143 has been considered to be essential in E. coli when cells are starved and its deletion in bacteria shows no obvious phenotype (Baba et al., 2006;Hauser et al., 2012).
Until now, few studies characterize the DUF domain-containing proteins in A. paragallinarum. In this study, most of the genes encoding for numerous DUF domain-containing proteins were identified in A. paragallinarum. Among these proteins, majority of them were upregulated and only a few were downregulated in bacteria under the condition of iron restriction. Remarkably, DUF454, DUF4198, and DUF2318 domain-containing proteins were significantly promoted, while DUF1345 domain-containing protein was significantly suppressed in the iron-restriction group than in the control group. According to the sRNA analysis, we speculate that the expression pattern of several DUF domaincontaining proteins could be regulated by sRNAs. This is the first research to characterize DUF domain-containing proteins in A. paragallinarum, implying that various DUF domaincontaining proteins may have positive roles in transcriptional regulation of physiological processes to cope with environmental iron restriction. Further experimentation will be needed to scrutinize the molecular and functional aspects of these crucial DUF domain-containing proteins in detail.

Stable Expression of Fic Protein Might
Be Attributed to the Regulatory of sRNAs in A. paragallinarum Under Iron Restriction Fic (filamentation induced by cAMP) protein is extensively present in bacteria, with a single member found in animals (Lu et al., 2016). It plays essential roles in drug tolerance, bacterial pathogenicity, and stress response. Previous studies have shown that Fic domain containing protein is able to post-translationally modify proteins (Mukherjee et al., 2011;Feng et al., 2012;Castro-Roa et al., 2013;Harms et al., 2016Harms et al., , 2017. In addition, Fic protein in non-pathogenic bacteria has been implicated in the regulation of bacterial cell division and programmed death (Worby et al., 2009). For instance, a mutation of the fic gene in E. coli can result in filamentation (Utsumi et al., 1982(Utsumi et al., , 1983. Although Fic protein has attracted considerable interest for its functions in bacterial infections and stress responses, their functions and related molecular mechanisms in iron homeostasis in A. paragallinarum had been unknown until recently. Here, we found that this protein was present in A. paragallinarum. However, there was no significance of gene expression of this protein between the iron-restriction group and the control group. Thus, it seems that iron-restriction may not result in the abnormal expression of Fic protein in A. paragallinarum. Identification and bioinformatic characterization of sRNAs further confirmed that Fic protein could be regulated not only by upregulated sRNA0022 but also by downregulated sRNA0004 and sRNA0036 in bacteria under the iron-restriction condition compared with that under the control condition. Consequently, we speculate that Fic protein is still stably expressed in A. paragallinarum under iron restriction, which might be largely attributed to the existence of these regulatory sRNAs. Herein, our results demonstrate for the first time about the effects of iron restriction on the expression of Fic protein as well as the tight relationship between this protein and regulatory sRNAs in A. paragallinarum.

CONCLUSION
We investigate the more profound transcriptomic profiles for the first time and illustrate changes of transcriptomes in A. paragallinarum in response to starvation using RNAseq. Compared with iron sufficiency, iron restriction results in many transcriptomic changes including changes of global gene expression and induction of specific genes. The bioinformatic analyses demonstrate that the Hut protein and DUF domain-containing proteins are preferentially activated in A. paragallinarum following iron restriction and have a positive role in iron acquisition and heme utilization in response to iron starvation. In addition, several novel sRNAs are identified in A. paragallinarum and have potential regulatory roles in iron homeostasis, which is likely to facilitate the stable expression of Fic protein. Researches on these genes will enrich our knowledge about the mechanisms of growth and pathogenesis of A. paragallinarum, which can be regarded as potential targets for the development of new prophylactic and therapeutic measures against A. paragallinarum.

DATA AVAILABILITY STATEMENT
The data presented in the study are deposited in the (NCBI SRA) repository, accession number (SRR13291129).

AUTHOR CONTRIBUTIONS
CH, FX, FL, DL, ZH, and HS carried out the experiments, analyzed the data, and wrote the manuscript. CH, JL, and HS designed the study and supervised the project. CH, XZ, FX, GL, YH, JL, and HS assisted in the data analysis and discussion. CH and HS drew the figures. All authors read and approved the final manuscript.