Antarctic microbial wars: observations on the antagonisms between Fungi, Archaea, and Bacteria in Livingston Island

1. Introduction

Antarctica, despite being isolated, remote, and difficult to access because of its geographic position and characteristics, can be considered as a major climate generator [1]. Microorganisms are few of ones who can survive in the Antarctic harsh terrestrial environments and form the major component of the biomass within the water ecosystems where they are also at the basis and at the end of every food chain. Antarctic microorganisms possess specific adaptations, allowing them to survive in these challenging environments [2,3]. However, these adaptations render them susceptible to environmental changes, thereby emphasizing the need to investigate them further.

Unfortunately, Antarctic microbial communities are still poorly studied, mostly because of the challenging environments and the remoteness from the needed adequate laboratory infrastructure. However, with the advances of the Next Generation Sequencing (NGS)-based metagenomics techniques, this task became more feasible, relying on in situ direct DNA isolation, followed by the transportation of the extracted DNA samples to the research centers and laboratories. Some good examples of such metagenomic studies are the works of Picazo et al. [4], Kim et al. [5], Coleine et al. [6], and Fernández et al. [7]. Yet, the focus of these studies was put on the characterization of the fungal, archaeal, and bacterial compositions within the soils and the aquatic environments but not on the existing relationships between these three domains.

In January – February 2022 during the 30th Bulgarian polar expedition to Livingston Island (South Shetland Islands, Maritime Antarctica), one of the research projects that members of the team worked on was focused on the taxonomical characterization of some soil and water microbial ecosystems. When we took a first look at the results, we were surprised to find low archaeal richness, as well as an inversely proportional dependence between the richesses of Archaea and Fungi, and to a lesser extent between bacteria and archaea, and fungi and bacteria. Despite an antagonism between Bacteria and Fungi has been reported for agricultural soils [8], no data was available for the Arctic and the Antarctic regions, so we decided to investigate these domains’ antagonisms which could be related to and/or to be more acute in the harsh and in most cases oligotrophic Antarctic environment. To do this we proposed to quantify this antagonism based on some of the metagenomic sequencing data parameters.

2. Materials and Methods

2.1. Sampling sites and sampling methods

All sampling sites are listed in Table 1 and are located around the Bulgarian Antarctic Base “St. Kliment Ohridski” on Livingston Island, Maritime Antarctica. For the solid samples, 150-250 mg of biomass or sediment were collected in as much as possible sterile conditions with a sterile spatula and put into sterile Eppendorf tubes. For the water samples, 2 liters of water were collected in new and unused plastic bottles which were washed 3 times with the water of the sampling site before the collection. All samples were transported immediately to the research station, and DNA was usually isolated within 45-60 minutes, except for the Hannah Point samples whose transportation to Bulgarian Antarctic Base required about 90-120 minutes.

2.2. DNA extraction

To reduce to a minimum possible shift in the microbial compositions of the collected samples, total DNA from them was isolated just after the collection. For the solid samples, the extraction was performed with ZR Soil Microbe DNA MiniPrep Kit (Zymo Research Corp., USA) according to the manufacturer’s Instructions manual. For the water samples, two liters of the water were filtered through a 0.2 µm filter, and then the biomass was washed off with 700 µl of the ZR Soil Microbe DNA MiniPrep Kit’s Bashing Beads buffer. The eluted biomass was transferred within the kit’s bashing beads tubes and further processed as the solid samples. The final elution step was performed in 50 µl of the provided Elution buffer. Two microliters of the eluted DNA were used for the determination of the concentration on a Quantus Fluorimeter (Promega Corp., USA). The eluted DNA samples were kept and transported at -20 °C.

2.3. Metagenomic sequencing

The extracted DNA was shipped on dry ice to the Novogene Company Ltd. (Cambridge, UK) for sequencing on the Illumina HiSeq 2 x 250 paired-end reads platform. The V3-V4 region of the genes encoding the 16S rRNA was chosen for the amplification of eubacterial 16S rRNA with primers 341F (5’-CCTAYGGGRBGCASCAG-3’) and 806R (5’-GGACTACNNGGGTATCTAAT-3’). For fungi, the ITS2 region was amplified using primers ITS3-2024F (5’-GCATCGATGAAGAACGCAGC-3’) and ITS4-2409R (5’-TCCTCCGCTTATTGATATGC-3’), while for archaea the V4-V5 16S rRNA region was amplified using primers Arch519F (5’- CAGCCGCCGCGGTAA-3’) and Arch915R (5’- GTGCTCCCCCGCCAATTCCT-3’). The sequencing was performed to generate 30,000 tags per sample. Novogene Company LTD performed the primary bioinformatic processing.

2.4. Data processing and operational taxonomic units (OTUs) analyses

