Sungeun Lee

Ella T Sieradzki

Christina Hazard
email: christina.hazard@ec-lyon.fr

Graeme W Nicol
email: graeme.nicol@ec-lyon.fr.

Viruses of soil ammonia oxidising archaea identified using a novel DNA stable isotope probing approach for low GC mol% genomes

Ammonia oxidising archaea (AOA) are a ubiquitous component of microbial communities which can dominate ammonia oxidation in some soils. While we are beginning to understand soil virus dynamics, we have no knowledge of the composition or activity of those infecting AOA or their potential to influence processes.

This study aimed to identify viruses infecting autotrophic AOA in two soils (pH 4.5 and 7.5) by following transfer of assimilated CO2-derived 13 C via DNA stable isotope probing in CsCl isopycnic buoyant density (BD) gradients and metagenomic analysis.

Incorporation of 13 C into low GC mol% AOA genomes (~35-40%) increased DNA BD and separation from unenriched DNA, but also resulted in co-migration with dominant non-13 C-enriched high GC mol% bacterial genomes (>57%) reducing AOA sequencing depth and contig assembly. We therefore developed a hybrid approach where AOA and virus genomes were assembled from metagenomes from low BD CsCl fractions from 12 C-CO2 incubations with subsequent mapping of high BD 13 C-vs 12 C-enriched non-assembled reads to identify activity. AOA metagenome assembled genomes (MAGs) were distinct between the two soils and represented a broad diversity of active lineages. Sixty-four viral OTUs infecting AOA were also distinct between soils, with 42% enriched in 13 C. Comparative analyses demonstrated that soil AOA viruses were distinct from characterised bacterial or archaeal viruses and auxiliary metagenomic viral genes included an AOA-specific multicopper oxidase involved in copper uptake essential for AOA metabolism. These observations in soils with distinct AOA communities indicate that virus infection of AOA is likely a frequent process during soil nitrification.

Main text

Autotrophic ammonia-oxidising archaea (AOA) of the class Nitrososphaeria are a ubiquitous component of soil microbial communities and often dominate the first step of nitrification and nitrification-associated N2O emissions when ammonia is supplied at low rates via organic matter mineralisation [1], slow-release fertilisers [2] or in acidic soils [3]. Prophages and other virus-associated protein-encoding genes are abundant in AOA genomes suggesting frequent interaction. While viruses infecting marine AOA have been characterised through metagenomic approaches [4] and cultivation [5],

those infecting soil AOA or other nitrifier groups are currently uncharacterised.

Recent advances have demonstrated that soil virus communities are dynamic in a wide range of soils [e.g. 6,7] and augmenting virus loads modulate C and N fluxes [START_REF] Li | MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph[END_REF][START_REF] Kang | MetaBAT 2: an adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies[END_REF]. Nevertheless, identifying active interactions with specific populations or functional groups remains challenging due to the vast diversity of both host and viruses in soil. Recent use of stable-isotope approaches has investigated whole community host-virus dynamics [START_REF] Wu | MaxBin 2.0: an automated binning algorithm to recover genomes from multiple metagenomic datasets[END_REF][START_REF] Alneberg | Binning metagenomic contigs by coverage and composition[END_REF] or interactions between individual host-virus populations specific to a functional process and substrate [START_REF] Parks | CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes[END_REF]. The aim of this study was to utilise the latter approach with 13 CO2-based DNA-SIP to focus on nitrification-associated interactions and to test the hypothesis that viruses are a dynamic component of soil AOA activity.

Hybrid analysis of GC mol% fractionation and 13 C-DNA-SIP for identifying active AOA hosts and viruses

Microcosms were established (see Supplementary Text) with two sandy loam soils of contrasting pH (4.5 or 7.5) which were characterised previously as containing distinct AOA communities [START_REF] Menzel | Fast and sensitive taxonomic classification for metagenomics with Kaiju[END_REF]. Nitrification was stimulated by addition of urea (two pulses of 100 µg N g -1 soil (dry weight)) and incubated at 25 o C for 30 days (Supplementary Figure 1) with a 5% 12 C-or 13 C-CO2 headspace to target viruses of autotrophic nitrifiers via transfer of assimilated carbon. DNA was extracted and subject to isopycnic centrifugation in CsCl gradients, in which DNA migrates as a function of GC mol% and isotopic enrichment, before fractionation and analysis of genomic DNA of different buoyant densities [START_REF] Parks | CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes[END_REF]. Quantification of AOA genomic DNA distribution across CsCl gradients demonstrated growth of AOA populations (Figure 1A; Supplementary Figure 2B) and genomic DNA in high buoyant density (HBD) fractions (>1.719 g ml -1 ) was then pooled for each replicate microcosm and metagenomes sequenced (Supplementary Table 1). While 13 C-incorporation clearly separated DNA from growing vs non-growing AOA populations, 13 C-enriched genomes (average AOA contig GC mol% content 38%) co-migrated within the CsCl gradient with community 12 C genomes with higher GC mol% (>57%), reducing AOA host and virus sequencing depth, contig assembly and binning. Consequently, only one medium quality AOA MAG and no AOA-associated virus contigs ³10 kb were identified (Figure 1B).

To facilitate analysis of low GC mol% genomes while also identifying activity, we developed a hybrid approach combining metagenomic analysis of low buoyant density (LBD) 12 C DNA with read-mapping of sequenced HBD DNA (Figure 1; Supplementary Figure 2A). AOA and virus genomes were assembled from metagenomic libraries of 12 C LBD DNA only (<1.699 g ml -1 ) representing both dormant and active populations. Individual sequence reads from HBD DNA libraries from both 12 C-and 13 C-incubations were then mapped to identify activity. In comparison to HBD, assemblies from LBD libraries increased the relative proportion of contigs from the AOA-containing Thaumarchaeota phylum (NCBI taxonomy (Thermoproteota GTDB taxonomy)) and Thaumarchaeota-infecting virus contigs from 0.3 to 2.8% and 0 to 22.3%, respectively (Supplementary Figure 3). In LBD metagenomes, 123 mediumand high-quality MAGs from both soils were identified (Supplementary Table 2) including 9 AOA MAGs representing 6 genera (GTDB classification [START_REF] Chaumeil | GTDB-Tk: a toolkit to classify genomes with the Genome Taxonomy Database[END_REF]) within the Nitrososphaerales and Nitrosopumilales orders and one nitrite oxidising bacterium (NOB) Nitrobacter MAG (MAG 15) (Figure 1B, Supplementary Figure 4). Mapping of HBD sequence reads revealed that 8 of the 10 MAGs were derived from active autotrophic populations (Figure 1C). Consistent with previous studies characterising the distribution of soil AOA from different taxa [START_REF] Roux | VirSorter: mining viral signal from microbial genomic data[END_REF], this included 3 Nitrosotalea genus MAGs from pH 4.5 soil and 4 Nitrososphaera genus MAGs from pH 7.5 soil. No evidence of autotrophic growth was observed for a further two AOA MAGs within the Nitrososphaerales. Both represented populations with no cultivated representatives at the genus level with MAG 112 representing TA-21 (amoA-gene lineage NS-b-2) [START_REF] Roux | VirSorter: mining viral signal from microbial genomic data[END_REF] and MAG 63 representing a novel genus (designated here as CS63) (amoA-gene lineage NS-d-2.1). The latter is a representative of the amoA gene-defined NS-dlineage which, while widely found in soil, is rarely implicated in contributing to ammonia oxidation. While the recovery of these two MAGs may simply be representative of AOA populations not active under the incubation conditions used, they are also consistent with organisms possessing physiologies distinct from cultivated AOA.

