# Insights on the particle-attached riverine archaeal community shifts linked to seasons and to multipollution during a Mediterranean extreme storm event

Mégane Noyer, Maria Bernard, Olivier Verneau, Carmen Palacios

► To cite this version:

Mégane Noyer, Maria Bernard, Olivier Verneau, Carmen Palacios. Insights on the particle-attached riverine archaeal community shifts linked to seasons and to multipollution during a Mediterranean extreme storm event. Environmental Science and Pollution Research, 2023, 10.1007/s11356-023-25637-x. hal-04028116

HAL Id: hal-04028116

https://hal.inrae.fr/hal-04028116

Submitted on 10 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 Insights on the particle-attached riverine archaeome at different
2 seasons and in response to multipollution during a Mediterranean
3 extreme storm event

4 Mégane Noyer¹,², Maria Bernard³,⁴, Olivier Verneau¹,²,⁵, Carmen Palacios¹,²,*

5 ¹Univ. Perpignan Via Domitia, CEFREM, UMR5110, F-66860, Perpignan, France.

6 ²CNRS, CEFREM, UMR5110, F-66860, Perpignan, France.

7 ³Univ. Paris-Saclay, INRAE, AgroParisTech, GABI, 78350, Jouy-en-Josas, France.

8 ⁴INRAE, SIGENAE, 78350, Jouy-en-Josas, France.

9 ⁵Unit. for Environmental Sciences and Management, North-West University, ZA-2520,
10 Potchefstroom, South Africa.

11 *Corresponding author: Carmen Palacios; carmen.palacios@univ-perp.fr. Centre de
12 Formation et de Recherche sur les Environnements Méditerranéens UMR 5110 CNRS-UPVD
13 Université de Perpignan Via Domitia 52 Avenue Paul Alduy 66860 Perpignan Cedex, France.
14 Phone: + 33468661791. Fax: + 33468662096.

15

16 Acknowledgments

