High-Resolution 3D Genome Map of Brucella Chromosomes in Exponential and Stationary Phases

The spatial structure of chromatin plays important roles in normal cell functions and in the regulation of gene expression. Three-dimensional genome sequencing has been performed in many mammals and plants, but the availability of such data for bacteria, especially intracellular pathogens, is still limited. ABSTRACT The three-dimensional (3D) genome structure of an organism or cell is highly relevant to its biological activities, but the availability of 3D genome information for bacteria, especially intracellular pathogens, is still limited. Here, we used Hi-C (high-throughput chromosome conformation capture) technology to determine the 3D chromosome structures of exponential- and stationary-phase Brucella melitensis at a 1-kb resolution. We observed that the contact heat maps of the two B. melitensis chromosomes contain a prominent diagonal and a secondary diagonal. Then, 79 chromatin interaction domains (CIDs) were detected at an optical density at 600 nm (OD600) of 0.4 (exponential phase), with the longest CID being 106 kb and the shortest being 12 kb. Moreover, we obtained 49,363 significant cis-interaction loci and 59,953 significant trans-interaction loci. Meanwhile, 82 CIDs of B. melitensis at an OD600 of 1.5 (stationary phase) were detected, with the longest CID being 94 kb and the shortest being 16 kb. In addition, 25,965 significant cis-interaction loci and 35,938 significant trans-interaction loci were obtained in this phase. Furthermore, we found that as the B. melitensis cells grew from the logarithmic to the plateau phase, the frequency of short-range interactions increased, while that of long-range interactions decreased. Finally, combined analysis of 3D genome and whole-genome transcriptome (RNA-seq) data revealed that the strength of short-range interactions in Chr1 is specifically and strongly correlated with gene expression. Overall, our study provides a global view of the chromatin interactions in the B. melitensis chromosomes, which will serve as a resource for further study of the spatial regulation of gene expression in Brucella. IMPORTANCE The spatial structure of chromatin plays important roles in normal cell functions and in the regulation of gene expression. Three-dimensional genome sequencing has been performed in many mammals and plants, but the availability of such data for bacteria, especially intracellular pathogens, is still limited. Approximately 10% of sequenced bacterial genomes contain more than one replicon. However, how multiple replicons are organized within bacterial cells, how they interact, and whether these interactions help to maintain or segregate these multipartite genomes are unresolved issues. Brucella is a Gram-negative, facultative intracellular, and zoonotic bacterium. Except for Brucella suis biovar 3, Brucella species have two chromosomes. Here, we applied Hi-C technology to determine the 3D genome structures of exponential- and stationary-phase Brucella melitensis chromosomes at a 1-kb resolution. Combined analysis of the 3D genome and RNA-seq data indicated that the strength of short-range interactions in B. melitensis Chr1 is specifically and strongly correlated with gene expression. Our study provides a resource to achieve a deeper understanding of the spatial regulation of gene expression in Brucella.

differences in the 3D genome interactions of the two B. melitensis chromosomes between the different growth states were identified. Combined analysis of the 3D genome and transcriptome sequencing (RNA-seq) data revealed that the strength of short-range interactions in B. melitensis Chr1 is specifically and strongly correlated with gene expression.