Active AOA viruses

After predicting contigs of virus origin using established bioinformatic tools, those from AOA-infecting viruses were identified with a custom database of hallmark genes from AOA prophages and analysis of other shared homologues (see Supplementary Text).

Sixty-four contigs (each representing a different viral operational taxonomic unit (vOTU)) were predicted as derived from AOA-infecting viruses in LBD DNA metagenomic libraries. In addition, 5 and 3 contigs were derived from viruses infecting ammonia oxidising bacteria (AOB) and NOB, respectively. Virus and host genomes typically have similar GC mol% [START_REF] Guo | VirSorter2: a multi-classifier, expert-guided approach to detect diverse DNA and RNA viruses[END_REF] and AOA-infecting virus contigs (average 36%) were similar to their hosts (38%). As with their hosts, AOA virus communities were distinct between the two pH soils (Figure 2A) as also observed previously for total virus communities [START_REF] Ren | Identifying viruses from metagenomic data using deep learning[END_REF] with 29 and 34 vOTUs significantly more relatively abundant in pH 4.5 and 7.5 LBD libraries compared to the other soil, respectively. Mapping of HDB sequence reads revealed that 42% of contigs were enriched in 13 C and thus confirming infection during incubation and transfer of C from AOA hosts to viruses.

Analysis of genome-wide sequence comparisons using translated open reading placed 60 of the 64 soil AOA viral contigs within one cluster but displayed no similarity to recently described viruses of marine AOA [5] or characterised non-AOA archaeal and bacterial viruses (Figure 2B). Eleven translated large terminase (TerL) sub-unit gene sequences found in contigs were similar to those recovered from other predicted AOA virus or prophage genomes only (e-value <10 -5 ) (Figure 2C). In addition to containing viral hallmark genes encoding portal, endonuclease, tail and capsid proteins, contig VC_053596 contained a type 1 multicopper oxidase (MCO1) homologue [START_REF] Nayfach | CheckV assesses the quality and completeness of metagenome-assembled viral genomes[END_REF] (Figure 2D) likely involved in the sequestration of copper from the environment which is essential for ammonia oxidation and electron transport in AOA [START_REF] Bushnell | A fast, accurate, splice-aware aligner[END_REF]. As MCO1 is found exclusively in AOA genomes [START_REF] Nayfach | CheckV assesses the quality and completeness of metagenome-assembled viral genomes[END_REF] it may represent a viral auxiliary metabolic gene contributing to a core process in AOA physiology.

Summary

These results from two soils with distinct AOA communities indicate that virus infection of AOA during soil nitrification is likely a common and dynamic process. Although viruses were only identified from contigs representing partial viral genomes, analysis of individual marker genes and genome-wide similarity comparisons revealed that those infecting AOA are highly divergent from known viruses. While one Nitrobacter MAG was obtained and both AOB and NOB viruses were also predicted, the hybrid approach developed in this study used both 12 C-and 13 C-enriched DNA to specifically target low GC mol% DNA and was tailored for the analysis of AOA populations and their viruses. A broader analysis of all viruses infecting nitrifying prokaryotic populations would therefore likely require a more comprehensive sequencing effort of all fractions throughout the CsCl gradient using quantitative SIP metagenomics [START_REF] Roux | Minimum Information about an Uncultivated Virus Genome (MIUViG)[END_REF] GTDB genus [START_REF] Chaumeil | GTDB-Tk: a toolkit to classify genomes with the Genome Taxonomy Database[END_REF] and amoA gene lineage-affiliation [START_REF] Roux | VirSorter: mining viral signal from microbial genomic data[END_REF] are given.

Supplementary text

Materials and methods

Soil microcosms, DNA extraction and stable-isotope probing

Soil microcosms were established in 144 ml serum vial bottles with 10 g soil (dry weight equivalent). Microcosms were amended urea solution (100 µg urea-N g -1 soil (dry weight)) at day 0 and 15 with the water content being 30% and 32% (w/w) after amendment, respectively. Headspace gas was amended with 5% (v/v) 12 C-CO2 (Air Liquide) or 13 C-CO2 (Sigma-Aldrich) (99% atom enriched) with microcosms opened every 3-4 days to maintain aerobic conditions before re-establishing CO2 concentrations. Microcosms were destructively sampled at day 0, 15 and 30 in triplicate for each treatment and soil stored at -20°C. Ammonium, nitrite and nitrate concentrations were determined immediately after sampling using standard colorimetric assays [1]. Total genomic DNA was extracted from 0.5 g soil using a CTAB-buffer phenol:chloroform:isoamyl alcohol protocol [2] and subject to isopycnic centrifugation in CsCl gradients as previously described [START_REF] Lee | Methane-derived carbon flows into host-virus networks at different trophic levels in soil[END_REF]. Briefly, 8 ml polyallomer tubes were filled with CsCl dissolved in Tris EDTA buffer (buoyant density 1.71 g ml -1 ) and 6 µg genomic DNA before sealing and centrifugation at 152,000 x g in an MLN80 rotor (Beckman-Coulter). CsCl gradients were then fractionated into ~350 µl aliquots before determining buoyant density via refractive index and recovery of DNA by PEG precipitation [2]. The predicted GC mol% of DNA in specific fractions was determined using the calculation of Schildkraut et al. [3].

