dataesr/scientific-paragraphs-categorization huggingface dataset:
cnf1-positive strains of E. coli ST131 segregate between two clade-specific virulence profiles We then defined strain contents in virulence factors (VF) and acquired antibiotic-resistance genes (RG) to perform an unbiased analysis of their distribution into clusters, using a latent block model approach. Acquired antibiotic-resistance genes in ST131 genomes were identified with ResFinder 46. Profiles of virulence factors were defined with the database published by Petty and colleagues  31. The unsupervised clustering procedure retained a total of 10 RG-clusters and 7 VF-clusters (Figure  3A ). Differences in number of VFs and RGs between clusters were all significant (Figure  3B ). We found that cnf1-positive strains were scattered among several RG clusters (Figure  3A , left panel). By contrast, most cnf1-positive strains segregated into the cluster VF4 (84% of cnf1-positive strains, n=609) with the remaining 16% strains being distributed between VF1 (15%) and other VF clusters (1%) (Figure  3A , right panel).
cnf1-positive strains display dominant expansion in ST131-H30Rx/C2 We next analyzed the temporal distribution of cnf1-positive strains within clades and subclades. Using a Generalized Linear Models (GLM) approach, we first verified within our dataset the increase of fimH30-positive isolates over time (clade C) in E. coli ST131 that was maximal in H30Rx/C2 (P<2 10 -16 ) (Figure  4A ). We also noted a significant increase in the proportion of cnf1-positive strains over time in E. coli ST131 (Figure  4B , top panel). The GLM was then fitted on years, clades, and subclades. We tested the significance of the year effect and P-values were corrected for multiple comparisons using Tukey's method. The year effect was not significant for clade A, B, or subclade H30R/C1 (Figure  4B ). Instead, we observed a significant increase of the proportion of cnf1-positive strains within H30Rx/C2 over time (P=1.
25 10 -11  ). In addition, the GLM fitted curves predicted that the prevalence of cnf1positive strains within H30Rx/C2 sublineage would be approximately 50% (confidence interval of 95% [43% to 58%] in 2018; [47% to 64%] in 2019). Predictive values were confronted to the prevalence of cnf1 in ST131 strains isolated in 2018 or 2019 in a second independent dataset up-loaded from EnteroBase in September 2020. This confirmed the rising prevalence of cnf1-positive strains within the sublineage H30Rx/C2 up to 45% in 2018 and 48% in 2019. In conclusion, we identified a dominant expansion of cnf1-positive strains within ST131-H30Rx/C2.
DISCUSSION Initially thought to be absent in the Escherichia coli ST131 lineage, the cnf1 gene was estimated to be found in approximately 15% of this lineage, among 99 isolates from distinct geographical locations across the world, in 2014  31, 52. Large-scale genetic analysis of more than five thousand isolates of E. coli ST131 from EnteroBase, a database widely used by clinicians, gives here sufficient statistical power to unveil a dominant expansion trend of cnf1-positive strains within clade H30Rx/C2. Our analysis supports the hypothesis of a recent expansion of a large phylogenetic subcluster of cnf1-positive ST131-H30Rx/C2 strains circulating between humans and dogs  53, 54. In addition, we document a stable population dynamic of cnf1-positive H30R/C1 strains within clade C1. This raises the question of whether cnf1 confers a fitness advantage at the population level.
Our compelling findings ascribed such a feature of cnf1 to specific genetic backgrounds, thereby enhancing the expansion and dissemination of a subpopulation of ST131-H30Rx/C2 within the ST131 lineage. Furthermore, we report the high prevalence of cnf1 gene in the three sequence types ST73, ST12 and ST127 of E. coli that have different antibiotic resistance profiles. Together, this points to a role of cnf1 in the dynamics of ExPEC that is independent from antibiotic resistance genetic backgrounds. The rising prevalence of cnf1-positive H30Rx/C2, and evidence of their mobilization between humans and dogs  53  , suggest that cnf1 enhances the dissemination of H30RxC2 within households with companion animals, which is likely driven by an increased ability to compete for GIT colonization. In further support of this conclusion, we found a prevalence of 24% of cnf1-positive strains in the group companion animals from the EnteroBase database.
ST131 dataset structure and phylogenomic analysis The database used for phylogenetic and statistical analyses consists of whole-genome sequences of E. coli ST131 isolates collected by mining EnteroBase from 1967 to 2018  41. Leaning on Find ST(s) tool from EnteroBase, we retained a total of 5,231 genome assemblies and associated metadata, including information of the isolation date, country and source of isolates (Sup. table  5 ). Phylogeny of ST131 isolates was resolved using core non-recombinant SNPs defined with Parsnp (in total 37,304 SNPs)  67  and Gubbins v2. 3. 4  68. A maximumlikelihood tree was then estimated with RAxML v8. 2. 8 applying a general time-reversible substitution-model with a gamma distribution rate across sites and with an ascertainment bias correction  69  and the resulting tree was edited with the interactive Tree of Life (iTol) v4 program  70.
In silico antimicrobial resistance and virulence-associated markers GyrA and ParC protein sequences were retrieved from the EnteroBase annotated genomes, and aligned with the mafft L-INS-I approach  71. After a visual inspection of the alignment, inhouse customized perl scripts (https://github. com/rpatinonavarrete/QRDR) were used to identify the amino acids at the quinolone resistance-determining region (QRDR) (positions 83 and 87, and 80 and 84 in GyrA and ParC, respectively). Search for cnf1 and hlyA alleles in ST131 genomes dataset was carried out by Blastn analysis. Sequences were next aligned with Muscle 72 and curated to remove incomplete sequences. SNPs were then extracted using SNP-sites  73. To determine strain specific VF profiles, annotated VFs from UPEC described in  31  were translated and pBLASTed against ST131 genomes dataset considering only hits with e-value < 10 -5 and identical matches > 95% (sup.
Co-clustering method Statistical analyses were performed using R software version 3. 6. 0. A total of 20 strains from the collection of 5,231 strains of E. coli ST131 were removed from the analysis due to incomplete associated metadata. The clustering of strains with specific virulence or acquired antibiotic-resistance gene profiles was performed with binary latent block model, implemented in the R package blockcluster  75. In this package, the model, a mixture of Bernoulli distributions proposed by  76  , is estimated using an efficient EM algorithm. As proposed by the authors, the number of clusters was estimated by maximizing the ICL criterion on a bidimensional grid of parameters making this unsupervised classification procedure automatic.
Pan-genome analysis The pangenome of E. coli ST131 was estimated using Roary, a high-speed pan genome pipeline analysis tool  77. Roary returns as output, the gene presence/absence matrix. The matrix was curated to retain genes present in at least 50 genomes and less than 3980 genomes (7678 sequences), that constituted our accessory genes pool dataset. Hierarchical clustering analysis was then conducted by using the pheatmap package in R (cran. rproject. org/web/packages/pheatmap/index. html). The gene presence/absence file generated by Roary was further analyzed using Scoary 45 with a significant Bonferroniadjusted P-value < 0. 05 for genes associated to cnf1-positive lineages (Sup. Table  8 ).
ACKNOWLEDGMENTS This work was supported by the " Fondation ARC" PJA 20191209650 , the  "Fondation pour la Recherche Médicale" (Equipe FRM   2016 ,  DEQ20161136698 ),  Ligue Nationale contre le Cancer Subvention de Recherche Scientifique ,  RS20/75-63  and the  French National Research Agency  ( ANR-10-LABX-62-IBEID ,  INCEPTION ) and  ANR-17-CE17-0014. The plasmid pKOBEG was kindly provided by  Jean-Marc Ghigo.
Acknowledgments We thank all our colleagues at  ENS Paris  and the  University of Chile  for fruitful discussions and useful information on the Chilean subduction. We also thank the comments from  František Gallovič , three other anonymous reviewers and associate editor, and Editor  Rachel Abercrombie  for checking the analysis process and improving the manuscript. Most of the computations are held at the French national computing center  GENCI/ CINES  under the Grants  A0050406700  ( 2018-2019 ),  A0070406700  ( 2019-2020 ), and  A0090406700  ( 2020-2021 ). A part of this work is funded by the  "Institut National des Sciences de l'Univers" (IN-SU-CNRS) , the  "Réseau Sismologique & Géodésique Français" (RESIF) , and the  "Agence Nationale pour la Recherche" (ANR) , under Grants  ANR-11-EQPX-0040  and  ANR-19-CE31-0003.
Data Availability Statement The seismicity catalog is from International Seismological Centre (http://www. isc. ac. uk/ and https://doi. org/10. 31905/d808b825). The interseismic coupling data  (Métois et al. , 2016)  can be obtained via http:// perso. univ-lyon1. fr/marianne. metois/. The used seismological data of Centro Sismológico Nacional of Universidad de Chile can be obtained via http://www. sismologia. cl/. The used continuous GPS data can be obtained via http://gps. csn. uchile. cl/ and the strong motion data from http://evtdb. csn. uchile. cl/. The Global CMT solution  (Dziewonski et al. , 1981; Ekström et al. , 2012)  is obtained via https://www. globalcmt. org/. The forward modeling codes of BIEM  (Aochi et al. , 2000)  and FDM  (Aochi & Madariaga, 2003 ) are available at http://doi.
org/10. 5281/zenodo. 1472238 and http://doi. org/10. 5281/zenodo. 3763864. GAUL is an open source programming library, accessible at http://gaul. sourceforge. net/index. php. Seismic Analysis Code  (Goldstein & Snoke, 2005)  and Genetic Mapping Tools (https://www. generic-mapping-tools. org) are used for signal processing and illustrations (last access October 24, 2019).
Our results show that most of the SS-DMPs (80. 5%) remained significantly associated at our threshold P < 6. 4 × 10 -08 while the reminder showed suggestive associations, of P < 1 × 10 -04 , indicating that therapy applied to the SS patients does influence the DNA methylation patterns associated with SS only in a minority of sites, with modest effects (Supplementary Table  3 ). In fact, a comparison of the effect size of the 118 significant SS-DMPs obtained in these 2 analyses, when including or not the therapy as covariate in the linear regression, showed high correlation (Pearson's correlation R = 0. 99, P = 2. 2 × 10 -6 ). Similar results were also observed in the replication cohort, for which we observed a high replication rate (81. 0%) (Supplementary Table  3 ). Differences in DNA methylation levels and variability of CpGs located in the X-chromosome were also evaluated in females of the discovery cohort, we could not find any significant differences in females at a significance level corrected for multiple testing (data not shown).
The top associations are labeled with gene names. (c) Plots showing high correlation between an average of DNA methylation quantified as β-values at the promoters of the most significant SS-associated DMRs and gene expression at the logarithmic scale. R software  73  and Adobe Illustrator (https:// www. adobe. com/) was used to create figures. ence of anti-Ro/SSA and/or anti-La/SSB autoantibodies previously  36. Stratified analyses based on anti-La/SSB and anti-Ro/SSA autoantibodies positivity show in our data that the epigenetic signature is observed only when positive anti-Ro/SSA patients were compared with healthy subjects (Supplementary Table  3 , Fig. 2a-b ). Indeed, all associations became far from genome-wide significant (P > 0. 006) when the model that includes positive and negative SS patients was adjusted by anti-Ro/SSA presence or when only patients negative for SSA were contrasted with controls (P > 0.
Interestingly, we discovered a large group of genes with differential methylated regions implicated in collagen metabolism such as COL11A2, COL18A1, COL27A1, COL13A1, COL23A1 and COL5A1 that were not previously identified. In addition, many DMRs, especially those that fell in CGI, were located within long non-coding RNAs (LncRNAs), as for example RP11-723C11. 2, RP11-326C3. 7, RP11-480D4. 2 or RP11-89K21. 1 (see Supplementary Table  5 ). We performed functional enrichment analyses for the set of genes showing lower DNA methylation in the SS group (64%) and separately those with increased DNA methylation (34%). On one hand, DMRs with negative signs showed an enrichment in well known pathways implicated in SS such as interferon and cytokine signaling (R-HSA-913531 and R-HSA-1280215, respectively, from Reactome database), as well as in NOD-like receptor signaling pathway (KEGG path:hsa04621), viral carcinogenesis (KEGG path:hsa05203), transcriptional misregulation in cancer (path:hsa05202), necroptosis (KEEF path:hsa04217).
R software  73  and Adobe Illustrator (https:// www. adobe. com/) was used to create figures. Table  1. Genetic associations between HLA variation, presence of Anti-La/SSA autoantibodies and epigenetic IFN signature. β reflects the additive effect of allele dosage for different HLA alleles and P is the associated significance level. The association between HLA genetic variation, SS and SSA was determined by means of logistic regression adjusted by sex and age. The association between HLA genetic variation and epigIFN was determined by linear regression models adjusted by sex, age, cell proportions and batch effects. epigIFN refers to the epigenetic IFN signature. DNA methylation at IFI44L gene (cg13452062) was used as a proxy for epigIFN. Patients exhibiting DNAm > 0. 8 were classified as negative epigIFN.
Moreover, we searched for additional evidence of their functional role in public databases of eQTL and GWAS studies. Importantly, genetic variants regulating DNA methylation at the intergenic region of the CCR cluster in chr3 (Fig. 4a ) and DNA methylation at the downstream region of IFI44 gene (Fig. 4b ) show convincing evidence of their link with other SADS and have been described as eQTL for the same genes (Supplementary Table  16 ). Likewise, genetic variants regulating gene expression of the genes TRIM27 (Fig. 4c ), BTN2A2, UNC119B, GBP5 (Fig. 4d ) and SLFN5 show significant associations with other SADS and a link with gene expression has been previously observed in eQTL studies for these genes (Supplementary Table  16 ). Altogether, these results support the scenario in which differences in DNA methylation and gene expression are intermediates of genetic risk for SS.
eQTL > 0. 05). (g) In SS patients, NUBI expression decreases with the dose of the minor A-allele of rs77466830 (ß SS. eQTL = -0. 15, P SS. eQTL = 3. 5 × 10 -04 ); however, in the healthy population it remains stable (P CTRL. eQTL > 0. 05). (h) In SS patients, PLSCR1 expression increased with the dose of the minor A-allele of rs56077428 (ß SS. eQTL = 0. 357, P SS. eQTL = 0. 0208); however, in the healthy population it remains stable (P CTRL. eQTL > 0. 05). Green boxplots represent SS population, while grey boxplots represent the healthy control population. DNA methylation is quantified with β-values, gene expression is at the logarithmic scale. R software  73  and Adobe Illustrator (https:// www. adobe. com/) was used to create figures. could generate antibodies that cross-react with autoantigens, such as Anti-Ro/SSA and Anti-La/SSB, by molecular mimicry  46.
We also could not evaluate specific epigenetic changes related to the development of some important disease phenotypes for SS, as for example lymphomas. Defective DNA methylation in SS patients with lymphomas and other cancers has been previously reported  69  and, recently a relationship between lymphoma and the DNA-methylation based interferon signature has been described  45. However we did not have enough samples in our dataset to assess the epigenome-wide contribution in SS. Future studies that target specific disease phenotypes, that incorporate higher sample sizes, and add other cell-specific layers of molecular information will add more to the understanding of the genetic regulation and transcriptional consequences of epigenetic variation associated with SS and other systemic autoimmune diseases.
Material and methods Participant recruitment. Samples included in this study were recruited from the PRECISESADS study, which is a European multi center, non-randomized, and observational clinical study with recruitment performed between December 2014 and December 2018 at 19 institutions in 9 countries. PRECISESADS included patients affected by systemic autoimmune (around 400 by disease or group of diseases: rheumatoid arthritis, scleroderma or systemic sclerosis, primary Sjögren's syndrome, systemic lupus erythematosus, and primary antiphospholipid syndrome, mixed connective tissue disease and undifferentiated connective tissue disease) and 554 healthy controls  70, 71. A total of 558 samples were included in the core analyses of the present study (Supplementary Table  1 ), that included SS patients and healthy controls for which DNA methylation and genotype data was available.
The protection of the confidentiality of records that could identify the included individuals is ensured as defined by the EU Directive 2001/20/EC and the applicable national and international requirements relating to data protection in each participating country. Genome-wide DNA methylation data and differential analyses. DNA methylation data of the discovery cohort was obtained using the Infinium Methylation EPIC BeadChip (Illumina, San Diego, CA, USA) that covers more than 800,000 CpG sites. For replication cohort, we had also available genome-wide DNA meth-ylation data profiled obtained by Infinium Methylation 450 K BeadChip (Illumina, San Diego, CA, USA), which covers more than 400,000 CpG sites, most of them included in the EPIC array. After DNA extraction from whole peripheral blood and bisulfite conversion, the genome for each sample was amplified, fragmented and hybridized to the corresponding Illumina arrays according to the manufacturer's protocol.
The quality control of samples and the normalization of the data was performed using the meffil R software  73, 74  Samples were excluded based on the detection P criteria > 99%, poor bisulfite conversion based on control dashboard check, and sex mismatches according to failed chromosome X and Y clustering. Probes were filtered out based on detection P > 0. 01 in > 95% of samples. Additionally, probes located at the X and Y chromosomes were separated in different datasets to avoid gender bias. Probes with genetic variants at their CpG sites were also excluded. After applying these filtering steps we obtained 776,284 and 433,337 autosomic probes in the discovery dataset and in the replication dataset, respectively. A total of 17,530 probes located in the X chromosome passed the QC in the discovery cohort and 9282 in the replication cohort.
A DMRs was considered as significant with a threshold FDR < 0. 05 and with a minimum number of 5 CpGs associated to the region. We used the mCSEA R package  52  to perform expression quantitative trait methylation (eQTM) analysis to discover significant associations between these DMRs and an expression alteration in the closer genes. In this case, the leading edge CpGs of each region is first defined and averaged for each region in each sample. Then, Pearson's correlation coefficient is calculated between each region's methylation and the proximal genes expression within 1500 base pairs upstream and downstream from the region.
RNA-seq data and differential expression analysis. We analysed gene expression data available from 179 SS patients and 247 healthy subjects. RNA-seq data was obtained as 50 bp single end reads using an Illumina HiSeq 2500 sequencer and its proprietary base caller. Fastq data was aligned against human genome GRCh19 using STAR aligner  76  with GENCODE v19 [URL: https:// www. genco degen es. org/ relea ses/ 19. html] annotation as single pass alignment. Read counts were obtained using RSEM  77. Samples were included only if they had more than 8 million reads, a RIN score > 8 and did not appear as outliers in a principal component analysis. After filtering, a total of 42,878 genes were included in the study. Differential expression analysis was performed using a linear model that regressed out SS status on readcounts as gene expression levels adjusted by age, sex, the first genetic principal component, the cell proportions and sequencing pool (POOL) and RNA integrity number (RIN) as batch effects.
A Bonferroni-adjusted P-value threshold < 1. 16 × 10 -6 was established as significant to identify the differentially expressed genes (DEGs) in a discovery cohort formed by 135 SS patients and 174 healthy controls. In the replication cohort (44 SS patients Genotype profiling and imputation methods. Genomic DNA from whole blood was obtained by standard methods. The samples were genotyped using HumanCore -12-v1-0-B, InfiniumCoreExome-24v1-2 and InfiniumCoreExome-24v1-3, all of them of Illumina (San Diego, CA, USA). Only genetic variants, not indels, present on all three platforms were considered for data cleaning and analysis. Quality controls (QC) were performed using PLINK v1. 9  79. Genetic markers were removed if these had a call rate < 90%, exhibited a significant differential missingness between cases and controls (P < 1 × 10 -4 ), and showed a significant deviation from Hardy-Weinberg equilibrium (P < 0.
inter natio nalge nome. org/ categ ory/ phase-3/] were included as reference panel. PCAs were also calculated using SMARTPCA from the Eigensoft software  84  and the same independent set of markers. In this analysis, only the European population from 1000 Genomes, excluding Finns, was included as reference population of our samples  85  (see supplementary Fig. 1 ). Six standard deviations from the centroid were used as a threshold to filter individuals that deviate from the main European clustering. Moreover, samples with less that 55% European ancestry were excluded from all analyses. After genotyping QC and filtering of individuals for European ancestry classical HLA alleles were imputed for the PRECISESADS dataset (cases and controls) in the extended MHC region in chromosome 6  86. The SNP2HLA software  87  was used for imputation using a reference panel consisting of 5,225 European individuals in the Type 1 Diabetes Genetic Consortium 88 containing data of 8,961 variants across the MHC region, and two and four digit-resolution allelic identities of the HLA class I (HLA-A, HLA-B, and HLA-C) and II genes (HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQB1, and HLA-DRB1).
Genotypes were extracted as allele dosages for 298 HLA classical alleles. Quantitative trait loci analysis and genetic associations. Genetic regulation of DNA methylation and gene expression was explored by the quantitative trait locus (QTL) approach using Matrix eQTL  89. A linear regression model was performed to evaluate the effect of closely located SNPs (no further than 1 Mb to CpG site or TSS) with a MAF > 0. 05 on DNA methylation levels of DMPs and on gene expression patterns of DEGs. In this case, the gene expression data were firstly normalized into variance stabilizing transformations (VST) using DESeq2  78. The linear model was corrected by age, sex, disease status, the first genetic principal component, cell proportions, and the respective batch effect. A P-value threshold of false discovery rate (FDR) < 0.