The raw files generated from the NGS platform were first subjected to demultiplexing and trimming. Next, the FLASH V1.2.7 software tool [9] was used to merge the paired-end reads. The raw tags quality filtering was performed according to Bokulich et al. [10] and Caporaso et al. [11]. The UCHIME algorithm was used to detect chimera sequences by comparing the tags with the reference database (http://drive5.com/uchime/uchime_download.html) [12]. The effective tags were obtained by removing the chimera sequences, accordingly to Haas et al. [13]. The Uparse v7.0.1001 software [14] was used to assign the effective tags with ≥97% similarity into OTUs. For species annotation at each taxonomic rank, each eubacterial and archaeal representative sequence was compared against the SSUrRNA Database of the SILVA138 Database with the Mothur software [15,16]. For fungi, the comparison was performed again with the Mothur software against the Unite V8.2 database [17]. The alpha diversity analysis was estimated based on the community richness index by the Chao1 estimator (http://scikit-bio.org/docs/latest/generated/skbio.diversity.alpha.chao1.html#skbio.diversity.alpha.chao1) [18] and the ACE estimator (http://scikit-bio.org/docs/latest/generated/skbio.diversity.alpha.ace.html#skbio.diversity.alpha.ace) [19,20]. The community diversity was assessed by the Shanon index (http://scikit-bio.org/docs/latest/generated/skbio.diversity.alpha.shannon.html#skbio.diversity.alpha.shannon) and the Simpson index (http://scikit-bio.org/docs/latest/generated/skbio.diversity.alpha.simpson.html#skbio.diversity.alpha.simpson) [21,22]. Rarefaction curves were constructed in order to assess the reliability of the obtained data [23].

2.5. Antagonism analyses and relative measuring

The mean values of the number of the total tags, the number of the OTUs, the Shannon, the Simpson, The Chao1, and the ACE indices were calculated, then each sample was attributed to “+” or “-“ depending on if its own value was greater or lesser than the mean value for the given type of microorganisms. Then the three types of binome antagonisms in each sample were relatively quantified as the sums of the cases where a discrepancy was observed between the two compared types of microorganisms.

2.6. Pearson and Spearman correlation analysis

Pearson and Spearman correlation analyses were performed between the three groups of microorganisms for each of the alpha diversity metrics, as well as the effective tags. For testing the statistical significance of the results regarding the community as a whole, a significance threshold of p=0.05 was chosen. Pearson correlation assesses the presence of a linear correlation between two variables with normal distribution, while Spearman correlation assesses the presence of a monotonic relationship. These values range from -1 for the presence of a high negative correlation, 0 for the absence of a correlation, to 1 for a presence of a high positive correlation.

3. Results

3.1. Sequencing statistics and alpha diversity indices

The numerical values of the amplicon-based sequencing statistics and some of the alpha diversity indices are presented in Table 2, where the average values of each index in the three types of microorganisms are calculated on the last row. The rarefaction curves of each sample are shown in Figure 1.

3.2. Antagonism study

The numerical values for each index (the number of the total tags, the number of the OTUs, the Shannon, the Simpson, the Chao1, and the ACE indices) in each sample are presented as bar histograms regarding the mean values in Figure 2. To ease the interpretation of the results, each sample was attributed to a “+” (greater) or “-“ (lesser) for each index according to its proper value in comparison to the mean value for this index (Table 3). A summarization of the cases of discrepancies is presented in Table 4 as total numbers as well as percentages.

Pearson linear correlation showed a statistically significant negative correlation between the effective tags of archaea and fungi, but not between any of the other metrics. (Table 5). On the other hand, a very high positive correlation was shown between archaea and bacteria in regard to both the Simpson and Shannon diversity indices.

4. Discussion

The total number of bacterial and archaeal cells on Earth is estimated to be around 1.2×10³⁰ cells, distributed in five major habitats: deep oceanic subsurface, upper oceanic sediment, deep continental subsurface, soil, and oceans [24]. There are no such quantitative estimations for the Fungi, probably because they can exist in both unicellular and multicellular form. However, they also represent a substantial percentage of the Earth’s biomass – 12.7% [25] and show a tremendous species diversity and functional roles in almost all ecosystems [26]. Still, with very few exceptions, the mutual ecological relationships between these three domains of microorganisms as a whole are poorly investigated. Even with the advent of NGS-based metagenomics, most of the studies are focused on research within a single domain in a given ecological niche.

During the Antarctic research season 2021-2022, four water and fourteen terrestrial samples were collected around the Bulgarian Polar Base on Livingston Island in Maritime Antarctica. The major goal of the research was the characterization of some microbiotas in different ecological niches with the means of amplicon-based metagenomics (submerged microbial surface contaminations of rocks and algae, sediments, soils, biomasses from ice melting ponds, lithotelms and glacial lakes (tarns), fresh waters, and marine waters) (Table 1).

The obtained data were considered reliable and informative due to the slopes of the rarefaction curves which were reaching a plateau (Figure 1), meaning that the data were representative [27]. However, before the in-depth analysis of the sequencing data, two observations came into our sight. First, the number of the Archaea (as total tags as well as the number of OTUs) was much lower than those of Fungi and Bacteria. Second, the archaeal alpha diversity indices (Shannon, Simpson Chao1, and ACE) were from the same order as those of the two other studied domains, meaning comparable species richnesses in all three super kingdoms (Table 2) [27]. These two observations suggest that the presence of archaeal cells within the samples was significantly lower in number compared to bacterial and fungal cells, despite all three kingdoms exhibiting similar levels of richness.

One of the most probable explanations of this disparity is the assumption that some antagonism between Archaea and the other two domains should exist, in turn probably caused by the harsh, and in some cases oligotrophic, Antarctic conditions. So, we decided to relatively measure this antagonistic relationship most simply - by comparing in binomes: in how many samples when one of the domains is more than the observed average, the other is less than the observed average. These results, graphically presented in Figure 2, were difficult to quantify, so, a simpler comparison matrix composed of only “+” and “-“ was constructed (Table 3) which in turn was used in our quantitative analysis of the antagonism relationships (Table 4).

The number of the total tags, which represents the number of effective tags after merging and filtering the sequencing reads, is the most primary parameter because it directly reflects the distribution of the number of the sequenced DNA molecules within the sample. This antagonism was observed in 10 of 18 cases (56%) in all three binomial comparison groups: Fungi vs Archaea; Fungi vs Bacteria and Archaea vs Bacteria.

The next parameter which was analyzed was the number of OTUs in each sample. Analyzing the OTUs numbers has the advantage to reflect the taxonomic groups [28]. In this case the discrepancy between Archaea and Fungi was more distinct – it was observed in 2/3 of all samples and even more in the water samples - 3/4 of the cases. Concerning the binome Fungi-Bacteria, the discrepancy was less distinct than in the case of the total number of the tags (in only 44% of the samples) while between Archaea and Bacteria, discrepancies were observed within the same percentage as within the total number of the tags.

Concerning the four observed alpha diversity indices, no such clear tendencies could be observed. Still, in our opinion, they should also be taken into account because they reflect the microbial communities’ structures. The observed discrepancies in all the three binomes concerning the Shanon and the Simposon indices [21,22,29] varied between 17% and 50%. Both are used as estimators of the species richness and evenness, so the low percentages could be explained by the high OTUs values for the three domains within the analyzed samples. On the other hand, more interesting results were obtained for the Chao1 and the ACE richness estimators which are abundance-based [18,19,20]. Regarding the Chao1 index, a discrepancy in more than 50% of the cases was observed only for the binome Archaea-Bacteria. This can be explained by the fact that this index is based on the singletons and the doubletons, thus estimating the “missed” species, Bacteria being largely the most species-abundant superkingdom. On the other hand, discrepancies in 56% and 61% of the cases were observed respectively for the binomes Fungi-Archaea and Archaea-Bacteria. This tendency corresponds more to those observed for the numbers of the total tags and OTUs because the ACE richness estimator takes into account both the abundant and the rare species, and thus should be more informative when domains’ antagonisms are investigated.

The presence of a moderate, negative correlation between the effective tags of fungi and archaea, coupled with the lack of such a distinct correlation in the rest of the metrics hints at a possible explanation for these results. Firstly, it is important to note that the effective tags represent the sequencing reads before they are grouped into OTUs based on a 97% similarity threshold. If the majority of the effective tags that are associated with the presence of a correlation between two groups are also assigned to a small number of OTUs, it is most likely that any potential interactions are limited within a subset of the species in the sample, and thus will not be as apparent in the alpha diversity metrics. On the other hand, if no correlation is observed between the number of effective tags, yet such a correlation is present within the alpha diversity metrics, this could point to a large number of species contributing to the interaction. We see a similar trend in some of the interactions between archaea and bacteria, where the lack of correlation in effective tags is followed by a very high correlation in diversity.

Additionally, it is crucial to outline that the given samples represent vastly different ecological niches, encompassing sediment, salt water, fresh water, biomass, and surface contamination samples. In fact, roughly half of the samples are from a type of surface contamination, and thus the correlation results will inevitably be skewed in favor of the specific interactions within those types of communities. The annotated taxonomy, which is not included in this paper, also points to the existence of microbial profiles that are specific to each sampling location. Thus, the interactions that take place on the scale of the sampling locations will inevitably overlap with the presence of any macro-scale interactions that could be attributed to the Antarctic climate as a whole. However, due to the lack of samples for a large number of the sampling sites, it is difficult to exclude them from the dataset, without compromising how representative the sample is of the Antarctic community on Livingston Island as a whole.

An interesting case is sample S22 which was obtained from biomass on a submerged rock in the new nameless tarn located some several hundred meters North-Est of the Bulgarian Polar Base “St. Kliment Ohridski”. This is the only sample where the numbers of the total tags and the OTUs in Archaea exceed those of Fungi, and approach those of Bacteria. More interestingly, the number of the archaeal total tags is some 16% greater than the fungal one. However, the archaeal OTUs are more than 6 times more numerous than the fungal ones. Not surprisingly, the greatest Chao1 and ACE community richness estimators among Archaea are observed in this sample in comparison to all other samples we analyzed. At this stage of research, it would be difficult and to a somewhat extent speculative to explain this observation, still, it is worth noting that this tarn was formed only very recently – during the last pre-Covid19 Bulgarian Antarctic expedition in 2018-2019 permanent glacier existed in its location, and the tarn discovered only during the 2020-2021 expedition when the scientists returned to the Bulgarian Polar Base (the 2019-2020 expedition being only logistical for base maintenance). So, we cannot say if the glacier retracted to form the tarn one or two years before the sample was taken off, but surely this tarn represents a very new ecosystem. This fact makes plausible the speculation of primary succession of the super kingdoms in newly colonized areas. It has been observed that microorganisms are the first colonizers after the receding glaciers in the mountains [30] so, it is logical that the same process occurs in Antarctica. Unfortunately, the primary microbial succession in barren lands and other ecological niches is poorly investigated [31]. According to our hypothesis, Archaea and Bacteria are the first, while for Fungi to colonize the new environment more time is needed. If this hypothesis is true, it would explain why in the “old” ecological niches in our study in general Fungi dominate over Archaea.

One of the reasons we consider our data informative and reflecting the actual situation is that we chose a non-discriminative method for DNA isolation based on the physical destruction of the microbial cells. Another one is that we used the same DNA samples for the three types of metagenomic analyses in combination with equal sequencing efforts – the generation of at least 30,000 tags per sample, which resulted in 3-4 times more tags except for Archaea (meaning that the archaeal DNA is indeed less represented within the sample). This means that comparisons between the three investigated domains are possible, especially in light of the plateaued refraction curves, prerequisites for correct enough quantitative analyses in amplicon-based metagenomics where primer annealing biases could be a concern.

Our study, as a pioneering one, has some drawbacks. First, we acknowledge that physicochemical parameters in the different ecological niches studied could impact the observed antagonisms. Still, we believe that the tendencies we observed are visible enough. Unfortunately, studying the exact mechanisms of the antagonistic interactions was far beyond the scope of the research project. At present, we can only suspect that both types of antagonistic interactions are involved – active, due to the synthesis of different types of antimicrobials, and passive, resulting from the concurrence for nutrients and energy between the different types of microorganisms. The lack of replicas should be considered as another drawback. Nevertheless, replicating these experiments is not possible because of the location of the area where the study was performed – Livingston Island in Maritime Antarctica. Additionally, the polar bases are seasonal, meaning that they are accessible only during the Antarctic summer. Still, because of the seasonal landscape erosion caused by the retreat of the glaciers every summer, the small glacial ponds where most of the biomass samples were taken off form in different locations every year. Finally, the third major drawback of the research is the lack of a well-established methodology for quantification of the antagonisms between the three domains. Unfortunately, this drawback could not be adequately addressed because of the very scarce publications of similar studies. Here, we only propose a very basic methodology that worked well for a field study of a location with limited seasonal accessibility. For these reasons, this manuscript should be regarded as a preliminary communication that could be the basis for future research projects.

5. Conclusions

We could only speculate about the causes of the observed antagonism between the three major microbial phyla in Maritime Antarctica. However, our data suggest that some antagonisms exist, especially between Fungi and Archaea and between Archaea and Bacteria, based on the values of the numbers of the total tags and the OTUs, as well as on the values of the Chao1 and the ACE community richness indices. Furthermore, we observed some indications about the order of colonization of new barren ecological niches, according to which Archaea and Bacteria are the first, followed by Fungi. Based on our findings the best indicators of antagonism are the numbers of the OTUs and the ACE estimator, giving the most discriminative results. Unfortunately, because of the general lack of similar studies, not only for Antarctica but also worldwide, we cannot either compare our findings or explain the causes of the antagonisms.

Our results of the studies of pioneer microbial communities are only preliminary thus further studies including inter-domains functional relationships undoubtedly are needed. Yet, we demonstrated the power of NGS-based metagenomics in research on inter-domain relationships, as well as in the investigation of the primary colonization of new ecological recesses.