Quantitative PCR and metagenomic sequencing of fractionated DNA samples

Quantitative PCR (qPCR) was performed to determine the distribution of prokaryotic and nitrifier genomes through CsCl gradients using a Corbett Rotor-Gene 6000 realtime PCR cycler. Prokaryotic 16S rRNA, bacterial and archaeal amoA genes were quantified across the entire buoyant density gradients using primer pairs P1(341f)/P2(534r) [4], amoA1F/amoA2R [5], and crenamoA23F/crenamoA616R [6]), respectively. Each 25 μl PCR contained 12.5 μl 2X QuantiFast SYBR Green Mix (Qiagen), 1 μM of each primer, 2 μl of standard DNA, recovered DNA from each fraction (1/10 dilution) or water (no template negative control). Thermal cycling programs consisted of an initial denaturation step of 15 min at 95°C, followed by 30 cycles of 15 s at 94°C, 30 s at 60°C, 30 s at 72°C for the 16S rRNA gene assay or 15 s at 94°C, 30 s at 55°C, 30 s at 72°C for bacterial and archaeal amoA genes assays followed by melt-curve analysis. All assays had an efficiency between 96 and 99% with an r2 value ≥0.99. After assessing the distribution of prokaryotic and nitrifier genomes, DNA in fractions with a buoyant density between 1.676 to 1.699 g ml -1 (low buoyant density) and 1.719 to 1.750 g ml -1 (high buoyant density) were pooled for metagenome sequencing. LBD and HBD metagenomes were sequenced at IntegraGen (Paris, France) and the Joint Genome Institute (JGI, Berkeley, CA, USA), respectively, both using the NovaSeq sequencing platform (Illumina) with NovaSeq XP version 1 reagent kits and a S4 flowcell 150 bp paired reads.

Assembly, binning, and annotation

Individual reads of HBD metagenomes were quality-trimmed with the MetaWrap read_qc module [7]. De novo co-assembly of the 135 to 245 million quality-controlled reads per metagenome was performed using MEGAHIT version 1.2.9 [START_REF] Li | MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph[END_REF]. Binning was performed using metaBAT 2 [START_REF] Kang | MetaBAT 2: an adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies[END_REF], MaxBin2 [START_REF] Wu | MaxBin 2.0: an automated binning algorithm to recover genomes from multiple metagenomic datasets[END_REF] and CONCOCT [START_REF] Alneberg | Binning metagenomic contigs by coverage and composition[END_REF] implemented in MetaWRAP version 1.2.1 [7] and bin completion and contamination estimated using CheckM version 1.0.12 [START_REF] Parks | CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes[END_REF]. Taxonomic annotation of contigs was performed using Kaiju [START_REF] Menzel | Fast and sensitive taxonomic classification for metagenomics with Kaiju[END_REF] with the NCBI nr database (2021-02-24). Taxonomy of MAGs was assigned using GTDB-Tk version 0.3.2 [START_REF] Chaumeil | GTDB-Tk: a toolkit to classify genomes with the Genome Taxonomy Database[END_REF] with the Genome Taxonomy Database (GTDB, release date 2022-03-23).

Identification of putative viruses infecting nitrifiers

Viral contigs (VCs) were predicted from contigs ≥10 kb using VirSorter [START_REF] Roux | VirSorter: mining viral signal from microbial genomic data[END_REF], VirSorter 2.0 [START_REF] Guo | VirSorter2: a multi-classifier, expert-guided approach to detect diverse DNA and RNA viruses[END_REF] and DeepVirFinder [START_REF] Ren | Identifying viruses from metagenomic data using deep learning[END_REF] (see Supplementary Figure 5 for comparison of predictions). CheckV [START_REF] Nayfach | CheckV assesses the quality and completeness of metagenome-assembled viral genomes[END_REF] and manual curation (i.e. identification of viral hallmark genes such as major capsid, portal, terminase and integrase proteins, enrichment in non-annotated genes) were performed to confirm a potential viral origin of each VC.

The relative abundance of each VC in HBD and LBD metagenomes was determined using the MetaWRAP-Quant-bins module [7] and BBMap [START_REF] Bushnell | A fast, accurate, splice-aware aligner[END_REF] was used to confirm a detection threshold of ≥75% of contig length with ≥1x read coverage recruited at ≥90% average nucleotide identity [START_REF] Roux | Minimum Information about an Uncultivated Virus Genome (MIUViG)[END_REF]. The relative abundance of each VC in each HBD and LBD metagenome was calculated based on the length of sequence size and coverage and abundance was expressed as normalized copies per million reads (CPM). Gene prediction of each VC was performed using Prodigal version 2.6.3 with meta option [START_REF] Hyatt | Prodigal: prokaryotic gene recognition and translation initiation site identification[END_REF], and homology search was performed against the NCBI nr database using Diamond blastp [START_REF] Buchfink | Fast and sensitive protein alignment using DIAMOND[END_REF] and Interproscan5 [START_REF] Jones | InterProScan 5: genome-scale protein function classification[END_REF].

To identify contigs from viruses infecting nitrifiers specifically, a database of virus genes present in proviruses of GTDB nitrifier reference genomes was constructed using VIBRANT [START_REF] Kieft | VIBRANT: automated recovery, annotation and curation of microbial viruses, and evaluation of viral community function from genomic sequences[END_REF], PhageBoost [START_REF] Sirén | Rapid discovery of novel prophages using biological feature engineering and machine learning[END_REF] and VirSorter [START_REF] Roux | VirSorter: mining viral signal from microbial genomic data[END_REF]. Major capsid, terminase, portal and integrase protein sequences were then used as search queries against each VC using Diamond blastp (e-value < 10 -5 ). Genes representing potential nitrifier-specific auxiliary metabolic genes were manually checked in annotation lists.

A homologue-based approach was also used where an AOA host was predicted when 'best hit' shared homologues with those in genomes of Thaumarchaeota (NCBI taxonomy) were ≥3x more abundant compared to the second most dominant phylum [START_REF] Al-Shayeb | Clades of huge phages from across Earth's ecosystems[END_REF]. Unlike our previous DNA-SIP analysis of methanotroph populations and associated viruses in these soils [START_REF] Lee | Methane-derived carbon flows into host-virus networks at different trophic levels in soil[END_REF], spacer sequences in AOA MAG CRISPR arrays only matched predicted AOA viruses from other studies with 2 mismatches between spacer and protospacer sequences and none from this study (data not shown).

