Antibiotic Susceptibility and General Genomic Features of the 112 H. Pylori-Shi Isolates
A total of 112 H. pylori isolates from patients with digestive diseases were categorized by antibiotic susceptibility phenotypes. Patients’ mean age was 52.2 ± 14.8 (mean ± SD) years, and 54.5% (61/112) were male (Table 1). The prevalence of resistance to MTZ was the highest (65.2%), with 33.9% MTZ mono-resistant, followed by LEV (34.8%), with 4.5% LEV mono-resistant, CLA (16.1%) with 2.7% CLA mono-resistant, and no isolate was resistant to AMX. For dual-resistant isolates, the majority (20.5%) were resistant both to MTZ and LEV. A total of eight isolates were MDR isolates that were simultaneously resistant to MTZ, CLA, and LEV, and 28 were susceptible to all four antibiotics (Table 1). Gender, age, and gastroscopy diagnosis were not significantly associated with resistance (P > 0.05).
Table 1
Demographic information and antibiotic resistance rates.
Group
|
Total, n (%)
|
Antibiotic-Resistance Phenotypes, n (%)
|
Amoxicillin-R/mono-R
|
Metronidazole-R/mono-R
|
Clarithromycin-R/mono-R
|
Levofloxacin-R/ mono-R
|
Metronidazole-Clarithromycin dual-R
|
Metronidazole-Levofloxacin dual-R
|
Clarithromycin-Levofloxacin dual-R
|
MDR
|
Susceptible to four antibiotics
|
Gender
|
male
|
61 (54.5)
|
0 (0)/0 (0)
|
40 (65.6)/23 (37.7)
|
8 (13.1)/3 (4.9)
|
19 (31.1)/4 (6.6)
|
2 (3.3)
|
12 (19.7)
|
0 (0)
|
3 (4.9)
|
14 (23.0)
|
female
|
51 (45.5)
|
0 (0)/0 (0)
|
33 (64.7)/15 (29.4)
|
10 (19.6)/0 (0)
|
20 (39.2)/1 (2.0)
|
2 (3.9)
|
11 (21.6)
|
3 (5.9)
|
5 (9.8)
|
14 (23.0)
|
Age
|
19-29
|
8 (7.1)
|
0 (0)/0 (0)
|
5 (62.5)/2 (25.0)
|
2 (25.0)/1 (12.5)
|
3 (37.5)/0 (0)
|
0 (0)
|
2 (25.0)
|
0 (0)
|
1 (12.5)
|
2 (25.0)
|
30-39
|
18 (16.1)
|
0 (0)/0 (0)
|
8 (44.4)/7 (38.9)
|
3 (16.7)/2 (11.1)
|
3 (16.7)/3 (16.7)
|
1 (5.6)
|
0 (0)
|
0 (0)
|
0 (0)
|
5 (27.8)
|
40-49
|
17 (15.2)
|
0 (0)/0 (0)
|
12 (70.6)/6 (35.3)
|
4 (23.5)/0 (0)
|
6 (35.3)/0 (0)
|
1 (5.9)
|
3 (17.6)
|
1 (5.9)
|
2 (11.8)
|
4 (23.5)
|
50-59
|
29 (25.9)
|
0 (0)/0 (0)
|
20 (69.0)/10 (34.5)
|
4 (13.8)/0 (0)
|
9 (31.0)/0 (0)
|
2 (6.9)
|
7 (24.1)
|
1 (3.4)
|
1 (3.4)
|
8 (27.6)
|
60-69
|
27 (24.1)
|
0 (0)/0 (0)
|
20 (74.1)/11 (40.7)
|
2 (7.41)/0 (0)
|
11 (40.7)/1 (3.7)
|
0 (0)
|
8 (29.6)
|
1 (3.7)
|
1 (3.7)
|
5 (18.5)
|
70-79
|
13 (11.6)
|
0 (0)/0 (0)
|
6 (46.2)/2 (15.4)
|
2 (15.4)/0 (0)
|
6 (46.2)/1 (7.7)
|
0 (0)
|
2 (15.4)
|
0 (0)
|
2 (15.4)
|
4 (30.8)
|
80<
|
2 (1.8)
|
0 (0)/0 (0)
|
2 (100)/0(0)
|
1 (50)/0(0)
|
1 (50)/0(0)
|
0 (0)
|
1 (50)
|
0 (0)
|
1 (50)
|
0 (0)
|
Gastroscopy diagnosis
|
Chronic superficial gastritis
|
19 (17.0)
|
0 (0)/0 (0)
|
8(42.1)/5(26.3)
|
2(10.5)/1(5.3)
|
5(26.3)/2(10.5)
|
0(0)
|
2(10.5)
|
0(0)
|
1(5.3)
|
8(42.1)
|
Chronic atrophic gastritis
|
66 (58.9)
|
0 (0)/0 (0)
|
46(69.7)/20(30.3)
|
13(19.7)/0(0)
|
29(43.9)/3(4.5)
|
3(4.5)
|
16(24.2)
|
3(4.5)
|
7(10.6)
|
14(21.2)
|
Gastric/Duodenal ulcer
|
26 (23.2)
|
0 (0)/0 (0)
|
18(69.2)/12(46.2)
|
3(11.5)/2(7.7)
|
5(19.2)/0(0)
|
1(3.8)
|
5(19.2)
|
0 (0)
|
0 (0)
|
6(23.1)
|
Gastric cancer
|
1 (0.9)
|
0 (0)/0 (0)
|
1(100)/1(100)
|
0 (0)/0 (0)
|
0 (0)/0 (0)
|
0 (0)/0 (0)
|
0 (0)/0 (0)
|
0 (0)/0 (0)
|
0 (0)/0 (0)
|
0 (0)/0 (0)
|
Total
|
112
|
0 (0)/0 (0)
|
73 (65.2)/38 (33.9)
|
18 (16.1)/3 (2.7)
|
39 (34.8)/5 (4.5)
|
4 (3.6)
|
23 (20.5)
|
3 (2.7)
|
8 (7.1)
|
28 (25.0)
|
The overall sequence read length, contigs, scaffolds, ORFs, and guanine-cytosine (GC) content for the genomes of the 112 H. pylori isolates grouped by antibiotic resistance phenotypes are shown in Supplementary Table S2. On average, the read length of the 112 H. pylori isolates’ genomes was 1.60 billion bp long and encoded a range of 1,511 to 1,624 genes per genome. A low average guanine-cytosine content, 38.7%, was identified, and the chromosome sizes ranged from 1.52 to 1.69 Mb. All sequenced genomes harbored 36 tRNA genes and an average of two rRNA genes.
Genome-wide Gene Family Analysis and Resistance-Unique Genes of H. pylori Isolates
We initially conducted gene family analysis to define the core genetic content and unique-gene pool of specific antibiotic resistance. We obtained 2,439 genes from all the genomes, and among these genes, 2,194 genes were found shared between MTZ-S and -R, 2,101 genes were shared between CLA-S and -R, and 2,193 genes were shared between LEV-S and -R, accounting for 89.95%, 86.14% and 89.91% of the total genes, respectively, which constituted the core genetic content of the corresponding resistant and susceptible categories (Figure 1A-C). Moreover, we obtained 183, 9, and 23 unique genes in MTZ-, CLA-, and LEV-resistant categories, respectively, containing 75, 5, and 13 newly found genes which were not included in neither NCBI-NR, SwissProt, nor COG database (Figure 1D-F, Supplementary Table S9). Among the remaining 108, 4 and 11 unique genes defined as resistance-unique-gene pools of the three antibiotics resistance categories (Supplementary Table S3), only 31, 2, and 3 genes could be assigned to various COG functional classes, respectively (Figure 1D-F).
Distributions of the OMP Family Genes in Resistant and Susceptible Categories
To investigate whether the occurrence of genes of the OMP families was related to specific antibiotic resistance, the presence and absence of genes predicted to encode OMPs in different categories were analyzed, including the currently identified five paralogous gene families major OMP family including Hop subfamily and Hop-related (Hor) subfamily, Hof family, Hom family, iron-regulated OMPs, and efflux pump OMPs, ranging in size from 3 to 36 members, as well as 9 other putative OMPs not in paralogous families (Supplementary Table S4). We first compared the average levels of OMP genes in the resistant and the susceptible categories irrespective of the number of alleles of the genes. There were no significant differences in the average number of total OMP genes, the major OMP family genes, Hof and Hom family genes, and efflux pump OMP genes between isolates resistant and susceptible to the three antibiotics, respectively (Figure 2A-F). We intriguingly found that the average frequency of the iron-regulated OMP genes (fecA and frpB) in the CLA-resistance category was significantly lower than that in CLA-susceptible category (P = 0.0103). Notably, the frequencies of these genes in the isolates varied greatly. Nearly all the isolates contained the majority of the Hop subfamily genes (not including hopE and hopN), the efflux pump OMP genes, horB, horG, horD, hofA belonging to the hor or hof subfamily, and the gene coding BamA and murein lipoprotein 20 (Lpp20) (Figure 2G). However, horA, horC and hofG genes were present in very few strains. Remarkably, difference analysis showed the proportion of isolates carrying hopE in both MTZ- and LEV-susceptible categories was significantly higher than that in corresponding resistant category (Figure 2G). Another gene with significantly lower frequency in resistant categories were hofF (LEV).
Distributions of the Virulence-Associated Genes in Resistant and Susceptible Categories
To investigate whether the occurrence of virulence factor genes H. pylori carried was related to specific antibiotic resistance, the presence and absence of 108 predicted virulence-associated genes against the VFDB database in resistant and susceptible categories was analyzed, including 26 cag pathogenicity island (cagPAI) genes, 36 flagellar genes, and 47 other virulence-associated genes (Supplemental Table S5). Similarly, we compared the average levels of virulence-associated genes in the resistant and susceptible categories irrespective of the number of alleles of the genes. There were no significant differences in the number of total virulence-associated genes, cagPAI genes, and flagellar genes between resistant and the susceptible categories of the three antibiotics, respectively (Figure 3A-C). The presence frequency of these genes in the resistant and susceptible categories were compared, and there was no significant difference in the distribution of most of the genes (Supplemental Table S5). Dramatically, the frequencies of cagY (HP0527) and pflA (HP1274) were significantly higher (P = 0.0240) and lower (P = 0.0488) in MTZ-R category than that in MTZ-S, respectively (Figure 3D-E).
Distributions of the Efflux Pump Genes in Resistant and Susceptible Categories
To investigate whether the occurrence of efflux pump genes was related to specific antibiotic resistance, the presence and absence of 18 literature-based efflux pump genes were analyzed, including 12 RND family genes, 3 ABC family genes, and 3 MFS family genes (Figure 4). Among the 18 genes, 16 were present in almost all the 112 H. pylori isolates, and the remaining two genes, spaB (HP0600) and HP1181 were widely deleted in these isolates, detected in 25 and 31 isolates, respectively. Nevertheless, the proportion of isolates carrying spaB in the MTZ-R (P = 0.0002) and LEV-R (P = 0.0408) categories was obviously higher than that in the MTZ-S and LEV-S categories, respectively (Figure 4).
Correlation Analysis of CRISPR and Its Components With H. pylori Antibiotic Resistance
Given the natural function of the CRISPR-Cas system as an adaptive defense mechanism against foreign DNA and a controversial effect of the CRISPR-Cas system on the transfer of the acquired antibiotic resistance genes (ARGs) in bacterial genomes, we conducted predictions in the CRISPRs-Cas systems and acquired ARGs where resistance is conferred by a complete gene rather than mutations, seeking to explore the relationship among CRISPR-Cas, acquired ARGs, and antibiotic resistance in H. pylori. Unfortunately, the results suggested neither ARGs nor cas genes were found in 112 H. pylori isolates. However, a total of 153 CRISPR loci were predicted in 83 isolates, ranging in length from 65 to 496 bp.
Because of the different CRISPR contents of the isolates, the proportions of isolates with 0, 1, 2, 3, and 4 CRISPRs in resistant and susceptible categories of MTZ, CLA, and LEV were analyzed, respectively (Figure 5A-C). Notably, the proportion of the isolates with no CRISPR in the LEV-R category was significantly higher than that in the LEV-S category (P = 0.0075, Figure 5C). Because of the different spacer contents of the CRISPRs, the proportion of CRISPRs with 1, 2, 3, 4, and 6 spacers in resistant and susceptible categories was also analyzed (Figure 5D-F). Interestingly, the proportion of the CRISPRs with three spacers in the MTZ-R category was significantly higher than that in the MTZ-S category (P = 0.0331, Figure 5D), and the proportion of CRISPRs with two spacers in the LEV-R category was significantly lower than that in the LEV-S category (P = 0.0152, Figure 5F). Similarly, the proportion of the CRISPRs with two or three spacers in the resistance category of the other two antibiotics was also higher or lower than that of the corresponding susceptible category, respectively, although there was no significant difference (Figure 5D-F). In addition, detailed analysis of the DR and spacer sequences showed that a CRISPR combination containing three most common CRISPRs in H. pylori coexisted in seven isolates, which intriguingly were all MTZ mono-R isolates (Figure 5G). Of note, one of the DR in the CRISPR combination occurred in CRISPRs of 9 isolates, which included CRISPRs in seven MTZ mono-R isolates and 2 CRISPRs (Supplementary Figure S1) in two other MTZ-R isolates (a MTZ mono-R isolate and a MTZ and CLA dual-R isolate), causing the significantly higher proportion of the isolates with the DR in the MTZ-R category compared with that in the MTZ-S category (P = 0.0259, Figure 5G, Supplementary Table S10).
Phylogenetic Analysis and Relationship of Antibiotic Resistance Patterns and specific genetic loci with strain relatedness in H. pylori Isolates.
When constructing the maximum likelihood tree of the 112 H. pylori-Shi genomes with respect to different resistance patterns and the profile of the observed differently present subsystem genomic loci between susceptible and resistant categories for the corresponding antibiotic, no strain-relatedness was found to be associated to neither each single antibiotic phenotype nor the presence/absence of the differently present loci including genes coding OMPs (hopE, hofF and iron-regulated OMPs), virulence-associated genes (virB10/cagY and pflA), efflux pump genes (spaB), and CRISPRs with 2 or 3 spacers or the indicated DR sequences (Figure 6). Alternatively, these strains had distinct susceptibility and indicated loci profiles regardless of genetic relatedness among strains.
Genome-Wide Genetic Variations and Prevalence of Variations in H. pylori Resistome
We culled SNPs and Indels within the coding sequences and noncoding rRNAs and tRNAs from the genomes of all the H. pylori isolates when compared with H. pylori 26695 reference genome (Figure 7A), which showed high densities of the variations that could be as high as 136 Indel density and 75 SNP density corresponding to reference genome (Figure 7A). Among these genetic variations, a total of 388 non-synonymous SNPs (nsSNPs) and 1,718 frameshift Indels (fsIndels) were identified (Supplementary Table S8).
We obtained a list of 80 genes defined as the H. pylori resistome, whose variation, presence, or absence were known to either confer H. pylori antibiotic resistance or subsequently found to be associated with H. pylori resistance to clinically used antibiotics not limited to antibiotics with known phenotypes in this study. To investigate the prevalence of the genetic variations in the H. pylori resistome, we first assayed the variations in six most common resistance genes, including 23S rRNA, gyrA and gyrB genes, and rdxA, frxA, and fdxB genes, in which single-point mutations in specific regions are known to be responsible for CLA, LEV, and MTZ resistance, respectively (Table 2-4). Of these SNPs, a subset was significantly associated with antibiotic resistance, including A2143G, A2142G and C148T in the 23S rRNA genegene, and N87T, N87I, and D91G, D91Y in gyrA, as well as I38V in fdxB. Additionally, the assay detected four Indels in 23S rRNA, three Indels in gyrA, and one in frxA. Unexpectedly, the prevalence of one of the Indels (762_763insT) in 23S rRNA was significantly higher in the CLA-S category (53.19%) compared with that in the CLA-R category (Table 2). Furthermore, we overviewed the distribution of the numbers of the SNPs and Indels present in the remaining genes of the H. pylori resistome of the 112 isolates along with the antibiotic resistant profile of these genes (Supplementary Figure S2). Of the remaining 74 genes of the H. pylori resistome, 174 variations were obtained corresponding to 54 genes, involving 29 genes detected containing nsSNPs, 36 containing fsIndels, and 10 containing both nsSNPs and fsIndels (Supplementary Table S7). Additionally, the maximum number of Indels (ranging from 0 to 15) contained in these genes was up to three times the maximum number of SNPs (ranging from 0 to 3), suggesting to a certain extent that Indel mutations make greater contributions to the high genetic diversity of the H. pylori resistome than SNP mutations do (Supplementary Figure S2).
Table 2
nsSNPs and fsIndels in genes conferring resistance to CLA.
Gene
|
Nucleotide change
|
No. among CLA-R (n=18)
|
No. among CLA-S (n=94)
|
P
|
23S rRNA
|
SNP
|
|
|
|
|
A2143G
|
17 (94.44%)
|
2 (2.13%)
|
<0.0001
|
|
A2142G
|
2 (11.11%)
|
0 (0.00%)
|
0.0250
|
|
T2182C
|
6 (33.33%)
|
31 (32.98%)
|
0.9770
|
|
C148T
|
2 (11.11%)
|
0 (0.00%)
|
0.0250
|
|
T1025C
|
0 (0.00%)
|
1 (1.06%)
|
1.0000
|
|
G1027A
|
14 (77.78%)
|
52 (55.32%)
|
0.0760
|
|
Indel
|
|
|
|
|
761delAA (GAA to G)
|
16 (88.89%)
|
84 (89.36%)
|
1.0000
|
|
762_763insT (A to AT)
|
5 (27.78%)
|
50 (53.19%)
|
0.0482
|
|
764_765insTTT (C to CTTT)
|
8 (44.44%)
|
54 (57.45%)
|
0.3093
|
|
1036_1037insC (A to AC)
|
18 (100.00%)
|
94 (100.00%)
|
1.0000
|
|
1515_1516insT (C to CT)
|
3 (16.67%)
|
9 (9.57%)
|
0.4054
|
|
No change
|
0 (0.00%)
|
13 (13.83%)
|
0.1232
|
Table 3
nsSNPs and fsIndels in genes conferring resistance to LEV.
Gene
|
Sequence change
|
No. among LEV-R (n=39)
|
No. among LEV-S (n=73)
|
P
|
Nucleotide position
|
Amino acid position
|
gyrA
|
SNP
|
|
|
|
|
|
C261A
|
N87T
|
8 (20.51%)
|
0 (0.00%)
|
0.0001
|
|
C261G
|
N87I
|
9 (23.08%)
|
0 (0.00%)
|
< 0.0001
|
|
G271A
|
D91G
|
12 (30.77%)
|
0 (0.00%)
|
< 0.0001
|
|
G271T
|
D91Y
|
5 (1.67%)
|
0 (0.00%)
|
0.0043
|
|
Indel
|
|
|
|
|
|
2455_2456insAT (AC to AATC)
|
T819fs
|
32 (82.05%)
|
57 (78.08%)
|
0.6203
|
|
2460_2463del (CTTCG to C)
|
S820fs
|
34 (87.18%)
|
62 (84.93%)
|
0.7460
|
|
2472delT (AT to A)
|
N824fs
|
38 (97.44%)
|
66 (90.41%)
|
0.2577
|
|
No change
|
No change
|
1 (2.56%)
|
7 (9.59%)
|
0.2577
|
gyrB
|
SNP
|
|
|
|
|
|
C1801T
|
L601F
|
37 (94.87%)
|
73 (100.00%)
|
0.3423
|
|
No change
|
No change
|
2 (5.13%)
|
0 (0.00%)
|
0.1192
|
Table 4
nsSNPs and fsIndels in genes conferring resistance to MTZ.
Gene
|
Sequence change
|
No. among MTZ-R (n=73)
|
No. among MTZ-S (n=39)
|
P
|
Nucleotide position
|
Amino acid position
|
rdxA
|
No change
|
No change
|
73 (100.00%)
|
39 (100.00%)
|
NA
|
frxA
|
SNP
|
|
|
|
|
|
G499A
|
D167N
|
4 (5.48%)
|
4 (10.26%)
|
0.3497
|
|
T577A
|
C193S
|
63 (86.30%)
|
39 (100.00%)
|
0.3238
|
|
Indel
|
|
|
|
|
|
653dupA (T to TA)
|
X218delinsX
|
70 (95.89%)
|
39 (100.00%)
|
0.5503
|
|
No change
|
No change
|
3 (4.11%)
|
0 (0.00%)
|
0.5503
|
fdxB
|
SNP
|
|
|
|
|
|
A112G
|
I38V
|
65 (89.04%)
|
39 (100.00%)
|
0.0488
|
|
No change
|
No change
|
8 (10.96%)
|
0 (0.00%)
|
0.0488
|
ins: insert, del: delete, fs: frameshift, NA: not available |
Resistome-Wide Association Study (RWAS) of Genetic Variations
In order to explore MTZ, CLA, and LEV resistance-related or resistance-induced compensation related potential genetic variants, we compared the prevalence of the resistome-wide nsSNPs and fsIndels between resistant and susceptible isolates. We observed 23 nsSNPs and 48 fsIndels corresponding to 31 genes that were significantly enriched in antibiotic- resistant or -susceptible isolates (P < 0.05) We found that four variations within the integrase/recombinase gene xerD showed a significant association with both MTZ and CLA resistance. Some of these variations had strong linkage disequilibrium (LD). Both nsSNPs (xerD:D353N** and xerD:S275N*) were linked together but in low LD with other two xerD fsIndels. Similarly, the two fsIndels (xerD:I40/120_121 delete GG and xerD:L39/116_117 insert G) were also linked together (Figure 7B). A deletion within the site-specific DNA-methyltransferase gene hsdM (HP0850) (hsdM:P333/998_1002 delete AGCCG) was found to be significantly associated with both MTZ and LEV resistance, which also were in high LD with gyrA:N87I. Furthermore, we observed 16 variations within 15 genes simultaneously showing association with MTZ susceptibility and LEV resistance, and many of these genes code membrane-related proteins (hetA, msbA, hofC, oppA, tetA(p), and babB) or proteins involved in recombination (recG and ruvC) (Figure 7B). Moreover, several deletions within the ATP-dependent serine protease gene lon (especially the deletions at lon:T582, lon:K584, and lon:E588) showing a significant association with MTZ susceptibility were found with high LD with more than half of the other variations, including the resistance-conferring nsSNPs 23S rRNA:A2143G, 23S rRNA:A2142G, 23S rRNA:C148T, gyrA:N87T, and some of the remaining variants, most of which were associated with MTZ resistance or susceptibility (Figure 7B) (Supplementary Table S11).
Resistance- or Susceptibility-Associated Variations Occurred in Functionally Important Domains of Genes Within the H. pylori Resistome
For the purpose of relating the resistance- or susceptibility-enriched variations to specific protein domains, we further characterized these variations in the genes within the H. pylori resistome with corresponding protein structure data. Localization of the nsSNPs and fsIndels in the four proteins (Lon, BabB, XerD, and TrpS) with the largest number of phenotype-associated variations were summarized (Figure 8), and the functional domains of all the genes with the phenotype -associated variations within the H. pylori resistome were identified. Eight variations in Lon and three partial in TrpS closely located variations were identified although they were not in any domains (Figure 8A and D, Supplementary Figure S4). In BabB, two clusters of locationally close variation sites at L59, N60, G129, S130, and N135 were all in the sialic acid binding adhesin (SabA) N-terminal extracellular adhesion domain (Figure 8B, Supplementary Figure S4). Among four variations in XerD, only one nsSNP S275N was located in the phage integrase family domain that enables the protein to covalently link to DNA. Of the remaining 27 genes, 14 contained variations totally or partially located in their functional domains, including three genes encoding efflux pump transporters with an S267P mutation in the extracellular solute-binding domain of dppA, an R390K mutation in the ABC transporter domain of hetA, and an L304F mutation in the ABC transporter transmembrane region of msbA; three genes required for lipopolysaccharide (LPS) biosynthesis with a V76A mutation in the mur ligase family catalytic domain of murC, a deletion at A350 and an insertion at I351 in the mur ligase family glutamate ligase domain of murE, and a deletion at P129 in the region with glycosyltransferase activity of rfaF; two genes related to protein translation with an S32A mutation in the elongation factor Tu GTP binding domain of lepA and a deletion at A297 in the tRNA synthetase domain of iles, as well as the metalloregulator gene fur with an insertion at S148 and a deletion at E149 in the ferric uptake regulator region; the OMP hofC with an I206V mutation and a deletion at M276 in the major functional domain; the stress-response protein coding gene grpE with H145N mutation in its main functional domain; the chromosomal replication initiator gene dnaA with A200T mutation in the whole functional region; and the restriction enzyme alleles hsdM_2 and hsdM_3 with a deletion at P333 in the domain the N-6 DNA methylase activity and a deletion at N49 in the N-terminal domain, respectively (Supplementary Figure S4).