Microbial Ecology (2023) 86:2161–2172 https://doi.org/10.1007/s00248-023-02220-y RESEARCH A Metagenomic and Amplicon Sequencing Combined Approach Reveals the Best Primers to Study Marine Aerobic Anoxygenic Phototrophs Carlota R. Gazulla1,2 · Ana María Cabello3 · Pablo Sánchez2 · Josep M. Gasol2 · Olga Sánchez1 · Isabel Ferrera3 Received: 12 January 2023 / Accepted: 7 April 2023 / Published online: 6 May 2023 © The Author(s) 2023 Abstract Studies based on protein-coding genes are essential to describe the diversity within bacterial functional groups. In the case of aerobic anoxygenic phototrophic (AAP) bacteria, the pufM gene has been established as the genetic marker for this particular functional group, although available primers are known to have amplification biases. We review here the existing primers for pufM gene amplification, design new ones, and evaluate their phylogenetic coverage. We then use samples from contrasting marine environments to evaluate their performance. By comparing the taxonomic composition of communities retrieved with metagenomics and with different amplicon approaches, we show that the commonly used PCR primers are biased towards the Gammaproteobacteria phylum and some Alphaproteobacteria clades. The metagenomic approach, as well as the use of other combinations of the existing and newly designed primers, show that these groups are in fact less abundant than previously observed, and that a great proportion of pufM sequences are affiliated to uncultured representatives, particularly in the open ocean. Altogether, the framework developed here becomes a better alternative for future studies based on the pufM gene and, additionally, serves as a reference for primer evaluation of other functional genes. Keywords pufM gene · Primer evaluation · Metagenomics · Amplicon sequencing · AAP bacteria Introduction study them. Indeed, the sequencing of the ribosomal 16S rRNA gene allowed the first studies on the biogeography of The open ocean microbiota consists of approximately 1­ 029 marine bacterial communities [2–4]. In the last decade, the organisms that perform key biogeochemical processes development of high-throughput sequencing (HTS) meth- essential for ecosystem functioning [1]. However, only a ods together with worldwide oceanic surveys [5, 6] has small portion can be isolated and culture-independent tech- generated a massive amount of sequencing data obtained niques based on their genetic content are fundamental to with standardized methodologies, which has facilitated studying the marine microbiome at an unprecedent scale and has elucidated patterns of prokaryotic diversity, inter- * Carlota R. Gazulla actions, and connectivity around the globe (e.g., [7–13]). carlota.ruiz@uab.cat While most of the efforts have been performed at the whole * Olga Sánchez bacterioplankton community level, the focus on specific olga.sanchez@uab.cat functional groups allows the identification of microorgan- * Isabel Ferrera isms involved in a wide range of functions, such as carbon isabel.ferrera@ieo.csic.es and nitrogen fixation, ammonia oxidation, or light harvest- 1 ing, which are key to understand global biogeochemical Departament de Genètica i de Microbiologia, Universitat Autònoma de Barcelona, 08193 Bellaterra, Catalunya, Spain cycles [12, 14]. Studies based on protein-coding genes are 2 essential for this endeavor, since they may have experi- Departament de Biologia Marina i Oceanografia, Institut de Ciències del Mar, ICM-CSIC, 08003 Barcelona, Catalunya, enced horizontal gene transfer (HGT) processes [15], and Spain their phylogeny differs from the observed with the canoni- 3 Centro Oceanográfico de Málaga, Instituto Español de cal ribosomal 16S rRNA gene. Oceanografía, IEO-CSIC, 29640 Fuengirola, Málaga, Spain 13 Vol.:(0123456789) 2162 C. R. Gazulla et al. A polyphyletic group that has been extensively studied in metagenomics overcomes some of the PCR limitations, it the last two decades is that of the aerobic anoxygenic photo- generally only retrieves the most abundant members of the trophic (AAP) bacteria. Their discovery in the ocean surface bacterial community. In the case of a functional group like [16] implied a change of paradigm in our understanding of the AAP bacteria, with a relative abundance range between carbon cycling since they are heterotrophic organisms that 0.1 and 10% of the total bacterioplankton [38], metagenom- can also obtain energy from light. Although they derive a ics can limit a comprehensive knowledge of AAP diversity. fraction of their energy needs harvesting light using bac- In this context, the aim of this study is to evaluate the per- teriochlorophyll a, AAP bacteria are thought to be unable formance of existing and newly designed primers of the pufM to fix inorganic carbon, relying thus on dissolved organic gene in marine environments. We employed several combina- matter. Studies of the diversity of AAP communities are tions of existing and novel pufM primers and determined their based on the sequencing of the pufM gene that encodes the phylogenetic coverage. Then, using a selection of these prim- M subunit of the AAP bacteria reaction center. The first ers, we compared the taxonomic composition of marine AAP versions of pufM primers were designed based on sequences communities based on amplicon sequencing vs. metagenom- from cultured bacteria [17–20]. A comparison of around ics. The combination of in silico tests and phylogenetic cov- 200 sequences from cultivated bacteria and environmental erage analyses, together with its application to environmen- samples carried out by Yutin et al. [21] indicated that the tal samples, allow us to propose the optimal combination of environmental sequences of the pufM gene had a greater primers for future studies targeting the pufM gene. In addition, variability than the ones from cultured bacteria and these the approach developed here can serve as reference for future authors proposed new universal primers: pufM_uniF (for- studies involving primer evaluation of functional genes. ward) and pufM_uniR (reverse), hereafter called UniF and UniR, and an additional reverse primer called pufM_WAW. Although they were originally designed for marine environ- Methods ments, primers UniF and UniR were discarded in subsequent studies due to PCR amplification problems (e.g., [22]) and Building a pufM Database have mostly been used in freshwater ecosystems [23–25]. The combination of primers pufMF (forward), designed We built a pufM gene database containing around 1300 by Béjà et al. [20], and pufM_WAW (reverse), designed by sequences from isolates and metagenomes from marine Yutin et al. [21], was first proposed by Lehours et al. [26] on environments. For that purpose, we downloaded 697 pufM the basis of their specificity and efficiency after testing mul- sequences from the Genome Taxonomy Database (GTDB, tiple primer combinations. Since then, this combination has https://​gtdb.​ecoge​nomic.​org/, release 202) using AnnoTree been used in most studies analyzing AAP communities in [39] and then added to the database sequences from metagen- the marine environment. Most of these studies pictured AAP omics datasets such as those from the Tara Oceans Expedi- communities as mainly composed by Gammaproteobacteria tion [6], the Malaspina Expedition [5], the GOS Expedition and Alphaproteobacteria clades [26–33]. The few studies [34, 35], and the Blanes Bay Microbial Observatory (BBMO) based on metagenomics showed, however, far more diversity [28]. The taxonomic assignation of the metagenomic and a large fraction of AAP assemblages composed of mem- sequences was based on the phylogenetic tree generated in bers with no cultured representatives. For example, Yutin Gazulla et al. [27]. We classified all sequences into different et al., [34], using the Global Ocean Sampling (GOS Expedi- groups based on their class or order ranks and the phylo- tion) metagenomic shotgun data, described several groups groups A to L, previously established by Yutin et al. [34]. of AAP bacteria that were abundant in some specific areas These phylogroups were defined based on the puf operon of the ocean and that had hardly been recovered with ampli- organization and on the pufM gene phylogeny. Around 100 con-based methodologies. Some of these groups were also sequences could not be assigned to any phylogroup and were described to be abundant in samples from a Brazilian coastal clustered together in the “Others” group. bay using metagenomics [35]. A new group of AAP bacteria named “Candidatus Luxescamonaceae” (class Alphapro- Phylogenetic coverage, primer design, and in vitro teobacteria), with a putative potential for carbon fixation, performance was described from the Tara Oceans metagenomic dataset [36]. These differences among diversity surveys based on To evaluate the phylogenetic coverage of the primers pufMF metagenomics and amplicon-sequencing approaches are [20], UniF, UniR, and pufM_WAW [21], we aligned them likely due to primer biases [37]. In fact, previous discus- against our pufM database using the AlignTranslation function sions regarding possible biases in pufM amplification argued in the Decipher R package [40]. We visualized the alignment that primers pufMF and pufM_WAW may overestimate with Geneious Prime® 2021.1.1. and calculated the percent- some groups to the detriment of others [21, 26–29]. While age of sequences showing between 0 and 7 mismatches for 13 A Metagenomic and Amplicon Sequencing Combined Approach Reveals the Best Primers to Study Marine… 2163 each primer region within the whole database. To further DNA Extraction, pufM Amplification, Sequencing, evaluate these primers, we calculated in silico parameters such and Sequence Processing as the mean melting temperature, GC content, and ∆G values for hairpin, self-dimer and hetero-dimer formation, with the To analyze the performance of these primers, we selected OligoAnalyzer tool from IDT (Integrated DNA Technologies, 17 environmental samples belonging to datasets in which https://​eu.​idtdna.​com/​calc/​analy​zer). Based on the results of AAP communities had previously been analyzed: 9 samples these analyses (Table 1), we decided to design forward prim- from a seasonal study of the BBMO [28] and 8 samples ers that would combine with the reverse primer pufM_WAW, from the surface global ocean Malaspina Expedition [27] which showed good performance. The designing was car- (Table S3). For this subset of samples, both amplicon (with ried out using the DesignPrimers function of the Decipher primers pufMF/pufM_WAW) and metagenomic data were R package [40] and Geneious Prime® 2021.1.1. All primer available and they represented a comprehensive picture proposals went through the same in silico tests as the existing of the diversity at a seasonal and at a spatial scale. DNA primers, and their phylogenetic coverage was calculated as extraction and amplification of the pufM gene with prim- explained above. We came up with five candidates (Table S1) ers pufMF/pufM_WAW was done as explained in Auladell that were synthesized by ©Metabion International AG (https://​ et al. [28] and in Gazulla et al. [27] in samples from the www.​metab​ion.​com/). We performed PCR amplifications of BBMO and the Malaspina Expedition, respectively. Prim- the pufM gene using the different combinations of the newly ers pufMF_Y/pufM_WAW and UniF/UniR were used to designed primers as well as the primers from Béjà et al. [20] amplify different size fragments of the pufM gene follow- and Yutin et al., [21]. DNA from the AAP strains Congregi- ing the conditions described in Supplementary Informa- bacter litoralis (Gammaproteobacteria), Sandarakinorhabdus tion 1. Sequencing was performed in an Illumina MiSeq limnophila (Alphaproteobacteria), and from the purple sulfur sequencer (2 × 250 bp) at the Research and Testing Labora- bacteria Allochromatium vinosum (Gammaproteobacteria), all tory (http://​r tlge​nomics.​com/). Noteworthy, amplification three containing the pufM gene, was used as positive control, with primers UniF/UniR was only possible after a cleaning as well as several marine environmental samples from the step performed in the sequencing house using the TaKaRa BBMO and the Alboran Sea, which were known to contain ExTaq DNA polymerase (TaKaRa Bio Inc., Shiga, Japan). AAP bacteria. After several attempts of amplification using Sequences of the pufM gene from the BBMO metagenomes various PCR conditions (see Table S2), only the following set were generated as described in Auladell et al. [28]. Those of primers showed positive results (forward/reverse): pufMF/ from the Malaspina Expedition were retrieved from the pufM_WAW, UniF/UniR, and pufMF_Y/pufM_WAW. The Malaspina gene catalog when annotated as any of the fol- forward primer pufMF_Y was the only candidate from our lowing: pufM (prokka 1.14.6 [41]), K08928 for the Kyoto designed primer proposals that successfully amplified a gene Encyclopedia of Genes and Genomes orthologs (KEGG fragment. We confirmed the amplification of pufM fragments [42]), and PF00124 (Protein Families [43]). Annotations using Sanger sequencing, performed at the Genomics Unit were manually curated to filter out possible false positives. of the University of Málaga, Spain (www.​scai.​uma.​es). The The generation of the Malaspina Gene Database and its following material and methods refer only to these three com- annotation is described in Supplementary Information 2 and binations of primers. in Sánchez et al. [44]. Table 1  In silico parameters of primers pufMF [20], pufM_WAW, percentage refers to the percentage of sequences from our database UniF, UniR [21], and pufMF_Y (this study). Abbreviations: Tm, that hybridize with 0 mismatches (0 mismt) and 0, 1, or 2 mismatches mean melting temperature; mismt, mismatches. The hybridization (0–2 mismt). Degenerate nucleotides are underlined Primer name Sequence Length %GC Tm (°C) Hybridization (%) Hairpin* Self- Hetero-dimer * dimer* 0 mismt 0–2 mismt pufMF TAC GGS AAC CTG TWC TAC​ 18 bp 50% 62.9ºC 21.59% 69.42% − 0.3 − 5.05 − 7.94*** pufM_WAW​ AYN GCR AAC CAC CAN GCC CA 20 bp 60% 73ºC 82.50% 91.01% 1.25 − 12.32 - UniF GGN AAY YTN TWY TAY AAY CCN 26 bp 36.5% 54.2ºC 94.86% 99.9% 2.56 − 12.67 − 7.22*** TTY CA UniR YCC ATN GTC CAN CKC CAR AA 20 bp 52.5% 56.9ºC 88.50% 95.75% ** − 10.66 - pufMF_Y GGS AAY CTS TWY TAY AAY C 19 bp 42.1% 47.5ºC 55.5% 94.62% 2.3 − 10.53 − 8.22*** * Maximum value of ∆G (kcal/mole) ** No structure found for this sequence *** Hetero-dimer values are calculated for the following pairs of primers: pufMF/pufM_WAW, UniF/UniR, and pufMF_Y/pufM_WAW​ 13 2164 C. R. Gazulla et al. Sequence Data Processing and Statistical Analyses the frequency of mismatches in different AAP assemblages, we classified our database into groups according to their Each amplicon dataset was processed separately with cutadapt taxonomic rank and into the phylogroups A to L previously v1.16 [45] to remove primers and spurious sequences, and established by Yutin et al. [34] and commonly used in AAP with DADA2 v1.10 [46] to differentiate exact sequence vari- diversity surveys (Fig. 1). The forward primer pufMF had the ants and remove chimeras (details in Supplementary Informa- highest number of mismatches, especially for phylogroups A, tion 1). In total, we obtained three amplicon sequence variant B, C, D, and G, for which most of the sequences had more (ASV) tables, each one corresponding to each primer set com- than three mismatches in that primer region. In total, the bination. Sample BL110412 from the BBMO dataset was dis- pufMF primer showed perfect matches to only 21.6% of the carded due to a low number of reads (103 reads in the UniF/ sequences in our database. This primer has been commonly UniR assay). In addition, to compare the performance of the paired with the reverse primer pufM_WAW, which shows a different primer pairs, we joined the three amplicon datasets, better performance, and for which most of the sequences pre- by cutting all sequences to the same length (145 bp, the size of sent zero mismatches (80.4%). On the other hand, primers the smallest amplicon, obtained with primers UniF/UniR) and UniF and UniR showed a higher coverage for all the taxo- using the mergeSequenceTables function in DADA2, to ana- nomic groups, with a hybridization ratio (zero mismatches) lyze them as a single dataset. The phylogeny of the sequences of 94.9% for UniF and 88.5% for UniR in our database. To was inferred using the phylogenetic tree from Gazulla et al. further evaluate these primers, we calculated in silico param- [27] and the Evolutionary Placement Algorithm v0.3.5 [47]. eters such as the mean melting temperature, GC content, and Community composition, statistical analyses and figures were ∆G values for hairpin, self-dimer, and hetero-dimer formation performed in R v4.2.0 (R Core Team 2022) using packages individually for each primer (Table 1). All primers have simi- phyloseq [48], tidyverse [49], vegan [50], and ggplot2 [51] lar characteristics in terms of length, GC content, or degener- (see details in Supplementary Information 3). ate nucleotides. However, the forward primer UniF stands out as the longest oligonucleotide, with a very low percentage of GC and a high number (ten) of degenerate nucleotides. Results Design of New Primers for the Amplification Phylogenetic Coverage and PCR Parameters of the pufM Gene of the Existing Primers for the Amplification of the pufM Gene The design of new primers for the pufM gene was addressed from two perspectives: On one hand, we attempted to design We examined the phylogenetic coverage of existing primers primers for an upstream region that, combined with the pufMF [20], UniF, UniR, and pufM_WAW [21] against an reverse primer pufM_WAW, would generate longer ampli- in-house built database of the pufM gene consisting of >1300 cons and improve the taxonomic resolution. The pufM_ sequences from isolates and metagenomes, by calculating the WAW primer region is located at the end of the pufM gene number of mismatches for each position. Then, to determine (Fig. 2A), and since the phylogenetic coverage is high (Fig. 1; pufMF pufM_WAW UniF UniR pufMF_Y A B C D H Taxonomic groups L Gamma- (Pseudomonadales) K proteobacteria (Burkholderiales) I (Rhizobiales) J Rhodobacterales (Rhodobacterales) E Alpha- (Rhodobacterales) F proteobacteria (Rhodobacterales) G Sphingomonadales Other Alphaproteobacteria Others 0 2 4 6 0 2 4 6 0 2 4 6 0 2 4 6 0 2 4 6 Mismatches Fig. 1  Mismatches of primers pufMF [20], pufM_WAW, UniF, and UniR, [21], and pufMF_Y (this study) for the different AAP groups defined in this study. Dashed lines separate sequences with three or more mismatches 13 A Metagenomic and Amplicon Sequencing Combined Approach Reveals the Best Primers to Study Marine… 2165 Table 1), it was a good candidate for a reverse primer. On the primers. Primers pufMF/pufM_WAW and UniF/UniR ampli- other hand, we intended to improve the phylogenetic cover- fied fragments of 207 and 145 bp, respectively. Although it age of the primer pufMF and the in vitro performance of has been used before in marine [52] and freshwater ecosys- primer UniF, both hybridizing in the same conserved region tems [24], the combination of primers UniF/pufM_WAW did of the gene, by revising their design and including modifica- not result in positive amplification of our marine samples tions. For this purpose, we used the nucleotide and amino after several attempts under different PCR conditions. acid alignment of the sequences, we analyzed the percentage of each nucleotide at each position, and we associated these Oceanographic Context of the Environmental changes to the different phylogroups when possible (Fig. 2B). Samples Used for Primer Comparison A total of five forward primer candidates (Table S1) were and Sequencing Results tested in vitro using cultures and environmental samples and by varying the PCR conditions (annealing temperature, ­Mg+2 To test the performance of the different primers in natu- concentration, and primer concentration; Table S2). Three ral samples, we compared the composition of AAP com- primers were designed in upstream regions of the pufM gene, munities retrieved by metagenomics and that retrieved by while two were improved versions of the existing pufMF amplicon sequencing using different primer combinations. primer. We successfully obtained one forward primer, named Based on the in vitro results explained above, we selected “pufMF_Y” that combined with the existing reverse primer primers pufMF/pufM_WAW, which have been routinely pufM_WAW, amplified a 203-bp fragment of the pufM gene, used in marine environments (e.g., [26–30, 53]), primers while the others were discarded due to their poor perfor- UniF/UniR, mainly used in freshwater ecosystems [23, 25], mance. This primer hybridizes in the same conserved region and the newly designed forward primer pufMF_Y with the as primers pufMF and UniF (Fig. 2A); it improves the phylo- reverse pufM_WAW. Illumina sequencing with the three genetic coverage of primer pufMF, and it has a lower number primer combinations was performed for 17 samples that of degenerate nucleotides that primer UniF. Although the covered both spatial and temporal variabilities: eight open hybridization ratio is relatively low (55.45%) compared to ocean samples from the Malaspina Expedition [5], from the primer UniF, when we consider sequences with zero, one, or Pacific, the Atlantic, and the Indian Oceans, and nine coastal two mismatches (that are likely to amplify), the hybridiza- samples from the BBMO in the Mediterranean Sea, col- tion ratio increases up to 95%. Additionally, we performed lected at different seasons during years 2011 and 2012 [54] PCR amplification with different combinations of the existing (Fig. 3A; Table S3). These samples were part of previous Fig. 2  A Sequence logo and phylogenetic coverage of primer pufMF number of mismatches (≥ 19% of sequences) that could be associ- [20]. The table shows the percentage of sequences with each nucleo- ated with specific taxonomic groups (see color legend) are in bold tide in the different positions of the primer, based on the alignment and underlined. B Schematic representation of the pufM gene and the of our pufM database. Nucleotides representing positions with a high primers used in this study 13 2166 C. R. Gazulla et al. A 50 BBMO samples Taxonomic groups MP1857 A MP1672 0 MP2243 MP0311 B MP1517 C MP0778 MP1176 MP1421 D −50 H L −100 0 100 K (Pseudomonadales) METAGENOMIC APPROACH I (Burkholderiales) B 100 100 J (Rhizobiales) E (Rhodobacterales) 75 75 F (Rhodobacterales) Abundance Abundance G (Rhodobacterales) 50 50 Sphingomonadales Other Alphaproteobacteria 25 25 Others 0 0 AMPLICON APPROACH Primers pufMF / pufM_WAW C 100 100 100 75 75 75 Metagenomes (%) 50 50 50 25 25 25 0 0 0 0 25 50 75 100 Primers pufMF/pufM_WAW (%) Primers pufMF_Y / pufM_WAW D 100 100 100 75 75 75 Metagenomes (%) Abundance Abundance 50 50 50 25 25 25 0 0 0 0 25 50 75 100 Primers pufMF_Y/pufM_WAW (%) Primers UniF / UniR E 100 100 100 75 75 75 Metagenomes (%) 50 50 50 25 25 25 0 0 0 0 25 50 75 100 MP1857 MP1176 MP1672 MP2243 MP0311 MP0778 MP1421 MP1517 BL110208 BL110615 BL110802 BL111129 BL120214 BL120411 BL120620 BL120807 Primers UniF/UniR (%) BBMO MALASPINA 13 A Metagenomic and Amplicon Sequencing Combined Approach Reveals the Best Primers to Study Marine… 2167 ◂Fig. 3  (A) Stations from the Malaspina Expedition used in this study and “Pseudomonadales”, with sequences from the family (“MP” code). Samples from Blanes Bay Microbial Observatory Halieaceae). Besides, we also used the taxonomic groups (BBMO) are all from the same coastal site, yet collected at differ- ent times of the year. The code of Blanes samples indicates Blanes- proposed by Yutin et al. [34] –A to L– for comparison with year-month-day (e.g., BL110208 is from the 8th of February 2011). previously published studies. Some of these phylogroups (B) pufM taxonomic composition of the BBMO left) and Malaspina can be associated to known groups: phylogroup K contains (right) samples with the metagenomic approach. Below, the commu- sequences affiliated to the NOR5/OM60 clade, from family nity composition at each station retrieved with the amplicon approach and with the different primers combinations: primers pufMF/pufM_ Halieaceae, order Pseudomonadales, and phylogroup I is WAW (C), pufMF_Y/pufM_WAW (D), and UniF/UniR (E). For each related to the Burkholderiales order, both belonging to the primer combination, and each sample, we have represented the rela- Gammaproteobacteria. Phylogroup J has been associated to tive abundances of the taxonomic groups retrieved with metagenom- the Rhizobiales order, while phylogroups E, F, and G are ics vs. those retrieved with the amplicon approach. The dashed lines represent the 1:1 lines in which both the metagenomic and the ampli- associated to the Rhodobacterales order. Phylogroups A, B, con approach would indicate the same taxonomic community compo- C, D, H, and L have no taxonomically described representa- sition tive. Finally, some of our ASVs could be associated to the recently described “Candidatus Luxescamonaceae” family. Since this family clustered within phylogroups C and D in studies analyzing AAP communities by means of amplicon our phylogenetic tree, we classified them as C or D depend- sequencing with primers pufMF/pufM_WAW (Malaspina ing on their position in the tree. Sequences that could not be samples in [27] and BBMO samples in [28]). Besides, we further classified, were categorized as “Others”. obtained 176 and 62 predicted genes from metagenomic assemblies associated to the pufM gene from the Malaspina Performance of the Primers as Compared and the BBMO metagenomic datasets, respectively. Pre- to Metagenomics in Marine Environmental Samples dicted genes from Malaspina were between 101 and 1040 bp (N50 = 804 bp) in length, while genes from the BBMO were The composition of AAP communities varied largely between 708 and 1044 bp long (N50 = 966 bp). Regard- depending on the methodology (amplicon sequencing vs. ing the amplicon analyses, for the primer combination metagenomics) and the primer combination in the amplicon pufMF_Y/pufM_WAW, we obtained a total of 1904 ASVs, approaches. In the metagenomes from BBMO, phylogroup for the UniF/UniR, we retrieved 1294 ASVs, and with prim- K was always present and abundant in some months, while ers pufMF/pufM_WAW, we obtained a total of 418 ASVs. other groups peaked at specific times of the year (Fig. 3B). There were almost no shared ASVs between oceanic and This was the case of phylogroup A, whose relative abun- coastal environments (~3%; see Table S4 for details). In dance increased from < 1% up to 30% during winter, or the terms of primer efficiency, we obtained 1.2 million reads Rhodobacterales, which increased after the spring bloom. with primers pufMF_Y/pufM_WAW, 0.5 million reads with In the case of the Malaspina metagenomes, the uncultured primers UniF/UniR, and less than 0.3 million reads with phylogroups A, B, C, and D dominated in all samples, while primers pufMF/pufM_WAW. The high efficiency of prim- Gammaproteobacteria and Alphaproteobacteria sequences ers pufMF_Y/pufM_WAW and the amplicon size of 203 bp were always present but scarce (Fig. 3B). (vs. 145 bp with UniF/UniR primers) led to higher values of The different amplicon approaches resulted in communi- alpha diversity for this primer combination (Fig. S1). Com- ties with very different taxonomic composition (Fig. 3C–E; munities amplified with primers pufMF_Y/pufM_WAW Adonis test, p < 0.05), especially for primers pufMF/ had the highest richness (mean 173.5 ± 0.27) values, fol- pufM_WAW. The communities described with this primer lowed by primers UniF/UniR, with a significantly lower combination were dominated by phylogroup K (Pseudomo- mean observed diversity (125.1 ± 2.5; Tukey test, p < 0.05; nadales), and only samples BL110208, MP0311, MP1176, Fig. S1). Instead, the mean Shannon index value (pufMF_Y/ and MP1672 were dominated by ASVs associated to the pufM_WAW, 3.58; UniF/UniR, 3.36 primer) was compara- Alphaproteobacteria class (Fig. 3C). Noteworthy, we only ble for both approaches. Primers pufMF/pufM_WAW failed retrieved three sequences with very low abundances of phy- to amplify many sequences as compared to the previous logroup C (in sample MP0778) and none from phylogroups primers, which resulted in a significantly lower observed A, B or D. Communities inferred using primers pufMF_Y/ diversity (mean 43.9 ± 0.3; Tukey test, p < 0.05) and Shan- pufM_WAW were more diverse (Tukey test, p < 0.05) and non index values (mean 2.47; Tukey test, p < 0.05; Fig. S1). contained more groups (Fig. 3D). Although Pseudomon- We classified taxonomically all ASVs into 14 broad adales were still prevalent, especially in coastal samples, taxonomic groups, according to their order within the other groups appeared as relevant, such as phylogroups A Alphaproteobacteria (“Rhizobiales,” “Rhodobacterales,” and B, and Rhodobacterales from phylogroups E and G. “Sphingomonadales” and “Other Alphaproteobacteria”) Likewise, primers UniF/UniR also succeeded in amplifying and the Gammaproteobacteria classes (“Burkholderiales” phylogroups A, B, C, and D. The communities retrieved with 13 2168 C. R. Gazulla et al. these primers were similar to those described with prim- depending on the primer, the community structure was con- ers pufMF_Y/pufM_WAW, albeit with some differences. In served in the three approaches, and general ecological pat- samples from Blanes Bay, phylogroup G was predominant terns could be observed indistinctively of the used primer. within the Rhodobacterales. In turn, Malaspina samples were mainly dominated by phylogroups A, B, and C, while the relative abundance of Pseudomonadales was much lower. Discussion To further compare the performance of the three prim- ers pairs, we plotted the relative abundance of the different The study of key marker genes together with the develop- taxonomic groups retrieved with each amplicon approach ment of the “-omics” techniques has increased our under- vs. the metagenomic approach (Fig. 3C–E) and summarized standing of marine diversity and biogeochemical cycles the phylogenetic coverage of the metagenomic and ampli- [14, 55]. Both metagenomics and amplicon approaches are con approaches (Fig. S2). All in all, primer pufMF failed at commonly used to target functional genes and to describe amplifying groups with no cultured representatives (phylo- their ecological patterns. Amplicon sequencing is easy to groups A, B, C, and D), while it overestimated the abundance implement, relatively cheap, and effective in capturing of Gammaproteobacteria and Rhizobiales. In turn, primers large numbers of sequence variants. However, due to the pufMF_Y/pufM_WAW and UniF/UniR better reflected high sequence variability of protein-coding genes, primer the composition of the communities as observed with the biases are common and can result in the misrepresentation metagenomic approach, assumed to be less biased than of the relative contribution of certain taxa. In contrast, amplicon approaches. Although primers pufMF_Y/pufM_ metagenomics is PCR-free and less biased for functional WAW seemed to overestimate the abundance of Gammapro- gene analysis, but it generally retrieves fewer copies of teobacteria, especially in open ocean samples (Fig. 3D), their marker genes and hinders the comprehensive study of performance was comparable to that observed with metage- functional groups with low abundances in the environ- nomes (no statistical differences after Tukey test, p > 0.05; ment. For example, the nifH gene, a genetic marker of Fig. S2). In turn, the performance of primers UniF/UniR was nitrogen-fixing populations, is usually studied using ampli- quite good for the most abundant phylogroups. con sequencing since the number of variants retrieved with Finally, given the differences in the taxonomic composi- metagenomic surveys is too low, if not undetectable [56, tion of samples analyzed with different primers combina- 57]. Likewise, most of the studies of AAP bacteria are tions, we wanted to test whether the community structure based on the partial amplification of the pufM gene (e.g. was conserved with primers UniF/UniR and pufMF_Y/ [18, 20, 27–29, 31, 33, 58, 59],), while only few have pufM_WAW. To do so, we computed Bray–Curtis dissimi- approached their study based solely on metagenomics larity matrices for each amplicon approach to compare the [34, 35]. Yet, the application of metagenomics allows the structure of communities and performed Mantel test cor- description of new metabolisms and new taxa in marine relations. Besides, we Hellinger-transformed the data and microbiomes, such as the discovery of new nitrogen fixa- applied a Procrustes test to assess the statistical significance tion pathways in surface ocean heterotrophs [56], and between each ordination. The matrices obtained with the the description of new AAP phylogroups previously not primers traditionally employed in AAP diversity surveys and reported [34]. In this study, we combined both approaches the primer combinations proposed in this study were highly to unveil the biases of existing primers for the pufM gene, correlated (Mantel tests; Table S5). Besides, the Procrustes design new ones, and test the performance of different tests showed no statistical differences between their spatial primer combinations in a variety of marine environments. ordination (Table S5). For further comparison of primer’s Debate regarding primer biases in the case of the pufM performance, we combined all the samples into non-metric gene is not new; it arose in studies that used the primer pair multidimensional scaling (NMDS) plots, first by merging pufMF/pufM_WAW and showed a large dominance of phy- the samples from the three amplicon approaches into one logroup K (Gammaproteobacteria) in AAP communities matrix (Fig. S3a) and then including the metagenomic data [26–29]. Lehours et al. [26] reported > 80% of relative abun- (Fig. S3b). Samples from coastal and open ocean environ- dance of this phylogroup in samples from the Mediterra- ments appeared in two clear distinct clusters (dispersion nean Sea and considered possible primer biases favoring the between samples is statistically significant, betadisper test, amplification of the Gammaproteobacteria but disregarded p < 0.001; Fig. S3c) for both plots. However, the clustering that problem because the same primer pair in Arctic waters of samples based on the approach (metagenomics and dif- led to low abundances of phylogroup K. Still, it is likely that ferent primer combinations) was less clear and not signifi- arctic bacterial communities have a very different composi- cant after testing the dispersion between samples (betadis- tion, due to temperature differences and the effect of ice melt per test, p = 0.0793; Fig. S3d). All in all, this indicates that in salinity. Gammaproteobacteria were also predominant in even though the diversity and taxonomic composition varied samples from the coastal Blanes Bay [28, 29]. In this area, 13 A Metagenomic and Amplicon Sequencing Combined Approach Reveals the Best Primers to Study Marine… 2169 although the seasonal trends of communities retrieved by other regions or combinations might work, nor that designs both metagenomics and amplicon sequencing had similar of primers that include both the pufL (upstream gene that trends for the predominant groups, primers pufMF/pufM_ encodes the L subunit of the bacterial reaction center) and WAW overestimated the contribution of phylogroup K and the pufM gene might produce better tools. failed to amplify sequences of some groups that appeared Finding universally conserved regions in a functional in the metagenomic approach [28]. In addition to the preva- gene is challenging, and sometimes, it is not possible to lence of Gammaproteobacteria in marine AAP communi- generate universal primers, as it happens, among others, ties, most studies also reported the presence of members of with genes nirS and nirK (NO-forming nitrite reductase the Alphaproteobacteria [26–29, 33, 53, 60, 61]. Our results genes) that rely on the use of clade-specific primers [62]. show that for the pufMF primer, the number of mismatches Although the use of universal primers is appropriate to of sequences affiliated to Gammaproteobacteria and some describe microbial communities, even perfectly matched Alphaproteobacteria orders is low in comparison to the mis- primers can exhibit preferential amplification; thus, beyond matches in phylogroups A, B, C, and D, which are almost the in silico testing, analyses with environmental samples absent in amplicon-based studies (Fig. 1). In fact, these phy- are also important for primer evaluation [63]. To test differ- logroups were described for the first time following the shot- ent primer combinations, we used samples from different gun metagenomic sequencing approach of Yutin et al. [34]. marine environments (coastal vs. open ocean) and from dif- They have no cultured representatives, and even though they ferent seasons, to include the spatial and seasonal variability were abundant in samples from the GOS Expedition [34], that exists in AAP assemblages, as reported previously [27, they have been barely retrieved in other studies, which might 28, 31, 34]. The metagenomic assay provided a bias-free be explained by the high number of mismatches within the representation of the most abundant components of AAP region of the commonly used pufMF primer (Figs. 1 and communities in samples from the BBMO and Malaspina 2B). On the contrary, the universal primers UniF, UniR, and (Fig. 3B). Even though the number of copies retrieved was pufM_WAW [21] have a very good phylogenetic coverage low, comparing these communities to those obtained through (Fig. 1). While the reverse primer pufM_WAW has been amplicon-sequencing was the best approach to analyze the extensively used in combination with pufMF, primers UniF/ biases of each primer combination. A previous analysis with UniR have barely been used in marine environments and samples from the BBMO already identified discrepancies in were discarded in previous studies as they repeatedly failed the taxonomic composition of communities with the differ- in amplifying under different conditions ([22, 26, 29] and ent methods, such as sequences from phylogroups A, B, and authors unpublished observations). Both primers, pufMF C that were only retrieved by metagenomics and were absent and UniF, hybridize in the same region of the pufM gene in the amplicon approach [28]. In this study, using the same (with 3-nucleotides shift between them), but UniF has 10 extracted DNA from BBMO samples and different primer degenerate nucleotides, a very low GC content, and low combinations, we were able to amplify sequences from these melting temperature (Table 1), which might explain why it phylogroups, which in fact constitute over 50% of the rela- is problematic in vitro. tive abundance in two samples (Fig. 3D, E). Likewise, we The primers designed in this study aimed at both improv- obtained a great proportion of these groups in samples from ing the coverage of the commonly used ones and producing Malaspina, which were completely missed with primers longer amplicons. While we did not succeed in the design pufMF/pufM_WAW, as shown in Fig. 3C and in Gazulla of primers in the upstream region, we were able to produce et al. [27]. These results are consistent with the reports an oligonucleotide with a higher hybridization ratio than from Yutin et al., [34] and Cuadrat et al. [35], in which they pufMF, and a lower number of degenerate nucleotides than describe a high proportion of sequences affiliated to phylo- UniF, while maintaining the amplicon size of ~ 200 bp. The groups A, B, C, and D in AAP communities from different detailed analysis of the nucleotide composition of primer marine environments [36]. The alpha diversity estimates pufMF (Fig. 2B) indicates that most of the mutations in the obtained with primers pufMF_Y/pufM_WAW and UniF/ primer region can be associated to different phylogroups and UniR (Fig. S1) surpassed by far the estimates described in generally represent silent mutations. The only exceptions are samples from the BBMO [28] and the Malaspina Expedi- sequences from phylogroup C and D that have a histidine tion [27] with primers pufMF/pufM_WAW. Interestingly, and a tryptophan respectively, instead of a tyrosine in the these previous studies had provided the most comprehensive first position (of the primer region), and sequences from datasets for AAP bacteria but were clearly underestimating phylogroup D that encode for a tyrosine instead of a phe- AAP diversity. nylalanine in position five. We identified 7 nucleotides with Overall, our results indicate that the taxonomic com- problematic mismatches, represented in bold and under- position of primers pufMF/pufM_WAW is biased towards lined in Fig. 2, that were considered when redesigning the phylogroup K (Pseudomonadales) (Fig. 3C), which is over- new primer pufMF_Y. Nevertheless, we cannot discard that estimated in almost all samples. The same happens with 13 2170 C. R. Gazulla et al. Sphingomonadales-like and Rhizobiales representatives and Author Contribution I.F. conceived the study and C.R.G. contributed with a small cluster of other Alphaproteobacteria. Previous to the design. Laboratory analyses were performed by C.R.G, A.M.C, and I.F; analyses of the metagenomic data were conducted by P.S., and studies reported high abundances of Gammaproteobacte- C.R.G. performed the rest of the analyses and wrote the first draft of the ria in the Mediterranean Sea [26, 28, 29], the Baltic Sea manuscript. O.S., J.M.G, and I.F. contributed substantially to manu- [53], the North Pacific Ocean [60], the Arctic Sea [61], the script drafting and funding acquisition. All authors read and approved east coast of Australia [33], or the tropical and subtropi- the final manuscript. cal global ocean [27]. Since these studies analyzed AAP Funding Open Access Funding provided by Universitat Auton- communities with the pufMF/pufM_WAW primers, it is oma de Barcelona. This work was supported by grant ECLIPSE likely that some of these results misrepresented the real (PID2019-110128RB-I00/AEI/10.13039/501100011033) funded by composition of AAP communities, just as we have shown the Agencia Estatal de Investigación from the Spanish Ministry of Sci- ence and Innovation. Authors affiliated to the Institut de Ciències del for samples of the Malaspina Expedition or the BBMO. Yet, Mar received the institutional support of the “Severo Ochoa Centre of albeit the exposed primer biases, the community structure Excellence” accreditation (CEX2019-000928-S). of the different approaches was conserved in different ordi- nation tests and the matrices strongly correlated, so pre- Data Availability Most of the analyses were performed in R version 4.2.0 (R Core Team 2022), and code is available in the following viously published results (e.g., the adaptation of different repository: https://​gitlab.​com/​crgaz​ulla/​analy​sing-​pufm-​prime​rs-​for-​ phylogenetic clades [31], their seasonal trends [28], or the marine-​aap-​studi​es. The pufM database is available in the Supplemen- ecological processes operating in the surface global ocean tary Material. Amplicon sequences have been deposited in the NCBI [27]) should not be discarded. Sequence Read Archive (SRA) under BioProject ID PRJNA919028. Sequences from the BBMO metagenomic dataset were published in To conclude, we used an extensive pufM database to Auladell et al., [28], and those from the Malaspina Expedition dataset show the limitations of the forward primer pufMF and pro- are deposited under BioProject ID PRJEB52452. pose some alternatives to determine the composition and diversity of AAP communities in marine environments. Declarations We revised existing primers for the pufM gene in the lit- Consent for Publication All the authors give their consent to allow erature, designed new ones and selected those with the best this manuscript to be published by the journal of Microbial Ecology. performance that were tested with environmental samples and benchmarked against metagenomics. We show that Competing Interests The authors declare no competing interests. the phylogenetic coverage of primer pufMF is very low for some taxonomic groups, and, as a result, amplification Open Access This article is licensed under a Creative Commons Attri- with this primer is biased towards phylogroup K (Gam- bution 4.0 International License, which permits use, sharing, adapta- tion, distribution and reproduction in any medium or format, as long maproteobacteria) and some orders of the Alphaproteo- as you give appropriate credit to the original author(s) and the source, bacteria class. Although Gammaproteobacteria are relevant provide a link to the Creative Commons licence, and indicate if changes components of AAP communities, several species with no were made. The images or other third party material in this article are cultured representatives, like phylogroups A, B, C, and D, included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in have been entirely underrepresented in the past and are in the article's Creative Commons licence and your intended use is not fact an important fraction of AAP assemblages. For future permitted by statutory regulation or exceeds the permitted use, you will studies analyzing marine AAP bacteria, we recommend need to obtain permission directly from the copyright holder. To view a using either primers pufMF_Y/pufM_WAW or UniF/UniR, copy of this licence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/. to guarantee an unbiased representation of their taxonomic composition. Supplementary Information The online version contains supplemen- References tary material available at https://​doi.​org/​10.​1007/​s00248-​023-​02220-y. 1. Whitman WB, Coleman DC, Wiebe WJ (1998) Prokaryotes: the Acknowledgements The authors thank the Group on Ecology of unseen majority. Proc Natl Acad Sci USA 95:6578–6583. https://​ Plankton and Environmental Challenges and the Laboratory of doi.​org/​10.​1038/​s41437-​017-​0034-1 Genetics and Molecular Biology for sharing their facilities at the 2. Hagström Å, Pommier T, Rohwer F, Simu K, Stolte W, Svensson “Centro Oceanográfico de Málaga” (COMA, IEO-CSIC), where the D, Zweifel UL (2002) Use of 16S ribosomal DNA for delinea- experimental work of this manuscript was performed. We also thank tion of marine bacterioplankton species. Appl Environ Microbiol Prof. Michal Koblížek for providing a culture of Congregibacter 68:3628–3633. https://​doi.​org/​10.​1128/​AEM.​68.7.​3628-​3633.2​ 002 litoralis and Dr. Frederic Gich for doing so with Sandarakinorhabdus 3. Fuhrman JA, Hewson I, Schwalbach MS, Steele JA, Brown MV, limnophila. We would also like to thank Dr. Adrià Auladell for help Naeem S (2006) Annually reoccurring bacterial communities with preliminary sequence analyses. Computing analyses were are predictable from ocean conditions. Proc Natl Acad Sci USA performed at the Marine Bioinformatics Core facility (MARBITS, 103:13104–13109. https://​doi.​org/​10.​1073/​pnas.​06023​99103 https://​marbi​ts.​icm.​csic.​es) of the Institut de Ciències del Mar (ICM- 4. Pommier T, Canbäck L, Riemman K, Boström H, Simu K, Lund- CSIC) in Barcelona, Spain. berg P, Tunlid A, Hagström Å (2007) Global patterns of diversity 13 A Metagenomic and Amplicon Sequencing Combined Approach Reveals the Best Primers to Study Marine… 2171 and community structure in marine bacterioplankton. Mol Ecol bacteria in a freshwater lake. Environ Microbiol Rep 15:60–71. 16:867–880. https://​doi.​org/​10.​1111/j.​1365-​294X.​2006.​03189.x https://​doi.​org/​10.​1111/​1758-​2229.​13131 5. Duarte CM (2015) Seafaring in the 21st century: the Malaspina 24. Piwosz K, Villena-Alemany C, Mujakić I (2022) Photoheterotro- 2010 circumnavigation expedition. Limnol Oceanogr Bull 24:11– phy by aerobic anoxygenic bacteria modulates carbon fluxes in 14. https://​doi.​org/​10.​1002/​lob.​10008 a freshwater lake. ISME 16:1046–1054. https://​doi.​org/​10.​1038/​ 6. Karsenti E, Acinas SG, Bork P et al (2011) A holistic approach to s41396-​021-​01142-2 marine eco-systems biology. PLoS Biol 9:7–11. https://​doi.​org/​ 25. Piwosz K, Vrdoljak A, Frenken T et al (2020) Light and primary 10.​1371/​journ​al.​pbio.​10011​77 production shape bacterial activity and community composition of 7. Lima-Mendez G, Faust K, Henry N et al (2015) Determinants of aerobic anoxygenic phototrophic bacteria in a microcosm experi- community structure in the global plankton interactome. Science ment. mSphere 5:4. https://​doi.​org/​10.​1128/​mSphe​re.​00354-​20 348(6237):1262073. https://​doi.​org/​10.​1126/​scien​ce.​12620​73 26. Lehours AC, Cottrell MT, Dahan O, Kirchman DL, Jeanthon C 8. Sunagawa S, Coelho LP, Chaffron S et al (2015) Structure and (2010) Summer distribution and diversity of aerobic anoxygenic function of the global ocean microbiome. Science 348:6237. phototrophic bacteria in the Mediterranean Sea in relation to envi- https://​doi.​org/​10.​1126/​scien​ce.​12613​59 ronmental variables. FEMS Microbiol Ecol 74:397–409. https://​ 9. Mestre M, Ruiz-González C, Logares R, Duarte CM, Gasol JM, doi.​org/​10.​1111/j.​1574-​6941.​2010.​00954.x Sala MM (2018) Sinking particles promote vertical connectivity 27. Gazulla CR, Auladell A, Ruiz-González C, Junger PC, Royo- in the ocean microbiome. Proc Natl Acad Sci USA 115:E6799– Llonch M, Duarte CM, Gasol JM, Sánchez O, Ferrera I (2022) E6807. https://​doi.​org/​10.​1073/​pnas.​18024​70115 Global diversity and distribution of aerobic anoxygenic photo- 10. Logares R, Deutschmann I, Junger PC et al (2020) Disentangling trophs in the tropical and subtropical oceans. Environ Microbiol the mechanisms shaping the surface ocean microbiota. Microbi- 24:2222–2238. https://​doi.​org/​10.​1111/​1462-​2920.​15835 ome 8:1–17. https://​doi.​org/​10.​1186/​s40168-​020-​00827-8 28. Auladell A, Sánchez P, Sánchez O, Gasol JM, Ferrera I (2019) 11. Ruiz-González C, Mestre M, Estrada M et al (2020) Major imprint Long-term seasonal and interannual variability of marine aerobic of surface plankton on deep ocean prokaryotic structure and activ- anoxygenic photoheterotrophic bacteria. ISME 13:1975–1987. ity. Mol Ecol 29:1820–1838. https://​doi.​org/​10.​1111/​mec.​15454 https://​doi.​org/​10.​1038/​s41396-​019-​0401-4 12. Louca S, Parfrey LW, Doebeli M (2016) Decoupling function and 29. Ferrera I, Borrego CM, Salazar G, Gasol JM (2014) Marked sea- taxonomy in the global ocean microbiome. Science 353:1272– sonality of aerobic anoxygenic phototrophic bacteria in the coastal 1277. https://​doi.​org/​10.​1126/​scien​ce.​aaf45​07 NW Mediterranean Sea as revealed by cell abundance, pigment 13. Herndl GJ, Bayer B, Baltar F, Reinthaler T (2022) Prokaryotic life concentration and pyrosequencing of pufM gene. Environ Micro- in the deep ocean’s water column. Ann Rev Mar Sci 15:1. https://​ biol 16:2953–2965. https://​doi.​org/​10.​1111/​1462-​2920.​12278 doi.​org/​10.​1146/​annur​ev-​marine-​032122-​115655 30. Jeanthon C, Boeuf D, Dahan O, Le Gall F, Garczarek L, Bendif 14. Ferrera I, Sebastián M, Acinas SG, Gasol JM (2015) Prokaryotic EM, Lehours AC (2011) Diversity of cultivated and metaboli- functional gene diversity in the sunlit ocean: stumbling in the cally active aerobic anoxygenic phototrophic bacteria along an dark. Curr Opin Microbiol 25:33–39. https://​doi.​org/​10.​1016/j.​ oligotrophic gradient in the Mediterranean Sea. Biogeosciences mib.​2015.​03.​007 8:1955–1970. https://​doi.​org/​10.​5194/​bg-8-​1955-​2011 15. Wiedenbeck J, Cohan FM (2011) Origins of bacterial diversity 31. Lehours AC, Enault F, Boeuf D, Jeanthon C (2018) Biogeo- through horizontal genetic transfer and adaptation to new ecologi- graphic patterns of aerobic anoxygenic phototrophic bacteria cal niches. FEMS Microbiol Rev 35:957–976. https://​doi.​org/​10.​ reveal an ecological consistency of phylogenetic clades in dif- 1111/j.​1574-​6976.​2011.​00292.x ferent oceanic biomes. Sci Rep 8:1–10. https://​doi.​org/​10.​1038/​ 16. Kolber ZS, van Dover CL, Niederman RA, Falkowski PG (2000) s41598-​018-​22413-7 Bacterial photosynthesis in surface waters of the open ocean. 32. Jiao N, Zhang Y, Zeng Y, Hong N, Liu R, Chen F, Wang P (2007) Nature 407:177–179. https://​doi.​org/​10.​1038/​35025​044 Distinct distribution pattern of abundance and diversity of aerobic 17. Nagashima KVP, Hiraishi A, Shimada K, Matsuura K (1997) anoxygenic phototrophic bacteria in the global ocean. Environ Horizontal transfer of genes coding for the photosynthetic reac- Microbiol 9:3091–3099. https://​doi.​org/​10.​1111/j.​1462-​2920.​ tion centers of purple bacteria. J Mol Evol 45:131–136. https://​ 2007.​01419.x doi.​org/​10.​1007/​PL000​06212 33. Bibiloni-Isaksson J, Seymour JR, Ingleton T, van de Kamp J, 18. Achenbach LA, Carey J, Madigan MT (2001) Photosynthetic and Bodrossy L, Brown MV (2016) Spatial and temporal variability phylogenetic primers for detection of anoxygenic phototrophs in of aerobic anoxygenic photoheterotrophic bacteria along the east natural environments. Appl Environ Microbiol 67:2922–2926. coast of Australia. Environ Microbiol 18:4485–4500. https://​doi.​ https://​doi.​org/​10.​1128/​AEM.​67.7.​2922-​2926.​2001 org/​10.​1111/​1462-​2920.​13436 19. Tank M, Thiel V, Imhoff JF (2009) Phylogenetic relationship of 34. Yutin N, Suzuki M, Teeling H, Weber M, Venter JC, Rusch DB, phototrophic purple sulfur bacteria according to pufL and pufM Béjà O (2007) Assessing diversity and biogeography of aerobic genes. Intern Microbiol 12:175–185. https://​doi.​org/​10.​2436/​20.​ anoxygenic phototrophic bacteria in surface waters of the Atlantic 1501.​01.​96 and Pacific Oceans using the Global Ocean Sampling expedition 20. Béjà O, Suzuki MT, Heidelberg JF et al (2002) Unsuspected metagenomes. Environ Microbiol 9:1464–1475. https://​doi.​org/​ diversity among marine aerobic anoxygenic phototrophs. Nature 10.​1111/j.​1462-​2920.​2007.​01265.x 415:630–633 35. Cuadrat R, Ferrera I, Grossart HP, Dávila AMR (2016) Picoplank- 21. Yutin N, Suzuki MT, Béjà O (2005) Novel primers reveal wider ton bloom in global south? A high fraction of aerobic anoxygenic diversity among marine aerobic anoxygenic phototrophs. Appl phototrophic bacteria in metagenomes from a coastal bay (Arraial Environ Microbiol 71:8958–8962. https://​doi.​org/​10.​1128/​AEM.​ do Cabo—Brazil). OMICS 20:76–87. https://d​ oi.o​ rg/1​ 0.1​ 089/o​ mi.​ 71.​12.​8958-​8962.​2005 2015.​0142 22. Koh EY, Phua W, Ryan KG (2011) Aerobic anoxygenic phototro- 36. Graham ED, Heidelberg JF, Tully BJ (2018) Potential for primary phic bacteria in Antarctic sea ice and seawater. Environ Micro- productivity in a globally-distributed bacterial phototroph. ISME biol Rep 3:710–716. https://​doi.​org/​10.​1111/j.​1758-​2229.​2011.​ 12:1861–1866. https://​doi.​org/​10.​1038/​s41396-​018-​0091-3 00286.x 37. Polz MF, Cavanaugh CM (1998) Bias in template-to-product 23. Villena-Alemany C, Mujakić I, Porcal P, Koblížek M, Piwosz K ratios in multitemplate PCR. Appl Environ Microbiol 64:3724– (2022) Diversity dynamics of aerobic anoxygenic phototrophic 3730. https://​doi.​org/​10.​1128/​aem.​64.​10.​3724-​3730.​1998 13 2172 C. R. Gazulla et al. 38. Koblížek M (2015) Ecology of aerobic anoxygenic phototrophs in growth curves document large differences in response of indi- aquatic environments. FEMS Microbiol Rev 39:854–870. https://​ vidual groups of marine bacteria to the top-down and bottom-up doi.​org/​10.​1093/​femsre/​fuv032 controls. mSystems 6:e00934-21. https://​doi.​org/​10.​1128/​mSyst​ 39. Mendler K, Chen H, Parks DH, Lobb B, Hug LA, Doxey AC ems.​00934-​21 (2019) Annotree: visualization and exploration of a functionally 53. Mašín M, Zdun A, Ston-Egiert J, Nausch M, Labrenz M, Moulis- annotated microbial tree of life. Nucleic Acids Res 47:4442–4448. ová V, Koblížek M (2006) Seasonal changes and diversity of aero- https://​doi.​org/​10.​1093/​nar/​gkz246 bic anoxygenic phototrophs in the Baltic Sea. Aquatic Micro Ecol 40. Wright ES (2016) Using DECIPHER v2.0 to analyze big biologi- 45:247–254. https://​doi.​org/​10.​3354/​ame04​5247 cal sequence data in R. R Journal 8:352–359. https://​doi.​org/​10.​ 54. Gasol JM, Cardelús C, Morán XAG et al (2016) Seasonal patterns 32614/​rj-​2016-​025 in phytoplankton photosynthetic parameters and primary produc- 41. Seemann T (2014) Prokka: rapid prokaryotic genome annotation. tion at a coastal NW Mediterranean site. Sci Mar 80:63–77 https://​ Bioinformatics 3:2068–2069. https://​doi.​org/​10.​1093/​bioin​forma​ doi.​org/​10.​3989/​scimar.​04480.​06E tics/​btu153 55. Grossart HP, Massana R, McMahon KD, Walsh DA (2020) Link- 42. Kanehisa M, Goto S, Sato Y, Kawashima M, Furumichi M, Tan- ing metagenomics to aquatic microbial ecology and biogeochemi- abe M (2014) Data, information, knowledge and principle: back cal cycles. Limnol Oceanogr 65:S2–S20. https://​doi.​org/​10.​1002/​ to metabolism in KEGG. Nucleic Acids Res 42:199–205. https://​ lno.​11382 doi.​org/​10.​1093/​nar/​gkt10​76 56. Delmont TO, Quince C, Shaiber A et al (2018) Nitrogen-fixing 43. El-Gebali S, Mistry J, Bateman A et al (2019) The Pfam protein populations of Planctomycetes and Proteobacteria are abundant families database in 2019. Nucleic Acids Res 47:D427–D432. in surface ocean metagenomes. Nat Microbiol 3:804–813. https://​ https://​doi.​org/​10.​1093/​nar/​gky995 doi.​org/​10.​1038/​s41564-​018-​0176-9 44. Sanchez P, Sebastián M, Pernice M et al (2023) Marine pico- 57. Cornejo-Castillo FM, Zehr JP (2021) Intriguing size distribution plankton metagenomes from eleven vertical profiles obtained by of the uncultured and globally widespread marine non-cyanobac- the Malaspina Expedition in the tropical and subtropical oceans. terial diazotroph Gamma-A. ISME 15:124–128. https://​doi.​org/​ bioRxiv. https://​doi.​org/​10.​1101/​2023.​02.​06.​526790 10.​1038/​s41396-​020-​00765-1 45. Martin M (2013) Cutadapt removes adapter sequences from high- 58. Oz A, Sabehi G, Koblízek M, Massana R, Béjà O (2005) Roseo- throughput sequencing reads. EMBnet J 17:10. https://​doi.​org/​10.​ bacter-like bacteria in Red and Mediterranean Sea aerobic 14806/​ej.​17.1.​200 anoxygenic photosynthetic populations. Appl Environ Microbiol 46. Callahan BJ, Mcmurdie PJ, Rosen MJ, Han AW, Johnson AJA, 71:344–353. https://​doi.​org/​10.​1128/​AEM.​71.1.​344-​353.​2005 Holmes SP (2016) DADA2: high-resolution sample inference 59. Karr EA, Sattley WM, Jung DO, Madigan MT, Achenbach LA from Illumina amplicon data. Nat Methods 13:7. https://​doi.​org/​ (2003) Remarkable diversity of phototrophic purple bacteria in 10.​1038/​nmeth.​3869 a permanently frozen Antarctic lake. Appl Environ Microbiol 47. Barbera P, Kozlov AM, Czech L, Morel B, Darriba D, Flouri T, 69:4910–4914. https://​doi.​org/​10.​1128/​AEM.​69.8.​4910-​4914.2​ 003 Stamakis A (2019) EPA-ng: massively parallel evolutionary place- 60. Boeuf D, Cottrell MT, Kirchman DL, Lebaron P, Jeanthon C ment of genetic sequences. Syst Biol 6:365–369. https://​doi.​org/​ (2013) Summer community structure of aerobic anoxygenic pho- 10.​1093/​sysbio/​syy054 totrophic bacteria in the western Arctic Ocean. FEMS Microbiol 48. McMurdie PJ, Holmes S (2013) Phyloseq: an R package for repro- Ecol 85:417–432. https://​doi.​org/​10.​1111/​1574-​6941.​12130 ducible interactive analysis and graphics of microbiome census 61. Lehours A, Jeanthon C (2015) The hydrological context deter- data. PLoS One 8:4. https://​doi.​org/​10.​1371/​journ​al.​pone.​00612​ mines the beta-diversity of aerobic anoxygenic phototrophic bac- 17 teria in European Arctic seas but does not favor endemism. Front 49. Wickham H, Averick M, Bryan J (2019) Welcome to the Microbiol 6:1–9. https://​doi.​org/​10.​3389/​fmicb.​2015.​00638 Tidyverse. J Open Source Softw 4:1686. https://d​ oi.o​ rg/1​ 0.2​ 1105/​ 62. Bonilla-Rosso G, Wittorf L, Jones CM, Hallin S (2018) Design joss.​01686 and evaluation of primers targeting genes encoding NO-forming 50. Oksanen J, Simpson G, Blanchet F et al (2022) vegan: Community nitrite reductases: implications for ecological inference of denitri- Ecology Package. R package version 2.6–2. https://​cran.r-​proje​ct.​ fying communities. Sci Rep 6:1–8. https://​doi.​org/​10.​1038/​srep3​ org/​packa​ge=​vegan. Accessed 25 Apr 2023 9208 51. Wickham H (2016) ggplot2: Elegant graphics for data analy- 63. Parada AE, Needham DM, Fuhrman JA (2016) Every base mat- sis. Springer-Verlag, New York. https://​ggplo​t2.​tidyv​erse.​org. ters: assessing small subunit rRNA primers for marine microbi- Accessed 25 Apr 2023 omes with mock communities, time series and global field sam- 52. Fecskeová LK, Piwosz K, Šantic D, Šestanovic S, Tomaš AV, ples. Environ Microbiol 18:1403–1414. https://​doi.​org/​10.​1111/​ Hasunová M, Šolic M, Koblížek M (2021) Lineage-specific 1462-​2920.​13023 13