Phylogenomic, phylogenetic and predicted protein analyses

Phylogenomic analysis of AOA MAGs and reference genomes from cultivated AOA and MAGs was performed using an alignment of single copy marker genes generated with GToTree [START_REF] Lee | GToTree: a user-friendly workflow for phylogenomics[END_REF] with a subsequent maximum likelihood phylogeny calculated using FastTree [START_REF] Price | FastTree 2--approximately maximum-likelihood trees for large alignments[END_REF]. Single gene phylogenetic analyses were performed using alignments generated using MUSCLE [START_REF] Edgar | MUSCLE: multiple sequence alignment with high accuracy and high throughput[END_REF] and manually refined. Maximum likelihood phylogenetic trees were constructed using PhyML [START_REF] Guindon | New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0[END_REF] with automatic model selection.

Affiliation of AOA MAGs to amoA gene-defined lineages was determined from phylogenetic analysis using the curated reference database of Alves et al. [START_REF] Alves | Unifying the global phylogeny and environmental distribution of ammonia-oxidising archaea based on amoA genes[END_REF]. Where amoA genes were absent in incomplete Nitrosotalea MAGs, a broad-level affiliation was inferred from close phylogenomic relationships to complete genomes and their amoA lineage designation. Comparison of genome-wide similarity between AOA virus contigs and curated virus reference sequences was performed using ViPTree [START_REF] Mihara | Linking Virus Genomes with Host Taxonomy[END_REF] and the Virus-Host DB [START_REF] Nishimura | ViPTree: the viral proteomic tree server[END_REF] to generate a protein tree based on normalised tBLASTx scores. AOA VC (n=64)

20 .Figure legends Figure 1 .

 20legends1 Figure legends

Figure 2 .Supplementary Figure 1 .Supplementary Figure 2 .Supplementary Figure 3 .

 2123 Figure 2. Identification and determination of activity of viruses infecting AOA in pH 4.5 and 7.5 soils using a hybrid analysis of 12 C-and 13 C-enriched DNA. A Heatmap

  Figure 1

2 )MAG 9 ,

 29 Metagenome sequencing of HBD 13 C vs 12 C DNA for read-mapping onto LBD MAG and virus contigs 1) Metagenome sequencing of LBD 12 C DNA for contig assembly and identification of MAGs and viruses 3) Identification of active AOA & viruses Nitrosotalea, NT-a-1 Nitrosotalea devanaterra ND1 (GCA_900065925) MAG 107, Nitrosotalea, NT-a-1 Nitrosotalea sinensis ND2 (GCF_900143675) MAG 77, Nitrosotalea, NT-a-1 Ca. Nitrosotalea okcheonensis CS (GCF_900177045) MAG 61, Nitrosotalea, NT-a-1 Ca. Nitrososphaera evergladensis SR1 (GCF_000730285) Nitrososphaera viennensis EN76 (GCA_000698785.1) MAG 12, Nitrososphaera, NS-a-3.2.3 Ca. Nitrososphaera gargensis Ga9.2 Nitrosocosmicus franklandianus C13 (GCF_900696045) Ca. Nitrosocosmicus oleophilus MY3 (GCF_008389435) Ca. Nitrosocosmicus hydrocola G61 (GCF_008389435) QCWD01 sp013114725 (GCA_013114725) MAG 96, TH1177, NS-e-2.2 TH1177 sp014523485 (GCA_014523485) Ca. Nitrosocaldus islandicus (GCA_013538715) Ca. Nitrosocaldus cavascurensis (GCF_900248165) MAG 52, Nitrososphaera, NS-a-3.2.

3

  Supplementary

