Drosophila stock maintenance<\/p>\n
PCNA::EGFP flies used for interphase-staged Pico-C were generously provided by S. A. Blythe and E. Wieschaus (Princeton University). The sfGFP-Zelda and GAF-sfGFP flies were generated as described previously11<\/a>,60<\/a>. Jabba-Trap was expressed from the maternal \u03b1-tubulin 67 promoter, kindly provided by S. A. Blythe61<\/a>. Flies were maintained on standard cornmeal-agar food and raised at 25\u2009\u00b0C.<\/p>\n Embryo fixation and collection<\/p>\n A full step-by-step protocol covering embryo fixation can be accessed at protocols.io (https:\/\/doi.org\/10.17504\/protocols.io.rm7vz9ro5gx1\/v1<\/a>). Briefly, PCNA::EGFP flies were allowed to lay eggs on 0.4% acetic acid agar plates at 25\u2009\u00b0C for 1\u2009h pre-collection (NC12\u2013NC14) or 2\u2009h pre-collection (NC9\u2013NC11). After egg laying, plates were incubated for 30\u2009min (NC9\u2013NC11) or 1\u20131.5\u2009h (NC12\u2013NC14) before fixation. Embryo fixation was performed following the published Micro-C protocol57<\/a>, with the exception that EGS was used as the sole secondary crosslinker. Embryos were first crosslinked in 1% formaldehyde, quenched with Tris-HCl (final 0.75\u2009M), washed in phosphate-buffered saline with Triton X-100 (0.5% final concentration) and then crosslinked in 3\u2009mM EGS for 45\u2009min at room temperature. After a second quench, embryos were washed, staged and sorted. Interphase embryos showed clear nuclear PCNA signals and NC was assigned based on nuclear density. Sorted embryos were snap-frozen in liquid nitrogen and stored at \u221280\u2009\u00b0C.<\/p>\n Pico-C<\/p>\n The complete step-by-step protocol, including the micrococcal nuclease digestion test, is available at protocols.io (https:\/\/doi.org\/10.17504\/protocols.io.kqdg31nm1l25\/v1<\/a>). Pico-C libraries were generated following the published Micro-C protocol57<\/a>,62<\/a> with several key modifications. Nuclei were immobilized on concanavalin A-coated magnetic beads to increase DNA recovery and eliminate centrifugation steps before DNA extraction. Approximately 10 embryos (~60,000 nuclei) were used for the NC14 libraries, with input doubled at each earlier stage (up to ~300 embryos for NC9). After nuclear extraction, samples were resuspended in 900\u2009\u00b5l of MB1 and 100\u2009\u00b5l of 10\u00d7 binding buffer (200\u2009mM Hepes-KOH, pH 8, 100\u2009mM KCl, 10\u2009mM CaCl2, 10\u2009mM MnCl2 and 5\u2009mM spermidine) along with 10\u2009\u00b5l of pre-washed concanavalin A-coated beads. The mixture was rotated for 10\u2009min at room temperature to allow nuclei to bind the beads.<\/p>\n Bound nuclei were digested with micrococcal nuclease at 37\u2009\u00b0C for 10\u2009min with shaking (950\u2009rpm). The reaction was stopped by adding EGTA ((ethylenebis(oxonitrilo))tetra-acetate, final 20\u2009mM concentration) and samples were washed 3\u00d7 with cold MB2 buffer (50\u2009mM NaCl, 10\u2009mM Tris-HCl, pH 7.5, and 10\u2009mM MgCl2) while kept on the magnet. End-repair and phosphorylation were performed as previously described62<\/a> using T4 PNK and Klenow fragment.<\/p>\n Biotinylation was extended to 4\u2009h to improve ligation efficiency and reduce downstream PCR requirements. At the 2-h mark, 1.5\u2009\u00b5l of 100\u2009mM ATP and 50\u2009U of Klenow fragment were added and incubation continued for the remaining 2\u2009h at 25\u2009\u00b0C. The reaction was stopped with 0.5\u2009M EDTA and samples were washed twice with 500\u2009\u00b5l of cold MB3 buffer (50\u2009mM Tris-HCl, pH 7.5, and 10\u2009mM MgCl2).<\/p>\n Proximity ligation was carried out using T4 DNA ligase (10,000\u2009U total, added in two steps) for 5\u2009h at room temperature. Unligated ends were treated with exonuclease III to remove free biotin-dNTPs. Crosslinks were reversed by overnight incubation at 65\u2009\u00b0C in the presence of proteinase K and sodium dodecylsulfate (SDS). DNA was purified by extraction with phenol and chloroform. To maximize yield from low-input samples, total DNA was carried forward into the final library construction phase and size selected for dinucleosomal fragments (350\u2013500\u2009bp) after PCR amplification using Ampure XP beads. Final libraries were eluted in 10\u2009mM Tris-HCl, pH 8, and quantified by Qubit before sequencing.<\/p>\n Embryo injections<\/p>\n Injections were performed as previously described18<\/a>. Briefly, \u03b1-amanitin (0.5\u2009mg\u2009ml\u22121 in water) or water alone was injected evenly into a PCNA::EGFP embryo. Embryos older than NC8 after injection were discarded. The remaining embryos were monitored until NC14, then either fixed and staged for Pico-C or pooled (10 embryos per batch) in 30\u2009\u00b5l of TRIzol and flash-frozen for RNA extraction. All embryos were hand sorted to ensure correct staging. A subset of embryos was allowed to develop, revealing consistent arrest in amanitin-injected embryos, whereas controls developed normally as shown before18<\/a>.<\/p>\n RNA-seq in embryos<\/p>\n For each replicate, 40 embryos were pooled and homogenized with a metal pestle. RNA was extracted by chloroform phase separation and purified from the aqueous layer using the NEB Monarch RNA Cleanup Kit (NEB, cat. no. T2030). All samples showed RNA integrity no. >7 on a Bioanalyzer. Libraries were prepared with the Watchmaker mRNA Enrichment and Library Prep Kits according to the manufacturer\u2019s protocol, dual indexed, pooled and sequenced on an Illumina NextSeq 2000 to a depth of \u226580 million paired-end, 60-bp reads per sample.<\/p>\n Single-embryo ATAC\u2013seq on embryos depleted for Zelda and GAF<\/p>\n Single stage 5 (NC14) embryos of the genotype sfGFP-Zld;Jabba-Trap\/CyO;GAF-sfGFP and sfGFP-Zld; Sp\/CyO;GAF-sfGFP controls, identified based on morphology, were harvested and processed for ATAC\u2013seq as described previously11<\/a>. Libraries were prepared from four replicates of both experimental and control embryos and paired-end sequencing was performed on a NovaSeq X Plus.<\/p>\n Immunohistochemistry of embryos depleted for Zelda and GAF<\/p>\n The sfGFP-Zelda;Jabba-Trap\/CyO;sfGFP-GAF and sfGFP-Zelda;Sp\/CyO;sfGFP-GAF control flies were put in cages, allowed to lay for 1.5\u2009h and the embryos were aged for an additional hour at 25\u2009\u00b0C. Drosophila embryo collection, fixation and staining were performed as previously described11<\/a>. The following antibodies were used: rabbit anti-GFP (Abcam, cat. no. ab290, 1:500); DyLight 488-conjugated goat anti-rabbit (Thermo Fisher Scientific, cat. no. 35552, 1:1,000); and DAPI (Thermo Fisher Scientific, cat. no. D1306, 1:1,000). Images were taken using a Nikon Ti-2e Epifluorescent microscope or Nikon A1R+ confocal and processed by Fiji or ImageJ63<\/a> and Adobe Photoshop software.<\/p>\n Pico-C data processing<\/p>\n Pico-C reads were mapped to the dm6 genome using BWA-MEM v0.7.17 (ref. 64<\/a>). Aligned reads were then processed with Pairtools v1.0.3 (ref. 65<\/a>) using the parse2 function to rescue complex walks and any reads below a mapping quality of 3 were removed. The resulting pairs files were sorted and converted to FAN-C v0.9.18 (ref. 66<\/a>) pairs format. Duplicates were removed and quality control metrics were extracted.<\/p>\n We then generated FAN-C .hic files. Principal component analysis (PCA) was performed on these .hic files using FAN-C, with a sample size of 100,000 interactions, focusing on distances between 100\u2009kb and 1\u2009Mb to assess variation and clustering across the samples. As biological replicates clustered well across all stages (Extended Data Fig. 1a<\/a>), the .hic files of replicates were merged. The merged files were binned at different resolutions and normalized using iterative correction. We then used a customized Python script to determine the bin size, where 80% of bins had >1,000 reads, as done previously67<\/a>, to estimate the optimal resolution of each stage. Finally, the merged FAN-C .hic files were converted to balanced Cooler .mcool files using the cooler zoomify utility from Cooler v0.9.3 (ref. 68<\/a>), with the options –balance –balance-args \u2018–mad-max 0 –max-iters 1000\u2019 to perform iterative matrix balancing.<\/p>\n Pico-C, Micro-C and Hi-C coverage tracks<\/p>\n Genome-wide coverage tracks from Pico-C data were generated using cooltools v0.5.3 (ref. 69<\/a>). Full or downsampled matrices (–count 100000000 using random sample) were used to compute coverage at multiple resolutions with cooltools coverage, producing raw TSV and BigWig outputs. Raw BigWigs were normalized to counts per million (CPM) using wiggletools and converted to the BigWig format via bedGraphToBigWig.<\/p>\n For quantifying signal across chromatin states, per-bin read counts were extracted from 1-kb cooler matrices, intersected with state annotations and scaled by overlap size to obtain total counts per state. Counts were then normalized by region length and total mapped reads to derive bins per million coverage, analogous to transcript per million normalization in RNA-seq.<\/p>\n Loop and boundary identification<\/p>\n Boundaries were detected from equally downsampled .hic files using FAN-C (500-bp bins, 16-kb window). Boundaries with insulation >0.45 and outside centromeric or low-mappability regions were retained. Stage-specific comparisons were performed in R using GenomicRanges v1.56 (ref. 70<\/a>) and ComplexUpset v1.3.3 (ref. 71<\/a>)with a 5-kb tolerance. All boundary calls are provided in Supplementary Data 1<\/a>.<\/p>\n Loops were identified using Mustache v1.3.2 (ref. 27<\/a>) at 250-bp, 1-kb, 2-kb, 4-kb and 16-kb resolutions to capture looping dynamics across distances. Loops with P\u2009>\u20090.01 or in low-mappability or low-centromeric regions were removed. Resolution-specific calls were merged iteratively with Bedtools v2.29.2 (ref. 72<\/a>) (pairToPair -slop 1000) to create a unified loop set. Stage-specific loop overlap (NC9\u2013NC14) was assessed with hicVennDiagram v1.2 (ref. 73<\/a>) applying the same overlap thresholds used for boundaries. All loop calls are provided in Supplementary Data 2<\/a>.<\/p>\n Compartment analysis<\/p>\n To identify chromatin compartments from Pico-C data, we used CALDER2 v0.7 (ref. 35<\/a>).<\/p>\n NC14 intrachromosomal contact matrices were generated with fanc dump at multiple resolutions (1\u201350\u2009kb) and converted to tab-delimited format for CALDER2. A log2(input-normalized H3K36me3)42<\/a>,74<\/a> BigWig track was used to define the correct A or B orientation. A 3-kb resolution was chosen because it yielded the strongest correlation (\u22650.4) with H3K36me3. Subcompartments were simplified to A.1, A.2, B.1 and B.2 for downstream analyses, excluding previously characterized centromeric regions16<\/a>. Compartment calls are provided in Supplementary Data 6<\/a>.<\/p>\n Aggregate analysis<\/p>\n Aggregate analysis for loops calls, boundaries and compartments was done by first calculating the expected interaction frequencies using our Pico-C matrices with cooltools expected-cis. Pileup analysis was then carried out using coolpup.py v1.1.0 (ref. 75<\/a>).<\/p>\n For boundary local pileups, the score in the right corner reflects the average insulation\u2014calculated by dividing the signals in opposing corners of the aggregate as described in the documentation. Compartment pileups were rescaled using the parameters –rescale –rescale_flank 1 –rescale_size 99 to normalize feature size and flanking regions.<\/p>\n Clustering of boundaries and loop anchors<\/p>\n The signal of various epigenetic marks and transcription was extracted for a 5-kb window around the center of individual boundaries and loop anchors (anchors within 1\u2009kb of each other were merged and treated as a single anchor). An enrichment score for each region was calculated by dividing the average signal of the center (\u2212500\u2009bp to +500\u2009bp) by the average signal from the flanking region (\u00b11\u2009kb to \u00b12.25\u2009kb). Hierarchical clustering of boundaries and loop anchors was performed on the enrichment scores using Euclidean distance and Ward\u2019s method. The number of clusters was selected empirically as the largest number that produced epigenetically distinct groups. For visualization purposes, UMAP dimensionality reduction was applied to the region by enrichment score matrix. Loops >500\u2009kb were excluded from this analysis, because they tend to show higher noise and weaker signal due to distance-dependent decay (16.1% of total loops). Cluster results are provided in Supplementary Data 3<\/a>.<\/p>\n ChromHMM input data processing<\/p>\n Public datasets were processed using standardized pipelines specific to each data type.<\/p>\n For ChIP\u2013seq data, adapter trimming and deduplication were performed with fastp v0.23.1, followed by alignment to dm6 using Bowtie2 v2.5.4 (ref. 76<\/a>) (–no-unal, –no-mixed and –no-discordant). SAM files were converted and sorted with Sambamba v1.0.1 (ref. 77<\/a>) and reads with mapping quality <20 or overlapping black-listed regions were removed. ATAC\u2013seq data were processed identically, with additional use of deepTools v3.5.678<\/a> alignmentSieve (–ATACshift). Open chromatin reads (<100\u2009bp) and nucleosomal fragments (180\u2013250\u2009bp) were separated. RNA-seq data were aligned using STAR v2.7.10a79<\/a>. The resulting BAM files were indexed using Sambamba.<\/p>\n Quality control checks for all datasets were performed using Samtools v1.18 (ref. 80<\/a>) and FastQC v0.12.1 (ref. 81<\/a>).<\/p>\n ChromHMM<\/p>\n ChromHMM v1.24 (ref. 29<\/a>) was used for chromatin state analysis. Datasets were grouped into five categories: pre-minZGA (NC1\u2013NC8), pre-majZGA(a) (NC9\u2013NC11), pre-majZGA(b) (NC12\u2013NC13), ZGA (NC14) and Zelda\u2212, which included Zld knockdown data from NC12\u2013NC14. For pre-majZGA stages, where data were available across both NC9\u2013NC11 and NC12\u2013NC13, shared datasets were used; otherwise, stage-specific data were included where available. Input controls were specified in the cell-marks file and, for RNA signal, flowthrough maternal RNA-seq data7<\/a> were used to minimize maternal transcript influence.<\/p>\n Deduplicated BAM files were used for binarization and model training (with LearnModel options -b 100 -r 400). Models with varying state numbers were tested and a 20-state model was selected based on log(likelihood). States were annotated with ChIPseeker v1.40 (ref. 82<\/a>) using enrichment across gene bodies, chromosomal distribution and ChromHMM emission probabilities, and ordered by developmental enrichment. Gene ontology enrichment was performed with clusterProfiler v4.12 (ref. 83<\/a>) using enrichGO (BP ontology, Padj\u2009<\u20090.05), with terms simplified and visualized as dot plots highlighting key biological processes. De novo motif enrichment was assessed using HOMER v4.11 (ref. 30<\/a>) (findMotifsGenome.pl) on regions showing peak enrichment for each state. ChromHMM tracks are provided in Supplementary Data 5<\/a>.<\/p>\n Metaloci analysis<\/p>\n ChIP\u2013seq and RNA-seq data were processed per replicate using deepTools bamCompare to generate log2(ratio BigWigs) when inputs are available. For RNA-seq, flowthrough was used as input to control for maternal effects. Replicates were merged using bigwigAverage. The dm6 genome was divided into 200-kb windows (100-kb step), excluding bins with >5% unmappability and\/or centromeric overlap.<\/p>\n Pico-C and genomic signals were integrated using an adapted METAloci framework40<\/a>. Spatial autocorrelation was applied to assess correlations between chromatin interactions and genomic features. Stage comparisons (NC12 versus NC14) used two-sided Wilcoxon\u2019s tests (*P\u2009<\u20090.05, **P\u2009<\u20090.01, ***P\u2009<\u20090.001) and violin plots show global Moran\u2019s distributions with mean values highlighted.<\/p>\n RNA-seq analysis (\u03b1-amanitin-treated versus water-treated embryos)<\/p>\n Adapters and low-quality bases were trimmed using fastp v0.24.0 (ref. 84<\/a>) with –detect_adapter_for_pe and reads were aligned to the dm6 genome using STAR v2.7.11b (–outFilterMismatchNoverLmax 0.04, –outFilterMultimapNmax 1 and –outSAMtype BAM SortedByCoordinate). Only uniquely mapping reads were retained. Gene-level quantification was performed with featureCounts v2.0.8 (ref. 85<\/a>) in paired-end mode (-p).<\/p>\n Differential expression analysis between \u03b1-amanitin-treated and water-treated embryos was performed using DESeq2 v1.44.0 (ref. 86<\/a>). Genes with \u226510 counts in at least 2 samples were retained.<\/p>\n A set of stable maternal transcripts based on prior annotations8<\/a> was used for normalization (\u22651,000 counts in \u22652 samples, |log2(fold-change)|\u2009<\u20091, log(counts per min)\u2009>\u20091). Variance-stabilized counts were used for visualization and significance was defined as Padj\u2009<\u20090.05.<\/p>\n Scale-factor-normalized BigWig tracks were generated with deepTools (bamCoverage, bin size\u2009=\u20091) and averaged across replicates using bigwigAverage.<\/p>\n Chromosight analysis<\/p>\n Chromosight v1.6.3 (ref. 44<\/a>) was used to quantify boundary and loop strengths using –subsample to have equally sampled contact maps from \u03b1-amanitin-injected, water-injected and Zelda or GAF DKD embryos and matched controls. For boundary analysis, scores were computed at 500-bp resolution using –pattern borders. Loop strengths were quantified using NC14 loop BEDPE coordinates across multiple resolutions (1\u20138\u2009kb) to capture interactions at varying genomic distances. The highest-resolution score was retained per loop. Loop scores were analyzed at the anchor level by assigning each loop anchor independently to predefined merged loop-anchor cluster regions.<\/p>\n Scores were compared between samples and controls in R using unpaired Wilcoxon\u2019s rank-sum tests (***P\u2009<\u20090.001, **P\u2009<\u20090.01, *P\u2009<\u20090.05). Scores were grouped by boundary or loop\u2013anchor cluster, with anchors participating in multiple loops linked to more than one score. To track feature score trends, normalized score differences (\u22121 to 1) were calculated per locus-positive values indicating stronger signals in controls and negative in experimental samples. Outliers were removed using the interquartile range method and mean and median differences were summarized per cluster.<\/p>\n Model training and evaluation<\/p>\n We adapted the second-stage Orca model46<\/a>, composed of a hierarchical sequence encoder and a multilevel cascading decoder, to predict the chromatin structures at 125-bp, 250-bp, 500-bp and 1-kb levels from 250-kb DNA sequences. NC12 and NC14 Pico-C data were used and genomic sequences were retrieved from the dm6 reference genome. All chromosomes except chr2L:0\u201319\u2009Mb were used for training, with this region reserved for validation.<\/p>\n We adopted the same training strategy as Orca46<\/a>. The training process with stochastic gradient descent took about 1,500,000 steps (250-kb sequence with batch size 8 and learning rate 0.01 with momentum 0.98) on a server with 4 NVIDIA Tesla v100 (32\u2009GB) graphics processing units. Model performance on the holdout region was assessed by concatenating and flattening 1-kb prediction matrices and calculating Pearson\u2019s correlation between predicted and observed Pico-C contact maps.<\/p>\n Multiplexed in silico mutagenesis<\/p>\n We employed the same in silico mutagenesis screening approach utilized in Orca46<\/a> at 1-kb resolution to identify individual motifs critical for genome interaction. The disruption impact on local genome interactions is measured by 250-kb structural impact score, which is the average absolute log(fold-change) of interactions between the disruption position and all other positions in the 250-kb window.<\/p>\n For motif enrichment analysis, Drosophila nonredundant motifs were downloaded from the JASPAR 2024 database87<\/a>. Motif matches for each 10-bp site were scanned for after extending by 10-bp flanking sequence on each side, and a maximum log(odds score) over the 30-bp window for each motif was calculated. To quantify the enrichment of motifs, two-sided Student\u2019s t-test (without assuming equal variance) was performed to compare the motif log(odds scores) of the sites with impact score >0.02 for NC14 and >0.01 for NC12 against the background of 100,000 sites randomly drawn among all 10-bp sites screened. Fold enrichment was also computed on the same sites with a motif log(odds threshold) of 10.<\/p>\n Motif-specific mutagenesis<\/p>\n To assess how specific motifs influence global genome architecture, we analyzed boundary regions following genome-wide mutagenesis. The validation region was tiled into 250-kb windows (125-kb step) and accessible motif instances were defined by intersecting ATAC\u2013seq peaks with motif matches. Each site was extended to 14\u2009bp and randomly mutated to ensure consistent sequence length.<\/p>\n Boundaries identified in the original Pico-C data were used for pileup analysis (log(fold-change) over background within 20\u2009kb). Boundary strength was quantified by the insulation score, defined as the difference between the sums of the diagonal versus off-diagonal quadrants of the interaction matrix centered on each boundary. The insulation change was computed as the difference in the score before and after mutation.<\/p>\n Five motifs\u2014Zelda, GAF, M1BP, Dref and sqz\u2014were tested individually and in combination. Pairwise effects were assessed by comparing insulation changes at shared boundaries and combined effects were evaluated by mutating all accessible instances of the five motifs within each window and comparing total insulation changes to the sum of individual perturbations.<\/p>\n Structural impact score distribution<\/p>\n Mean structural impact scores were calculated across ChromHMM state regions by averaging scores within annotated regions. Promoter-like states (7 and 19) were further subdivided based on overlap with TSSs. Differences across states were tested using ANOVA (P\u2009<\u20092\u2009\u00d7\u200910\u2212\u00b9\u2076) and confirmed with a nonparametric Kruskal\u2013Wallis test (P\u2009<\u20092.2\u2009\u00d7\u200910\u2212\u00b9\u2076). Pairwise comparisons were performed using Tukey\u2019s honest significant difference (HSD) test and the results were visualized with mean scores\u2009\u00b1\u200995% confidence intervals and a compact letter display indicating statistically distinct groups.<\/p>\n Effect size distribution<\/p>\n To quantify the effect size of factors identified as significant in our machine learning analysis, we first gathered the predicted effect data from in silico mutagenesis for each factor and a corresponding set of random regions as controls. The effect size for each factor was determined by calculating the mean difference between the factor-specific predictions and the control values.<\/p>\n State local enrichment analysis at loop anchors and boundaries<\/p>\n Chromatin state enrichment within loop anchors and boundaries was assessed across developmental stages. Loop and boundary coordinates were intersected with stage-specific ChromHMM annotations and local background regions were defined by extending each feature \u00b1150\u2009kb. Enrichment was calculated as the observed\/expected ratio of each state within loops or boundaries relative to their local background, log2(transformed).<\/p>\n ATAC\u2013seq analysis<\/p>\n Paired-end ATAC\u2013seq reads were processed as above. BAM files were filtered (MAPQ\u2009\u2265\u200910), deduplicated and merged with sambamba. Tn5 insertion bias was corrected using deepTools alignmentSieve –ATACshift. Coverage tracks were generated with bamCoverage (RPGC normalization, 5-bp bins, blacklist excluded).<\/p>\n Accessible regions were identified with HMMRATAC v1.2.10 (ref. 88<\/a>). We note that, in DKD embryos, HMMRATAC detected an increased number of peaks, predominantly corresponding to accessible regions in controls. HOMER de novo motif enrichment was performed on regions that lost accessibility in DKD embryos. Called peaks are provided in Supplementary Data 7<\/a>.<\/p>\n Data visualization<\/p>\n Chromatin contact data were visualized using HiGlass v0.4.7 (ref. 89<\/a>). pyGenomeTracks v3.8 (ref. 90<\/a>) was used to display chromatin contact data in Figs. 1<\/a> and 3<\/a> and Extended Data Figs. 1<\/a> and 6<\/a>, using a natural log with scale factor set to 100. All other panels displaying chromatin contact data were plotted with CoolBox v0.3.8 (ref. 91<\/a>) using a log10 scale.<\/p>\n Statistics and reproducibility<\/p>\n Statistical analyses were performed in R (v4.4.0) and data visualization was carried out using the ggplot2 package (v3.5.1). Box plots were defined with boxes spanning from the first quartile (25th percentile) to the third quartile (75th percentile), with the median represented by a horizontal line within the box. Whiskers extend to the most extreme data points within 1.5\u00d7 the interquartile range from the box. No collected data were excluded and data collection and analysis blinding were not applicable to this study. All statistical tests were two sided.<\/p>\n Reporting summary<\/p>\n