The order of macrophage 4D genome coordinates gene transcription during differentiation and infection


 The highly organized three-dimensional genome is crucial for gene transcription. However, it remains elusive how the order of the genome architecture related to its function. Here, we developed a single-cell Hi-C method and proposed TAD “degree of disorder” as a measure of genome organizational patterns, which is correlated with the chromatin epigenetic states, gene expression and co-regulation, and chromatin structure variability in individual cells. Upon Mycobacterium tuberculosis infection, NF-κB enters into the nucleus, binds to the target genome regions and initiates systematic chromatin conformation reorganization. Furthermore, we identified a remote NF-κB enriched enhancer promotes the expression of PD-L1 through chromatin loop, which could be a potential anti-tuberculosis and even anti-tumor therapeutic target. The integrated Hi-C, eQTL, and GWAS analysis depicted the atlas of the long-range target genes of tuberculosis susceptible loci. Among which SNP rs1873613 is located in the anchor of a dynamic chromatin loop with LRRK2, whose inhibitor AdoCbl could be an anti-tuberculosis drug candidate. Our study provides comprehensive resources for the 4D genome of immunocytes and sheds insights into the genome organization order and the coordinated gene transcription.


INTRODUCTION
It has been many decades since Schrödinger posed the intriguing assumptions about the order 55 or entropy in organisms and the chromosomes containing code-script of the individual's 56 development and function in 1944 1 . While the hereditary code-script in the chromosomes has been 57 well deciphered by sequencing now, there is still a big gap to understand the order or entropy of 58 the chromosomes in the cell nucleus. It has been well-established that the folding and organization 59 of the chromosomes in the nucleus is a process with a high degree of non-randomness and order 2,  94 differentiation and infection 95 To investigate the 4D genome and epigenetic dynamics during differentiation into different 96 lineages of macrophages, we used THP-1 cells, a highly plastic human monocytic cell line 97 that can be differentiated into macrophage-like cells and activated by M.tb, as the cell model 18,19 . 98 In our study, we named the THP-1 cells before differentiation, after differentiation, and infection 99 by M.tb as Thp1-mono, Thp1-macro, and Thp1-M.tb respectively. After PMA (phorbol 12-100 myristate-13-acetate) treatment, the nuclei of differentiated cells were transformed from ellipsoid 101 to kidney-shaped (Fig. 1a), implying a possible 3D genome conformational change. The distinct  During these processes, we observed heterogeneous phenotypes between individual cells, 110 such as some cells have typical macrophages morphology including amoeboid shape, largest 111 surface area and cell volume, while other cells remain as spherical shape during differentiation. 112 Moreover, highly diverse invaded intracellular bacteria numbers and invasion times were observed 113 between different macrophages, suggesting distinct degrees of host immune system activation. To 114 investigate the heterogeneous chromatin architecture of these cells at single cell level, we 115 developed a single-cell-indexed DLO Hi-C (sciDLO Hi-C) based on DLO Hi-C 20 and two rounds 116 of molecular barcoding to capture the chromatin conformation of individual cells (Fig. 1b). The  Fig. 1g, h). Next, the single cell data were merged to generate a 127 pooled Hi-C matrix, which is similar to the bulk Hi-C data (Extended Data Fig. 1i). To investigate 128 the heterogeneity of these three type of cells, we employed scHiCTools 23 to classify all the 129 individual cells based on their 3D genome structures. As shown in Fig. 1c, three distinct clusters 130 of cells were identified, representing Thp1-mono, Thp1-macro, and Thp1-M.tb, respectively.

131
Among which, the Thp1-mono has obvious boundaries between the other two types of cells, 132 whereas Thp1-macro and Thp1-M.tb are partially overlapped. This data reflected that PMA 133 treatment could uniformly reprogram monocyte to a distinct cell type macrophage, whereas 134 macrophage activation caused by M.tb infection was much more heterogeneous. Notably, 135 compared to Thp1-macro and Thp1-mono cells, Thp1-M.tb were enriched with significantly more 136 chromatin contacts around the immune genes (Fig. 1d). For example, the simulated chromatin 137 structures of individual cells showed that the chromatin interaction between innate immune-related 138 genes NOD2 and BRD7 were more contacted in Thp1-M.tb cells (Fig. 1e,f). 139 141 gene transcription 142 To explore the order and stochasticity of genome organization, we further analyzed chromatin 143 contact patterns in individual cells. We observed that while the individual chromatin contacts 144 displayed a certain extent of heterogeneity, the chromatin structures in specific regions of a cell 145 in certain process were uniformly folded. Take the innate immune-related STAT1 gene locus as an 146 example ( Fig. 2a-d). Most of the individual Thp1-macro cells had stochastic chromatin interaction 147 7 pattern between STAT1 promoter and an potential enhancer (Fig. 2a,c and Extended Data Fig.   148 2a), while the promoter-enhancer contact patterns of Thp1-M.tb cells in this region were much 149 more consistent (Fig. 2b,d). These data suggested that in the highly organized genome regions, 150 there was an intrinsic non-randomness or certainty for chromatin folding at both single cell and 151 bulk cell level.  Table 4). We retained all the significant interaction 158 spots in the Hi-C matrix to filter the stochastic genome loci with random chromatin interactions.