Table 2 .

 2 Summary of 123 medium-and high-quality metagenome assembled genomes from low buoyant density DNA ranked by phylum (GTDB taxonomy).

	MAG	Completeness	Contamination	GC	N50	Size
	ID	(%)	(%)	mol%	(bp)	(bp)	Phylum	Class	Order	Family	Genus
		83.7	0.2	66.9	24,140	2,864,512 Acidobacteriota	Acidobacteriae	Acidobacteriales
		82.5	0.0	60.3	15,391	3,027,569 Acidobacteriota	Acidobacteriae	Acidobacteriales	Acidobacteriaceae	PALSA-350
		80.7	5.1	54.0	21,626	3,847,244 Acidobacteriota	Acidobacteriae	Acidobacteriales	Koribacteraceae	Gp1-AA145
		77.1	3.0	67.5	17,321	3,322,481 Acidobacteriota	Vicinamibacteria	Vicinamibacterales	UBA2999
		73.9	0.9	58.4	14,749	2,220,118 Acidobacteriota	Acidobacteriae	UBA7541	UBA7541
		67.3	3.4	63.9	11,528	2,512,459 Acidobacteriota	Acidobacteriae	Acidobacteriales
		66.0	5.5	59.3	15,405	2,887,037 Acidobacteriota	Acidobacteriae	UBA7541	UBA7541	Palsa-295
		53.3	0.9	56.3	9,893	2,452,695 Acidobacteriota	Acidobacteriae	UBA7541	UBA7541
		50.8	3.4	61.9	9,900	2,500,689 Acidobacteriota	Acidobacteriae	Acidobacteriales	Acidobacteriaceae	KBS-83
		50.4	2.0	52.4	8,779	1,812,072 Acidobacteriota	Blastocatellia	Pyrinomonadales	Pyrinomonadaceae	OLB17
		97.3	1.9	66.6	36,117	2,871,784 Actinobacteriota	Thermoleophilia	Solirubrobacterales	Solirubrobacteraceae
		97.0	2.6	66.0	142,875 3,274,599 Actinobacteriota	Thermoleophilia	Solirubrobacterales	Solirubrobacteraceae
		97.0	3.1	65.7	71,608	2,181,555 Actinobacteriota	Thermoleophilia	Solirubrobacterales	70-9
		82.0	0.1	63.9	17,436	1,359,101 Actinobacteriota	Acidimicrobiia	Acidimicrobiales	RAAP-2	RAAP-2
		79.2	2.2	68.5	54,707	2,098,665 Actinobacteriota	Thermoleophilia	Solirubrobacterales	70-9
		70.8	1.2	69.6	12,884	2,594,437 Actinobacteriota	Actinobacteria	Nanopelagicales
		63.0	3.0	70.1	11,121	2,066,628 Actinobacteriota	Acidimicrobiia	Acidimicrobiales	RAAP-2	Bog-756
		60.3	0.9	62.8	10,343	1,801,421 Actinobacteriota	Acidimicrobiia	Microtrichales	Ilumatobacteraceae	UBA668
		67.5	0.9	61.7	11,258	2,006,440 Armatimonadota	Fimbriimonadia	Fimbriimonadales	Fimbriimonadaceae
		96.5	0.7	44.9	52,178	3,501,364 Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	JJ008
		94.8	1.5	39.6	298,057 4,276,415 Bacteroidota	Bacteroidia	AKYH767
		92.0	1.7	44.7	26,881	3,591,767 Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	UBA8621
		86.9	7.5	38.5	28,454	3,315,570 Bacteroidota	Bacteroidia	AKYH767-A	OLB10
		84.0	0.7	49.6	26,115	3,041,112 Bacteroidota	Bacteroidia	AKYH767	2-12-FULL-35-15
		78.6	2.1	39.6	15,282	3,408,942 Bacteroidota	Bacteroidia	AKYH767
		77.5	0.2	54.5	15,405	2,291,910 Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae
		76.6	3.4	44.4	12,939	2,749,566 Bacteroidota	Bacteroidia	AKYH767	b-17BO	PALSA-968
		75.2	3.8	40.6	13,274	4,299,582 Bacteroidota	Bacteroidia	AKYH767-A	OLB10
		73.6	4.2	34.8	11,826	2,228,512 Bacteroidota	Ignavibacteria	Ignavibacteriales	Ignavibacteriaceae	BMS3ABIN03
		73.2	1.9	51.5	13,585	2,332,193 Bacteroidota	Kapabacteria	Palsa-1295	Palsa-1295
		71.5	0.0	42.1	13,719	2,460,739 Bacteroidota	Bacteroidia	NS11-12g	UBA955
		67.2	0.5	37.7	10,543	1,634,027 Bacteroidota	Bacteroidia	Sphingobacteriales	Sphingobacteriaceae
		66.5	1.5	42.0	11,489	3,546,402 Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae
		65.7	3.4	34.7	9,330	2,102,202 Bacteroidota	Bacteroidia	Flavobacteriales	Flavobacteriaceae	Flavobacterium
		65.4	0.0	39.1	18,038	2,305,465 Bacteroidota	Ignavibacteria	SJA-28	OLB5
		64.5	7.3	39.0	8,876	2,679,712 Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Taibaiella_B
		52.0	3.2	38.0	9,613	3,001,274 Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Palsa-955
		51.3	4.4	39.4	8,197	2,857,311 Bacteroidota	Bacteroidia	Chitinophagales	BACL12	UBA7236
		50.8	1.0	41.9	11,602	2,718,158 Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae
		52.6	1.8	42.5	10,353	1,916,736 Bdellovibrionota	Bdellovibrionia	Bdellovibrionales	Bdellovibrionaceae	Bdellovibrio
		91.4	0.0	61.2	14,550	3,678,188 Chloroflexota	UBA5177	UBA5177	UBA5177
		75.9	2.2	62.3	12,810	2,934,925 Chloroflexota	UBA4733	UBA4733	UBA4733
		75.1	2.8	70.3	18,220	2,481,858 Chloroflexota	Ellin6529	CSP1-4	CSP1-4	Palsa-1033
		64.2	1.4	55.1	23,141	1,640,525 Chloroflexota	Anaerolineae	Anaerolineales	UBA11579
		63.5	8.2	48.7	8,645	3,693,765 Chloroflexota	Anaerolineae	Anaerolineales	envOPS12	OLB14
		62.9	1.0	68.7	11,300	3,380,230 Chloroflexota	Ktedonobacteria	Ktedonobacterales

* * * * * * * * * * * * * *

Acknowledgements

This work was funded by the AXA Research Fund awarded to GN, ANR grant 'CONSERVE' awarded to CH and GWN, and an EC MSCA Fellowship 'DIVOBIS' awarded to ETS. Sequencing data generated under JGI Community Science Program Proposal No. 503702 awarded to GWN and CH was conducted by the US DOE JGI, a DOE Office of Science User Facility, and supported by the Office of Science of the US DOE under Contract No. DE-AC02-05CH11231. The authors would like to thank Dr. Willm Martens-Habbena (University of Florida) for continued discussions on soil AOA ecophysiology and taxonomy, Dr. Melina Kerou (University of Vienna) for helpful

Data availability

Metagenome sequence data are deposited in NCBI's GenBank under BioProject accession nos. PRJNA621418-PRJNA621429 ( 13 C data) and PRJNA868779 ( 12 C data).

Supplementary text references


Bertran Auvert
email: bertran.auvert@uvsq.fr

Joelle Sobngwi-Tambekou

Ewalde Cutler
email: ewaldecutlerewaldec@nicd.ac.za

Marthi Nieuwoudt
email: marthinieuwoudtmarthin@nicd.ac.za

Pascale Lissouba
email: pascale.lissouba@inserm.fr

Adrian Puren
email: adrianp@nicd.ac.za

Dirk Taljaard

Prof. Bertran Auvert

Effect of male circumcision on the prevalence of high-risk human papillomavirus in young men: results of a randomized controlled trial conducted

Keywords: HPV, male circumcision, Africa, randomized controlled trial

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