17 We would like to thank Cristiana Cravo-Laureau for her suggestions on an advance version of
18 this manuscript. We are thankful to B. Reoyo-Prats, D. Aubert, A. Sanchez-Garcia, J. Sola, N.
19 Delsaut, S. Kunesch and C. Menniti (CEFREM laboratory) for helping us during sampling. We
20 thank Research and Testing Laboratories for kindly providing NGS primer assay list. We
21 would like to acknowledge the Genotoul bioinformatics platform Toulouse Midi-Pyrénées and
22 Sigenae group (http://bioinfo.genotoul.fr; Toulouse, France) for computing and storage
23 resources. We also thank D. Ning and N. Xiao for MENAP support
24 (http://ieg4.rccc.ou.edu/mena/main.cgi).
25

**Abstract**

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

Even if Archaea deliver important ecosystem services and are major players in global biogeochemical cycles, they remain poorly understood in freshwater ecosystems. To our knowledge, no studies specifically address the direct impact of xenobiotics on the riverine archaeome. Using environmental DNA metabarcoding of the 16S ribosomal gene, we previously demonstrated bacteriome significant responses to pollutant mixtures during an extreme flood in a typical Mediterranean coastal watercourse. Here, using the same methodology, we sought to determine whether archaeal community shifts were also driven by environmental stressors during the same flood. Further, we wanted to determine how archaea taxa compared at different seasons. In contrast to the bacteriome, the archaeome showed a specific community in summer compared to winter and autumn. We also identified a significant relationship between *in situ* archaeome shifts and changes in physicochemical parameters along the flood, but a less marked response to river hydrodynamics than bacteria. New urban-specific archaeal taxa, which were significantly related to multiple stressors, were identified. Through statistical modeling of both domains, our results demonstrate that Archaea, seldom considered as bioindicators of water quality, could provide a rapid assessment of microbial pollution risk and thus have the potential to improve monitoring methods of watersheds.

**Keywords:** microbial ecotoxicology; water quality; extreme storm event; coastal Mediterranean rivers; sewer overflow; environmental DNA; metabarcoding.
45

# 1. Introduction

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

In the 1970s, Woese and collaborators highlighted a new domain of life, the Archaea, distinct from Bacteria and Eukaryota domains (Woese and Fox 1977; Woese et al. 1978). With the advent of Next Generation Sequencing (NGS), archaeal research and knowledge are expanding (Adam et al. 2017; Bang and Schmitz 2018). Archaea are recognized as major players in the global biogeochemical cycles of carbon, nitrogen, hydrogen and sulfur (Casamayor and Borrego 2009; Offre et al. 2013; Castelle et al. 2015) and there are at least two metabolisms, essential for nutrient cycling, which are carried out exclusively by archaea: methanogenesis and anaerobic methane oxidation (Joye 2012). Despite their importance for ecosystem functioning and their ubiquitous presence (Chaban et al. 2006; Casamayor and Borrego 2009; Herfort et al. 2009; Auguet et al. 2010), the environmental archaeome has been overlooked for decades. Most studies have concentrated on marine archaea (Zinger et al. 2012; Zeglin 2015), and the diversity and importance of archaea in other environments have been largely disregarded (Casamayor and Borrego 2009; Zinger et al. 2012; Zeglin 2015; Adam et al. 2017).

Only recently, continental freshwater habitats have emerged as one of the largest reservoirs of archaeal diversity (Auguet et al. 2010; Zinger et al. 2012). Some studies have linked archaeal communities to water physicochemical properties such as pH, temperature and nutrients (Herfort et al. 2009; Wang et al. 2018; Lei et al. 2020). *Methanobrevibacter smithii*, for example, was proposed as a potential indicator of human-specific sewage pollution (Johnston et al. 2010; McLellan and Eren 2014). To assess water quality, several studies have focused on freshwater archaea (Cannon et al. 2017; Dila et al. 2018), notably for their potential remediation role in heavily contaminated urban rivers (Samson et al. 2019; Lei et al. 2020). Cannon *et al.* (2017) found that a rain event induced changes in the structure of microbial communities, including archaea, from environmental DNA (eDNA) and stressed the importance of considering hydrological conditions when studying riverine microbiomes. Both chemical and microbial pollutants reach surface waters via point sources (such as urban sewage) or diffuse sources (linked to runoff) of pollution. It is well known that during rainfall events, particles from terrestrial soils and river basin sediments remobilize, and affect water quality (Garcia-Esteves et al. 2007; Dumas et al. 2015; Faure et al. 2015) because soils and sediments store pathogens, nutrients and pollutants. Thus, suspended particles have a crucial role in the transfer of contaminants to surface waters through runoff and in the resuspension of river sediments during storms (Turner and Millward 2002; Amalfitano et al. 2017). Moreover, some studies have reported that riverine archaea tend to be associated with particles and derive only from
78 allochthonous inputs (Crump and Baross 2000; Casamayor and Borrego 2009; Hu et al. 2018).
79 This is all the more important in Mediterranean regions, where extreme hydrological events
80 are expected to become more intense and frequent due to climate change (Cowling et al. 2005;
81 Blanchet et al. 2016). Furthermore, first-flush events during rainfalls lead to Combined Sewer
82 Overflows (CSOs), which carry large loads of contaminant mixtures over surface waters
83 through in-sewer solids resuspension (Ashley et al. 1992; Osorio et al. 2012; Oursel et al. 2014;
84 Reoyo-Prats et al. 2017, 2018). We previously demonstrated that an extreme Mediterranean
85 flood produced shifts in the particle-attached bacterial compartment from eDNA, severely
86 affecting resident riverine communities (Noyer et al. 2020). In this study, we went a step further
87 by using metabarcoding of the 16S ribosomal RNA gene (rDNA) to explore how the diversity
88 of the particle-attached archaeome changed within the same environmental samples, and to
89 compare the responses between the bacteriome and the archaeome. Next, we modeled shifts in
90 both communities using physical parameters as well as several families of chemical parameters
91 as environmental forces, which is a first in environmental microbiology. This study addressed
92 the following questions: how did the fluviatile particle-attached archaeal community change
93 between seasons and how did it evolve over the course of a heavy rain event? How our findings
94 compare to other studies of archaeal alpha and beta diversity in lotic ecosystems? Did
95 environmental parameters drive structural shifts in the archaeome as they did in the
96 bacteriome? Was there a strong relationship between key taxa and environmental dynamics?
97 How did the seasonal and temporal succession of bacterial and archaeal communities compare?
98 The answer to these questions could, in general, provide insights on the use of microorganisms
99 in water quality assessment and, in particular, help on the rapid risk assessment of multiple
100 pollutants in aquatic ecosystems.

## 2. Materials and Methods

### 2.1. Study site, sampling information

103 Sampling took place in the Têt River, a watercourse representative of Mediterranean coastal
104 watercourses with a torrential regime that discharges into the Gulf of Lion (Southeast of
105 France) (Dumas et al. 2015; Reoyo-Prats et al. 2017) downstream from the Perpignan city
106 wastewater treatment plant, the main threat to the water quality of this river (Fig. 1a, Conseil
107 Général des Pyrénées Orientales 2009, 2012; Reoyo-Prats et al. 2017). The Vinça dam, situated
108 40 km upstream from the sampling station, controls the downstream river flow, particularly
during floods. For sampling methodology and sampling site see Reoyo-Prats *et al.* (2017). In short, ten liters of river water were collected in summer, winter, and autumn, before, during and after a 5-year flood that we followed 24h/24h during more than 4 days (see Fig. 1b for further details).

First flushes during the flood brought the highest levels of *E. coli* and *Enterococci* ever detected in this river (Reoyo-Prats et al. 2017) as well as relative higher abundance of other typical sewer bacteria (Noyer et al. 2020) and higher levels of dissolved pharmaceuticals (Reoyo-Prats et al. 2018), which indicated the moment at which sewer overflows occurred (Fig. 1b). Environmental DNA sampling was previously described by Noyer *et al.* (2020). In short, a liter of mixed-water sample was entirely (or until clogged) filtered through cellulose acetate MF-Millipore membrane filters with 3 µm porosity (Merck Darmstadt, Germany) and repeated three times to obtain three replicates per sample.

## 2.2. Nucleic acid extraction, 16S rRNA gene sequencing, and sequence analyses to obtain Operational Taxonomic Units (OTUs) contingency table

Nucleic acid extraction and bacterial sequencing from eDNA were described in a previous study (Noyer et al. 2020) and were performed in triplicate per sample. For archaea sequencing, a single set of DNA replicates from each sample was sent to the Research and Testing Laboratory (RTL, Texas, USA). Two more replicates were later sent to the Genome Quebec laboratory (GQ, Montreal, Canada), except for samples t61 and t68 which were only sequenced once again for technical reasons. DNA from samples t19 and t23 from the first replicate, which was already sequenced at RTL, was also sent to GQ to be sequenced again for comparative purposes. Sequencing targeted V3 and V4 hypervariable regions of the 16S rDNA by using 519wF (5'-CAGCMGCCGCGGTAA-3') and 1017R (5'-GGCCATGCACCWCCTCTC-3') universal archaeal primers (Borrego et al. 2020) and was performed on an Illumina MiSeq sequencer using a 2x300bp paired-end protocol. Libraries were generated by pooling equimolar ratios of amplicons before sequencing. In contrast to RTL sequences, those provided by GQ contained primers. In order to pool sequences together, we eliminated primers from GQ sequences using Cutadapt version 1.8.3 for Unix (Martin, 2011). Default options were used with the exception of sequences treatment as paired with -g and -G options for forward and reverse primers respectively, taking wildcards into account, discarding untrimmed sequences and setting an overlap of 14 bp, a quality-base of 33 and an error-rate of 0.1. Harmonized sequences were pooled together and archived before being uploaded to the
Galaxy platform (Afgan et al., 2016): http://sigenae-workbench.toulouse.inra.fr. FROGS pipeline (v. r3.0-3.0, Escudié et al. 2018) was used to process sequences to form OTUs and to taxonomically affiliate them as described in Noyer et al. (2020) except that VSEARCH (v2.6.2, Rognes et al. 2016) was used as a read pair assembler, which allows a higher number of sequences to be conserved when amplicon sizes are highly variable, as was the case for our dataset. Read length was set to 300pb and the amplicon size was set to 450pb for the minimum and optimized to 545pb for the maximum. We assigned each OTU up to the species taxonomic level based on blast alignment using the affiliation tool (v. r3.0-2.0) and the Silva 138.1 database with a pintail score of 100, which allowed for the most accurate affiliation possible. No archaea hits were removed hereafter. To resolve the taxonomic ambiguity of OTUs that were multi-affiliated within the FROGS pipeline when using blast against the curated pintail 100 score silva database, we blasted these OTUs against the NCBI 16S nucleotide collection database using the megablast algorithm.
<figure><img src="image_2.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_3.png" /></figure>

<table><thead><tr><th>Summer drought</th><th>Winter drought</th><th>Autumn drought</th><th>Autumn flood</th></tr></thead><tbody><tr><td>SD</td><td>WD</td><td>t0</td><td>12 sampling times: t17 → t109</td></tr></tbody></table>
154

**Fig. 1** Têt River archaeome sampling sites and environmental parameters measured in the same samples. (a) Watershed of the river with sampling site (black arrow), located after wastewater treatment plan (WWTP) of the city of Perpignan, combined sewers (black crosses) and water reservoirs indicated as grey rectangles (adapted from Reoyo-Prats et al. 2017). (b) Environmental parameter dynamics in the Têt River at different seasons and along an autumn flood. For sample names, see table below figure. Autumn sample names are followed by a number that indicates the sampling time in hours after t0, which was sampled at autumn basal level water discharges. Sampling took part at crucial moments of the flood that occurred thereafter: at first flushes (t17-t19-t23), before the flow peak (t32-t37), during the flow peak (t38-t41), following the release of water from the upstream Vinça reservoir (t44-t61) and during the return to basal level (t68-t86-t109). ProxyDyn1 corresponds to the dynamics of particulate organic carbon (/20 mg/l), which represented the dynamics of water flow, also represented in figure, total suspended solids, total organic carbon, total nitrogen, and terbuthylazine parameters. ProxyDyn2 corresponds to aminomethyl phosphonic acid (AMPA, µg/l), which represented glyphosate, phosphate, copper, temperature, E. coli, enterococci, diclofenac, sulfamethoxazole and carbamazepine parameters. ProxyDyn3 corresponds to lead (/150 µg/g) in the representation of the dynamics of cadmium, zinc, and conductivity parameters. ProxyDyn4 corresponds to pH (/70), which represented cobalt, nickel, and chrome parameters. Three parameters, Diuron, 2,4D and NO<sup>3</sup>-, had a unique dynamic. For further details on statistical analyses for environmental parameters, see Noyer et al. (2020).

155

156

157

158

159

160

161

162

163

164

165

166

167

168
169
170
171

**2.3. Archaeal diversity analyses**

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
194

visualize replicate dissimilarities. Once outliers checked, beta dissimilarities were recalculated
on the averaged OTU abundances. To this end, the number of OTUs decreased to less than
10,000 after the dataset was normalized, thus allowing qualitative Unifrac and quantitative
Weighted-Unifrac (W-U) dissimilarities, which also consider phylogeny of OTUs, to be
included within FROGS. PCoA and hierarchical clustering Ward.D2 methods were used to
visualize archaeal community dissimilarities among samples. A one-way analysis of similarity
(ANOSIM, Clarke 1993) was performed to test significant differences between sample groups
resulting from hierarchical clustering. To further check for significant differences in archaeal
community shifts at the class level, the Mann-Whitney test (MW) was implemented with R
software.
203

195

196

197

198

199

200

201

202

204

## 2.4. Statistical analyses for inference

205

2.4.1. Constrained (canonical) ordination analyses by environmental parameters dynamics

206

Physicochemical environmental parameters were previously measured (Reoyo-Prats et al. 2017) in the same samples in which nucleic acid extractions were performed. Measured parameters included pH, temperature, conductivity, flow, total suspended solids, particulate organic carbon, total organic carbon, nitrogen and 250 pesticide molecules, 90 pharmaceutically active compounds, polycyclic aromatic hydrocarbons and polychlorinated biphenyls, nutrients, trace metals in the particulate fraction and fecal indicators (load in *Escherichia coli* and Enterococci). Collinearity issues resolution and the seven major environmental dynamics retained for further analyses are described in Noyer *et al.* (2020). For clarity, the retained variables that will be used for further analyses as proxies of correlated environmental variables are summarized in Fig. 1b. Constraint-based ordination analyses were then used to evaluate the relationships between the normalized OTU contingency table and retained environmental parameter dynamics using *vegan* R package version 2.5-3 (Oksanen et al. 2018). Detrended correspondence analysis (DCA) performed on the OTU dataset rendered a first axis gradient length of 3.7, so both canonical correspondence (CCA) and redundancy analyses (RDA) were performed (ter Braak 1988). We also performed a Hellinger transformed-based RDA (tbRDA) (Legendre and Gallagher 2001) and a distance-based RDA (dbRDA) using all beta dissimilarities from the previous section. Permutation analyses of variance were used to evaluate the significance of constraint-based models, axes, and variables. Variables were tested by adding each of them independently and the number of permutations was set to 10,000. To further determine which OTUs best responded to environmental variables, we reduced the OTU matrix to a percentage of abundance so that CCA/RDA modeling became

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226
227

significant. Instead of 0.005%, as reported by Bokulich et al. (2012) and Noyer et al. (2020) for bacteria, we found a percentage of 0.05% for the archaeal dataset. RDA modeling of OTUs with ≥ 0.05% of the total read number (i.e. keeping 2,688,328 reads), allowed for the OTU goodness of fit (GOF) to be calculated. As for bacteria, OTUs retained for further analyses had a GOF ≥ GOFaverage, which is considered a conservative approach to OTU selection.

228

229

230

231

232

### 2.4.2. Network construction via module eigengene analysis

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

We used the Molecular Ecological Network Analyses Pipeline (MENAP: http://ieg4.rccc.ou.edu/MENA/) to build the relationships among OTUs following the developer's recommendations (Deng et al. 2012; Tao et al. 2018) but sample-specific OTUs were kept for network construction by changing "OTUs present at least in one sample". An automatically generated similarity threshold value (0.32) was obtained with Random Matrix Theory (RMT)-based method, which allowed network construction ensuring that the connections between microorganisms were non-random (R² of power-law = 0.28). The network was separated into modules via the short random walk method (Pons and Latapy 2005), which had the highest modularity (0.078) (Newman 2004). Module higher-order organization was then performed via eigengene analysis (Langfelder and Horvath 2007) using default parameters to obtain the correlation significance between modules and environmental parameters dynamics. Cytoscape software (v. 3.7.1, Shannon et al. 2003) was used to visualize this constraint-based network. We also used MENAP to explore the relationship between OTUs from Bacteria (from the previous study by Noyer et al., 2020) and Archaea domains and flood environmental dynamics. The network was constructed with the same parameters as before except zero counts were replaced "by 0.01" instead of "on paired values", to avoid losing an important number of OTUs when using the matrix of joined OTUs from both domains. The highest modularity was obtained via the leading eigenvector method (0.520). The RMT threshold was 0.940 and the R² of power law = 0.863.

252

## 3. Results

253

### 3.1. Changes in particle-attached archaeal diversity along a storm event in the Têt River

254

255

256

257

To study changes in an urban archaeome along a storm event, a high-resolution environmental DNA sampling from river water was performed during a major flood that occurred in autumn 2013. Two other seasons were also sampled during the summer and winter drought periods.
258

All reported sequence data passed standard quality controls of certified sequencing companies.

259

A total of 1,804,302 reads corresponding to 350,106 operational taxonomic units (OTUs) were identified. Most OTUs (98%) were assigned to the archaeal domain, the rest were eliminated from further analysis.

260

261

262

Alpha diversity statistical analyses showed significant differences among samples for Shannon (KW = 0.03, Fig. 2) and Pielou (KW = 0.03, Table S1) indices only. These two indices changed similarly along seasons and during the flood, with equivalent significant differences between samples except for the Pielou index for the summer drought, which was significantly higher than the flow peak (t38-t41) indicating a greater evenness in OTU abundances in this sample.

263

In general, both indices decreased significantly at t41, i.e. right after the flood peak, with respect to t17 and t19. Three hours later, a significant increase in diversity was noticed (t44).

264

Even if the observed OTU number, Chao1 and Fisher indices did not show significant differences (Table S1), they had the same change pattern as the Shannon index (Fig. 2), except for a much lower value in the summer drought sample. This result contrast with the Pielou and Simpson indices (Table S1), which emphasize the evenness component of diversity, as opposed to Fisher and Shannon indices, which appear more related here to the richness component of diversity (Magurran 2004).

265

266

267

268

269

<figure><img src="image_4.png" /></figure>

275

276

<caption>**Fig. 2** Alpha diversity of the archaeome of the Têt river along time. Changes in observed OTU number (blue) and Shannon index (red) along the flood (tX) and at summer and winter droughts (SD and WD respectively). For the</caption>

277
Shannon index, different letters indicate a significant difference between samples (dunn.test<0.025) and red arrows show major significant differences. Observed OTU number was not significantly different (KW=0.14). The dotted profile is the flow level at each sampling point (see Fig. 1, also for sample names). Even though the absence of replicates for t61 and t68 samples impeded statistical testing, they are represented through time for comparison.

After normalization of OTU matrix once singletons were excluded, a total of 1,377 OTUs (3,777 sequences/sample) were retained for further analyses. We noticed a great difference in the taxonomic composition of the summer drought (SD) with respect to all other samples. SD had a significantly lower amount of Nitrososphaeria and Thermoplasmata classes when compared to all samples (Fig. 3, MW p=0.012 and p=0.021 respectively) and a higher amount of Nanoarchaeia (Nanoarchaeota phylum, Fig. 3) and Bathyarchaeia (Crenarchaeota phylum) classes when compared to most samples. These large differences in the summer sample community were confirmed with the beta diversity analyses (data not shown), as all samples significantly separated from the SD sample with all beta dissimilarity indices, contrarily to WD, which always clustered with t17, t19, and t32. Additionally, the taxonomic composition of the WD sample was very similar to t17, t19, and t32 samples, mainly composed of Nitrososphaeria class, followed by Nanoarchaeia, Thermoplasmata, and then three classes, Methanobacteria (Euryarchaeota phylum), Methanomicrobia (Halobacterota phylum) and Bathyarchaeia at variable smaller abundances (Fig. 3).

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

**Fig. 3** Composition of archaeal communities averaged across replicates. Histogram of relative abundances (a) of the six major phyla and (b) of the ten major classes. Samples are organized according to sampling time from left to right: summer drought (SD), autumn flood (sample names are followed by a number that indicates the sampling time in hours after the beginning of the flood at t0), and winter drought (WD). The dotted profile is the flow level at each sampling point (see Fig. 1 for further details).
Given these results, only flood samples were used to further explore structural diversity
changes using beta dissimilarities. These samples comprised a total of 1,739 OTUs (19,355
sequences/sample) after eliminating singletons and normalization. Using different beta
dissimilarities allowed for a better assessment of which differences are responsible for the
community structure (either presence/absence or abundance and/or phylogeny of OTUs).
Using only OTU presence/absence with Jaccard qualitative dissimilarity, three communities
were significantly differentiated (Fig. S2.1a). The first community group included samples
collected at t0, t17, t19 and t32, the second included samples t23, t37, and t38, and the third
included the rest of the samples (t41-t109), which was the most differentiated cluster. But when
phylogenetic relationships were considered using Unifrac qualitative dissimilarity, t41 sample
grouped with t23, t37 and t38 instead (Fig. S2.1b). With Unifrac, this last cluster was the most
differentiated. When Bray-Curtis quantitative dissimilarity, which considers OTU abundance,
was used, three community groups of samples were distinguished (Fig. S2.1c). The first
included from t17 to t38 samples, the second included t0, t41, and t68 samples, and the third
t44, t61, t86, and t109 samples. When considering phylogenetic relationships using W-U,
observed community groups coincided with those of Bray-Curtis dissimilarity (Fig. 4, bold
black line), but the most significantly differentiated communities were those from the group of
t23, t37 and t38 (axis 1) and then t41 and t0 samples (axis 2). Even though the first taxonomic
changes occurred from the first flood sample at t17, samples t23, t37 and t38 had a particular
significant increase in Methanobacteria (Euryarchaeota phylum) (Fig. 3, MW, p=0.001).
Sample t0 differentiated significatively from samples from the end of the flood (from t44 to
t109) only when considering OTUs presence/absence but not when abundance alone or with
phylogeny were considered. This sample had a significantly lower abundance of taxa from the
Thermoplasmata class with respect to that group of samples (Fig. 3b, MW p=0.003). Regarding
t41 sample community, it had a particular taxonomy, significantly dominated by
Nitrososphaeria class (MW p=0.002) and W-U dissimilarity significantly separated t41 from
the rest of the samples (Fig. 4). With regard to the group of samples at the end of the flood
event, from t44 to t109, we could notice a significant increase of Thermoplasmata class with
respect to all other flood samples (Thermoplasmatota phylum, MW p=0.002, Fig. 3b).
<figure><img src="image_7.png" /><figcaption>Fig. 4 Structure of archaeal communities averaged across replicates. (a) Principal Coordinate Analysis (PCoA) and (b) hierarchical clustering with Ward D2 linkage method using Weighted-Unifrac dissimilarity computed on OTU average abundance. Lines indicate ANOSIM significant groups. Sample names are followed by a number that indicates the sampling time in hours after t0 (see Fig. 1 for details). Significant codes ** and *** indicate p-value < 0.01 and < 0.001, respectively.</figcaption></figure>

<table><caption>**Table 1** Summary of constraint-based multivariate statistical models on archaea OTU matrix averaged over replicates and without singletons. (a) Permanova significance of the five models tested and the percentage of biological variance that is explained by each model using permutation test with anova.cca function. (b) Axes and modeled variables significance after permanova using anova.cca of significant models in (a). Axes not shown were not significant. p-values significance codes: (*** < 0.001 < ** < 0.01 < * < 0.05).</caption><thead><tr><th>a</th><th>OTUs matrix<br>transformation</th><th>Model<br>significance</th><th>Variance (%)</th></tr></thead><tbody><tr><td rowspan="4">dbRDA</td><td>Jaccard</td><td>0.003**</td><td>63.38%</td></tr><tr><td>Unifrac</td><td>0.008**</td><td>66.35%</td></tr><tr><td>Bray-Curtis</td><td>0.326</td><td>60.04%</td></tr><tr><td>Weighted-Unifrac</td><td>0.004**</td><td>76.58%</td></tr><tr><td colspan="2">CCA</td><td>0.765</td><td>56.44%</td></tr><tr><td colspan="2">RDA</td><td>0.134</td><td>61.4%</td></tr></tbody></table>

<table><thead><tr><th rowspan="2">b</th><th rowspan="2">OTUs matrix<br>transformation</th><th colspan="8">Statistical significance of modeled variables</th></tr><tr><th>CAP1</th><th>ProxyDyn1</th><th>ProxyDyn2</th><th>ProxyDyn3</th><th>ProxyDyn4</th><th>2.4D</th><th>Diuron</th><th>NO<sub>3</sub><sup>-</sup></th></tr></thead><tbody><tr><td rowspan="3">dbRDA</td><td>Jaccard</td><td>0.002** (22.74)</td><td>0.001***</td><td>0.002**</td><td>0.001***</td><td>0.105</td><td>0.095</td><td>0.001***</td><td>0.096</td></tr><tr><td>Unifrac</td><td>0.040* (30.24)</td><td>0.001***</td><td>0.021*</td><td>0.046*</td><td>0.458</td><td>0.383</td><td>0.036*</td><td>0.014*</td></tr><tr><td>Weighted-Unifrac</td><td>0.006** (42.89)</td><td>0.006**</td><td>0.009**</td><td>0.039*</td><td>0.050*</td><td>0.004**</td><td>0.002**</td><td>0.166</td></tr></tbody></table>
346

## 3.2. Modeling archaeome diversity according to multiple contaminant dynamics

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

Constrained multivariate analyses were performed to determine if the retained environmental dynamics (Fig. 1b) could statistically explain the observed community structure and diversity shifts through time. Significance and percentage of variance explained by all models tested were summarized in Table 1a. Jaccard, Unifrac and Weighted-Unifrac dbRDAs supported significantly (p-value < 0.01) a link between environmental parameters included in the model and our community data (Fig. S2.2). Models using Jaccard and Unifrac dissimilarities explained between 63 and 66% of the total variance, respectively. Only the first canonical axis was significant in both models, with a higher proportion of the variance in the dissimilarity matrix explained when phylogeny was considered using Unifrac (33% vs 23% for Jaccard, Fig. S2.2a-b). In both models, the same four environmental dynamics were significant, ProxyDyn1, ProxyDyn2, ProxyDyn3, and Diuron (Table 1b). W-U dbRDA model performed best, explaining 78% of variance in the OTU matrix. The first axis was significant and explained 44% of the total variance. Four environmental dynamics were significant according to this model, ProxyDyn1, ProxyDyn2, 2,4-Dichlorophenoxyacetic acid (2,4D), and Diuron (Table 1b, Fig. S2.2c). Finally, only when the raw matrix was reduced to OTUs ≥ 0.05% of total read number (see section 2.4.1 for details) tbRDA became significant (p-value = 0.046, DCA first axis length < 3), with one significant axis and three significant dynamics, ProxyDyn1 (p-value = 0.032), ProxyDyn2 (p-value = 0.007) and Diuron (p-value = 0.001) (Fig. 5). This matrix included 53 OTUs, and 284,506 reads, i.e. 3% of total OTUs, representing 69% of the total reads) which were conserved for further analyses. This model explained 66.32% of the total variance, and the first axis was significant and explained 37.6% of the variance. Samples well projected to ProxyDyn2 and Diuron were t23 and t37, followed by t17 and t19, and to a lesser extent t38 (see Fig. 5), which was well projected to ProxyDyn1. Notice however that axis 2 was not significant.
<figure><img src="image_8.png" /></figure>

**Fig. 5** Redundancy analysis (RDA) biplot with scaling by sites on the normalized matrix of OTUs with an abundance ≥ 0.05%. The model explained 66.32% of the variance (p<0.05). Significance for axes and environmental dynamics after permanova analyses are indicated, p-value significance codes: ***<0.001<**<0.01<*<0.05. Sample names are followed by a number that indicates the sampling time in hours after the beginning of the flood at t0. For further details on sample names and retained environmental variable dynamics, see Fig. 1b. Perpendicular grey lines represent the projection of the corresponding samples onto the corresponding dynamics and approximate the value of that sample along the variable (Legendre and Legendre 2012).

**3.3. Constrained molecular ecological network analyses by environmental parameters dynamics**

**3.3.1. Archaeal network**

Network analysis using MENAP resulted in four modules containing all OTUs and 774 links. Module eigengene analysis allowed correlation of modules with environmental variables. Only one module was significantly correlated to the environmental dynamics of the flood, particularly with ProxyDyn2 and Diuron (Fig. 1b). This module included 28 OTUs of which 17 had a significant module membership (Fig. 6a), comprised of 11 with positive (52,149 reads) and six with negative (18,573 reads) correlations (Table S3.1). Worth noting is that OTUs with positive module membership (Fig. 6b) had positive interactions with each other (Fig. 6a) and negative interactions with OTUs with negative module membership (Fig. 6c), and OTUs with negative module membership had positive interactions with each other and negative interactions with OTUs with positive module membership. All positively correlated OTUs were abundant in samples t23 (20,246 reads, which represent 39% of reads along the flood
from these OTUs, Table S3.1), t37 (20%), t19 (16%) and t17 (13%). A total of 60% of reads from these OTUs belonged to the Methanobacteria class (Euryarchaeota phylum). Two OTUs were particularly abundant, OTU3, matching *Methanobrevibacter (Mbr.) smithii* at 99% similarity after blastn search against NCBI database, and OTU11 matching *Methanobacterium (Mba.) palustre* (100%). The next most abundant OTUs also matched methane-related taxa: OTU41, which matched *Methanosaeta (Msa.) concilii*, 100%, from Halobacterota phylum), OTU19 (*Methanogranum* sp. 98.04%. from Thermoplamatota phylum), OTU26 (*Mbr. acidurans*, 99%, from Euryarchaeota phylum) and OTU16 (Methanomethylophilacea, 99.61%, from Thermoplamatota). The other positively correlated OTUs (not highlighted in bold in Table S3.1) did not match any further than the family level after a blast search of the NCBI database. Finally, negatively correlated OTUs were abundant in samples from the end of the flood (t44-t109) and belonged to Nitrososphaeria (Crenarchaeota phylum, representing 85% of the sequences of these OTUs) and Thermoplasmata (15%) classes (Thermoplasmatota phylum).

<figure><img src="image_9.png" /><figcaption>**Fig. 6** (a) Molecular ecological network of the unique module significantly positively correlated with environmental dynamics, particularly ProxyDyn2 and Diuron (see Fig. 1b for details and sample names). (b) Histogram of total reads of OTUs with positive module membership in function of samples along the flood. (c) Same for OTUs with negative module membership. Environmental dynamics are followed by module correlation value between parenthesis and the p-value significant code as follows: `***<0.001<**<0.01<*<0.05`. OTUs are colored according to their phylum. The positive and negative connectivities between OTUs are indicated by blue and red lines, respectively. Only OTUs with a significantly correlated abundance profile with module are represented in this figure.</figcaption></figure>
418

### 3.3.2. Archaeal and bacterial joining network

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

To further explore the microbial relationships with environmental parameters, we performed a second eigengene analyses by first constructing a network using the 53 retained archaeal OTUs from this study, and the 260 retained bacterial OTUs from the previous study using the same samples and methodology (Noyer et al. 2020). The network we obtained before considering environmental parameters consisted in 200 nodes (OTUs) and 967 links and was composed of 21 modules. When considering environmental parameters, the network was composed of 136 OTUs, 442 links and 13 modules (Fig. 7, Table S3.2). All retained environmental dynamics (Fig. 1b) were positively or negatively correlated with one or more modules and there were 127 bacterial and 9 archaeal significant OTUs. More than half of the total reads (53%) were represented in module 1 (Fig. 7). Archaeal OTUs were present in five modules (Fig. 7). Six archaeal OTUs represented 41% of the total reads among significant OTUs in the network and were linked to ProxyDyn2 and/or Diuron. Notably, OTU3 and OTU41 in module 1, OTU16, OTU26, and OTU61 (*Mba. acidurans*, 98.42%) in module 2, and OTU47 (*Mba. formicicum*, 100%) in module 5 (Table S3.2). They were mainly present in t17, t19, t23, and t37 samples. Two archaeal OTUs, 94 and 44, in module 4 correlated to $NO_3^-$ (4% of the total reads, Fig. 7), and belonged to Nitrosotaleaceae family (Crenarchaeota phylum) and were mainly present towards the end of the flood event. Of the 127 bacterial OTUs present in this joined network, 49 were present in the bacterial only network from the previous study (Noyer et al. 2020, in purple in Table S3.2) and were significantly linked to the same environmental dynamics. Worth noting is that the bacteria in modules which were positively correlated with ProxyDyn2 and Diuron were also the most relatively abundant in t17, t19, t23, and t37 samples. The most abundant of these bacterial OTUs was OTU5, which matched *Arcobacter cryaerophilus* and represented 20% of the total reads included in the network.
<figure><img src="image_10.png" /></figure>

Fig. 7 Summary of joined bacterial and archaeal network analysis. Each colored rectangle represents a module whose number appears at the bottom left corner of each rectangle together with the percentage of reads in the module out of the total number of reads analyzed within the network. The number of OTUs and the proportion of reads within each module for each domain: archaea (A) and bacteria (B) are also indicated in each rectangle. At the top left of each rectangle the environmental dynamics are indicated with module eigengene correlation value between parentheses and the p-value significant code as follows: ***<0.001<**<0.01<*<0.05. The rectangle on the bottom right corner represents the six modules whose percentage of reads is less than 1% of the total number of reads analyzed in the network.

# 4. Discussion

## 4.1. Particle-attached archaeome diversity at three different seasons was higher and more even than in all other lotic ecosystems studied so far

Archaea have been studied in only a few river environments (Herfort et al. 2009; Hu et al. 2016, 2018; Cannon et al. 2017; Samson et al. 2019; Lei et al. 2020; Cao et al. 2020; Pinto et al. 2020; Shen et al. 2021) and to the best of our knowledge, no study has characterized their diversity over different seasons in lotic ecosystems. Here, eDNA was extracted from riverine water samples from summer, autumn, and winter and characterized using high-throughput Illumina sequencing of the archaeal rDNA in a coastal Mediterranean watercourse, the Têt River (Southeast France). Shannon alpha diversity (Fig. 2) from summer, winter and autumn drought samples were not significantly different (SD = 6.22; t0 = 6.25 and WD = 7.39). This was also the case for all alpha diversity indices (Table S1). Lei et al. (2020) obtained lower values of Shannon index (ranging 4.07-5.72) through Illumina sequencing of a highly polluted
river in China. Along a river in India sampled during summer and sequenced by metagenomic analyses, Samson *et al.* (2019) found a Shannon diversity varying from 3.12 to 4.79, depending on the sampling site. Hu *et al.* (2018) studied microbial diversity in a high-elevation river in China using universal primers and found even lower Shannon values for archaea along the river (2-2.7) and a Pielou evenness index very high (0.77-0.82) compared to ours (ranging 0.44-0.51 in Têt drought samples along seasons). These authors used 97% identity to define OTUs, while the other two studies did not specify how they defined the OTUs in their studies. Nevertheless, given the lower number of OTUs retained in these studies, the differences in alpha diversity with respect to our study could be due to less throughput sequencing. Borrego *et al.* (2020) studied particle-attached archaea in a Mediterranean water reservoir in summer with the same specific primers as in the present study, and they also found a lower Shannon alpha diversity than ours (4.9 ± 0.15). Gobet *et al.* (2014), in a comparative study of different alpha diversity indices between ARISA versus pyrosequencing methodologies, determined significant correlation results among methods, particularly with Shannon index, without a need for correction for sequencing depth. In summary, although the Têt River archaeal diversity was higher and less even (lower Pielou index) when compared to other riverine environments studied so far, additional studies with higher throughput and standardized sequencing pipeline analyses are needed to better understand the ecology of fluviatile archaea.

Studies on lotic ecosystems have revealed the presence of three major phyla, Crenarchaeota, Euryarchaeota and Thaumarchaeota using the previous silva database (Hu et al. 2016, 2018; Samson et al. 2019; Lei et al. 2020). When considering archaeal taxonomy changes since the last silva database release, the major phyla observed in the Têt river corresponded to those same phyla, but we also found the Nanoarchaeota phylum (Fig. 3a). In a recent study, which used the same primers as those used in our study, Woesarchaeia class (now called Nanoarchaeia) was found in a small water reservoir sampled in summer (Borrego et al. 2020). In our study, the summer drought community showed, through beta diversity analyses, the greatest difference with respect to communities present in autumn and winter and was also dominated by this class of the Nanoarchaeota phylum, followed by Bathyarchaeia class (Crenarchaeota phylum). Primer bias could thus be responsible for the absence of this class in previous riverine studies. Nevertheless, winter drought and dry autumn weather (t0) samples had a smaller proportion of these classes at the expense of Nitrososphaeria class (Crenarchaeota phylum, Fig. 3). On the other hand, Herfort *et al.* (2009) sequenced five DGGE fragments amplified from surface waters of the Têt River in June 2005 that matched two classes, a marine benthic group Crenarchaeota and LDS/RCV from Halobacteriales order
499

(currently at Halobacteria class from Halobacterota phylum). In this study, we also found taxa
from the genus *Methanomicrobia*, which were present in all samples independently of season,
500
501 and are classified within Halobacterota and thus closely related to Halobacteria. While we also
502 found taxa from the marine benthic group A from Crenarchaeota phylum (not shown), they
503 were absent in the summer drought. Given the technical limitations of the DGGE fingerprinting
504 method, a biased result cannot be discarded.

500

501

502

503

504

505

## 4.2. Particle-attached archaeome shifts gave evidence of allochthonous inputs into the Têt River along a Mediterranean extreme flood

506

507

In this manuscript, we characterized the temporal shift of the particle-attached riverine archaeal
508
509 community throughout an autumn storm event that included a 5-year flood and led to an
instantaneous peak discharge of 270 m³/s. Cannon *et al.* (2017) is the only study that
510 characterized river archaeal composition in water samples from a contaminated urban river
511 before and after transient rainfall using Illumina sequencing of the rDNA. They found a weak
512 resolution in archaeal sequences of samples from before and after the rain event, which they
513 concluded might be due to a primer bias. In this study, we demonstrated major changes in
514 archaeal taxonomy, and therefore community shifts observed through beta diversity, occurred
515 at specific moments along the flood event and not when comparing samples from before and
516 after the flow peak. It is therefore possible that the absence of differences in Cannon *et al.*
517 (2017) was due to the resilience of archaeal communities at the end of the rainfall as was the
518 case at the Têt River. Communities from the end of the flood, from t44 to t109 were indeed
519 more similar to the t0 sample (sampled before the rainfall started, see Fig. 4). A first shift in
520 the structural diversity occurred after t17, particularly in samples t23 and t37. These samples
521 corresponded to moments of multipollution peaks that occurred at storm first flushes CSOs and
522 inputs from in-sewer sediment resuspension (Fig. 1b, see also Reoyo-Prats et al. 2017).
523 Methane-related classes such as Methanobacteria, Methanomicrobia and Methanosarcinia
524 became significantly abundant at these moments, and our results support these classes as major
525 components of riverine communities related to pollutants, organic matter inputs and/or hypoxic
526 sediments, as reported by Casamayor and Borrego (2009) with respect to the Euryarchaeota
527 phylum. All these classes belonged indeed to the Euryarchaeota before the last silva database
528 release. Lei *et al.* (2020) also found a predominance of different Euryarchaeota methanogens
529 in water samples from a black odor, highly polluted Chinese river and more particularly the
530 large presence of *Methanobacterium* genus (Methanobacteria class) in water as well as in hotter

508

509

510

511

512

513

514

515

516

517

518

519

520

521

522

523

524

525

526

527

528

529

530
and lower-oxygen, downstream sediments, what corroborates our hypothesis on the origin of this river community shift. A second major community shift during the storm event occurred at t41, just after the maximum water discharges (Fig. 3). Compared to Lei *et al.* (2020), who did not observe a significant variation in alpha diversity along vertical and horizontal river samples, we noticed a significant decrease in the archaea alpha diversity at t41 (Table S1). This can be explained by the presence of three OTUs representing 90% of all reads in this sample (not shown), which belonged to Nitrosopumilaceae and Nitrososphaeraceae families (Nitrososphaeria class). Nitrososphaeria is the only known ammonia-oxidizing archaea (AOA) class which accomplishes nitrification in all kinds of environments (Pinto et al. 2020) but has been specifically found attached to terrestrial soil and freshwater sediments (Sonthiphand et al. 2013; Li et al. 2018). The AOA predominated assemblage at t41 is therefore derived from soil runoff or sediments from upstream environments or from the resuspension of deep river sediments, as they became predominant over polluted in-sewer sediments resuspension and urban runoff-related assemblages that predominated until t38. The last community shift included a significant enrichment of Thermoplamata class in samples t44 to t109, that coincided with the start of the second peak of flow, which occurred due to the release of waters from the upstream Vinça reservoir (Reoyo-Prats et al. 2017). This class has been associated with anoxic, sulphide-rich lentic sediments (Fillol et al. 2015; Compte-Port et al. 2017). The release of water from this reservoir during regular floods is performed via a valve situated at the bottom of the dam. As the bottom reservoir waters are anoxic (Fovet et al. 2020) and flow increase at the bottom level of the dam could potentially lead to sediment release as well, both are potential causes for the increase of Themoplasmata when dam discharges became predominant at t44.

### 4.3. Major environmental forces were linked to particle-attached archaeome shifts during an extreme event through comprehensive modelling analyses

Among the few papers in the literature that have linked archaeal diversity to environmental parameters, all addressed the effect of one family of parameters on changes in archaeal communities, such as nutrients (Herfort et al. 2009; Hu et al. 2016, 2018; Lei et al. 2020; Cao et al. 2020) or metals (Mahamoud Ahmed et al. 2020; Shen et al. 2021). This is, therefore, the first study on the Archaea domain where both diversity and a large panoply of physicochemical parameters (notably nutrients, trace metals, pharmaceuticals, and pesticides) have been analyzed for the same samples. To analyze such complex datasets, we used a comprehensive
analysis as such performed on bacteria (Noyer et al. 2020). On one side, Bokulich et al. (2012) suggested quality-filtering strategies to eliminate artifacts before interpretation of results and recommended a conservative OTU threshold of 0.005% for bacteria. On the other side, it has been suggested to use multivariate cut-offs instead of arbitrary thresholds to delineate abundant versus rare OTUs, as the latest are largely dependent on sequence coverage (Jia et al. 2018). Furthermore, multivariate cut-offs can be set considering environmental parameters when available (Gobet et al. 2010). Based on these studies, we designed a strategy with the purpose to define this cut-off at the point where the constraint-based model which considers OTU relative abundance (tbRDA model) became significant. The environmental response of the archaeal community was obtained with a tenfold higher threshold (0.05%) than bacteria (0.005%, Noyer et al. 2020). This result indicated that archaeal taxa that responded to environmental pollutants were ten times more abundant than bacteria. Less abundant archaeal taxa were therefore less crucial than those same bacterial taxa in the response of their respective communities to environmental parameters. We believe that was the reason why the constraint-based model obtained when using Bray-Curtis dissimilarity i.e., that considers OTU abundance only, turned out non-significant for the archaeome (Table 1a) in contrast to the bacteriome. The other constrained multivariate analyses models that were significant, as well as the constrained network analyses, gave evidence of two dynamics, ProxyDyn2 and Diuron, as responsible for major structural diversity shifts observed on riverine particle-attached archaea during the extreme flood event (Fig. 5 and 6b). These specific dynamics were discharged in the watershed by CSOs as well as by urban runoff (Reoyo-Prats et al. 2017) from t17 to t38 and represented several environmental pollutants including pesticides, as well as copper and dissolved pharmaceutical products and a contaminant, phosphate (see Fig. 1b for further details). These parameters have been found to affect bacterial communities in freshwater ecosystems, and metals and nutrients also affect archaeal communities in freshwaters (Hu et al. 2016; Lei et al. 2020; Shen et al. 2021). Nevertheless, to our knowledge, the effect of xenobiotics such as pesticides or pharmaceutical products on archaeal communities from freshwater ecosystems has not yet been specifically reported. But, the pesticide glyphosate, for instance, interferes with the aromatic-acids pathway in microorganisms, including archaea, and alters microbial communities (van Bruggen et al. 2021). On the other hand, eutrophication by phosphate has been found to decrease glyphosate degradation by biofilms and increase AMPA accumulation in surface waters (Carles et al. 2019), what could therefore indirectly affect microbial communities. The third significant dynamics according to the tbRDA, ProxyDyn1, projected on axis 2, was mainly linked to the flow peak sample at t38, but axis 2 turned out not significant
597

(Fig. 5). ProxyDyn1 was neither correlated to the significant module of network analyses (Fig. 6).

598

599

## 4.4. Key players in the response of the riverine archaeome to multiple stressors

600

601

602

603

604

605

606

607

608

609

610

611

612

613

614

615

616

617

618

619

620

621

622

623

624

625

626

627

628

One of the major objectives of this study was to identify OTUs which could play a key role in the archaeal community response to environmental changes derived from the delivery of xenobiotics into the river by different pollutant sources. OTUs with positive membership to the significant module from constraint-based network analysis were affiliated with Euryarchaeota, Halobacterota and Thermoplasmatota (Fig. 6a, Table S3.1). These major OTUs played a significant role in the response of the archaeome to multiple stressors derived from point sources of pollution, as they were most abundant in samples t17-t23, t37, and t38 (Fig. 6b, Table S3.1) and were linked to ProxyDyn2 and Diuron but not to ProxyDyn1, which is a proxy for diffuse sources of pollution. The predominant OTU11, affiliated with *Mba. palustre*, is a species isolated for the first time in peat bogs that has the ability to use secondary alcohols to produce methane (Zellner et al. 1988; Chaban et al. 2006). *Mbt. smithii* (OTU3) was the second most abundant key player, which is known to be derived from the human gastrointestinal tract (Miller et al. 1982; Oliveira et al. 2016) and identified as a potential indicator of sewage (Johnston et al. 2010; McLellan and Eren 2014). The third most abundant OTU41 was affiliated with *Msa. concilii*, which is involved in methane production from acetate (Zwain et al. 2017) and is abundant in-sewer biofilms (Sun et al. 2014). What was striking in the present study was that other OTUs, which had not yet been linked to CSOs and/or pollutant inputs in the literature, acted as major OTUs in the response of the archaeome to pollutant mixtures (Fig. 6a, Table S3.1). Three of them matched Methanobacteria (OTU26) or Thermoplasmata (OTUs 16 and 19) classes and were also related to methane; and the other three major OTUs were affiliated to Nanoarchaeia and Bathyarchaeia classes (OTUs 39, 66 and 34). Lastly, several OTUs affiliated to Thermoplasmata and Nitrososphaeria classes, had a significant negative membership to the module of network analysis linked to ProxyDyn2 and Diuron (Fig. 6c). These key players were instead linked to samples from the end of the flood event, from t44 to t109. As evidenced above, they were taxa from allochthonous origins. All these urban taxa, which were identified as either positively or negatively correlated with major stressors, could be used as alternative bioindicators for rapid risk assessment of the impact of multiple stressors on aquatic ecosystems, as proposed for bacteria (McLellan and Eren 2014; Dila et al. 2018; Noyer et al. 2020).
## 4.5. Comparison between particle-attached archaeal and bacterial diversity at seasons and along the flood event

Few studies so far have compared archaeal and bacterial structural diversity in lotic ecosystems (Cannon et al. 2017; Hu et al. 2018) and, to the best of our knowledge, the differences in the responses of both domains to seasons have not yet been studied. Fluviatile archaeal response to environmental dynamics differed from the bacterial response in three ways. First, the archaeome showed a specific community in summer compared to winter and autumn samples, while bacterial communities from these three seasons clustered together when compared to communities along the flood (Noyer et al. 2020). Second, environmental forces structured the archaeome diversity differently depending on beta diversity distance-based statistical models, particularly when considering qualitative vs. quantitative dissimilarities. These models were not always significant (Table 1), while for bacteria all models were significant. Third, while the bacterial community structured differently at the two multipollution events of this extreme flood, which were derived from CSOs and the flow peak (Reoyo-Prats et al. 2017), archaeal structural shifts could only be interpreted after further statistical modeling of OTU abundances and extreme event dynamics. On the one hand, the archaeal riverine resident communities shifted significantly according to two environmental dynamics only, ProxyDyn2 and Diuron (Fig. 5 and 6). These dynamics were derived from point sources of pollution, such as CSOs, in-sewer sediments resuspension and urban runoff conducted through CSOs. On the other hand, bacteria also responded significantly to other environmental parameters linked to diffuse sources of pollution. Notably, at the highest water discharge at t38 and t41, bacteria showed a specific community that was correlated to ProxyDyn1 and was thus attributed to allochthonous inputs from runoff (Reoyo-Prats et al. 2017). The archaeome shifted at t41 into an anoxic community, most likely related to the resuspension of deep sediments, and was not significantly related to ProxyDyn1 (Fig. 3 and 5) and therefore the flow peak. In conclusion, archaea can help to better understand the origin of watershed sediments and can thus be of interest for quality assessments of suspended matter.

To further understand the different responses of archaea and bacteria domains of life to river contaminants and hydrodynamics, we used a network analysis combined with a module eigengene analysis of bacterial and archaeal taxa best fitted to significant beta diversity models. One remarkable finding was the joint presence of archaea *Mbr. smithii* and *Msa. concilii* with bacteria *Arcobacter cryaerophilus*, *Bacteroides graminisolvens*, *Cloacibacterium normanense* and *Macellibacteroides fermentans* in Module 1, which was significantly correlated with
ProxyDyn2 and Diuron. All of these species have already been identified as potential indicators of human-specific sewage pollution (Dick and Field 2004; Johnston et al. 2010; McLellan and Eren 2014; McLellan and Roguet 2019) and *A. cryaerophilus* is a known human pathogen (Collado et al. 2010). While no archaea has yet shown pathogenic effects on humans (Cavicchioli et al. 2003; Bang and Schmitz 2018), archaea and bacteria can share genes through horizontal gene transfer, which is particularly enhanced in anaerobic environments (Fuchsman et al. 2017). Biofilms that develop in urban pipes are anaerobic (Guisasola et al. 2008), and the co-habitation of archaea with bacterial pathogens in urban systems can increase the risk of antibiotic resistance gene transfer. Sewers have indeed been identified as reservoirs of antibiotic-resistance bacteria carried by human pathogens (Millar and Raghavan 2017; Auguet et al. 2017). Similarly, resistance to pollutants could be enhanced between both domains through the same mechanisms. Furthermore, the presence of both domains in the same habitat pointed out to potentially common as well as complementary metabolic and physiological functions. Another interesting finding that emerged from this analysis was the presence of archaea together with bacteria OTUs in module 4, which significantly correlated with NO₃⁻ (Fig. 7, Table S3.2). Notably, two Nitrososphaeria class OTUs, known to be dominant in particle- attached ammonia-oxidizing archaeal communities (Cai et al. 2019; Pinto et al. 2020) were linked to module 4. Wang *et al.* (2018) also demonstrated a significant influence of dissolved inorganic nitrogen on the composition of bacterial and archaeal communities along an urban river. In the present study, the addition of bacterial OTUs in the network analysis strengthened the importance of nitrates to drive shifts in archaeal assemblages, otherwise mitigated when archaea were modeled alone, which enhance the interest of studying different domains of life to better understand environmental drivers of community structure in natural ecosystems.

## 5. Conclusion

Our study provides the first overview of archaeal community shifts along an extreme storm event that led to multiple pollutions in a typical coastal Mediterranean watercourse. Shifts were also compared with changes in archaeal diversity during three seasons. Archaea from the Têt river showed a specific community in summer compared to winter and autumn samples, and a higher alpha diversity and lower evenness could be observed in this river compared to other, yet less-thorough studies on riverine archaeal communities. Further studies on the spatio-temporal shifts in archaeal assemblages in these ecosystems are therefore urgently needed to
better understand seasonal shifts as well as their ecological diversity. For the first time, a comparison of the response of archaea alone and together with bacteria in a fluviatile ecosystem has shed light on the similarities and differences in their responses to seasons and when facing multiple stressors derived from an extreme event. The fluviatile bacteriome and archaeome did not respond in the same way to environmental forces. Extreme events were stronger at structuring bacterial communities than seasons, while the opposite was observed in archaeal communities. In contrast to bacteria, which responded quickly and significantly to both sewage overflow and river hydrodynamics and associated environmental parameter changes, the archaeal community shifted in response to multipollution derived from point sources and from the resuspension of deep anoxic sediments but not so clearly at the river flow peak. Archaeal taxa already known as urban-specific, as well as new archaea, mainly methane-related and never identified as urban-specific taxa, predominated assemblages during multiple stress events and were confirmed through statistical modeling of archaeal alone and together with bacteria. These taxa could be used as bioindicators of point sources of pollution. Our results highlight fluviatile archaea, seldom considered as bioindicators of water quality, could provide a rapid risk assessment of multiple pollutants in aquatic ecosystems, as is the case for bacteria. Furthermore, a better understanding of parallel shifts in assemblages from both domains of life when confronted with multiple stressors could help to improve how urban watersheds are monitored and would thus be extremely helpful in the management of pollution risk.

713

## Data Availability

Sequencing data are deposited on NCBI under BioProject ID PRJNA602803.

716

## Supplementary Information

Additional Supporting Information may be found in the online version of this article at the publisher's website.

720

### Supplementary Information S1

**Table S1.** Alpha diversity indices at the Têt river with Kruskal-Wallis test results and graphical representation for each index.

724

### Supplementary Information S2

**Fig. S2.1.** Beta diversity based on additional dissimilarities.

**Fig. S2.2.** Jaccard, Unifrac and Weighted-Unifrac distance-based RDA triplot.
728 Supplementary Information S3

729 Table S3.1. Key player significant archaeal OTUs in module eigengene analysis significantly
730 related to environmental dynamics and graphical representation of their total reads along time.

731 Table S3.2. Key player significant bacterial or archaea OTUs in modules derived from module
732 eigengene analysis significantly related to environmental dynamics.

733

734 Statements & Declarations

735 Funding

736 This work was supported by scholarships from Ecole Doctorale Energie et Environnement (E2-

737 UPVD) and from Region Occitanie to MN (n. 2017-19-ED.305). Supporting projects were:

738 DEBiMicro (2013 BQR CEFREM to CP), StepBiodiv (2015-2017 VEOLIA-Eau Perpignan to

739 OV); and DEBi2Micro (2016-17 EC2CO CNRS INSU to CP).

740

741 Competing Interests

742 The authors have no relevant financial or non-financial interests to disclose.

743

744 Contributions

745 Megane Noyer and Carmen Palacios contributed to conception and design of the study,

746 organized the database, performed the statistical analyses and wrote the first draft of the

747 manuscript. Maria Bernard helped with metabarcoding analyses. All authors contributed to

748 manuscript revision, read and approved the submitted version.

749

750 Ethical Approval

751 Not applicable.

752

753 Consent to Participate

754 All authors are informed and agree to the study.

755

756 Consent to Publish

757 The authors declare no competing interests.
758 **References**

759 Adam PS, Borrel G, Brochier-Armanet C, Gribaldo S (2017) The growing tree of Archaea: new perspectives on their diversity, evolution and ecology. The ISME Journal 11:2407–2425. https://doi.org/10.1038/ismej.2017.122

762 Amalfitano S, Corno G, Eckert E, et al (2017) Tracing particulate matter and associated microorganisms in freshwaters. Hydrobiologia 800:145–154. https://doi.org/10.1007/s10750-017-3260-x

765 Ashley RM, Wotherspoon DJJ, Coghlan BP, McGregor I (1992) The Erosion and Movement of Sediments and Associated Pollutants in Combined Sewers. Water Science and Technology 25:101–114. https://doi.org/10.2166/wst.1992.0184

768 Auguet J-C, Barberan A, Casamayor EO (2010) Global ecological patterns in uncultured Archaea. The ISME Journal 4:182–190. https://doi.org/10.1038/ismej.2009.109

770 Auguet O, Pijuan M, Borrego CM, et al (2017) Sewers as potential reservoirs of antibiotic resistance. Science of The Total Environment 605–606:1047–1054. https://doi.org/10.1016/j.scitotenv.2017.06.153

773 Bang C, Schmitz RA (2018) Archaea: forgotten players in the microbiome. Emerging Topics in Life Sciences 2:459–468. https://doi.org/10.1042/ETLS20180035

775 Blanchet J, Molinié G, Touati J (2016) Spatial analysis of trend in extreme daily rainfall in southern France. Climate Dynamics 51:799–812. https://doi.org/10.1007/s00382-016-3122-7

778 Bokulich NA, Subramanian S, Faith JJ, et al (2012) Quality-filtering vastly improves diversity estimates from Illumina amplicon sequencing. Nature Methods 10:57–59. https://doi.org/10.1038/nmeth.2276

781 Borrego C, Sabater S, Proia L (2020) Lifestyle preferences drive the structure and diversity of bacterial and archaeal communities in a small riverine reservoir. Scientific Reports 10: https://doi.org/10.1038/s41598-020-67774-0

784 Cai X, Yao L, Hu Y, et al (2019) Particle-attached microorganism oxidation of ammonia in a hypereutrophic urban river. Journal of Basic Microbiology 59:511–524. https://doi.org/10.1002/jobm.201800599

787 Cameron ES, Schmidt PJ, Tremblay BJ-M, et al (2021) Enhancing diversity analysis by repeatedly rarefying next generation sequencing data describing microbial communities. Sci Rep 11:22302. https://doi.org/10.1038/s41598-021-01636-1

790 Cannon MV, Craine J, Hester J, et al (2017) Dynamic microbial populations along the Cuyahoga River. PLOS ONE 12:e0186290. https://doi.org/10.1371/journal.pone.0186290

793 Cao X, Wang Y, Xu Y, et al (2020) Adaptive variations of sediment microbial communities and indication of fecal-associated bacteria to nutrients in a regulated urban river. Water 12:1344. https://doi.org/10.3390/w12051344
796

797

798

Carles L, Gardon H, Joseph L, et al (2019) Meta-analysis of glyphosate contamination in surface waters and dissipation by biofilms. Environment International 124:284-293. https://doi.org/10.1016/j.envint.2018.12.064

799

800

801

Casamayor EO, Borrego CM (2009) Archaea in inland waters. In Encyclopedia of Inland Waters Likens, G (ed) Oxford, UK: Academic Press, Elsevier

802

803

804

Castelle CJ, Wrighton KC, Thomas BC, et al (2015) Genomic expansion of domain archaea highlights roles for organisms from new phyla in anaerobic carbon cycling. Current Biology 25:690-701. https://doi.org/10.1016/j.cub.2015.01.014

805

806

Cavicchioli R, Curmi PMG, Saunders N, Thomas T (2003) Pathogenic archaea: do they exist? BioEssays 25:1119-1128. https://doi.org/10.1002/bies.10354

807

808

Chaban B, Ng SYM, Jarrell KF (2006) Archaeal habitats — from the extreme to the ordinary. Canadian Journal of Microbiology 52:73-116. https://doi.org/10.1139/w05-147

809

810

Clarke KR (1993) Non-parametric multivariate analyses of changes in community structure. Austral Ecology 18:117-143. https://doi.org/10.1111/j.1442-9993.1993.tb00438.x

811

812

813

Collado L, Kasimir G, Perez U, et al (2010) Occurrence and diversity of Arcobacter spp. along the Llobregat River catchment, at sewage effluents and in a drinking water treatment plant. Water Research 44:3696-3702. https://doi.org/10.1016/j.watres.2010.04.002

814

815

816

817

Compte-Port S, Subirats J, Fillol M, et al (2017) Abundance and co-distribution of widespread marine archaeal lineages in surface sediments of freshwater water bodies across the iberian peninsula. Microb Ecol 74:776-787. https://doi.org/10.1007/s00248-017-0989-8

818

819

Conseil Général des Pyrénées Orientales (CG66). Suivi de la qualité des cours d'eau du bassin versant de la Têt Année - Rapport final année 2012.

820

821

Conseil Général des Pyrénées Orientales (CG66). Suivi de la qualité des cours d'eau du bassin versant de la Têt Année - Rapport final année 2009.

822

823

824

825

826

Cowling RM, Ojeda F, Lamont BB, et al (2005) Rainfall reliability, a neglected factor in explaining convergence and divergence of plant traits in fire-prone mediterranean-climate ecosystems: Rainfall reliability in mediterranean-climate ecosystems. Global Ecology and Biogeography 14:509-519. https://doi.org/10.1111/j.1466-822X.2005.00166.x

827

828

829

Crump BC, Baross JA (2000) Archaeaplankton in the Columbia River, its estuary and the adjacent coastal ocean, USA. FEMS Microbiology Ecology 31:231-239. https://doi.org/10.1111/j.1574-6941.2000.tb00688.x

830

831

Deng Y, Jiang Y-H, Yang Y, et al (2012) Molecular ecological network analyses. BMC Bioinformatics 13:113. https://doi.org/10.1186/1471-2105-13-113

832

833

834

Dick LK, Field KG (2004) Rapid estimation of numbers of fecal Bacteroidetes by use of a quantitative PCR assay for 16S rRNA genes. Applied and Environmental Microbiology 70:5695-5697. https://doi.org/10.1128/AEM.70.9.5695-5697.2004
835

836

837

Dila DK, Corsi SR, Lenaker PL, et al (2018) Patterns of host-associated fecal indicators driven by hydrology, precipitation, and land use attributes in great lakes watersheds. Environmental Science & Technology. https://doi.org/10.1021/acs.est.8b01945

838

839

840

Dumas C, Ludwig W, Aubert D, et al (2015) Riverine transfer of anthropogenic and natural trace metals to the Gulf of Lions (NW Mediterranean Sea). Applied Geochemistry 58:14-25. https://doi.org/10.1016/j.apgeochem.2015.02.017

841

842

Escudié F, Auer L, Bernard M, et al (2018) FROGS: Find, Rapidly, OTUs with Galaxy Solution. Bioinformatics 34:1287-1294. https://doi.org/10.1093/bioinformatics/btx791

843

844

845

Faure D, Bonin P, Duran R, the Microbial Ecology EC2CO consortium (2015) Environmental microbiology as a mosaic of explored ecosystems and issues. Environ Sci Pollut Res 22:13577-13598. https://doi.org/10.1007/s11356-015-5164-5

846

847

848

849

Fillol M, Sànchez-Melsió A, Gich F, M. Borrego C (2015) Diversity of Miscellaneous Crenarchaeotic Group archaea in freshwater karstic lakes and their segregation between planktonic and sediment habitats. FEMS Microbiology Ecology 91: https://doi.org/10.1093/femsec/fiv020

850

851

852

Fovet O, Ndom M, Crave A, Pannard A (2020) Influence of dams on river water- quality signatures at event and seasonal scales: The Sélune River (France) case study. River Research and Applications 36:1267-1278

853

854

855

Fuchsman CA, Collins RE, Rocap G, Brazelton WJ (2017) Effect of the environment on horizontal gene transfer between bacteria and archaea. PeerJ 5:e3865. https://doi.org/10.7717/peerj.3865

856

857

858

859

Garcia-Esteves J, Ludwig W, Kerhervé P, et al (2007) Predicting the impact of land use on the major element and nutrient fluxes in coastal Mediterranean rivers: The case of the Têt River (Southern France). Applied Geochemistry 22:230-248. https://doi.org/10.1016/j.apgeochem.2006.09.013

860

861

862

Gobet A, Boetius A, Ramette A (2014) Ecological coherence of diversity patterns derived from classical fingerprinting and N ext G eneration S equencing techniques. Environmental Microbiology 16:2672-2681. https://doi.org/10.1111/1462-2920.12308

863

864

865

Gobet A, Quince C, Ramette A (2010) Multivariate Cutoff Level Analysis (MultiCoLA) of large community data sets. Nucleic Acids Research 38:e155-e155. https://doi.org/10.1093/nar/gkq545

866

867

Guisasola A, de Haas D, Keller J, Yuan Z (2008) Methane formation in sewer systems. Water Research 42:1421-1430. https://doi.org/10.1016/j.watres.2007.10.014

868

869

870

Herfort L, Kim J, Coolen M, et al (2009) Diversity of Archaea and detection of crenarchaeotal amoA genes in the rivers Rhine and Têt. Aquatic Microbial Ecology 55:189-201. https://doi.org/10.3354/ame01294

871

872

873

Hu A, Wang H, Li J, et al (2016) Archaeal community in a human-disturbed watershed in southeast China: diversity, distribution, and responses to environmental changes. Appl Microbiol Biotechnol 100:4685-4698. https://doi.org/10.1007/s00253-016-7318-x
874

875

876

Hu Y, Cai J, Bai C, et al (2018) Contrasting patterns of the bacterial and archaeal communities in a high-elevation river in northwestern China. Journal of Microbiology 56:104–112. https://doi.org/10.1007/s12275-018-7244-y

877

878

879

Jia X, Dini-Andreote F, Falcão Salles J (2018) Community Assembly Processes of the Microbial Rare Biosphere. Trends in Microbiology 26:738–747. https://doi.org/10.1016/j.tim.2018.02.011

880

881

882

883

Johnston C, Ufnar JA, Griffith JF, et al (2010) A real-time qPCR assay for the detection of the nifH gene of Methanobrevibacter smithii, a potential indicator of sewage pollution: A qPCR assay for detecting Methanobrevibacter smithii. Journal of Applied Microbiology 109:1946–1956. https://doi.org/10.1111/j.1365-2672.2010.04824.x

884

885

886

887

888

889

Joye SB (2012) A piece of the methane puzzle. Nature 491:538–539. https://doi.org/10.1038/nature11749

Langfelder P, Horvath S (2007) Eigengene networks for studying the relationships between co-expression modules. BMC Systems Biology 1, 54:. https://doi.org/10.1186/1752-0509-1-54

890

891

Legendre P, Gallagher ED (2001) Ecologically meaningful transformations for ordination of species data. Oecologia 129:271–280. https://doi.org/10.1007/s004420100716

892

893

894

895

896

Legendre P, Legendre L (2012) Numerical Ecology. Elsevier. ISBN: 978-0-444-53869-7

Lei M, Li Y, Zhang W, et al (2020) Identifying ecological processes driving vertical and horizontal archaeal community assemblages in a contaminated urban river. Chemosphere 245:125615. https://doi.org/10.1016/j.chemosphere.2019.125615

897

898

899

Li M, Wei G, Shi W, et al (2018) Distinct distribution patterns of ammonia-oxidizing archaea and bacteria in sediment and water column of the Yellow River estuary. Sci Rep 8:1584. https://doi.org/10.1038/s41598-018-20044-6

900

Magurran AE (2004) Measuring biological diversity. Blackwell Pub, Malden, Ma

901

902

903

904

Mahamoud Ahmed A, Tardy V, Bonnineau C, et al (2020) Changes in sediment microbial diversity following chronic copper-exposure induce community copper-tolerance without increasing sensitivity to arsenic. Journal of Hazardous Materials 391:122197. https://doi.org/10.1016/j.jhazmat.2020.122197

905

906

McLellan SL, Eren AM (2014) Discovering new indicators of fecal pollution. Trends in Microbiology 22:697–706. https://doi.org/10.1016/j.tim.2014.08.002

907

908

909

McLellan SL, Roguet A (2019) The unexpected habitat in sewer pipes for the propagation of microbial communities and their imprint on urban waters. Current Opinion in Biotechnology 57:34–41. https://doi.org/10.1016/j.copbio.2018.12.010

910

911

912

McMurdie PJ, Holmes S (2013) phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE 8(4):e61217. https://doi.org/10.1371/journal.pone.0061217
913 McMurdie PJ, Holmes S (2014) Waste not, want not: why rarefying microbiome data is 914 inadmissible. PLoS Comput Biol 10:e1003531. 915 https://doi.org/10.1371/journal.pcbi.1003531

916 Millar JA, Raghavan R (2017) Accumulation and expression of multiple antibiotic resistance 917 genes in *Arcobacter cryaerophilus* that thrives in sewage. PeerJ 5:e3269. 918 https://doi.org/10.7717/peerj.3269

919 Miller TL, Wolin MJ, Macario AJ (1982) Isolation of Methanobrevibacter smithii from Human 920 Feces. Applied and Environmental Microbiology 43:227-232

921 Newman MEJ (2004) Fast algorithm for detecting community structure in networks. Physical 922 Review E 69: https://doi.org/10.1103/PhysRevE.69.066133

923 Noyer M, Reoyo-Prats B, Aubert D, et al (2020) Particle-attached riverine bacteriome shifts in 924 a pollutant-resistant and pathogenic community during a Mediterranean extreme storm 925 event. Science of The Total Environment 139047. 926 https://doi.org/10.1016/j.scitotenv.2020.139047

927 Offre P, Spang A, Schleper C (2013) Archaea in biogeochemical cycles. Annual Review of 928 Microbiology 67:437-457. https://doi.org/10.1146/annurev-micro-092412-155614

929 Oksanen J, Blanchet FG, Friendly M, et al (2018) Vegan: community ecology package. R 930 package version 2.5-3.

931 Oliveira SS, Sorgine MHF, Bianco K, et al (2016) Detection of human fecal contamination by 932 nifH gene quantification of marine waters in the coastal beaches of Rio de Janeiro, 933 Brazil. Environ Sci Pollut Res 23:25210-25217. https://doi.org/10.1007/s11356-016- 934 7737-3

935 Osorio V, Marcé R, Pérez S, et al (2012) Occurrence and modeling of pharmaceuticals on a 936 sewage-impacted Mediterranean river and their dynamics under different hydrological 937 conditions. Science of The Total Environment 440:3-13. 938 https://doi.org/10.1016/j.scitotenv.2012.08.040

939 Oursel B, Garnier C, Zebracki M, et al (2014) Flood inputs in a Mediterranean coastal zone 940 impacted by a large urban area: Dynamic and fate of trace metals. Marine Chemistry 941 167:44-56. https://doi.org/10.1016/j.marchem.2014.08.005

942 Pinto OHB, Silva TF, Vizzotto CS, et al (2020) Genome-resolved metagenomics analysis 943 provides insights into the ecological role of Thaumarchaeota in the Amazon River and 944 its plume. BMC Microbiology 20: https://doi.org/10.1186/s12866-020-1698-x

945 Pons P, Latapy M (2005) Computing Communities in Large Networks Using Random Walks. 946 In: Yolum pInar, Güngör T, Gürgen F, Özturan C (eds) Computer and Information 947 Sciences - ISCIS 2005. Springer Berlin Heidelberg, Berlin, Heidelberg, pp 284-293

948 R Core Team (2018) R: A language and environment for statistical computing. R Foundation 949 for Statistical Computing, Vienna

950 Reoyo-Prats B, Aubert D, Menniti C, et al (2017) Multicontamination phenomena occur more 951 often than expected in Mediterranean coastal watercourses: Study case of the Têt River
990 Wang L, Zhang J, Li H, et al (2018) Shift in the microbial community composition of surface
991 water and sediment along an urban river. Science of The Total Environment 627:600-
992 612. https://doi.org/10.1016/j.scitotenv.2018.01.203

993 Woese CR, Fox GE (1977) Phylogenetic structure of the prokaryotic domain: The primary
994 kingdoms. Proceedings of the National Academy of Sciences 74:5088-5090.
995 https://doi.org/10.1073/pnas.74.11.5088

996 Woese CR, Magrum LJ, Fox GE (1978) Archaebacteria. Journal of Molecular Evolution
997 11:245-252. https://doi.org/10.1007/BF01734485

998 Zeglin LH (2015) Stream microbial diversity in response to environmental changes: review
999 and synthesis of existing research. Frontiers in Microbiology 6:454:
1000 https://doi.org/10.3389/fmicb.2015.00454

1001 Zellner G, Bleicher K, Braun E, et al (1988) Characterization of a new mesophilic, secondary
1002 alcohol-utilizing methanogen, Methanobacterium palustre spec. nov. from a peat bog.
1003 Archives of Microbiology 151:1-9. https://doi.org/10.1007/BF00444660

1004 Zinger L, Gobet A, Pommier T (2012) Two decades of describing the unseen majority of
1005 aquatic microbial diversity: Sequencing aquatic microbial diversity. Molecular Ecology
1006 21:1878-1896. https://doi.org/10.1111/j.1365-294X.2011.05362.x

1007 Zwain HM, Aziz HA, Ng WJ, Dahlan I (2017) Performance and microbial community analysis
1008 in a modified anaerobic inclining-baffled reactor treating recycled paper mill effluent.
1009 Environmental Science and Pollution Research 24:13012-13024.
1010 https://doi.org/10.1007/s11356-017-8804-0

1011
a

<figure><img src="image_11.png" /></figure>

b

<figure><img src="image_12.png" /></figure>

<table><thead><tr><th>Summer drought</th><th>Winter drought</th><th>Autumn drought</th><th>Autumn flood</th></tr></thead><tbody><tr><td>SD</td><td>WD</td><td>t0</td><td>12 sampling times: t17 → t109</td></tr></tbody></table>
<table><tr><td>Observed</td><td>5061</td><td>20357</td><td>38462</td><td>39025</td><td>37930</td><td>36234</td><td>35920</td><td>30821</td><td>22451</td><td>42347</td><td>59266</td><td>36544</td><td>40449</td><td>42965</td><td>31284</td></tr><tr><td>Shannon</td><td>6.22</td><td>6.25</td><td>7.44</td><td>7.52</td><td>6.81</td><td>7.04</td><td>6.65</td><td>6.11</td><td>5.24</td><td>7.63</td><td>7.65</td><td>6.30</td><td>7.76</td><td>8.03</td><td>7.39</td></tr></table>

<figure><img src="image_13.png" /></figure>
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /></figure>
<figure><img src="image_17.png" /><figcaption>Figure5</figcaption></figure>
<figure><img src="image_18.png" /></figure>

a

<figure><img src="image_19.png" /></figure>

Dyn2 (0.67)**

Diuron (0.8)***

b

<figure><img src="image_20.png" /></figure>

c

<figure><img src="image_21.png" /></figure>
Dyn2 (0.8)***
Diuron (0.77)**

Dyn2 (0.74)**
Diuron (0.77)**
Dyn4 (-0.66)*

Dyn3 (0.82)***

2

14.3%

3

<figure><img src="image_22.png" /></figure>

13.7%

4

<figure><img src="image_23.png" /></figure>

0.6%

14

6.6%

25

<figure><img src="image_24.png" /></figure>

1

53.2%

OTUs number

2

<figure><img src="image_25.png" /></figure>

24.7%

10

<figure><img src="image_26.png" /></figure>

28.5%

Read proportion

Dyn1 (-0.71)**
NO₃⁻ (0.64)*

12

4.6%

8

<figure><img src="image_27.png" /></figure>

Dyn1 (-0.60)*
NO₃⁻ (0.73)**

13

3.8%

13

<figure><img src="image_28.png" /></figure>

Diuron
(0.63)*

<figure><img src="image_29.png" /></figure>

2.5% 0.2%

5

2.7%

Dyn4
(-0.72)*

<figure><img src="image_30.png" /></figure>

2.1% 0.2%

3

2.3%

NO₃⁻ (0.58)*

2

<figure><img src="image_31.png" /></figure>

4.2%

2

<figure><img src="image_32.png" /></figure>

0.2%

4

4.4%

Dyn1 (0.66)*

21

<figure><img src="image_33.png" /></figure>

15

3.5%

Dyn3 (0.90)***

11

<figure><img src="image_34.png" /></figure>

1.3%

2.4D (0.97)***

10

<figure><img src="image_35.png" /></figure>

1.1%

16: 0.8%

17: 0.4%

6: 0.4%

8: 0.4%

9: 0.2%

7: 0.2%
**Fig. 1** Têt River archaeome sampling sites and environmental parameters measured in the same samples. (a) Watershed of the river with sampling site (black arrow), located after wastewater treatment plan (WWTP) of the city of Perpignan, combined sewers (black crosses) and water reservoirs indicated as grey rectangles (adapted from Reoyo-Prats et al. 2017). (b) Environmental parameter dynamics in the Têt River at different seasons and along an autumn flood. For sample names, see table below figure. Autumn sample names are followed by a number that indicates the sampling time in hours after t0, which was sampled at autumn basal level water discharges. Sampling took part at crucial moments of the flood that occurred thereafter: at first flushes (t17-t19-t23), before the flow peak (t32-t37), during the flow peak (t38-t41), following the release of water from the upstream Vinça reservoir (t44-t61) and during the return to basal level (t68-t86-t109). ProxyDyn1 corresponds to the dynamics of particulate organic carbon (20 mg/l), which represented the dynamics of water flow, also represented in figure, total suspended solids, total organic carbon, total nitrogen, and terbuthylazine parameters. ProxyDyn2 corresponds to aminomethyl phosphonic acid (AMPA, µg/l), which represented glyphosate, phosphate, copper, temperature, E. coli, enterococci, diclofenac, sulfamethoxazole and carbamazepine parameters. ProxyDyn3 corresponds to lead (150 µg/g) in the representation of the dynamics of cadmium, zinc, and conductivity parameters. ProxyDyn4 corresponds to pH (70), which represented cobalt, nickel, and chrome parameters. Three parameters, Diuron, 2,4D and NO3-, had a unique dynamic. For further details on statistical analyses for environmental parameters, see Noyer et al. (2020).

**Fig. 2** Alpha diversity of the archaeome of the Têt river along time. Changes in observed OTU number (blue) and Shannon index (red) along the flood (tX) and at summer and winter droughts (SD and WD respectively). For the Shannon index, different letters indicate a significant difference between samples (dunn.test<0.025) and red arrows show major significant differences. Observed OTU number was not significantly different (KW=0.14). The dotted profile is the flow level at each sampling point (see Fig. 1, also for sample names). Even though the absence of replicates for t61 and t68 samples impeded statistical testing, they are represented through time for comparison.

**Fig. 3** Composition of archaeal communities averaged across replicates. Histogram of relative abundances (a) of the six major phyla and (b) of the ten major classes. Samples are organized according to sampling time from left to right: summer drought (SD), autumn flood (sample names are followed by a number that indicates the sampling time in hours after the beginning of the flood at t0), and winter drought (WD). The dotted profile is the flow level at each sampling point (see Fig. 1 for further details).

**Fig. 4** Structure of archaeal communities averaged across replicates. (a) Principal Coordinate Analysis (PCoA) and (b) hierarchical clustering with Ward D2 linkage method using Weighted-Unifrac dissimilarity computed on OTU average abundance. Lines indicate ANOSIM significant groups. Sample names are followed by a number that indicates the sampling time in hours after t0 (see Fig. 1 for details). Significant codes ** and *** indicate p-value < 0.01 and < 0.001, respectively.

**Fig. 5** Redundancy analysis (RDA) biplot with scaling by sites on the normalized matrix of OTUs with an abundance ≥ 0.05%. The model explained 66.32% of the variance (p<0.05). Significance for axes and environmental dynamics after permanova analyses are indicated, p-value significance codes: ***<0.001<**<0.01<*<0.05. Sample names are followed by a number that indicates the sampling time in hours after the beginning of the flood at t0. For further details on sample names and retained environmental variable dynamics, see Fig. 1b. Perpendicular grey lines represent the projection of the corresponding samples onto the corresponding dynamics and approximate the value of that sample along the variable (Legendre and Legendre 2012).

**Fig. 6** (a) Molecular ecological network of the unique module significantly positively correlated with environmental dynamics, particularly ProxyDyn2 and Diuron (see Fig. 1b for details and sample names). (b) Histogram of total reads of OTUs with positive module membership in function of samples along the flood. (c) Same for OTUs with negative module membership. Environmental dynamics are followed by module correlation value between parenthesis and the p-value significant code as follows: ***<0.001<**<0.01<*<0.05. OTUs are colored according to their phylum. The positive and negative connectivities between OTUs are indicated by blue and red lines, respectively. Only OTUs with a significantly correlated abundance profile with module are represented in this figure.

**Fig. 7** Summary of joined bacterial and archaeal network analysis. Each colored rectangle represents a module whose number appears at the bottom left corner of each rectangle together with the percentage of reads in the module out of the total number of reads analyzed within the network. The number of OTUs and the proportion of
reads within each module for each domain: archaea (A) and bacteria (B) are also indicated in each rectangle. At the top left of each rectangle the environmental dynamics are indicated with module eigengene correlation value between parentheses and the p-value significant code as follows: ***<0.001<**<0.01<*<0.05. The rectangle on the bottom right corner represents the six modules whose percentage of reads is less than 1% of the total number of reads analyzed in the network.
<table><caption>**Table 1** Summary of constraint-based multivariate statistical models on archaea OTU matrix averaged over replicates and without singletons. (a) Permanova significance of the five models tested and the percentage of biological variance that is explained by each model using permutation test with anova.cca function. (b) Axes and modeled variables significance after permanova using anova.cca of significant models in (a). p-values significance codes: (*** < 0.001 < ** < 0.01 < * < 0.05).</caption><thead><tr><th>a</th><th>OTUs matrix<br>transformation</th><th>Model<br>significance</th><th>Variance (%)</th></tr></thead><tbody><tr><td rowspan="4">dbRDA</td><td>Jaccard</td><td>0.003**</td><td>63.38%</td></tr><tr><td>Unifrac</td><td>0.008**</td><td>66.35%</td></tr><tr><td>Bray-Curtis</td><td>0.326</td><td>60.04%</td></tr><tr><td>Weighted-Unifrac</td><td>0.004**</td><td>76.58%</td></tr><tr><td>CCA</td><td></td><td>0.765</td><td>56.44%</td></tr><tr><td>RDA</td><td>Hellinger</td><td>0.134</td><td>61.4%</td></tr></tbody></table>

<table><thead><tr><th rowspan="2">b</th><th rowspan="2">OTUs matrix<br>transformation</th><th colspan="8">Statistical significance of modeled variables</th></tr><tr><th>CAP1</th><th>ProxyDyn1</th><th>ProxyDyn2</th><th>ProxyDyn3</th><th>ProxyDyn4</th><th>2.4D</th><th>Diuron</th><th>NO<sub>3</sub><sup>-</sup></th></tr></thead><tbody><tr><td rowspan="3">dbRDA</td><td>Jaccard</td><td>0.002** (22.74)</td><td>0.001***</td><td>0.002**</td><td>0.001***</td><td>0.105</td><td>0.095</td><td>0.001***</td><td>0.096</td></tr><tr><td>Unifrac</td><td>0.040* (30.24)</td><td>0.001***</td><td>0.021*</td><td>0.046*</td><td>0.458</td><td>0.383</td><td>0.036*</td><td>0.014*</td></tr><tr><td>Weighted-Unifrac</td><td>0.006** (42.89)</td><td>0.006**</td><td>0.009**</td><td>0.039*</td><td>0.050*</td><td>0.004**</td><td>0.002**</td><td>0.166</td></tr></tbody></table>
Click here to access/download

**Supplementary Material**

Supplementary_Material_S1.xlsx
Click here to access/download

**Supplementary Material**

Supplementary_Material_S3.xlsx
Click here to access/download

**Supplementary Material**

Supplementary_Material_S2.pdf
Click here to access/download

**Supplementary Material**

ResponseReviewersComments.pdf