159
As the distance between the spots might reflect the similarity of the chromatin folding patterns, 160 the overall mean distance between the spots was then calculated to represent the whole DoD value  Next, the TAD around MYC gene was used an example to illustrate the detailed DoD 171 dynamics. As this genome locus with abundant enhancer and suppressor elements around the MYC 172 gene is crucial for cell pluripotency and stemness 24 . In Thp1-mono, the chromatin contacts in this 173 region were randomly distributed and had no obvious interaction patterns (DoD=2.74). Notably, 174 8 the Hi-C matrix of this TAD displayed more highly organized chromatin contact patterns 175 (DoD=2.42) after differentiation (Fig. 2f), suggesting a systematic non-random promoter-enhancer 176 or promoter-repressor 24 interaction occurred in the individual cells during differentiation. Notably, 177 we observed a stripe structure right above MYC gene after differentiation (Fig. 2f, marked with 178 dashed box), which may be attributed to the highly organized chromatin contacts of MYC gene to 179 a series of cis-elements or genes in this TAD.  (Fig. 2h), suggesting the stochasticity of TAD 187 organization is highly related to gene transcription.

189
As low DoD reflects that the individual cells have similar folding patterns and the chromatin 190 organization is less stochastic, the genes in these TADs may undergo similar transcriptional 191 regulation. In this scenario, the genes in highly organized (low DoD) TADs may be more likely to 192 be synchronously coregulated due to the cooperation within the same transcriptional regulatory 193 complex. To test this hypothesis, we defined the average gene coregulation score (CRS) within a 194 TAD during differentiation and infection, as shown in Extended Data Fig. 2g. Notably, we found 195 that TADs with lower DoDs had overall higher CRSs (Fig. 2i), supporting that highly organized 196 chromatin tends to be more synchronously (or less randomly) coregulated.    low DoD values. In contrast, TADs with high DoD contain significantly more transcription 228 10 repressive signals and heterochromatin signals, such as H3K27me3 and H3K9me3 (Fig. 3b). Our 229 data implied that the genome architecture of TADs with low DoDs were more sophisticatedly 230 organized, which might be mediated by the elegant cooperation with a greater number of 231 transcriptional regulatory factors and more accessible chromatin within TADs. Consistent with this 232 hypothesis, we observed that the DoD was indeed highly negatively correlated with ATAC signals 233 within the TADs (Fig. 3c), suggesting that the DoD is correlated to the chromatin accessibility and 234 the subsequent binding of transcriptional regulatory proteins.  (Fig. 3d,e). Moreover, the ATAC peak and active chromatin marks such as H3K4me3 240 and H3K27ac were enriched in the promoter region. Consistently, the expression of HERC5 and 241 HERC6 were also synchronously co-upregulated (Fig. 3f,g). By simulating the spatial structure of 242 this region, we found that in the low DOD state, the spatial distance between HERC5 and HERC6 243 was more adjacent (Fig. 3h,  values, more active enhancer/promoter enrichment, higher gene transcription levels and more 252 concerted coregulation. In contrast, the disordered TADs with higher "TAD entropy" levels 253 generally exhibit lower transcription and more random coregulation (Fig. 3j). Together, these 254 results suggested that the order of chromatin structure is highly correlated with its epigenetic     we found that GBP family genome region were enriched with these two phase separation marker 296 proteins (Fig. 4d). Our ATA C -seq and H3K4me3 ChIP-seq data indicate that the chromatin state 297 is activated and may bind with more regulatory proteins after infection (Fig. 4a). Furthermore, the   Table 6. Notably, these genes are significantly enriched in immunity pathways, which is not immunofluorescence assay revealed that NF-κB (p65) was indeed translocated into the nucleus 328 upon infection (Fig. 5d).

330
The chromatin accessibility and RNA expression analysis showed that, upon M.tb infection, 331 the NF-κB target loci turned to be more open and the expression level of target genes was 332 significantly up-regulated compared to the random genes (Fig. 5e). Moreover, the chromatin loops 333 related to the NF-κB target genes were more strengthened compared to the random genes ( Fig. 5f). 334 We further investigated the chromatin remodeling of the typical NF-κB target gene loci, such as   The single cell Hi-C data also demonstrated that the local chromatin interaction pattern became 345 more ordered during infection (Fig. 5h). The chromatin loop analysis revealed that these IFIT 346 genes were regulated by a single upstream cis element enriched with NF-κB binding motif ( Fig.   347 5i,j). Before infection, the cis-element was linked to IFIT1B and IFIT5 (Fig. 5i), which were 348 dynamically reshuffled to link with IFIT1, IFIT2, and IFIT3 after infection (Fig. 5j). Interestingly, 349 it has been demonstrated that IFIT1, IFIT2, and IFIT3 proteins can interact with each other and 350 form a complex to perform antipathogenic functions 31 . Consistently, we observed synchronously 351 upregulated expression of these three genes in the newly established loops (Fig. 5j). Notably, 352 IFIT1B and IFIT5, which are not located in this spatially adjacent hub, were not upregulated. These 353 data demonstrated a high order chromatin structure mediated elegant spatial and temporal co-354 regulation of the NF-κB target genes transcription during immunoresponse (Fig. 5k).  (Fig. 6a-c). 364 Notably, t h e ATA C -seq data showed that these enhancer and promoter regions were more 365 accessible and enrichment active chromatin modifications after M.tb infection (Fig. 6b,c), 366 indicates that more regulatory proteins are enriched in the enhancer region. Transcription factor 367 binding motif analysis revealed that both enhancer and promoter (Fig. 6b,c and Extended Data 368 Fig. 6a,b) regions can be bound by NF-κB. What's more, by ChIP-qPCR, we proved that NF-κB 369 (p65) was significantly enriched in this enhancer region after infection (Fig. 6d).  Fig. 7c,d). Importantly, the target genes of TB susceptibility loci 394 discovered by this integrated omics analysis, such as ASAP1 and LRRK2, have been reported to be 395 involved in TB pathogenesis 17, 34 , supporting the integrity of our multiomics analysis.

397
We further analyzed LRRK2, which has a strong loop and eQTL correlation with the TB 398 susceptibility SNP rs1873613. The DoD in the TAD was decreased after M.tb infection (Fig. 7b). 399 Moreover, the simulated chromatin structure revealed that SNP rs1873613 and LRRK2 tended to 400 be adjacent in space upon infection (Fig. 7c). The ENCODE ChIP-Seq data showed that this SNP 401 is located in the enhancer region of LRRK2 with the binding site of STAT3, SPI1, and FOS 402 (Extended Data Fig. 7e). Moreover, this SNP is located right in the anchor of an LRRK2 chromatin 403 loop that was significantly strengthened after infection (Fig. 7b). This data suggests that this TB 404 susceptibility SNP locus was dynamically reshuffled upon infection to spatially link with LRRK2 405 to regulate its gene transcription. Consistently, the eQTL analysis data showed that a T:C mutation 406 can indeed increase the expression of LRRK2 in the lungs (Fig. 7d). This SNP in the patient may showed that AdoCbl could significantly inhibit the proliferation of M.tb in the lungs (Fig. 7h,i). 416 This result was further supported by hematoxylin-eosin staining of tissue sections, showing that 417 AdoCbl significantly inhibited the initiation of lung and spleen lesions (Fig. 7j). These data            623 The nuclei were stained by DAPI (Thermo) and diluted by nuclei wash buffer. The density of

785
The reference genome was divided into fragments according to the restriction enzyme sites, and 786 the uniquely mapped sequence pairs were aligned to the restriction enzyme fragments. If two ends 787 from a paired-read were mapped to the same restriction enzyme fragment, the paired-read was 788 considered a self-ligation product. If two ends from a paired-read were mapped to two adjacent 789 restriction enzyme fragments, the paired-read was considered a religation product. Both the self-

853
The hiclib library was applied to distinguish A and B compartments as previously described 61 .

854
If the values in the first eigenvector were higher than 0, the corresponding bins were marked as A 855 compartments. If the values in the first eigenvector were smaller than 0, the corresponding bins 856 were marked as B compartments. The consistency of the overall transcriptional change direction (upregulation or 881 downregulation) of the genes within a TAD during differentiation and activation was defined as 882 the TAD CRS. The regulation direction score (DS) of a TAD t was defined as follows: where sgn is the sign function. The TAD CRS was defined as the absolute value of the DS :   Thp1-mono and Thp1-macro around MYC gene region by using bulk cell DLO Hi-C data sets. g, Relationship between TAD DoD and average gene expression level (normalized read count) in Thp1mono, Thp1-macro, and Thp1-M.tb, respectively. The DoD value higher than median was de ned as "High DoD"; lower than the median was de ned as "Low DoD". Signi cance of differences was measured by Kolmogorov-Smirnov test. h, Correlation between the TAD DoD and A, B compartments. In x axis, B/A means the TADs which contains both A and B compartments. Signi cance of differences was measured by unpaired one-sided t test. i, Correlation between the TAD DoD and coregulation score (CRS). TADs with lower TAD DoD preferentially have higher CRS, underlying more coordinated gene regulation in such TADs. x axis is the average value of the quantile groups with the lowest, low, high and highest TAD DoDs, respectively. j, Absolute values of the log2 fold changes in gene expression (y-axis) in the boundary not changed but DoD value signi cant changed TADs (dynamic DoD), boundary and DoD unchanged TADs (static DoD), boundary shifted TADs (shifted), boundary fused TADs (fused), and boundary divided TADs (divided) during differentiation and infection.

Figure 3
Low DoD is associated with active epigenetic modi cation and high gene expression level. a,b, Correlation between TAD DoD and chromatin epigenetic pro le. The y axis shows the normalized ChIPseq peaks per TAD. The DoD value higher than median was de ned as "High DoD"; lower than the median was de ned as "Low DoD". Signi cant differences were measured by unpaired one-sided t test. c, Correlation between TAD DoD and chromatin accessibility. Each dot represents one TAD. x-axis is the number of ATAC peaks per Mb. d,e, Chromatin contact matrix and DoD value changes around HERC5 and HERC6 gene region after M.tb infection. f,g, RNA expression level, chromatin accessibility, H3K27ac, and H3K4me4 histone modi cation changes after M.tb infection. h,i, Simulated three-dimensional structure of HERC5 and HERC6 gene region during M.tb infection by using bulk cell DLO Hi-C data sets. j, The relationship between DoD value, chromatin folding order, chromatin accessibility, gene transcription coregulation, and histone modi cation.

Figure 4
Remodeling of the chromatin con guration of guanylate-binding protein (GBP) gene clusters. a, Comparison of TAD DoD value, chromatin interaction matrices, chromatin loops, RNA expression levels and epigenetic modi cations around the GBP loci in bulk cells before and after M.tb infection. b, Comparing the single cell chromatin contact differences using merged single cell chromatin contact matrix. Each color of the dots in the matrix represents the chromatin contact from the same cell. c, Single cell chromatin contacts between GBP2, GBP4, GBP5, and GBP7 and RNA expression before and after M.tb infection. d, Depiction Hi-C chromatin loops before and after infection and BRD4 and MED1 ChIPseq peaks around GBP gene clusters. e, Co-localization between BRD4 and the GBP gene cluster by IF and DNA-FISH in Thp1-macro and Thp1-M.tb cells. IF, DNA-FISH, and merged channels (overlapping signal in white) are shown in separate images. The dashed line highlights the nuclear periphery, determined by DAPI staining. The rightmost column shows the area in the yellow box in greater detail. Each cell type, we counted 60 GBP gene loci, 18.3% (11/60) GBP loci were co-located with BRD4 in Thp1-macro cells and 65.0% (39/60) GBP loci were co-located with BRD4 in Thp1-M.tb cells.

Figure 5
Reorganization of chromatin architecture around NF-κB target gene sites. a, Overall TAD DoD of THP-1 cells during differentiation and infection. The y-axis represents the average DoD value of TADs in each state. b, The relationship between TAD DoD and chromatin loop during differentiation and activation. The x-axis represents log2 fold change of TAD DoD value, and the y-axis represent the number of loop changes in each TAD. c, Compare NF-κB enrichment around transcription start site (± 3 Kb) between strengthened loop anchor gene and random control gene. The gene in the heatmap were sorted according to the enrichment intensity of NF-κB. ChIP-seq data is from ENCODE (GM12891 cell line). d, Immuno uorescence analysis of NF-κB (p65) subcellular location before and after M.tb infection. e, Log2 fold changes of ATAC peaks and RNA expression level of NF-κB target genes (n=435) and random control genes. Signi cant differences were calculated by unpaired one-sided t test.    Identi cation of long range regulatory target genes of TB susceptibility loci and novel anti-TB drug targets by integrated multi-omics analysis. a, Identi cation of long range regulatory target genes of TB susceptibility loci, LRRK2 by integrated eQTL, Hi-C and GWAS multi-omics analysis. b, The rs1873613 SNP is located in the anchor of the LRRK2 gene chromatin loop, which is signi cantly strengthened after M.tb infection, as shown in the Hi-C matrix. c, Simulated typical single cell chromatin structure of SNP rs1873613 and LRRK2 gene region before and after M.tb infection. d, Regulation of LRRK2 gene expression by the rs1873613 SNP. The x axis indicates three genotypes of individuals from the Genotype-Tissue Expression (GTEx) dataset. e,f, Immuno uorescence analysis and quanti cation of the colocalization of the M.tb (H37Ra-RFP) and Rab7 in bone marrow-derived macrophages treated with the LRRK2 inhibitor AdoCbl and PBS (control). g, CFU assays of M.tb in bone-marrow-derived macrophage cells treated with AdoCbl and PBS (control