A recent meta-analysis estimated worldwide human papillomavirus (HPV) prevalence among women at 10•4% [START_REF] De Sanjose | Worldwide prevalence and genotype distribution of cervical human papillomavirus DNA in women with normal cytology: a metaanalysis[END_REF]. HPV genotypes are divided into "high-risk" and "low-risk" genotypes, on the basis of their association with cervical lesions. The high-risk human papillomavirus (HR-HPV) types are more frequently found in pre-malignant or malignant lesions and are associated with cancers of the cervix, vulva, vagina, anus, and penis [START_REF]Human papillomavirus HPV vaccines: Technical information for policy-makers and health professionals[END_REF][START_REF] Baldwin | Condom use and other factors affecting penile human papillomavirus detection in men attending a sexually transmitted disease clinic[END_REF][START_REF] Castellsague | Male circumcision, penile human papillomavirus infection, and cervical cancer in female partners[END_REF]. A causal association between cervical cancer and HR-HPV has now been established [START_REF] Castellsague | Male circumcision, penile human papillomavirus infection, and cervical cancer in female partners[END_REF][START_REF] Bosch | Prevalence of human papillomavirus in cervical cancer: a worldwide perspective. International biological study on cervical cancer (IBSCC) Study Group[END_REF][START_REF] Koutsky | A cohort study of the risk of cervical intraepithelial neoplasia grade 2 or 3 in relation to papillomavirus infection[END_REF][START_REF] Lehtinen | Human papillomavirus infection, risk for subsequent development of cervical neoplasia and associated population attributable fraction[END_REF][START_REF] Rozendaal | PCR based high risk HPV testing is superior to neural network based screening for predicting incident CIN III in women with normal cytology and borderline changes[END_REF][START_REF] Morris | Why circumcision is a biomedical imperative for the 21(st) century[END_REF] and worldwide HR-HPV prevalence of cervical carcinomas has been estimated at 99.7% [START_REF] Walboomers | Human papillomavirus is a necessary cause of invasive cervical cancer worldwide[END_REF]. Cervical cancer is the most common cancer affecting women in developing countries, with over 70% of cases occurring in Africa being attributed to HR-HPV genotypes 16 and 18 [START_REF] De Sanjose | Worldwide prevalence and genotype distribution of cervical human papillomavirus DNA in women with normal cytology: a metaanalysis[END_REF][START_REF] Schiffman | Human papillomavirus and cervical cancer[END_REF][START_REF] Munoz | Against which human papillomavirus types shall we vaccinate and screen? The international perspective[END_REF]. Thus, any factor reducing the probability of acquiring or transmitting HPV will also considerably reduce the burden of disease, especially in the developing world [START_REF] Castellsague | Male circumcision, penile human papillomavirus infection, and cervical cancer in female partners[END_REF].

Observational studies have suggested that HPV prevalence is reduced among circumcised men compared with uncircumcised men [START_REF] Baldwin | Condom use and other factors affecting penile human papillomavirus detection in men attending a sexually transmitted disease clinic[END_REF][START_REF] Castellsague | Male circumcision, penile human papillomavirus infection, and cervical cancer in female partners[END_REF][START_REF] Castellsague | HPV and circumcision: a biased, inaccurate and misleading meta-analysis[END_REF][START_REF] Hernandez | Circumcision and human papillomavirus infection in men: a site-specific comparison[END_REF][START_REF] Vaccarella | Prevalence and determinants of human papillomavirus infection in men attending vasectomy clinics in Mexico[END_REF]. Nevertheless, such an association has not yet been proven using a randomized controlled trial (RCT). The objective of this study was to analyze the effect of male circumcision (MC) on HR-HPV prevalence using data collected during a male circumcision RCT conducted in Orange Farm (South Africa), which demonstrated a partial protective effect of MC on the acquisition of HIV by young males.

Methods

Collection of data

The technical details of the trial (ANRS-1265) have been published elsewhere [START_REF] Auvert | Randomized, controlled intervention trial of male circumcision for reduction of HIV infection the ANRS 1265 Trial[END_REF] and only a summary will be presented in this article. Between February 2002 and July 2004, 3274 uncircumcised males, aged 18 to 24, were recruited, randomized into two groups and followed up. MC was offered immediately after randomization to the intervention group and after the end of the follow-up period to control group participants. During each follow-up visit at 3, 12 and 21 months, circumcision status was assessed by a nurse through genital examination. In addition, information about sexual behavior was collected, including number of partners as a function of time, number of sexual contacts with each partner, condom use and age of partners.

During 262 consecutive days, from March 7 th , 2005 to November 24 th , 2005, a urethral swab was collected by a nurse from all participants coming for the 21-month visit. All participants signed a written consent form for this test to be performed. Due to limited funding, the collection of swabs for HR-HPV testing was not started earlier. The urethra was chosen because the detection of HPV in this anatomical site is probably not affected by circumcision status. These swabs were analyzed to assess the association between the prevalence of HR-HPV strains and male circumcision. A urethral swab was also collected at a follow-up visit about 6 weeks after MC from all control group participants who took part in a nested study designed to compare two circumcision methods. To ascertain that the detection of HR-HPV was not affected by circumcision status, we used these swabs to compare the prevalence of HR-HPV among the nested study participants before and after circumcision. Lastly, study participants were asked to give a first-void urine sample in order to test for urogenital Neisseria gonorrhea (NG), which presence was used as a biological marker of sexual behavior.

Laboratory methods

Specimens were frozen at -20°C immediately after collection and kept frozen until processing. DNA was extracted from the urethral swabs using the MagNA Pure LC (Roche) instrument, with the Roche MagNA Pure LC DNA I Isolation Kit. Swabs were lysed in 500 µl of the kit lysis buffer for 30 minutes at room temperature. The MagNa Pure external lysis protocol was used to extract DNA from the lysis buffer into a 100 µl eluate. 50 µl of the eluate was used for screening (Roche Amplicor HPV test) and 50 µl eluate was used for genotyping (Roche Linear Array Genotyping test). This standardized PCR-based method can detect 13 HR genotypes of HPV (i.e., genotypes [START_REF] Auvert | Randomized, controlled intervention trial of male circumcision for reduction of HIV infection the ANRS 1265 Trial[END_REF][START_REF] Rosenbaum | The Central Role of the Propensity Score in Observational Studies for Causal Effects[END_REF]31,33,35,39,45,51,52,56,58,59, and 68). Because of the combined probe of the assay for HPV-52 and in order to be conservative, samples were classified as HPV-52 positive only when they were negative for genotypes 33, 35 and 58. Negative results with a negative internal beta globin PCR control were excluded. All positives were genotyped. An HPV-positive sample was defined by a sample where at least one HR-HPV was detected. In some analyses we also considered multiple HPV samples defined for samples where at least two HR-HPV genotypes were detected. Urine specimens were tested for NG by PCR (Roche Cobas Amplicor PCR).

Data analysis

The intention-to-treat and as-treated prevalence rate ratios (PRRs) of HR-HPV-positivity and NG-positivity were estimated using univariate log-Poisson regression These analyses were repeated multivariately by controlling for ethnic group, age, education, number of lifetime partners, marital status, number of non-spousal partners in the past 12 months, condom use in the past 12 months, number of sex acts in the past 12 months and HIV status.

To assess the potential impact of HIV acquisition, which is reduced by MC and is associated with HPV infection [START_REF] Ng | Epidemiology of human papillomavirus infection among fishermen along Lake Victoria Shore in the Kisumu District, Kenya[END_REF], these analyses were repeated after excluding those who seroconverted for HIV during the follow-up period (n=25).

