Seasonal Dynamics and Diversity of Haemosporidians in a Natural Woodland Bird Community in Slovakia

Abstract

Despite the ubiquity of disease seasonality, mechanisms behind the fluctuations in seasonal diseases are still poorly understood. Avian hemosporidiosis is increasingly used as a model for ecological and evolutionary studies on disease dynamics, but the results are complex, depending on the focus (hosts, parasites, vectors) and scale (individuals, community, populations) of the study. Here, we examine the local diversity of haemosporidian parasites and the seasonal patterns of infections, parasite richness, and diversity in a natural woodland bird community in Slovakia. In 35 avian species, we detected 111, including 19 novel, haemosporidian cytochrome b lineages. The highest numbers of lineages were detected during spring and autumn, corresponding with higher avian species richness and infection prevalence in the avian community during these periods of time. Nevertheless, the haemosporidian community in the local breeders in summer was relatively stable, Haemoproteus lineages dominated in the local avian haemosporidian community, and only few parasite lineages were abundant within each genus. While prevailing Leucocytozoon infections in spring suggest that the majority of sampled birds wintered in the Mediterranean region, Plasmodium infections in spring can be due to relapses in reproductively active short-distance migrants. Multiple haemosporidian infections, both intra- and inter-generic ones, were common in the local avian community. Infection intensity peaked during summer and tended to be higher in older birds, pointing to the role of supressed immunity in reproductively active birds.

1. Introduction

Despite the ubiquity of seasonal disease incidences, mechanisms behind the fluctuations in seasonal diseases are poorly understood [1,2,3]. Disease seasonality is thought to be driven by seasonal changes in: (1) abiotic conditions affecting pathogen survival outside hosts, (2) host behaviour affecting host abundance and contact rate, (3) host immune function affecting host susceptibility to pathogens, and (4) the abundance of vectors and heterospecific hosts affecting pathogen availability [3]. However, the identification of the underlining mechanisms behind the seasonality is difficult due to the complexity of interactions between the individual drivers [2]. For example, disease seasonality can exhibit different patterns depending on whether the focus is on a single host species or a wider host community [4]. Additionally, disease detection may strongly depend on infection intensity [5], which can importantly confound the patterns of disease prevalence and seasonality.

Avian hemosporidiosis is increasingly used as a model for ecological and evolutionary studies on disease dynamics [6]. Seasonality in avian parasitaemia was first described by Danilewsky [7], suggesting that parasitaemia in wild birds in Ukraine was higher in warm seasons [6]. Beaudoin et al. [8] developed a functional model of avian malaria prevalence for songbirds of the temperate climate zone in North America, identifying prevalence peaks in spring and late summer/early autumn. The spring peak is usually explained by infection relapses, which could be induced by changing hormone levels, internal clock of parasites, or through stimulation by vectors [9,10,11,12]. Strenuous physical activity during spring migration can also contribute to the higher spring infection prevalence [12,13]. Yet, long-distance migratory birds might show higher parasite prevalence and diversity than short-distance migrants or resident birds, as the former encounter parasites across more diverse ecological settings [14,15,16,17,18]. High bird densities at stopover sites may also increase the risk of avian hemosporidiosis in migratory birds [19].

Infection seasonality may be related to fluctuations in absolute and relative numbers of birds of different age cohorts. Most birds contract infection as fledglings or during the first weeks after fledging, likely due to weaker immune function [20,21]. In contrast, older birds often show higher prevalence than juveniles due to the accumulation of haemosporidian infections throughout their lifetime [22,23]. For example, Emmenegger et al. [24] studied the effect of age in European bee-eaters Merops apiaster, reporting that haemosporidian prevalence for birds in their second calendar year was lower by 10% compared to that in older birds. Nonetheless, infection seasonality in birds can occur even within age cohorts. While high infection intensity can decline soon with the development of the immune system in fledglings [20], breeding activity can temporarily increase infection intensity through the supressed immune system in adult birds, e.g., [21]. Consequently, in addition to a relatively high abundance of immunologically naïve juveniles, the infection peak at the end of summer can also be due to higher numbers of immunologically supressed adult breeders towards the end of the breeding season.

The local structure of haemosporidians and avian hosts can also affect seasonal dynamics in infection prevalence and intensity in the local avian community. Haemosporidian lineages that show higher host specificity perform better in related host species with a similar immune system [25,26]. The relative abundance of host species also is important because abundant host species are thought to reach higher parasite prevalence and diversity than less common host species [4,27]. In turn, generalist parasites may infect a wider range of host species than specialists, but they can reach lower intensities in their hosts compared to parasites specialised on fewer host species [26,28]. Therefore, seasonal fluctuations in infection prevalence and intensity can be related to intra- and inter-seasonal fluctuations in local host and parasite communities, both in terms of diversity and species composition.

Numerous studies reported that the same parasite lineages can perform differently between host species [4,16,19,26,29,30,31,32]. For example, Lynton-Jenkins et al. [31] reported a decreased prevalence of one group of Leucocytozoon lineages in great tits between autumn and spring, while the prevalence of the same lineages increased in blue tits. Furthermore, a study by Hunag et al. [21] detected seasonally distinct infection patterns of Haemoproteus lineage WW2 in garden warbler Sylvia borin and Phylloscopus spp. Interestingly, these authors also recorded differences in the intensity of infection of lineage PARUS1 between juvenile and adult birds. Furthermore, differences in haemosporidian lineage prevalence were also shown between populations within host species [33]. In fact, local adaptations were suggested to play role in certain host-parasitic associations. In Sardinia, endemic haemosporidian lineages were dominant in sedentary birds compared to short-distance and long-distance migrants [34]. The community-wide outcome of these associations for disease seasonality, therefore, can be strongly dependent on the local structure and strength of host-parasite associations [26].

Here, we examine the local diversity of haemosporidian parasites in a natural wetland bird community in Slovakia. Our aim is to describe the seasonal patterns of infection prevalence, infection intensity, and parasite richness and diversity. Moreover, associations between avian species and haemosporidian cytochrome b lineages are examined to establish their local structure, abundance, and seasonal fluctuations, as well as the levels of structural and phylogenetic host specificity of the lineages detected.

2. Materials and Methods

2.1. Study Area, Field Methods and Study Species

Wild birds were sampled in south-east Slovakia near the Drienovec village at the Drienovec Bird Ringing Station (48°1′11.303″ N; 17°4′2.255″ E; see [35]). The study site (ca. 7.7 ha) is represented by a mosaic of woody wetland and forest meadow ecotones at 190 m a.s.l. Birds were captured, banded, and sampled under the permits of the Ministry of the Environment of the Slovak Republic No. 269/132/05-5.1_p and 9830/2017-6.3. Birds were captured using mist nets for three years (2017, 2018, 2019) between April and November. The blood sample was taken from a brachial vein and ring code, species and, if possible, age and sex were recorded for each bird sampled, and the birds were subsequently released. Blood samples were stored in 70% ethanol at 4 °C until DNA extraction (within 7 months). Totally, 1851 birds of 61 species were blood-sampled (spring: April, n = 444; summer: June–July, n = 474; autumn: September–November, n = 933; Figure 1). Birds were sampled each year in the second half of April, between mid-June and mid-July, and between mid-September and the beginning of November. These three sampling periods were chosen for our study area and climatic zone to obtain representative samples for haemosporidian community composition in birds during spring migration, breeding, and autumn migration periods, respectively.