RESULTS AND DISCUSSION
Generation of high-resolution whole-genome contact maps for B. melitensis. The 3D chromosome structure reveals the intrinsic function of genome organization in biological processes. Chromosome conformation capture combined with deep sequencing has been used previously to explore bacterial genome organization (3). To reveal the spatial structure and interaction of two chromosomes in Brucella, we applied Hi-C technology (2), which utilizes cross-linking with formaldehyde in conjunction with spatially constrained ligation to assess the average spatial proximity of genomic loci of early-exponential-phase (optical density at 600 nm [OD 600 ] of 0.4 [AOD0.4]) (Fig. 1A) and stationary-phase (AOD1.5) (Fig. 1B) B. melitensis cultures. Each Hi-C experiment yields a genome-wide library of 1-kb ligation products whose frequencies reflect the 3D structure of the genome. As expected, the two chromosomes in the resulting contact maps exhibited a strong diagonal signal, reflecting frequent contacts between adjacent loci. This main diagonal represents the high-frequency interactions between DNA loci on the same arm of the chromosome (i.e., intra-arm interactions). A second, less prominent diagonal reflects lower-frequency interactions between loci on one arm of the chromosome and those on the opposite arm (i.e., interarm interactions) (23). Thus, the contact heat maps of the two B. melitensis chromosomes contain a prominent diagonal and a fuzzy secondary diagonal ( Fig. 1A and B). These interactions resemble those of unichromosomal bacteria, such as C. crescentus, and of Chr2 of V. cholerae, a bacterium with a multipartite genome (21). According to the interaction matrix, the sum of interactions within a certain range upstream and downstream of each bin on the chromosome can be calculated, and the proportion of the interactions between the bin and the entire chromosome reflects the compactness of the local region on the genome. This method has been used to study the effect of a structural protein in E. coli on the degree of local spatial compactness of the genome (1). Figure 1C and D show the change in the degree of genome compaction of the two B. melitensis chromosomes at the early exponential and stationary phases. Then, contact decay curves of the genome-wide 2-kb resolution of AOD0.4 and AOD1.5 were integrated, and the differences in distances and interaction frequencies between the samples were compared and analyzed (Fig. 1E). The conversion of contact maps into 3D structures at OD 600 s of 0.4 (left) and 1.5 (right) ( Fig. 1F; see Movies S1 and S2 in the supplemental material) shows that Chr1 and Chr2 fold into a helicoidal shape with its two replichores tightly interlaced. The positions of origins (red balls) and termini (purple balls) are annotated in the figures. These contact maps provide insights into the interplay of all loci in the two B. melitensis chromosomes.
Analysis of chromatin interaction domains (CIDs) in the 3D genome of B. melitensis. Chromatin interaction domains (CIDs) refer to the high-level structures in genome space with relatively strong local interactions found on prokaryotic interaction maps (14). These are similar to the eukaryotic topologically associated domains (TADs) (42). CIDs are related to important structural and regulatory roles in bacteria. Figure 2A to D show the distribution of CIDs on chromosomes at a 2-kb resolution, including interaction heat maps, insulation scores, and CID boundaries. In B. melitensis at an OD 600 of 0.4, 79 CID regions were detected, with the longest being 106 kb and the shortest being 12 kb. In B. melitensis at an OD 600 of 1.5, 82 CID regions were detected, with the longest being 94 kb and the shortest being 16 kb (Table S1). This may suggest that as bacteria transition to stationary phase, the number of CIDs increases and the range becomes smaller. In order to compare the stability of the CID boundary structure between the two growth stages, we compared the CID boundaries between the samples according to the 2-kb resolution CID analysis in each sample and identified the CID boundaries that were shared between the samples and unique to each sample. Figure S1 shows that AOD0.4 (logarithmic-phase sample) shared 41 CIDs with AOD1.5 (plateau sample), while AOD0.4 and AOD1.5 had 36 and 39 unique CIDs, respectively. This result shows that the variation in CIDs during bacterial growth is relatively high.
CIDs are usually divided into border and interior regions. According to the results of CID analysis, the gene density distribution and GC content distribution of CID in B. melitensis were calculated. The results showed that there was no significant difference  Table S3) at an OD 600 of 0.4 or 1.5. In eukaryotes, TAD boundaries are enriched in highly expressed genes, and structural proteins such as CTCF (CCCTC-binding factor) and cohesin are usually enriched at TAD boundaries in the human genome (43). The characterization of CID or TAD boundaries is of great significance in understanding the mechanism by which the genomic topology of this species is established (44). We used FIMO software (45) to scan the motifs at the CID boundary and found that certain motifs were enriched in the boundary region. We also visualized these motif-binding sequences. The top four CID boundary enrichment motif sequences are shown in Fig. 3A to D; detailed information is shown in Table S4. These motifs are associated with transcription factors, indicating that CID boundaries are dependent on transcription.
Subsequently, the gene sets located within the CID boundary in B. melitensis were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis (Fig. S3). GO analysis showed that the biological functions of the CID boundary genes were mainly enriched in metabolism and cellular processes, and their subcellular localization was mainly concentrated in the membrane and cell parts. Molecular functions were mainly enriched in catalytic activity, binding, and transport function ( Fig. S3A and B). KEGG analysis showed that B. melitensis CID boundary genes were mainly involved in metabolism-related pathways ( Fig. S3C and S3D). These results suggest that CID boundary genes may play important roles in B. melitensis growth and metabolism. Furthermore, to explore the functions of genes within the boundaries of the CIDs that were unique to each sample, the gene sets located within the boundaries of differential CIDs were functionally annotated (Fig. 4). Figure 4A and B show that the differential CIDs of AOD0.4 and AOD1.5 are mainly enriched in metabolic processes, and the subcellular localization of these genes is mainly enriched in the membrane and cell parts. The enriched molecular function terms are mainly related to catalytic activity, binding, and transport function. Figure 4C and D show that differential CID boundary genes are mainly involved in metabolism-related pathways.
Analysis of interacting loci in B. melitensis at the early exponential and stationary phases. DNA loci that are distant within the genome may become spatially close via the formation of loops, bringing different biological domains closer, which impacts transcriptional regulation (16,46). Fit-Hi-C software (47) was employed to statistically analyze the significance of interactions between every two loci in the whole genome at a 1-kb resolution to detect the chromatin loop structure. We obtained 49,363 significant cisinteraction loci and 59,953 significant trans-interaction loci at the logarithmic growth phase (Table S5). We also obtained 25,965 significant cis-interaction loci and 35,938 significant trans-interaction loci at the stable growth phase (Table S5). Figure 5A and B show the significant cis-interaction loci of B. melitensis at OD 600 s of 0.4 and 1.5. Figure 5C and D show the significant trans-interaction loci of B. melitensis at OD 600 s of 0.4 and 1.5.
Hi-C interaction studies have shown that most interactions usually occur within chromosomes, and interchromosome interactions were long considered to represent in-matrix noise. However, in olfactory neurons, the genome forms a large number of specific transinteractions in 3D structure, regulating the expression of different olfactory receptor genes (48), suggesting that trans-interactions may have undiscovered biological significance. In addition, recent studies have reported a number of trans-interactions in filamentous fungi (7). In our study, the interaction between replication origins and termini was particularly strong in the trans-interaction Circos plot (Fig. 5C and D), indicating that these trans-interactions might play important roles in the replication process of B. melitensis.
Significant cis-interaction loci can also be classified according to different distance ranges. B. melitensis has two chromosomes. Considering that the difference in chromosome length will affect the statistical results, we drew histograms based on distance calculations in the circular genome. Figure S4 shows that when the interaction distance on the genome is less than 200 kb, the percentage of cis-interactions is remarkably high regardless of the bacterial growth state. The percentage of cis-interactions decreased dramatically for interaction distances between 200 and 600 kb. However, there was a clear increase in the frequency of cis-interactions as the distance continued to increase (600 to 800 kb). When the interaction distance was longer than 800 kb, the frequency of cis-interactions decreased. In addition, a comparison of Fig. S4A and B shows that as the bacteria grew from logarithmic to plateau phase, the proportion of short-range interactions (,200 kb) increased, while that of long-range interactions (.200 kb) decreased.
To visualize the differences in interactions between the early exponential and stationary phases of Brucella cultures more intuitively, an interaction subtraction matrix was utilized (49). Two heat maps representing the single chromosome 1-kb resolution interaction subtraction matrixes were obtained ( Fig. 6A and B). Figure 6C and D show heat maps of the significant difference interaction subtraction matrix at 5-kb resolution within the red region in Fig. 6A and B. The main differential interactions involve B. melitensis replication origins and termini. At an OD 600 of 0.4, the interactions between the origin of replication and other genes are strong, and these interactions are significantly reduced after the bacteria enter the stable stage. This indicates that the genes at these positions may be involved in the growth and replication of B. melitensis. Transcription is involved in local chromatin structure in Chr1. A previous study showed that the boundary genes of CIDs are often highly expressed (1, 14). In our study, we first constructed a gene expression pattern cluster diagram for the differential genes screened by RNA-seq (Fig. 7A) of exponential-and stationary-phase B. melitensis. Then, the gene expression characteristics of the CID border and interior were detected. Figure 7B and C show that gene expression of CID border and interior regions are significantly different at both OD 600 s of 0.4 and 1.5, with CID border regions having higher gene expression levels than interior regions.
Moreover, it was reported that transcription level, transcription rate, and transcript length drive the formation of chromosomal interaction domain boundaries (14,15). In our study, to verify the correlation of transcription and the 3D genome structure in B. melitensis, scatterplots between 3C signal and the gene transcription level of AOD0.4 and AOD1.5 were created as shown in Fig. S5A and B, respectively. The results showed a positive correlation. Then, we analyzed the curves of the short-distance interaction strength and gene expression of the two chromosomes under different culture conditions. Figure 8A and B show the strong correlation between the transcription levels Remarkably, these correlations were also observed in contact maps of E. coli and B. subtilis generated by different laboratories using different enzymes and cross-linking conditions (1). However, Chr2 exhibits a weak correlation (PC of ,0.21) between transcription level and short-range interaction strength ( Fig. 8C and D), similar to the pattern observed in M. pneumoniae and C. crescentus (1). These results imply the existence of transcription-induced constraints in B. melitensis Chr1 that favor interactions between neighboring loci.
The underlying mechanisms between transcription and short-range contacts remain to be deciphered. Nucleoid-associated proteins (NAPs) are a class of nucleic acid-binding basic proteins present in prokaryotic cells that can bind to bacterial chromosomes to form compact structures-that is, nucleoids. Local binding of NAP H-NS in E. coli inhibits short-range interactions, whereas HU and Fis promote long-range contacts in different ways (1). A recent study found that the transcription factor Rok forms large nucleoprotein complexes in B. subtilis, and the complexes interact robustly with each other over long distances. These Rok-dependent long-range interactions lead to the formation of anchored chromosome loops that spatially isolate large sections of DNA (50). However, whether these NAPs play important roles in the changes in the 3D genome during the growth of B. melitensis and whether this mechanism is applicable to B. melitensis remain to be studied further.
Conclusions. This study provides whole-genome contact heat maps of the exponential-and stationary-phase B. melitensis chromosomes at a 1-kb resolution. In total, 79 CID regions, 49,363 significant cis-interaction loci, and 59,953 significant trans-interaction loci were detected in B. melitensis at the logarithmic growth phase. Meanwhile, a total of 82 CID regions, 25,965 significant cis-interaction loci, and 35,938 significant trans-interaction loci were obtained for B. melitensis in the stable growth phase. Furthermore, combined analysis of 3D genome and RNA-seq data revealed that the strength of short-range interactions in B. melitensis Chr1 is specifically and strongly correlated with gene expression. Taken together, the results of our study provide a model to explain certain differences of Brucella in different growth periods from a 3D genome perspective and will serve as a resource for deeper study of the spatial regulation of gene expression in Brucella, which will contribute to further understanding of Brucella's growth and pathogenesis.