To evaluate a possible imbalance between the groups, analyses were repeated when controlling for propensity score coded in quintiles [START_REF] Rosenbaum | The Central Role of the Propensity Score in Observational Studies for Causal Effects[END_REF].

Statistical analyses were performed using the statistical package SPSS for Windows version 8 (SPSS, Chicago, Illinois, United States) and R programming language (version 2.6.1) [START_REF] Team | R: A language and environment for statistical computing[END_REF].

Ethics

Results

Baseline characteristics of the 1264 participants from whom a urethral swab was collected at the 21-month visit are reported by randomization group in Table 1. These characteristics are similar and not statistically different except for HIV status. The mean (median) duration in days of follow-up among the intervention and control group was 644 (637) vs. 649 (637), respectively.

Table 2 presents the intention-to-treat univariate association between HR-HPV prevalence and male circumcision at the scheduled 21-month visit. HR-HPV prevalence was significantly lower among men of the intervention group. As indicated in Figure 1, the percentage of each of the 13 HR-HPV genotypes was always lower among men of the intervention group than among men of the control group. In the intention-to-treat comparison, the difference was significant for genotypes 18, 31, 45, 52, 56, 58, and 68.

Table 2 shows that the protective effect of MC on HR-HPV is higher in the as-treated analysis than in the intention-to-treat analysis. The protective effect is also higher in both analyses when controlling for potential confounders, including HIV status and reported sexual behavior cofactors. HR-HPV was associated with HIV status in both analyses with aPRRs of 2.2 (95%CI 1.5-3.3) and 2.2 (95%CI 1.5-3.2), respectively. When those who seroconverted for HIV during follow-up were excluded from the analysis, the results indicated in Table 2 remained practically unchanged with P-values less than 0.009 and a relative variation of the PRRs and aPRRs of less than 5.2%. This suggests that the effect of MC on HR-HPV is independent of the effect of MC on HIV. The aPRRs were almost identical when the analyses were adjusted for the propensity score in addition to the other covariates. Multiple HR-HPV prevalence was 7.0% (89/1267; 95% CI: 5.7% -8.6%). It was significantly lower among men of the intervention group compared with men of the control group (4.2% vs. 9.9%; PRR=0.43; 95% CI: 0.28 -0.66; Among men with at least one HR-HPV, multiple HR-HPV prevalence was also lower among men of the intervention group (44.3% vs. 28.7%; PRR=0.64; 95% CI: 0.45 -0.94; P=0.020).

As indicated in Table 3, NG prevalence was similar in the two groups. Among men of the control and intervention groups, median number of lifetime partners was 4.1 and 4.2, respectively, (P=0.49; Kruskal-Wallis test) and the proportion of consistent condom users was 17.4% and 19.7%, respectively, (P=0.45; Fisher exact test). These findings suggest that the protective effect of MC on HR-HPV cannot be attributed to a difference of sexual behavior between the two groups.

During the study period, 371 men of the control group were circumcised and had a urethral swab taken before and after male circumcision. The average (median) duration between the two swab collections was 59 (43) days. As expected, the HR-HPV prevalence was the same the two samplings (23.7% vs. 23.9%; Sign test). The proportion of males with multiple HR-HPV genotype infections was not significantly different (10.2% vs. 12.1%; P=0.40; Sign test). These results indicate that the as-treated effect of male circumcision on HR-HPV prevalence shown in Table 2 cannot be attributed to an easier detection of HR-HPV by urethral swabbing in uncircumcised men compared with circumcised men.

Discussion

Using data collected during the MC trial conducted in Orange Farm (South Africa), we demonstrated an independent and partial protective effect of male circumcision on HR-HPV prevalence. The effect was shown on HR-HPV prevalence and not incidence due to the available biological samples in this MC trial. This effect remained unchanged when the analysis was adjusted for possible confounding factors such as sexual behavior and condom use. Due to a) the randomization, b) the results of the propensity analysis, and c) the absence of obvious differences in the gonorrheal prevalence and sexual behavioral characteristics between the intervention and control groups, the HR-HPV prevalence difference between the two groups is likely attributable to male circumcision. In this view, the difference observed is probably the consequence of a difference in HR-HPV incidence between circumcised and uncircumcised men. Indeed, in this study HR-HPV prevalence is likely a proxy for HR-HPV incidence because among young men HPV prevalence is rising as a function of age [START_REF] Okesola | Prevalence of human papilloma virus genital infections in transmitted diseases clinic attendees in Ibadan[END_REF].

This study has some limitations. Biological samples were not collected throughout the followup period, so the HR-HPV status at inclusion is unknown. This information would have allowed us to compare HR-HPV incidence as a function of MC status and HR-HPV prevalence between intervention groups at inclusion. Because some participants were certainly already infected by HR-HPV at inclusion, the effect measured on prevalence is an underestimation of the true effect of MC. Secondly, the fact that participants were not blind to the intervention could have led to sexual behavior change and bias. Lastly, HR-HPV detection was performed on urethral swabs which are likely to miss infections [START_REF] Giuliano | The optimal anatomic sites for sampling heterosexual men for human papillomavirus (HPV) detection: the HPV detection in men study[END_REF]. Thus, the prevalence of HR-HPV infections in our cohort is likely to be underestimated since detection in the urethra is significantly lower than detection in the glans, corona sulcus or the penis shaft [START_REF] Giuliano | The optimal anatomic sites for sampling heterosexual men for human papillomavirus (HPV) detection: the HPV detection in men study[END_REF][START_REF] Aguilar | Human papillomavirus in men: comparison of different genital sites[END_REF]. However, we believe that there is no risk of non-differential misclassification since we did not find any difference when we compared the urethral HR-HPV prevalences before and after circumcision among a sub-sample of participants. Hence we believe that HRinfections will be underestimated equally among the two arms and will have no effect on PRRs. Despite this loss of power, our study evidenced a significant protective effect of MC against HPV infection.

In this study we could not detect an effect of MC on some HR-HPV genotypes such as genotypes 16 and 33. The apparent variation of the effect of MC according to genotype can be due to a true variation of the effect of MC according to genotypes or to random variation. This possible variation of the effect of MC according to genotypes should be further investigated for example by combining the results of this study with the results of the other MC trials conducted in Kenya and Uganda [START_REF] Gray | Male circumcision for HIV prevention in men in Rakai, Uganda: a randomised trial[END_REF][START_REF] Bailey | Male circumcision for HIV prevention in young men in Kisumu, Kenya: a randomised controlled trial[END_REF].