2.2. DNA Extraction

DNA extractions for samples collected in 2017 and 2019 were performed using the QIAamp DNA Blood Kit (Qiagen, Hilden, Germany) following the manufacturer’s recommendation. DNA from samples taken in 2018 was extracted using a standard phenol-chloroform extraction with ethanol precipitation [36]. Extracted DNA was resolved to a final concentration of ca. 100 ng/µL and stored at −20 °C until subsequent analyses. The quantity and quality of DNA samples was assessed by NanoPhotometer Pearl (Implen, Munich, Germany).

2.3. PCR Analyses

The DNA samples were examined for haemosporidian infection intensity using quantitative real-time PCR (qPCR) targeting 182 bp fragment of cytochrome b (cyt b) gene [37]. All reactions were carried out using GoTaq qPCR Master Mix (Promega, Madison, WI, USA) on a CFX96 real-time thermocycler (Bio-Rad, Hercules, CA, USA). The total volume of the reactions was 20 μL, containing 10 μL of GoTaq qPCR Master Mix 2×, 0.5 μL of each primer (10 μM concentration), 6 μL of molecular grade water, and 3 μL of DNA template (ca. 300 ng). The following cycling conditions were used: 95 °C for 2 min, followed by 40 cycles of 95 °C for 30 s, and 64 °C for 35 s (with a plate read) followed by a final melt curve analysis using instrument default settings. The samples were run in duplicate for all of the samples, together with two non-template controls to check for non-specific amplifications. At the end of the reactions, the amplification curves and melting curves were inspected to obtain values of threshold cycles (Ct) for each sample and determine false positives [21]. For quantification of parasites, seven samples with known quantities [37] were included in each reaction to establish amplification efficiencies and standard curves.