MATERIALS AND METHODS
Bacterial strains, media, and growth conditions. Bacterial culture was performed as reported previously (51,52). Briefly, B. melitensis 16M was grown at 37°C in either tryptic soy broth (TSB) (BD; no. 211825) or tryptic soy agar (TSA) (BD; no. 236950). The experiments on B. melitensis were conducted in a biosafety level 3 laboratory.
Chromosome conformation capture with deep sequencing. Hi-C (chromosome conformation capture with deep sequencing) experiments were performed as described previously, with slight modifications (7,14). B. melitensis cells were grown in TSB, and bacterial growth was measured via optical density at 600 nm (OD 600 ). Bacterial cultures were collected when the OD 600 value reached 0.4 (exponential phase) or 1.5 (stationary phase). To fix the long-range DNA interactions, bacterial DNA was cross-linked with 3% formaldehyde (methanol free; Pierce) for 30 min at room temperature and another 30 min at 4°C. The reaction was stopped with 375 mM glycine for 15 min at 4°C, and the cells were washed with SET buffer (20 mM Tris [pH 8.0], 25 mM EDTANa 2 Á2H 2 O [pH 8.0], 75 mM NaCl) by centrifugation at 2,000 Â g for 5 min at 4°C. The supernatant was discarded, and the cells were resuspended thoroughly in 3 mL SET buffer. Aliquots of bacteria were frozen at 280°C.
Hi-C library construction and sequencing were carried out as described previously (7). Briefly, the genomic DNA from the bacteria was cross-linked for 1 h with 3% formaldehyde at 4°C and quenched with glycine at a final concentration of 0.375 M for 15 min. The cross-linked cells were lysed by incubation with lysozyme at 37°C for 10 min. Endogenous nucleases were inactivated with 0.5% SDS, and chromatin DNA was digested with 50 U Sau3AI (New England BioLabs, Inc.). Restriction fragment ends were labeled with biotinylated cytosine nucleotides by biotin-14-dCTP (TriLINK), followed by ligation using 50 U T4 DNA ligase (New England BioLabs, Inc., Ipswich, MA). After reversing the cross-links, the ligated DNA was extracted through using a QIAamp DNA minikit (Qiagen) according to the manufacturers' instructions. Purified DNA was sheared to a length of 300 to 500 bp and was further blunt-end repaired and A-tailed, with adaptor added, using the NEBNext Ultra II DNA library prep kit for Illumina (NEB, catalog no. E7645S). Point ligation junctions were pulled down with Dynabeads MyOne conjugated with streptavidin C1 (Thermo Fisher), followed by PCR amplification. Then, the PCR product was subjected to paired-end sequencing by BGI on the MGISEQ-2000 platform with PE150. Two biological replicates were used for each sample.
RNA isolation. RNA isolation was performed as described previously (51,52). Briefly, B. melitensis 16M was grown at 37°C and 180 rpm to an OD 600 of 0.6 to ;0.8. The growth of the bacterial culture (5 mL per sample) was stopped by adding 0.625 mL of prechilled stop solution (5% phenol-95% ethanol). Cells were harvested by centrifugation for 10 min at 6,000 Â g, and the resulting pellet was stored at 280°C. Frozen cell pellets were thawed on ice and resuspended in 100 mL per pellet of TE buffer (10 mM Tris-HCl, 1 mM EDTA [pH 8.0]). Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). The RNA samples were then subjected to treatment with DNase I (TaKaRa, Dalian, China) to remove genomic DNA. The concentration and quality of each sample were evaluated using a NanoDrop instrument (Thermo Fisher Scientific, Waltham, MA, USA) and agarose gel electrophoresis.
RNA-seq. RNA library construction and sequencing were performed as described (1). Briefly, total RNA was stabilized in vivo with RNAprotect bacterial reagent (Qiagen) and purified with an RNeasy Protect bacterial minikit (Qiagen). rRNA was depleted with a RiboZero-rRNA removal kit for bacteria (Epicenter). Libraries were prepared according to the instructions accompanying the Illumina TruSeq SBS Transcription kit.
Hi-C data analysis. The sequenced raw reads were filtered and quality evaluated using Trimmomatic (53) and FastQC software (54) with the default parameters to generate clean reads. The clean reads obtained were mapped onto the genome of the B. melitensis 16M strain, and valid paired reads were generated by the hiclib pipeline. The reliability of the Hi-C data from the two replicates was assessed with Genome DISCO. Chromatin interaction domains (CIDs) were identified by calculating the insulation score values using the software cworld-dekker. The intra-and interchromosomal interactions at the resolution of 1-kb bins were selected by Fit-Hi-C software (47), with thresholds of a P value of #0.01, false-discovery rate (q) value of #0.01, and contact count of .2. Motif scanning was carried out in the JASPAR database (http://jaspar.genereg.net) using MEME Suite (http://meme-suite.org).
Establishment of 3D genomic model of B. melitensis. The 3D model of the B. melitensis genome was predicted from Hi-C data by the PASTIS software using a 1-kb resolution matrix with the multidimensional scaling (MDS) module (55). In the 3D model of B. melitensis based on the obtained Hi-C data, the origin and terminus regions were annotated with Centurion software (56).
Data availability. The Hi-C and RNA-seq data of the B. melitensis 16M strain generated in this study have been submitted to the NCBI BioProject database (https://www.ncbi.nlm.nih.gov/bioproject/) under accession no. PRJNA918778 and PRJNA916855, respectively.