The protective effect has a magnitude corresponding to what could have been expected from observational studies. Castellsague and colleagues reported in their meta-analysis an odds ratio of 0.56 (95% CI: 0.39-0.82) [START_REF] Castellsague | HPV and circumcision: a biased, inaccurate and misleading meta-analysis[END_REF], while Baldwin and colleagues found an aRR of 0.44 (95% CI: 0.23-0.81) [START_REF] Baldwin | Condom use and other factors affecting penile human papillomavirus detection in men attending a sexually transmitted disease clinic[END_REF]. Similarly, Hernandez and colleagues found that uncircumcised men had nearly a two-fold (RR=1.96; 95% CI: 1.02-3.75) increased risk of oncogenic HPV infections [START_REF] Hernandez | Circumcision and human papillomavirus infection in men: a site-specific comparison[END_REF]. Based on the results of the RTC, there is now clear evidence that male circumcision decreases the heterosexual risk of HR-HPV acquisition by males. HR-HPV is a major public health problem because of its causal association with malignancies, especially cervical cancer in women. Hence this study illustrates why MC has long been thought to be protective against cervical cancer [START_REF] Morris | Why circumcision is a biomedical imperative for the 21(st) century[END_REF]. Indeed, as shown in our study, MC reduced the risk of HR-HPV infection among men and consequently reduced the exposure of women to HR-HPV. Thus, the risk of cervical cancer is lowered because of the causal link between HR-HPV and cervical cancer among women [START_REF] Castellsague | Male circumcision, penile human papillomavirus infection, and cervical cancer in female partners[END_REF][START_REF] Bosch | Prevalence of human papillomavirus in cervical cancer: a worldwide perspective. International biological study on cervical cancer (IBSCC) Study Group[END_REF][START_REF] Koutsky | A cohort study of the risk of cervical intraepithelial neoplasia grade 2 or 3 in relation to papillomavirus infection[END_REF][START_REF] Lehtinen | Human papillomavirus infection, risk for subsequent development of cervical neoplasia and associated population attributable fraction[END_REF][START_REF] Rozendaal | PCR based high risk HPV testing is superior to neural network based screening for predicting incident CIN III in women with normal cytology and borderline changes[END_REF][START_REF] Morris | Why circumcision is a biomedical imperative for the 21(st) century[END_REF].

Since three randomized controlled trials have shown a partial protective effect of MC on the acquisition of HIV by males in Africa [START_REF] Auvert | Randomized, controlled intervention trial of male circumcision for reduction of HIV infection the ANRS 1265 Trial[END_REF][START_REF] Gray | Male circumcision for HIV prevention in men in Rakai, Uganda: a randomised trial[END_REF][START_REF] Bailey | Male circumcision for HIV prevention in young men in Kisumu, Kenya: a randomised controlled trial[END_REF], the effect of MC on HR-HPV reinforces the WHO-UNAIDS recommendation for the implementation of MC programs in countries with a high HIV prevalence, a low MC prevalence and a high MC acceptability [START_REF]New Data on Male Circumcision and HIV Prevention: Policy and Programme UNAIDS-WHO[END_REF]. These countries, mainly in Southern and Eastern Africa, are those where the affordability of the HPV vaccine remains a problem. Moreover, the protective effect of MC may supplant HPV vaccines in terms of genotype coverage and target group age-range.

  The research protocol was reviewed and approved by the University of Witwatersrand Human Research Ethics Committee (Medical) on February 22 nd , 2002 (protocol study no. M020104). The trial was also approved by the Scientific Commission of the French National Agency for AIDS Research (ANRS; protocol study no. 1265; 2002, decision No. 50) and authorization was obtained from the City of Johannesburg, Region 11, on 25 February 2002. This trial has been registered at http://www.clinicaltrials.gov under the number NCT00122525.

Table 1 . Background characteristics, reported sexual behavior and HIV prevalence at the 21-month visit (

 1 1) Chi-square or Fisher exact test when appropriate (2) Kruskal-Wallis test (3) During the past 12 months (4) At some time during the past 12 months (5) Among those having had sexual intercourse during the past 12 months

	Control	Intervention	P
	n=627	n=637

Table 2 . Association between High-risk Human Papillomavirus (HR-HPV) prevalence and male circumcision

 2 . PRR= prevalence rate ratio * PRR adjusted for ethnic group, age, education, number of lifetime partners, marital status, number of non-spousal partners in the past 12 months, condom use in the past 12 months, number of sex acts in the past 12 months and HIV status. Distribution of the high-risk HPV genotypes as a function of randomization group.

		4%
	Percentage of cases	1% 2% 3%	Control Intervention
		0%
		16 18 31 33 35 39 45 51 52 56 58 59 68
		HPV genotypes
	Figure 1.

Acknowledgments

The authors would like to thank all those who agreed to take part in this study, answer the questions put to them, and provide swab samples. The authors would like to thank Gaph Sipho Phatedi for his management of the recruitment process and Yvon de La Soudière for his help in the management of the data set. The authors would like to thank the general practitioners who performed the male circumcisions for the study (Dr Bhekuyise Gwala, Dr George Shilaluke and Dr Dumiso Zulu). The authors would like to thank Goliath Gumede for the clinical investigation and Zodwa for interviewing all the respondents. They would also like to thank Bongiwe Klaas for the data capture, Mabel Hunter and the recruitment staff and all the assistants (Cynthia Dlamini; Sidwell Dumisi; Benjamin Masitenyane; Robert Tsietsi Mbuso; Anthony Motha; Sibongiseni Mpetsheni; Jabulani Nhlapo; Joseph Ntsele; Male Chakela; Audrey Tshabalala; Donald Mashamba; Nkululeko Nhlapo) for their cooperation and support. Lesley Short, Moses Mashiloane, Beulah Miller, Beverley Singh, Sarah Hloma and the HIV serology laboratory of the National Institute for Communicable Diseases, Johannesburg, South Africa provided technical assistance in regard to the laboratory testing and administration.

Financial support ANRS grant 1265 (France), Gates Foundation (USA) grant 33759, NICD (South Africa), INSERM (France).

NG prevalence % (positive/total) PRR (95% CI; P) aPRR* (95% CI; P)