{'fraction_non_alphanumeric': 0.038736979166666664, 'fraction_numerical': 0.051310221354166664, 'mean_word_length': 4.913618864292589, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-coding germline GATA3 variants alter chromatin topology and contribute to pathogenesis of acute lymphoblastic leukemia

Inherited non-coding genetic variants confer significant disease susceptibility in many cancers. However, the molecular processes of by which germline variants contribute to somatic lesions are poorly understood. We performed targeted sequencing in 5,008 patients and identified a key regulatory germline variant in GATA3 strongly associated with Philadelphia chromosome-like acute lymphoblastic leukemia (Ph-like ALL). By creating an isogenic cellular model with CRISPR-Cas9 system, we showed that this variant activated a strong enhancer that significantly upregulated GATA3 transcription, which in turn reshaped the global chromatin accessibility and 3D genome organization. Remarkably, this genotype switch induced a chromatin loop between the CRLF2 oncogene and a distal enhancer, similar to the somatically acquired super-enhancer hijacking event in patients. GATA3 genotype-related alterations in transcriptional control and 3D chromatin organization were further validated in Ph-like ALL patients. Finally, we showed that GATA3 directly regulates CRLF2 and potentiates the oncogenic effects of JAK-STAT signaling in leukemogenesis. Altogether, our results provide evidence for a novel mechanism by which a germline non-coding variant contributes to oncogene activation epigenetic regulation and 3D genome reprogramming.


Introduction 51
Acute lymphoblastic leukemia (ALL) is the most common cancer in children and there is 52 growing evidence of inherited susceptibility to this hematological malignancy (Hunger and 53 Mullighan, 2015;Moriyama et al., 2015;Pui et al., 2015). In particular, genome-wide 54 association studies (GWAS) have identified at least 9 genomic loci (i.e., CDKN2A/2B, 55 IKZF1, ARID5B,CEBPE,GATA3,TP63,LHPP,and ELK3) with common 56 variants that influence ALL risk (Papaemmanuil et al., 2009;Perez-Andreu et al., 2013;57 Sherborne et al., 2010;Trevino et al., 2009;Xu et al., 2013;Xu et al., 2015). These 58 variants cumulatively confer a substantial increase of ALL risk (Xu et al., 2015), and 59 explain a large proportion of the estimated heritability of this leukemia (Enciso-Mora et al., 60 2012). Interestingly, some ALL germline risk variants also co-segregate with specific 61 acquired genomic abnormalities in leukemia (Perez-Andreu et al., 2013;Trevino et al., 62 2009;Walsh et al., 2013), suggesting intricate interactions between somatic and germline 63 mutations during leukemogenesis. In particular, we have previously reported germline 64 intronic variants in the GATA3 gene associated with the risk of developing Philadelphia 65 chromosome (Ph)-like ALL (Perez-Andreu et al., 2013), a subtype characterized by a 66 leukemia gene expression profile resembling that of Ph-positive ALL with BCR-ABL1 67 fusion (Den Boer et al., 2009;Roberts et al., 2014). Each copy of the GATA3 risk allele 68 increased the risk of Ph-like ALL by 3.25-fold. Because Ph-like ALL is associated with 69 distinctive genomic lesions in the cytokine signaling pathway genes (e.g., 50% of cases 70 harbor CRLF2-rearrangements) (Roberts et al., 2014), it raises the question whether 71 germline genetic variation in GATA3 is directly or indirectly involved in the deregulation 72 of this pathway in ALL. So far, the exact variant that determines GWAS signal at the 73

Discussion 275
Both inherited germline and somatic genetic variations contribute to the pathogenesis of 276 different malignancies, including leukemias. Somatic genomic aberrations, i.e., mutations, 277 rearrangements, insertion/deletion, have been shown to drive overt leukemogenesis by 278 promoting the survival and proliferation of pre-leukemia hematopoietic cells. However, 279 the roles of inherited leukemia risk variants, especially those in intronic/intergenic loci, 280 remain largely unclear. For example, GWAS studies have identified 9 genomic loci with 281 common SNPs associated with susceptibility to childhood ALL, but there has been little 282 progress to move from descriptive association studies to identifying causative 283 mechanisms relating these variants to ALL pathogenesis. 284 Here we define the regulatory function of a non-coding SNP rs3824462 associated 285 with Ph-like ALL (Perez-Andreu et al., 2013). This variant strongly influences the 286 susceptibility to high-risk ALL and also prognosis, consistently across different ALL 287 treatment regimens (Migliorini et al., 2013;Perez-Andreu et al., 2013). In this work, we 288 first reported that the rs3824662 variant is located inside an enhancer element and the 289 risk allele showed significantly increased enhancer activity. Introducing the risk A allele 290 at rs3824662 by CRSPR/Cas9 editing in the wildtype GM12878 cells directly confirmed 291 its enhancer effects on GATA3 transcription. Using a variety of chromatin conformation 292 capturing techniques, we further demonstrated that this variant significantly reshaped 293 chromatin interactions both locally and also in a global fashion. A recent study showed 294 that GATA3 can act as a pioneer factor in the course of cellular reprogramming, making 295 previously condensed chromatin more accessible by recruiting BRG1, a chromatin 296 remodeling factor (Takaku et al., 2016). Similarly, our ATAC-seq data also suggested that 297 the C-to-A allele substitution at rs3824662 resulted in many newly-gained open chromatin 298 regions enriched for GATA3 binding sites, coupled with global 3D geome re-orginazation. 299 In particular, we observed hundreds of regions switched from the active and open 300 compartment to the repressive and compacted compartment. Among them are many 301 essential genes whose expression are altered in Ph-like ALL, likely due to the change of 302 chromatin environment. We also performed ATAC-seq, GATA3 ChIP-seq and Hi-C in a 303 panel of seven ALL patient samples with different genotypes at rs3824462. In these 304 analyses, we identified similar B-to-A switching in PON2 genes and novel looping events 305 betweeen CRLF2 and P2RY8 locus, indicating that these transcriptional regulation 306 mechanisms are indeed operative in Ph-like ALL patients. However, these human 307 leukemia samples harbor a plethora of somatic genomic abnoralities which likely 308 confounded the effects from germline GATA3 polymorphisms. 309 More interestingly, we found many GATA3 binding sites are located near the 310 breakpoints of translocation events observed in Ph-like ALL, suggesting its over-311 expression might be related with chromosomal instability and susceptiblity to 312 translocations. Therefore, we hypothesize that GATA3 over-expression might facilitate 313 enhancer hijacking, where a distal enhancer is rearranged to the proximity of oncogenes 314 and leads to oncogenesis without gene fusions (Groschel et al., 2014;Hnisz et al., 2016;315 Northcott et al., 2014;Weischenfeldt et al., 2017). To further explore the role of GATA3 316 in genome instability, we also explored its binding profile in breast cancer, as GATA3 317 abnormal expression has also been reported in certain human breast cancer subtypes. 318 We confirmed that GATA3 binding is also enriched in a breast cancer cell line (T47D) 319 translocation breakpoints region as well ( Figure S16A). Moreover, we also observed 320 GATA3 and BRG1 co-localize at these translocation breakpoints ( Figure S16B and 321 S16C), suggesting potential intricate interactions between GATA3, BRG1 and genome 322 instability. 323 Aberrantly high GATA3 expression has been also identified in other B cell 324 malignancies, such as classical Hodgkin lymphoma. Constitutive activation of NFkB and 325 Notch-1 leads to higher GATA3 expression in Reed Sternberg cells, which then 326 contributes to cytokine secretion (especially IL13) and signaling typical in Hodgkin 327 lymphoma (Stanelle et al., 2010). In contrast, GATA3 is not expressed in normal B cells 328 and in fact functions as a key regulator of lymphoid cell lineage commitment (B vs T 329 cells) (Banerjee et al., 2013). The data we present in the current study points to novel 330 roles of GATA3 in global cellular reprogramming and pathogenesis of B-cell malignancies. 331 In conclusion, we report here that the inherited genetic variant rs3824662 is a cis-332 acting enhancer variant associated with GATA3 transcription activation, which contributes 333 to Ph-like ALL leukemogenesis through regulating CRLF2 signaling. Our results suggest 334 that transcription factor-mediated epigenomic reprograming can directly influence 335 oncogene activity, and may be an important mechanism by which germline genetic 336 were mapped and aligned using the Burrows-Wheeler Aligner (BWA) (Li and Durbin,360 2009), and genetic variants were called using the GATK pipeline (version 3.1) (Poplin et 361 al., 2017), as previously described, and annotated using the ANNOVAR (Wang et al., 362 2010) program with the annotation databases including RefSeq (O'Leary et al., 2016). All 363 the GATA3 non-silent variants were manually reviewed in the Integrative Genomics 364 Viewer (Robinson et al., 2011). Association of genotype with Ph-like ALL status was 365 examined following our stablished statistical procedure (Perez-Andreu et al., 2013)

3D chromatin structure mapping by Hi-C 378
Hi-C in GM12878 cells and PDX samples were performed using the Arima-HiC kit as per 379 the manufacturer's instructions. Briefly, 1 million GM12878 WT, A/A cells and PDX 380 sample were fixed with 1% formaldehyde, digested with restriction enzyme, end-labeled 381 with Biotin-14-dATP, and then followed by ligation. The ligated chromatin was reverse-382 crosslinked and sonicated by Covaris E220 to produce 300-500 bp fragments . Biotin 383 labeled DNA fragments were isolated using dynabeads Streptavidin C1 beads and 384 followed by end-repair, adenylation, adaptor ligation and PCR amplification. The quantity 385 of the library was measured by both BioAnalyzer (Agilent) and Kapa Library Quantification 386 Kit (Kapa Biosystems). Finally, the library was performed pair-end 2x100bp high-387 throughput sequencing using HiSeq 2500 and Nova-seq (Illumina). 388

Cytokine-dependent growth assay in Ba/F3 cells and Ba/F7 cells 389
The full-length GATA3 and CRLF2 coding sequence were purchased from GE Healthcare 390 and cloned into the cL20c-IRES-GFP lentiviral vector. cl20c-CRLF2-IRES-GFP was 391 modified into cl20c-CRLF2-IRES-CFP, and lentiviral supernatants were produced by 392 A, Genetic variants associated with Ph-like ALL at the GATA3 locus discovered by targeted-sequencing. The purple dot indicates the variant with strongest association signal (rs3824662), and the blue box represents the 2.7 kb region flanking this variant. B, Chromatin state annotations from the Roadmap Epigenomics Project. Enhancer marks (H3K4me1 and H3K27ac) and repressive mark (H3K27me3) are plotted across 42 cell and tissue types for the GATA3 genomic region. These epigenomic data suggest that rs3824662 is located inside a hematopoietic cellspecific enhancer element (red boxes). The upper panel shows the H3K27ac profiling in this region. The H3K27ac signals are averaged by different tissue-types and plotted in 100bp bins. C, Luciferase reporter assay comparing the enhancer activities of the fragments containing either the rs3824662 risk A allele or wildtype C allele in GM12878, an immortalized B lymphoblastoid cell line with normal karyotype. T bars indicate standard deviations. D, Allelic analysis of H3K4me1 ChIP-seq data in three lymphoblastoid cell lines with heterozygous genotype at rs3824662 (GM12119, GM19200, and GM19219). Orange and blue bars indicate the percentage of ChIP-seq reads from the A allele and the C allele, respectively. E, Open chromatin status at the rs3824662 locus in 7 ALL PDX samples of different genotypes, as determined using ATAC-seq. The bottom panel represents a 2.8 kb region flanking rs3824662. F, Allelic analysis of ATAC-seq data in three Ph-like ALL samples with heterozygous genotype at rs3824662. Orange and blue bars indicate the percentage of ATAC-seq reads from the A allele and the C allele, respectively (P3: 31 A vs 1 C reads, P4: 8 A vs 0 C reads, P5 : 6 A vs 1 C reads).      K 2 R 6 8 3 G )

Luciferase reporter gene assay
A 1,120-bp region encompassing the GATA3 SNP rs3824662 was amplified using CloneAmp HiFi PCR Premix (Clontech) (primer sequences in Supplementary Table 4) and then cloned into the pGL4.23-mini/P vector with a minimal SV40 promoter upstream of the firefly luciferase gene sequence. For reporter assays, 2x10 6 SUP-B15, GM12878, Ba/F3 cells were resuspended in 100 μl of Nucleofector Solution Kit V (Lonza) with the addition of 1.9 μg of pGL4.23 constructs and 100ng of renilla pTK plasmid. Cells were electroporated and then incubated at 37°C for 24 hours with 5% CO2. Similarly, HEK293T cells (6x10 4 ) were plated on 96-well plate (flat bottom), and co-transduced with 95 ng pGL4.23 constructs and 5ng renilla pTK, and then incubated for 24 hrs. Luciferase activity was measured using the Dual-Glo Luciferase Assay system (Promega). Experiments were performed in triplicate. To control for cell number and transfection efficiency, firefly luciferase activity was normalized to renilla luciferase. Measurements were presented as a ratio relative to the activity of the pGL4.23-mini/P empty vector.

Inducible GATA3 overexpression and GATA3 knockdown
Full-length GATA3 cDNA was cloned into the lentiviral vector pLV-tetON (a gift from Dr. Chunliang Li at St. Jude Children's Research Hospital). Nalm6 and SUP-B15 cells were transduced with pLV-GATA3-tetON for 48 hours and then subjected to bleomycin selection (0.5 mg/ml). Single clones were established in which doxycycline-induced GATA3 overexpression was confirmed by RT-qPCR and immuno-blotting.
The lentiviral pLKO.1 constructs with GATA3 shRNA and scrambled shRNA were purchased from Sigma-Aldrich. Nalm6 and SUP-B15 cells were lentivirally transduced with pLKO.1-GATA3-shRNA or scrambled shRNA for 48 hours and then subjected to puromycin selection (1.0 ug/ml). The degree of GATA3 and CRLF2 knockdown was evaluated by RT-qPCR.

RNA-seq
Total RNAs were extracted from 5 million engineered GM12878 cells using Trizol (Invitrogen). cDNA libraries were prepared using SureSelect Strand Specific RNA Library Preparation Kit (Agilent). Briefly, polyA RNA was purified from 1000 ng of total RNA using oligo(dT) beads (Invitrogen) and then fragmented, followed by reverse transcription, end repair, adenylation, adaptor ligation and subsequent PCR amplification. The final product was checked by size distribution and concentration using BioAnalyzer High Sensitivity DNA Kit (Agilent) and Kapa Library Quantification Kit (Kapa Biosystems). Pair-end 2x50bp high-throughput sequencing was performed using HiSeq 2500 (Illumina).

ATAC-seq
ATAC-seq was performed as previously described  . A total of 30,000-50,000 live cells were collected, washed once in PBS and resuspended in 50 ul ATAC-seq lysis buffer (10 mM Tris-HCl pH 7.4, 10 mM NaCl, 3 mM MgCl2 and 0.1% IGEPAL CA-630) followed by immediate centrifugation for 10 min at 4 °C. Cell pellets were then resuspended in 50 ul reaction solution which contained 1XTD buffer and 2.5 ul Tn5 transpose (Illumina, FC-121-1030) and incubated at 37 °C for 60 mins. The fragmented DNA was purified by MinElute kit (Qiagen) and amplified by PCR. The final product was checked by size distribution and concentration using BioAnalyzer High Sensitivity DNA Kit (Agilent). Pair-end 2x50bp high-throughput sequencing was performed using HiSeq 2500 (Illumina).

ChIP-seq and ChIP-qPCR
ChIP-seq and ChIP-qPCR were performed as previously described (Shen et al., 2012).

Capture-C
Capture-C was performed as previously described (Huang et al., 2017). Briefly, 10 million engineered or wildtype GM12878 cells were fixed with 1% formaldehyde and digested with DpnII (NEB), followed by DNA ligation. The ligated chromatin was then reversecrosslinked, purified by phenol-chloroform and followed by sonication to produce 200-300 bp fragments using Covaris E220. Fragmented DNA was used to make libraries with the NEBNext DNA Library Prep Master Mix Set (NEB). Hybridization with 60 bp biotinylated capture probes (Supplementary Table 5) was performed with the xGen® Lockdown® Reagents (Integrated DNA Technologies). In brief, 3C libraries were dried and resuspended with hybridization reagents. 3 pmol pooled capture probes was mixed with the resuspended libraries and incubated for 72 hr at 47°C. After streptavidin beads purification and PCR, the pulldown material was treated with a second round of 24-hr incubation to improve specificity. The capture probes, ordered from Integrated DNA Technologies, flank DpnII sites proximal to rs3824662.

Capture C process
Capture C data were processed by Hi-C-Pro (Servant et al., 2015), we required one of the pairs should be mapped to the rs3824662 anchor regions. CHiCAGO (Cairns et al., 2016) was used to assign significant interactions linked to the captured rs3824662 fragments.

ChIP-seq data process
Pair-end sequencing data were mapped to the female reference genome (hg19) using bowtie2(version 2.3.4.3) (Langmead and Salzberg, 2012). Non-uniquely mapping reads (MAPQ<30) were removed, and PCR duplicate reads were removed by Picard (http://broadinstitute.github.io/picard/). ENCODE-chip-seq-pipeline2 were followed to call the peaks, the narrow peaks with Poisson p-value greater than 0.001 were removed to ensure good quality peaks for further analysis. To further qualify the predicted peaks, Reads Per Million (RPM) of IP data and input data in each peak region were calculated and the qualified peaks should pass the threshold of two-fold enrichment (RPMIP/RPMinput>2) and RPMIP-RPMinput>1. To check repeatability between biological replicates, firstly we divided the reference genomes into 10kb bins and computed the number of reads within each bin. The Pearson correlation coefficient between each biological replication was calculated using above-normalized 10kb bins reads. IDR with a threshold of 0.05 was used to measure the reproducibility of peaks from replicates . The peaks were re-centered and set to a fixed width of 250 bp and identified differential GATA3 peaks using the DiffBind R package (Ross-Innes et al., 2012), The genome-wide ChIP-seq signal tracks were generated by MACS2(V2.2.4) (Zhang et al., 2008) for TFs and histone marks.

The GWAS hints enrichment analysis.
The GWAS hints were download from https://www.ebi.ac.uk/gwas/. 10,279 diseaseassociated and at least identified by 2 articles single-nucleotide polymorphisms were selected. Investigated their distribution at differential GATA3 peaks 10kb flanking regions.

ATAC-seq data process
Pair-end sequencing data were mapped to the female reference genome (hg19) using bowtie2(version 2.3.4.3) (Langmead and Salzberg, 2012) with -X 2000 parameter. Nonuniquely mapping reads (MAPQ<30) were removed, and PCR duplicate reads were removed by Picard ( http://broadinstitute.github.io/picard/). ENCODE-atac-seq-pipeline were followed to call the peaks, the narrow peaks with Poisson p-value greater than 0.001 were removed to ensure good quality peaks for further analysis. To check repeatability between biological replicates, firstly we divided the reference genomes into 10kb bins and computed the number of reads within each bin. The Pearson correlation coefficient between each biological replication were calculated using above-normalized 10kb bins reads. IDR with a threshold of 0.05 was used to measure the reproducibility of peaks from replicates . The genome-wide ATAC-seq signal tracks were generated by MACS2(V2.2.4) (Zhang et al., 2008)

Hi-C data process
The Hi-C data were aligned to the female reference genome(hg19) by bwa mem model (Li and Durbin, 2009) with -SP5M parameters. The PCR duplications and lowquality aligned pairs were removed by pairtools (https: //github.com/ mirnylab/pairtools), the "UU", "UR" and "RU" types pair were kept for further analysis. We generated 5kb, 10kb, 25kb, 40kb,50kb,100kb muti-resolutions balanced cool file and hic file for visualization. Correlations between Hi-C replicates were calculated HiCRep (Yang et al., 2017). We combined biological replicates of Hi-c data from each engineered and GM12878 clone.
A and B compartments were identified using previously described (Lieberman-Aiden et al., 2009) with modifications. We construct raw 10kb Hi-C contact matrix without normalization of each cell type and patient, then calculated the expected interaction frequency between any two bins given the distance separating them in the genome. The observed/expected matrix was generated and then converted to a Pearson correlation matrix. Principal component analysis is applied to the correlation matrix similar as previously described. The value on first principal component for each bin was used to correlate with ATAC-seq signal to assign a genomic region to A or B compartment. If the sign of PC1 value changed between engineered GM12878 cell lines with different genotypes at rs3824662, we considered it as the A/B switch region.
The insulation score was calculated by the Perl script matrix2insulation.pl (Record Owner) at 40kb resolution matrix with "-ss 80000 --im iqrMean --is 480000 --ids 320000" parameters. The topologically associated domains were identified by the Perl script insulation2tads.pl, the 0.3 of min boundary strength was set as threshold.
The interaction loops were identified by Peakachu (Salameh et al., 2019) in 10kb resolutions, the models for predicting loops using H3K27ac 10% and CTCF 10% model.
The predicted loops pass the probability score threshold great than 0.8. For the differential loops, we first calculated the probability score of each pair cross the genome, then we compare the probability score of predicted loops in sample A and probability score of pairs in sample B, we set 2-fold change as the cut off.
The Virtual 4C plot were used the embed method from the 3D genome browser (Wang et al., 2018). Briefly, a bait (for example, CRLF2 locus) and flanking region were chosen, then the row overlapping the bait and flanking regions were extracted from the Hi-C matrix. The number of observed contacts was plotted with a smoothing window to obtain virtual 4C profiles. To ensure the interaction frequency from different library comparable, the interactions in chromosome X were normalized by the number of interactions in viewpoints.

GATA3 and BRG1 binding in cancer breakpoints
Ph-like patient breakpoints and T47D cancer cell line breakpoints were collected from previous study (Dixon et al., 2017;Reshmi et al., 2017) and expanded to +/-1kb region.
Motif enrichment analysis was performed using HOMER version 4.8 findMotifsGenome.pl function for exploring potential transcription factor binding within these expanded breakpoint regions. GATA3 binding signal in these expanded breakpoint regions was generated by deeptools using GATA3 ChIP-seq in the following cells: GM12878 C/C clone, GM12878 A/A clone, Nalm6, Nalm6_gata3overexpressing and T47D breast cancer cell (Adomas et al., 2014). BRG1 ChIP-seq in T47D breast cancer cell were collected from GSE112491 are plotted the same way.    and A allele at rs3824662 are from the same allele. Therefore, allelic expression derived from rs2229359 would directly inform the differential transcription activation effects of the A vs. C allele at rs3824662 in engineered GM12878 clones. B. Sanger sequencing of PCR products of GATA3 3 rd exon cDNA shows allelic expression of GATA3 in two GM12878 heterozygous clone cells by rs2229359 genotyping. C. Shown the transcription level associated with rs3824662-A allele vs. the transcript associated with wild type C allele (p value =0.0066 by student t-test).      ) in Ph-like ALL patient translocation breakpoints region by ChIP-seq (+/-1kb) in Nalm6 GATA3 ov (yellow) and Nalm6 GATA3 wt (blue) cells . Upright panel is GATA3 binding signal in 1000 random genomic regions in Nalm6 GATA3 ov (yellow) and Nalm6 GATA3 wt (blue) cells. GATA3 alone, JAK2 R683G alone, GATA3 with JAK2 R683G , or empty vector control. All the experiments were performed in triplicates (p value < 0.001 by 2way ANOVA). Ba/F7 cells with GATA3 and JAK2 R683G were treated with or without 10 ng/ml TSLP. All the experiments were performed in triplicate for three times independently (p value < 0.001 by 2way ANOVA).