A synthetic double-stranded DNA product (Eurofins Genomics, Ebersberg, Germany), designed from a 220 bp fragment of the conserved rDNA region of Plasmodium relictum (Accession #NC012426; [37]), was used as a positive control. The DNA was diluted to a starting concentration of 10⁶ copies/µL using online calculator on www.thermofisher.com (accessed on 1 May 2021). This starting solution was then serially diluted by 10-fold to prepare a series of solutions from 10⁶ copies of genomic DNA per µL down to 1 copy/µL (that is, there were seven dilutions: 10⁶, 10⁵, 10⁴, 10³, 10², 10¹, and 10⁰). To determine the microarray limit of detection, 3 µL of these diluted DNA samples were used as templates for the amplifications of the cyt b gene.

We applied two conservative criteria for accepting qPCR data as valid, following Friedl & Groscurth [38]. Firstly, we only accepted PCR runs with the correlation coefficient for the standard curve higher than 0.99. Secondly, if a run was accepted, we only accepted samples if SD values of parasite intensities (i.e., relative numbers of DNA copies; see Section 2.4) for the two replicates were lower than 50. This second criterion roughly corresponds to the criterion by Friedl & Groscurth [38] that SD of Ct values for replicates should be lower than 1. Finally, qPCR data for five of the accepted samples were conservatively excluded from analyses on parasite intensity because qPCR data suggested very high parasite intensities (>366 DNA copies per 100 avian red blood cells) for these samples.

In addition to qPCR analysis, all of the samples were analysed by nested PCR assays targeting a mtDNA cyt b gene fragment [39,40]. In the first step, HaemNF1 and HaemNR3 primers were used to amplify the 617 bp fragment. PCR reactions were performed in 20 cycles, which included initial denaturation at 94 °C for 2 min, followed by 20 cycles of denaturation at 94 °C for 30 s, annealing at 50 °C for 30 s and elongation at 72 °C for 45 s and consequently final elongation at 72 °C for 10 min. Subsequently, the amplicon was added as a template to the second step of nested PCR with primers HaemF and HaemR2, targeting a 480 bp amplicon of Haemoproteus sp. and Plasmodium sp., and the set of HaemFL and HaemR2L primers amplifying a 478 bp fragment of Leucocytozoon sp. The PCR protocol was adjusted so that the number of cycles was raised to 35 and the annealing temperature was lowered to 48 °C for Leucocytozoon. Each reaction with the total volume of 10 µL contained 1 µL of template DNA (ca. 20 ng), 5 µL of SuperHot Master Mix (2X) (Bioron, Ludwigshafen, Germany), 0.25 µL of each primer (with concentrations 10 pmol/µL) and 3.25 µL of miliQ water. PCR products from the second round were separated and visualized on 2% agarose gel with SYBR^® Safe DNA gel stain (Invitrogen, Carlsbad, CA, USA). Positive PCR products were purified using a Qiagen purification kit (Qiagen, Hilden, Germany). The purified PCR fragments were sent for sequencing in both directions (Macrogen Europe, Amsterdam, The Netherlands).

2.4. Sequence and Statistical Analysis

The sequences were edited and aligned using Unipro UGENE software v1.32 [41]. The resulting contigs of 478 and 476 bp for Haemoproteus/Plasmodium and Leucocytozoon parasites, respectively, were examined with the BLAST algorithm in GenBank (http://www.blast.ncbi.nlm.nih.gov/Blast.cgi (accessed on 1 July 2021)) and MalAvi [42] (http://130.235.244.92/Malavi/blast.html (accessed on 1 July 2021)) databases. The sequences without double peaks (n = 252/310, 225/255 and 95/131 sequences for Haemoproteus, Plasmodium and Leucocytozoon parasites, respectively) were considered as single-lineage infections and their identity was assessed by a BLAST search. The majority of these sequences matched (100%) with sequences published in GenBank and Malavi. The remaining single-infection sequences were submitted as new lineages to both databases (n = 10, 1 and 8 sequences for Haemoproteus, Plasmodium and Leucocytozoon parasites, respectively). GenBank accession codes (MalAvi lineage names) for the sequences determined in this study are: Haemoproteus spp.—MZ571097 (ERIRUB02), MZ571099 (PRUMOD04), MZ571100 (TURMER08), MZ571101 (ERIRUB03), MZ571102 (CYACAE09), MZ571103 (CYACAE08), MZ571106 (SYAT50), MZ571107 (SYAT51), MZ571108 (SYAT52), MZ571109 (SYCUR03); Plasmodium spp.—MZ571110 (TURMER09); Leucocytozoon spp.—MZ571111 (TURMER10), MZ571112 (TURMER11), MZ571113 (TUPHI14), MZ571115 (PARUS93), MZ571117 (TURMER12), MZ571118 (TURMER13), MZ571119 (TURMER14), MZ571120 (PARUS94).

The unresolved sequences with double peaks or positions containing two or more different bases were subjected to a BLAST search to identify the most closely matching lineage. If such sequences matched completely known lineages for some of the suggested bases at undefined positions, we considered such sequences as multiple infections involving one or more of the known lineages.

The Exp (H′) index, reflecting the effective number of species [43], is based on the Shannon diversity index H′ [44]. The Exp (H′) was calculated for avian hosts and parasite lineages for spring, summer and autumn.

Host specificity was examined considering two contexts: (1) structural host specificity, reflecting the relative ecological importance of each host species in the host community and (2) phylogenetic host specificity, reflecting the importance of relatedness among host species in the host community [45]. Structural host specificity was described using two indices: H₂ and d′ [46]. These metrics were developed for the analysis of bipartite networks to quantify the level of ecological specialization [47,48]. The H₂ index is a community-level index of specificity, reflecting a measure of structural specificity of the entire network; it ranges from 0 (the most generalist community) to 1 (the most specialist community; [46]). The d′ index is a measure of structural specificity at a species level, and ranges from 0 (the most generalist species) to 1 (most specialist species). Phylogenetic host specificity was described using the SES_MPD index [45], assuming that the more closely related host species have a greater tendency to share the same parasite. The index represents the standardized effect size of the mean phylogenetic distance between each pair of parasitised host species. Positive index values indicate greater phylogenetic distance among parasitised host species (low phylogenetic host specificity), while negative values indicate high phylogenetic host specificity [45].

Seasonal fluctuations in the prevalence and intensity of haemosporidian infection were examined at the level of the entire host community with Bayesian multivariate multilevel models. The data for infection prevalence represented prevalence for each species per season (spring, summer, and autumn) and year (2017–2019). Since we could not determine the infection status for single samples for 3 avian species, infection prevalence was examined for 58 of 61 avian species sampled. Parasite intensity, referring to the number of parasite DNA copies per 100 avian red blood cells, was determined for each sample based on qPCR data and sample DNA concentration following Friedl & Groscurth [38] and assuming that the average genome size of a passerine bird is 2.8 pg [49]. Parasite intensity data was examined with a Bayesian model for each individual bird of known age (2 age categories: birds hatched that year—HY, and birds hatched in previous years—AHY) for which we obtained valid parasite intensity data (n = 359 birds of 35 avian species). We could not determine age for 13 of the 372 birds (the 13 birds comprised 7 species) with valid parasite intensity data and, thus, this data was not included in the Bayesian model. The inclusion of all data, but excluding age as a categorical predictor, would not change the result on the association between date and parasite intensity.

As the data obtained from different species were not independent [50], we accounted for the phylogenetic dependency in both Bayesian models on parasite prevalence and intensity, constructing a covariance matrix of avian species following Hadfield & Nakagawa [51]. Avian phylogenies were obtained from Birdtree.org (http://birdtree.org (accessed on 1 May 2021)), which contains a Bayesian posterior distribution of phylogenies for 6670 avian species/operational taxonomic units [52]. We downloaded 1000 trees from the “Hackett sequenced species” data set for 58 avian species for parasite prevalence data and 35 avian species for parasite intensity data. The maximum credibility tree was selected from the downloaded trees using TreeAnnotator in BEAST v2.6.3 [53].

In addition to modelling phylogenetic covariance as a random effect, the random effect of avian species was included in both Bayesian models to account for the species effects that were unrelated to phylogenetic relationships [54]. The effect of species was nested in year that was used as an additional random effect in both models to account for inter-annual variation in response variables.

The Bayesian model on parasite intensity used log-transformed parasite intensity data (n = 359 intensity values) to satisfy the Gaussian distribution. Due to obtaining multiple observations of parasite intensity per avian species, for this response variable, we ran a phylogenetic model with repeated measurements following Bürkner [54]. Consequently, fixed effects included (1) age of birds (HY and AHY), (2) mean date of sampling per species (continuous variable where April 1st = 1), (3) second power of mean date of sampling per species to account for a potential quadratic relationship between date and parasite intensity, and (4) within species variability of date of sampling within species, which was calculated as a difference between date of sampling and mean date of sampling per species.

For parasite prevalence, we ran a Bayesian model assuming the binomial distribution. Prevalence was examined as proportions (infected birds/all birds) for each species per season and year (n = 220 prevalence values). The proportions between the number of all birds of a specific species and the number of all birds per season and year were used as weights to account for variation in the number of sampled birds per species. Season (spring, summer, and autumn) as a categorical predictor was used as a single fixed effect, as it was not possible to examine prevalence at the scale of date due to sample size limitations. Similarly, sample size limitations did not allow for examining the role of age in infection prevalence at the host community scale.

All statistical analyses and plots were performed with R software [55]. The vegan package [56] was used to obtain values of the Exp (H′) diversity index. The structural host specificity indices were calculated using the bipartite package [57,58] with “H2fun” and “dfun” functions. The bipartite package was used to prepare plots for host-parasite networks for each parasite genus. The SES_MPD index was calculated using the picante package [59]. To test the significance of the values of this index, we used a null model implemented in the picante package as “taxa.labels” that randomizes the names of the host species but preserves the distribution of phylogenetic branches [59]. Bayesian multilevel models were fitted using the brms package and “brm” function [60,61], following Bürkner [54], for the implementation of the phylogenetic relationships between species as well as the repeated measurements within species. The contingency tables were examined with the Chi² test using the “chisq.test” function in the stats package [55].

3. Results

3.1. Avian and Haemosporidian Communities

Our analysis involved 1851 birds of 61 species that were sampled over 3 years during spring, summer and autumn (Figure 2,

Table 1. List of avian species sampled and haemosporidian parasites detected in Drienovec, Slovakia, during 2017, 2018 and 2019. For each bird species we report migratory strategy (MS: SD = short-distance migrant, LD = long-distance migrant), number of sampled individuals during three years and for two age classes (AHY: after-hatch-year, HY: hatch-year), haemosporidian prevalence in percent (proportion), values of structural parasite specificity in terms of the d′ index (d′ H, P and L: specificity for Haemoproteus, Plasmodium and Leucocytozoon lineages, respectively), and haemosporidian cytochrome b lineages (H: Haemoproteus, P: Plasmodium and L: Leucocytozoon). New lineages are in bold. Note that it was not possible to identify lineage for each positive sample while multiple lineages were often identified in individual blood samples. Furthermore, the d′ index was not possible to calculate for each parasite genus due to sample size limitations. Sample sizes for prevalence may not correspond with the number of sampled birds per species as it was not possible to determine infection status for some (n = 76) of the samples.

Bird Species	MS	n (’17/’18/’19)	HY	AHY	Prevalence	d′ H	d′ P	d′ L	Haemosporidian Lineage (n)
Caprimulgiformes
Caprimulgus europaeus	LD	1 (0/1/0)	1	0	0 (0/1)	-	-	-
Gruiformes
Crex crex	LD	1 (0/1/0)	1	0	0 (0/1)	-	-	-
Rallus aquaticus	SD	1 (0/1/0)	0	1	100 (1/1)	-	-	-
Accipitriformes
Buteo buteo	SD	2 (1/0/1)	0	2	0 (0/2)	-	-	-
Coraciiformes
Alcedo atthis	SD	1 (0/1/0)	1	0	0 (0/1)	-	-	-
Piciformes
Dendrocopus major	SD	7 (1/1/5)	3	4	17 (1/6)	-	-	0.47	L: SFC8 (1)
Dryobates minor	SD	1 (1/0/0)	1	0	- (0/0)	-	-	-
Jynx torquilla	LD	2 (1/0/1)	1	1	0 (0/2)	-	-	-
Passeriformes
Laniidae
Lanius collurio	LD	5 (4/0/1)	2	3	40 (2/5)	1.00	-	-	H: RBS2 (2)
Corvidae
Garrulus glandarius	SD	11 (3/4/4)	4	7	30 (3/10)	1.00	-	1.00	H: CIRCUM05 (1), GAGLA02 (1), L: COCOR02 (1), GAGLA06 (2)
Remizidae
Remiz pendulinus	SD	6 (3/0/3)	5	1	0 (0/5)	-	-	-
Paridae
Cyanistes caeruleus	SD	121 (42/48/31)	94	28	64 (76/118)	0.56	0.17	0.52	H: CYACAE08 (1), CYACAE09 (1), PARUS1 (11), ROBIN1 (1), SYAT03 (1), P: BT7 (4), COLL1 (1), GRW11 (3), LINN1 (1), PADOM02 (1), SGS1 (7), TUPHI08 (1), TURDUS1 (25), L: PARUS4 (8), PARUS14 (1), PARUS15 (1), PARUS16 (1), PARUS18 (1), PARUS20 (1)
Parus ater	SD	10 (0/7/3)	9	1	70 (7/10)	0.06	0.04	1.00	H: PARUS1 (1), P: GRW11 (2), SGS1 (1), TURDUS1 (1), L: PERATE06 (1)
Parus major	SD	175 (35/86/54)	126	49	53 (92/172)	0.46	0.23	0.58	H: PARUS1 (11), SYAT03 (5), TURDUS2 (1), P: BT7 (4), GRW06 (1), GRW11 (5), SGS1 (27), TURDUS1 (20), L: BT1 (1), PARUS4 (12), PARUS7 (1), PARUS16 (4), PARUS19 (1), PARUS20 (1), PARUS22 (7), PARUS25 (1), PARUS81 (1), PARUS93 (1), PARUS94 (1)
Lophophanes cristatus	SD	1 (0/1/0)	1	0	0 (0/1)	-	-	-
Poecile montanus	SD	1 (1/0/0)	0	1	100 (1/1)	-	-	0.28	L: PARUS22 (1)
Poecile palustris	SD	23 (4/8/11)	14	9	30 (7/23)	0.17	0.17	0.32	H: PARUS1 (2), ROBIN1 (1), P: GRW11 (1), SGS1 (1), L: PARUS20 (1), PARUS22 (1)
Acrocephalidae
Hippolais icterina	LD	2 (0/1/1)	1	1	0 (0/1)	-	-	-
Acrocephalus arundinaceus	LD	1 (0/1/0)	1	0	0 (0/1)	-	-	-
A. palustris	LD	13 (2/5/6)	3	10	31 (4/13)	0.89	-	-	H: ARW1 (2),
A. scirpaceus	LD	5 (0/5/0)	3	2	40 (2/5)	0.70	1.00	-	H: ARW1 (1), P: SW5 (1)
A. schoenobaenus	LD	5 (0/5/0)	5	0	0 (0/5)	-	-	-
Locustellidae
Locustella naevia	LD	4 (0/4/0)	4	0	0 (0/4)	-	-	-
Hirundinidae
Hirundo rustica	LD	7 (4/3/0)	5	2	0 (0/7)	-	-	-
Sylvidae
Sylvia atricapilla	SD	372 (109/158/85)	158	214	60 (213/355)	0.71	0.14	0.50	H: PARUS1 (2), ROBIN1 (3), SYAT01 (29), SYAT02 (37), SYAT03 (25), SYAT07 (12), SYAT10 (1), SYAT11 (2), SYAT12 (2), SYAT13 (9), SYAT14 (1), SYAT16 (16), SYAT17 (1), SYAT28 (3), SYAT33 (1), SYAT44 (4), SYAT50 (1), SYAT51 (1), SYAT52 (1), SYBOR35 (1), TURDUS2 (1), WW2 (2), P: GRW11 (1), LINN1 (1), SGS1 (3), SYAT05 (2), L: PARUS4 (1), PARUS16 (1), SYAT22 (2), SYCON05 (1), SFC8 (1)
S. borin	LD	7 (1/6/0)	4	3	14 (1/7)	-	-	-
S. communis	LD	12 (3/9/0)	4	8	45 (5/11)	0.58	-	-	H: CWT2 (1), SYAT01 (1), WW2 (1)
S. curruca	LD	35 (14/10/11)	2	33	41 (13/32)	1.00	0.38	1.00	LWT1 (1), SYCUR03 (1), P: GRW11 (4), SGS1 (2), L: REB11 (1), RECOB3 (1), SYBOR23 (2)
Phylloscopidae
Phylloscopus collybita	SD	73 (14/16/43)	9	53	9 (6/65)	-	0.00	-	P: TURDUS1 (1)
P. sibilatrix	LD	4 (1/1/2)	1	3	50 (2/4)	0.81	-	-	H: PHSIB2 (1)
P. trochilus	LD	28 (7/16/5)	1	27	32 (9/28)	0.81	-	0.56	CWT4 (2), PARUS1 (1), WW1 (2), L: BT1 (1), SFC8 (1)
Aegithalidae
Aegithalos caudatus	SD	8 (2/3/3)	0	8	13 (1/8)	-	-	-
Regulidae
Regulus ignicapilla	SD	4 (0/1/3)	0	4	0 (0/3)	-	-	-
Regulus regulus	SD	20 (3/6/11)	10	10	20 (4/20)	-	0.00	0.47	P: TURDUS1 (1), L: BT2 (1)
Sittidae
Sitta europaea	SD	5 (2/2/1)	4	1	20 (1/5)	-	0.32	-	P: GRW11 (1)
Certhiidae
Certhia familiaris	SD	9 (2/3/3)	2	7	22 (2/9)	-	0.07	-	P: SGS1 (1)
Troglodytidae
T. troglodytes	SD	29 (12/4/13)	10	19	59 (16/27)	0.00	0.20	0.54	H: SYAT03 (1), P: BT7 (2), TURDUS1 (8), SGS1 (1), L: PARUS20 (1)
Turdidae
Turdus iliacus	SD	1 (0/0/1)	1	0	100 (1/1)	0.05	-	-	H: ROBIN1 (1)
Turdus merula	SD	96 (26/33/37)	45	51	67 (64/95)	0.94	0.50	1.00	TURDUS2 (17), TURMER08 (1), P: LINN1 (9), SYAT05 (16), TURDUS1 (4), TURMER09 (1), L: NEVE01 (2), TUMER03 (2), TUMER10 (1), TURMER10 (1), TURMER11 (1), TURMER12 (1), TURMER13 (1), TURMER14 (1)
Turdus philomelos	SD	40 (8/10/22)	22	18	78 (31/40)	0.74	0.30	1.00	H: TUPHI01 (6), TURDUS2 (1), P: LINN1 (8), SYAT05 (3), TURDUS1 (3), L: EUSE2 (1), SANG02 (1), STUR1 (2), TUPHI06 (7), TUPHI13 (4), TUPHI14 (1)
Muscicapidae
Muscicapa striata	LD	1 (0/0/1)	1	0	- (0/0)	-	-	-
Ficedula albicollis	LD	29 (6/7/16)	15	13	14 (4/29)	0.74	-	-	H: COLL2 (2), COLL3 (1), SYAT03 (1)
Ficedula hypoleuca	LD	3 (1/1/1)	1	1	33 (1/3)		-	-
Luscinia megarhynchos	LD	17 (6/6/5)	8	11	24 (4/17)	0.00	0.53	-	H: SYAT02 (1), P: AFTRU4 (1), LINN1 (2)
Erithacus rubecula	SD	354 (78/171/135)	225	125	25 (85/336)	0.65	0.13	0.68	H: ERIRUB02 (1), ERIRUB03 (1), PARUS1 (2), PHSIB1 (2), ROBIN1 (22), SYAT03 (4), P: BT7 (1), LINN1 (10), SGS1 (5), SYAT05 (2), SYBOR2 (1), TURDUS1 (7), L: BT2 (4), PARUS25 (1), SFC8 (2)
Phoenicurus phoenicurus	LD	6 (1/3/2)	3	3	17 (1/6)	-	-	-
Phoenicurus ochruros	SD	5 (1/1/3)	2	3	20 (1/5)	-	0.07	-	P: SGS1 (1)
Prunellidae
Prunella modularis	SD	75 (33/21/21)	54	21	38 (26/69)	0.60	0.18	-	H: DUNNO01 (1), PRUMOD04 (1), ROBIN1 (4), P: BT7 (4), COLL1 (1), SGS1 (3), SYAT05 (1), SYBOR2 (1), TURDUS1 (8)
Passeridae
Passer montanus	SD	5 (1/4/0)	2	3	60 (3/5)	-	0.07	-	P: SGS1 (1)
Motacillidae
Anthus trivialis	LD	8 (0/7/1)	5	3	25 (2/8)	1.00	0.16	1.00	H: YWT2 (1), P: BT7 (1), TURDUS1 (1), L: TRPIP1 (1)
Fringillidae
Fringilla coelebs	SD	31 (8/7/16)	5	26	61 (19/31)	0.97	0.08	1.00	H: CCF2 (1), CCF3 (1), CCF6 (3), CCF23 (1), PHSIB2 (1), P: SGS1 (2), TURDUS1 (3), L: BRAM3 (2)
F. montifringilla	SD	2 (0/0/2)	0	2	0 (0/2)	-	-	-
Coccothraustes coccothraustes	SD	44 (16/12/16)	19	25	60 (25/42)	1.00	0.07	1.00	H: HAWF1 (2), HAWF2 (6), HAWF6 (10), P: SGS1 (1), L: HAWF7 (8)
Linaria cannabina	SD	3 (2/1/0)	2	1	33 (1/3)	-	-	-
Carduelis carduelis	SD	47 (27/20/0)	29	18	13 (6/47)	-	0.13	-	P: GRW11 (1), LINN1 (1), SGS1 (2)
Chloris chloris	SD	12 (10/2/0)	7	5	0 (0/11)	-	-	-
Serinus serinus	SD	1 (1/0/0)	1	0	- (0/0)	-	-	-
Pyrrhula pyrrhula	SD	39 (3/11/25)	20	8	42 (16/38)	1.00	0.09	0.91	H: SISKIN1 (2), P: GRW11 (1), LINN1 (1), TURDUS1 (1), L: PARUS25 (5), PICVIR01 (1), PRUMOD01 (2), PYRPYR01 (1)
Emberizidae
Emberiza citrinella	SD	17 (4/9/4)	3	14	53 (8/15)	1.00	0.13	-	H: EMCIR01 (1), P: GRW11 (2), SGS1 (4), SYAT05 (1), TURDUS1 (1)
Emberiza cia	SD	1 (1/0/0)	0	1	0 (0/1)	-	-	-
Emberiza schoeniclus	SD	3 (0/1/2)	1	2	0 (0/2)	-	-	-

). The majority of the birds sampled comprised short-distance migrants and most of the birds were sampled in autumn rather than spring and summer (Figure 2). Both the richness and diversity of avian species sampled were relatively higher in autumn than spring and summer (

Table 2. Richness (S) and diversity of avian host species and haemosporidian cytochrome b lineages detected in Drienovec, Slovakia, during 2017, 2018 and 2019. Taxonomic diversity was examined in terms of the Exp (H′) index, reflecting the effective number of species based on the Shannon diversity index H’ for the given taxonomic group and time of year.

	Time of Year	S	Exp (H′)
Avian hosts (n = 61)	spring	38	12.6
	summer	30	10.9
	autumn	50	17
Haemosporidian lineages	spring	60	38.0
	summer	44	28.4
	autumn	65	24.4
Haemoproteus (n = 53)	spring	27	17.2
	summer	30	18.4
	autumn	26	13
Plasmodium (n = 14)	spring	8	6
	summer	5	4.2
	autumn	13	5.8
Leucocytozoon (n = 44)	spring	25	22.3
	summer	9	7.6
	autumn	26	15.2

).

Combining the results of quantitative RT-PCR (qPCR) and nested PCR, 767 of the 1775 (43%) bird samples were found to be positive for haemosporidian parasites; 76 samples provided ambiguous results and were not involved in further analyses. At least 1 positive sample was detected in birds of 71% (41/58) avian species (Table 1). Of the 767 positive samples, it was possible to determine the genus of haemosporidians for 700 samples: 675 (96%) samples were positive for Haemoproteus/Plasmodium and 258 (37%) for Leucocytozoon; 233 (33%) samples were positive for both Haemoproteus/Plasmodium and Leucocytozoon. Further, of the 700 samples for which we were able to determine the parasite’s genus, we were able to unambiguously determine at least one cytochrome b lineage for 587 samples: 300, 248 and 119 samples for Haemoproteus, Plasmodium and Leucocytozoon lineages, respectively, including 84 samples positive for both Haemoproteus/Plasmodium and Leucocytozoon parasites.

While the absolute number of positive samples was comparable between spring, summer and autumn for Haemoproteus, the number of positive samples for Plasmodium and Leucocytozoon was markedly higher in autumn than spring and summer (Figure 3). There was a significant relationship between time of year and haemosporidian genus for haemosporidian prevalence (Haemoproteus, Plasmodium, and Leucocytozoon: spring = 25, 9, and 9%, summer = 21, 5, and 4%, autumn = 13, 23, and 8%; Chi²-test: χ² = 16.77, df = 4, p < 0.01). This result on the relative number of positive samples was mainly due to the higher-than-expected prevalence of Plasmodium-positive samples in autumn and the lower-than-expected prevalence of Haemoproteus-positive samples in autumn.

We identified 53, 14 and 44 cytochrome b lineages for Haemoproteus, Plasmodium and Leucocytozoon, respectively, in birds of 35 avian species (Table 1). A total of 19 lineages (10, 1, and 8 for Haemoproteus, Plasmodium, and Leucocytozoon, respectively) were detected for the first time by this study (Table 1). While the richness of haemosporidian lineages was relatively higher in spring and autumn than summer, lineage diversity was markedly higher in spring than summer and autumn (Table 2). The latter result was largely due to the high diversity of Leucocytozoon lineages in spring (Table 2).

A similar number (5–6) of lineages dominated in Haemoproteus and Plasmodium, but only a single lineage dominated in Leucocytozoon (Figure 4). While the dominant Haemoproteus lineages, with the exception of TURDUS2, were regularly detected throughout the whole season, the dominant Plasmodium and Leucocytozoon lineages occurred predominantly in autumn (Figure 4). The most frequently detected haemosporidian lineages during spring and summer comprised five Haemoproteus (spring: SYAT03 and SYAT16; summer: SYAT03, SYAT02, PARUS1, and TURDUS2) and five Plasmodium (spring: SGS1, TURDUS1, LINN1, SYAT05, and GRW11) lineages (Figure 4).

Intra-generic multiple infections (infections involving different lineages of the same haemosporidian genus) were revealed for all the three haemosporidian genera, but multiple infections with more than two lineages were only detected in a single sample for Plasmodium (

Table 3. Relationship between the type of intra-generic infection and time of year for the frequency of haemosporidian cytochrome b lineages detected in Drienovec, Slovakia, during 2017, 2018 and 2019. Note that haemosporidian lineage identity was not known for all parasites involved in multiple infections.

	Infection Type
Time of Year	Single	Double	>Double	Total
Haemoproteus
spring	85	12	5	102
summer	70	19	8	97
autumn	97	10	4	111
Plasmodium
spring	29	6	0	35
summer	16	6	0	22
autumn	180	17	1	198
Leucocytozoon
spring	28	6	4	38
summer	11	2	7	20
autumn	56	10	7	73

). The proportion of multiple (≥double) infections was relatively higher in summer than spring and autumn for all the genera (the relationship between the type of infection, single vs. multiple, and time of year, Haemoproteus: χ² = 7.7, df = 2, p = 0.021; Plasmodium: χ² = 11.18, df = 2, p < 0.01; Leucocytozoon: χ² = 13.23, df = 2, p < 0.01).

We identified the following lineage combination for intra-generic multiple infections: SYAT02 + SYAT28 (n = 2), SYAT02 + SYAT07 (1), SYAT07 + SYAT12 (1), SYAT16 + SYATSYAT44 (1), SYAT16 + SYBOR35 (1), SYAT01 + SYAT17 (1), and SYAT03 + TURDUS2 (1) for Haemoproteus spp. and TUPHI08 + LINN1 (1) and BT7 + TURDUS1 (1) for Plasmodium spp. We were not able to determine more than one lineage involved in intra-generic multiple infections for Leucocytozoon.

Inter-generic multiple infections involving Haemoproteus-Leucocytozoon and Plasmodium-Leucocytozoon lineages also were detected in our study (

Table 4. Relationship between infection type, considering intra- vs. inter-generic infections, and time of year for the frequency of haemosporidian cytochrome b lineages detected in Drienovec, Slovakia, during 2017, 2018 and 2019. Note that haemosporidian lineage identity was not known for all parasites involved in multiple infections.

	Infection Type
Time of Year	Intra-Generic	with Leucocytozoon
Haemoproteus
spring	90	12
summer	87	10
autumn	100	11
Plasmodium
spring	30	5
summer	17	5
autumn	155	43
	Intra-generic	with Haemoproteus
Leucocytozoon
spring	26	12
summer	10	10
autumn	62	11
	Intra-generic	with Plasmodium
Leucocytozoon
spring	33	5
summer	15	5
autumn	30	43

). The relationship between infection type (intra-generic vs. inter-generic) and time of year was insignificant for Haemoproteus and Plasmodium lineages (proportions of Haemoproteus infections with or without Leucocytozoon lineages: χ² = 0.28, df = 2, p = 0.87; proportions of Plasmodium infections with or without Leucocytozoon lineages: χ² = 3.08, df = 2, p = 0.21). However, Leucocytozoon infections coinfected with Haemoproteus lineages occurred more often than expected in summer (proportions of Leucocytozoon infections with or without Haemoproteus lineages: χ² = 28, df = 2, p < 0.001) and with Plasmodium lineages more often than expected in autumn (proportions of Leucocytozoon infections with or without Plasmodium lineages: χ² = 52.04, df = 2, p < 0.001).

3.2. Seasonal Variation in Infection Intensity and Prevalence

Following the conservative criteria to accept the data of qPCR (see Methods), we were able to determine the intensity of haemosporidian infection for 372 (49%) of the 762 haemosporidian-positive samples. The intensity of infection ranged between 0.002 and 199 DNA copies/100 red blood cells (RBC), but it was highly skewed towards zero (median and 25–75% quantiles: 1.36, 0.09–9.04). Controlling for the random effects of avian species identity, year, and the phylogenetic relationships between avian species, we found that infection intensity at the level of the local avian community had a quadratic relationship with date (

Table 5. Bayesian multivariate multilevel model on the predictors of haemosporidian infection intensity in a local avian community in Drienovec, Slovakia, during 2017, 2018 and 2019.

Parameter	Estimate	Est. Error	L-95% CI	U-95% CI	Rhat	Bulk ESS	Tail ESS	pMCMC
Random effects (SD)
Phylogeny (35 levels)	0.05	0.03	0.01	0.11	1	1126	1665
Year (3 levels)	0.55	0.62	0.03	2.2	1	784	633
Year:Species (67 levels)	0.24	0.15	0.01	0.55	1	812	3163
Fixed effects
Intercept	−0.7	0.68	−2.05	0.65	1	2553	3620	0.134
Date	0.02	0.01	0	0.03	1	7369	21,192	0.013
Date^2	~0	~0	~0	~0	1	6810	20,008	0.006
Date/species	~0	~0	~0	~0	1	72,435	103,455	0.466
Age_HY	−0.25	0.16	−0.56	0.06	1	19,968	51,277	0.057
Residual effect (SD)
Sigma	1.18	0.05	1.09	1.28	1	15,783	48,731

). Infection intensity peaked during summer around mid-July (Figure 5A). Infection intensity in older (AHY) birds was almost twice as high as in younger (HY) birds (back-transformed estimated marginal means of the number of DNA copies/100 RBC: AHY—1.75; HY—0.98), and the analysis detected a strong trend for the effect of age (Table 5). Although the random effects were small, they were important based on 95% credible intervals (CI; Table 5).

In contrast to infection intensity, infection prevalence at the level of the local avian community reached its lowest values in summer (Figure 5B). Although infection prevalence in spring and autumn was comparable (back-transformed estimated marginal means: spring—42.8%, summer—32.1%, autumn—48.8%), prevalence was only significantly different between summer and autumn (one-sided hypothesis tests: spring vs. summer: posterior p = 0.09; summer vs. autumn: posterior p = 0.01). Random effects of avian species identity, year, and avian phylogeny affected infection prevalence based on 95% CI (

Table 6. Bayesian multivariate multilevel model on the predictors of haemosporidian infection prevalence in a local avian community in Drienovec, Slovakia, during 2017, 2018 and 2019.

Parameter	Estimate	Est. Error	L-95% CI	U-95% CI	Rhat	Bulk ESS	Tail ESS	pMCMC
Random effects (SD)
Phylogeny (58 levels)	0.12	0.05	0.04	0.25	1	38,921	35,894
Year (3 levels)	0.72	0.72	0.03	2.69	1	27,292	18,348
Year:Species (130 levels)	0.23	0.19	0.01	0.71	1	48,274	57,250
Fixed effects
Intercept	−0.06	0.95	−1.99	1.9	1	34,464	26,486	0.469
Season_Spring	−0.29	0.41	−1.11	0.51	1	80,281	68,270	0.242
Season_Summer	−0.83	0.38	−1.59	−0.09	1	103,349	80,439	0.014

).

3.3. Host-Parasitic Associations

The highest (>50%) infection prevalence in frequently sampled (n > 30) bird species was observed in song thrush Turdus philomelos (78%), Eurasian blackbird T. merula (67%), Eurasian blue tit Cyanistes caeruleus (64%), common chaffinch Fringilla coelebs (61%), Eurasian blackcap Sylvia atricapilla (60%), hawfinch C. coccothraustes (60%), and great tit Parus major (53%; Table 1).

Examining the structural parasite specificity in terms of the d′ index, and considering only species with >10 positive samples, hawfinch, song thrush, Eurasian bullfinch P. pyrrhula, lesser whitethroat Sylvia curruca, common chaffinch, and Eurasian blackbird can be considered to be primarily infected by Haemoproteus and Leucocytozoon, while the Eurasian blackbird also shows a tendency to be relatively highly infected by Plasmodium lineages (Table 1).

In line with the results for avian hosts, the structural host specificity in terms of the H₂ index indicated that at the parasite community level, Haemoproteus and Leucocytozoon lineages detected in our study in general represented intermediate host specialists, while Plasmodium lineages in general represented host generalists (

Table 7. List of haemosporidian cytochrome b lineages detected in Drienovec, Slovakia, during 2017, 2018 and 2019. For each lineage, we report values of structural host specificity in terms of the d′ and H₂ indices, phylogenetic specificity in terms of the SES_MPD index. The indices were calculated at the level of haemosporidian lineage (d′, SES_MPD) or genus (H₂). Significant values (p ≤ 0.05) for the SES_MPD index are denoted by asterisk. New lineages are in bold.

	d′	SESMPD	H2		d′	SESMPD	H2		d′	SESMPD	H2
Haemoproteus			0.70	Plasmodium			0.31	Leucocytozoon			0.69
ARW1	0.57	−1.67		AFTRU4	0.72	-		BRAM3	0.48	-
CIRCUM05	0.59	-		BT7	0.19	0.51		BT1	0.25	−0.25
CCF23	0.42	-		COLL1	0.27	0.54		BT2	0.16	−0.15
CCF2	0.42			GRW06	0.18	-		COCOR02	0.59	-
CCF3	0.42	-		GRW11	0.20	−2.64 *		EUSE2	0.38	-
CCF6	0.49	-		LINN1	0.18	−1.87 *		GAGLA06	0.64	-
COLL2	0.49	-		PADOM02	0.26	-		HAWF7	0.46	-
COLL3	0.43	-		SGS1	0.12	−1.00		NEVE01	0.31	-
CWT2	0.58	-		SW5	1.00	-		PARUS14	0.19	-
CWT4	0.50	-		SYAT05	0.26	−1.00		PARUS15	0.19	-
DUNNO01	0.27	-		SYBOR2	0.16	0.40		PARUS16	0.12	−1.06
EMCIR01	0.52	-		TUPHI08	0.26	-		PARUS18	0.19	-
GAGLA02	0.59	-		TURDUS1	0.13	0.26		PARUS19	0.13	-
HAWF1	0.43	-		TURMER09	0.31	-		PARUS20	0.19	−1.69 *
HAWF2	0.45	-						PARUS22	0.28	−2.99 *
HAWF6	0.46	-						PARUS25	0.29	−0.09
LWT1	0.40	-						PARUS4	0.17	−1.06
PARUS1	0.13	−3.00 *						PARUS7	0.13	-
PHSIB1	0.11	-						PARUS81	0.13	-
PHSIB2	0.53	0.34						PERATE06	0.61	-
RBS2	0.76	-						PYRPYR01	0.38	-
ROBIN1	0.11	−0.18						PRUMOD01	0.44	-
SISKIN1	0.44	-						PICVIR01	0.38	-
SYAT01	0.18	−3.35 *						REB11	0.40	-
SYAT02	0.18	0.37						RECOB3	0.40	-
SYAT03	0.07	−1.09						SANG02	0.38	-
SYAT07	0.16	-						SFC8	0.14	1.77
SYAT10	0.00	-						STUR1	0.44	-
SYAT11	0.10	-						SYAT22	0.10	-
SYAT12	0.10	-						SYBOR23	0.46	-
SYAT13	0.16	-						SYCON05	0.00	-
SYAT14	0.00	-						TRPIP1	0.65
SYAT16	0.17	-						TUMER03	0.31
SYAT17	0.00	-						TUPHI06	0.47
SYAT28	0.12	-						TUPHI13	0.46
SYAT33	0.00	-						TUMER10	0.23
SYAT44	0.13	-						TURMER10	0.23
SYBOR35	0.00	-						TURMER11	0.23
TUPHI01	0.47	-						TUPHI14	0.38
TURDUS2	0.27	−0.73						PARUS93	0.13
WW1	0.50	-						TURMER12	0.23
WW2	0.20	−3.35 *						TURMER13	0.23
YWT2	0.65	-						TURMER14	0.23
ERIRUB02	0.01							PARUS94	0.13
PRUMOD04	0.27
TURMER08	0.23
ERIRUB03	0.01
CYACAE09	0.19
CYACAE08	0.19
SYAT50	0.00
SYAT51	0.00
SYAT52	0.00
SYCUR03	0.40

). When focusing on individual lineages, we revealed a considerable variation in structural host specificity between lineages within all the three genera (Table 7). Interestingly, while some lineages such as Haemoproteus lineages SYAT01 and WW2 (associated with Sylviidae), Plasmodium lineage LINN1 (associated mainly with Turdidae and Muscicapidae), and Leucocytozoon lineages PARUS20 and PARUS22 (associated mainly with Paridae) showed a low structural specificity, they showed a significant phylogenetic host specificity (Table 1; Figure 6).

4. Discussion

In this study, we examined seasonal variation in local diversity and infection patterns of haemosporidian parasites for an avian community in a natural woody wetland in Slovakia. For this community, represented almost entirely (>99%) by passerines, with dynamic temporal changes in its structure and size, we detected 112 haemosporidian cytochrome b lineages, including 19 novel ones, in 35 avian species. The highest numbers of lineages were detected during spring and autumn, which corresponded with higher avian species richness and infection prevalence in the avian community during these periods. Conversely, infection intensity peaked during summer and tended to be higher in older birds.

The highest haemosporidian infection prevalence for commonly sampled bird species was observed in song thrush, Eurasian blackbird, Eurasian blue tit, common chaffinch, Eurasian blackcap, hawfinch, and great Tit. These thrushes, tits, finches, and a Sylvia warbler represent both locally and regionally abundant bird species. Our study, therefore, accords with the expectation that infection prevalence is a positive function of host abundance [4,62,63]. Nevertheless, despite their higher abundances and sampling rates, haemosporidian prevalence was relatively low in European robin Erithacus rubecula and common chiffchaff Phylloscopus collybita (see [64,65] for similar results), and their associated haemosporidians mostly represented structural and phylogenetic generalist lineages. Finally, even though multiple haemosporidian lineages showed similar prevalence in different hosts within the local avian community, exhibiting low structural specificity, some of these lineages were primarily found in phylogenetically related hosts: Haemoproteus lineages SYAT01 and WW2 in Sylviid warblers, Plasmodium lineage LINN1 in thrushes and old world flycatchers, and Leucocytozoon lineages PARUS20 and PARUS22 in tits.

4.1. Haemosporidian Diversity

Overall, Haemoproteus lineages dominated in the local avian haemosporidian community. Although the composition of Haemoproteus lineages fluctuated between spring and autumn, this genus showed similar species richness and diversity during spring, summer, and autumn periods. In turn, while Leucocytozoon lineages showed a comparable overall richness as those of the genus Haemoproteus, the richness and diversity of Leucocytozoon lineages dropped in summer, and Leucocytozoon parasites did not represent abundant avian haemosporidians in the locally breeding bird species. Lineages of the genus Plasmodium showed the lowest richness and diversity, even though the richness of lineages of this genus increased in autumn.

The highest haemosporidian lineage diversity was observed in spring, which was mainly due to the high diversity of Leucocytozoon lineages. The most common Leucocytozoon lineages detected in spring were HAWF7 and SFC8 associated with hawfinch and European robin, respectively. Valkiūnas [66] found that in spring, after arrival to the breeding grounds in North Europe, Leucocytozoon infections predominate in middle-distance migratory birds wintering in the Mediterranean region, while Haemoproteus infections are typical for long-distance (trans-Saharan) migrants wintering in Africa. Our study supports these findings, because bird species with a high infection prevalence in spring comprise short-distance migrants wintering in West and South Europe [67,68,69]. Furthermore, the seasonal dynamics of infections involving Leucocytozoon lineages suggests that the local parasite community for this genus is mainly sustained by birds breeding in North Europe and wintering in South Europe.

The lineages of the genus Plasmodium were least represented in the local avian community, though their abundance, but not diversity, increased in autumn. This result accords with previous work indicating an environmental barrier for more intensive transmission of avian malaria parasites in northern latitudes, with successful parasite transmission in these regions being mostly restricted to the warmest period of the year during late summer/early autumn [70]. Yet, generalist Plasmodium lineages TURDUS1 and SGS1 were fairly common in spring in several short-distance migrants: Eurasian wren T. troglodytes, song thrush, European robin, yellowhammer Emberiza citrinella, and common chaffinch. As these species belong to the typical early-season singing birds and TURDUS1 and SGS1 lineages were also most abundant in autumn, it is possible that these spring infections are caused by spring relapses in immunologically supressed, reproductively active, birds, e.g., [71,72]. It is important to note that Plasmodium infections could be underestimated in our study if the nested PCR procedure favoured Haemoproteus parasites.

4.2. Haemosporidian Multiple Infections

We found that intra-generic multiple infections were relatively more common in summer than spring or autumn, but this pattern was significant only for Leucocytozoon parasites. Additionally, samples infected with Leucocytozoon parasites were more likely found to be coinfected with Haemoproteus parasites in summer and with Plasmodium parasites in autumn. These results suggest that birds infected with Leucocytozoon parasites have a higher chance of being coinfected with other Leucocytozoon lineages and with Haemoproteus parasites during the summer period, and with Plasmodium parasites during autumn. Experimental inoculations of sporozoites of different genera, e.g., using sporozoites administration in a different sequence, could be performed to test this hypothesis. Importantly, additional information about dipteran vectors is needed to see how the ecology of vectors and bird hosts contribute to specific types of multiple infections. Finally, different sets of primers should be used in the future to study multiple infections, because many samples interpreted as single infections may readily involve multiple infections [73].

4.3. Haemosporidian Infection Prevalence and Intensity

Our results on infection prevalence at the community level support the functional model proposed by Beaudoin et al. [8]. Namely, we found elevated parasite prevalence in spring and autumn, even though the peak in autumn was more pronounced than that in spring. This result may indicate that the period of summer (June–July), when most of the bird species are still reproductively active at the study site, is relatively less important for parasite circulation and avian host health. This picture, however, is deceiving because infection intensity was found to peak in summer. Therefore, even though relatively fewer birds were infected with haemosporidians in summer at the community level, infected birds showed higher parasite intensity in this period. Moreover, older (after-hatch-year) birds tended to show higher parasite intensities than younger (hatch-year) birds. These results suggest that elevated parasite intensities in summer can be related to supressed immune functions in certain reproductively active birds [21]. As specific multiple infections occur more likely in summer (see Section 4.2), dominantly involving bird hosts primarily infected by Haemoproteus and Leucocytozoon parasites, it is possible that these multiple infections can also be responsible for increased parasite intensities in the avian community in summer. Importantly, these results on parasite prevalence and intensity imply that haemosporidian parasite transmission and circulation can locally be maintained by relatively few super-hosts, i.e., by birds with high infection intensities in sparsely infected avian communities.

In conclusion, this is one of the few studies examining temporal changes in the haemosporidian population structure at the avian community level. We found that a local avian community in Slovakia was characterised by Haemoproteus lineages, with Leucocytozoon lineages enriching the haemosporidian community mainly in spring after winter migration and, together with Plasmodium lineages, in autumn during autumn migration. In accordance with previous studies, despite the mixture of the locally breeding birds with spring and autumn migrants, the haemosporidian community in the local breeders in summer was relatively stable [66] and only few parasite lineages were found to dominate within each genus (Figure 4). We found that while several passerines in the local avian community were readily infected by various Haemoproteus and Leucocytozoon lineages, only the Eurasian blackbird was readily infected by parasite lineages of any of the three genera. These results deserve attention because the Eurasian blackbird is a suitable model system for diverse questions on the ecology and evolution of vector-borne pathogens, e.g., [35,74], but baseline data about infection patterns in this species are largely lacking for natural habitats (but see [26,75] for haemosporidians).