Comparative demography elucidates the longevity of parasitic and symbiotic relationships

Parasitic and symbiotic relationships govern vast nutrient and energy flows, yet controversy surrounds their longevity. Enduring relationships may engender parallel phylogenies among hosts and parasites, but so may ephemeral relationships when parasites colonize related hosts. An understanding of whether symbiont and host populations have grown and contracted in concert would be useful when considering the temporal durability of these relationships. Here, we devised methods to compare demographic histories derived from genomic data. We compared the historical growth of the agent of severe human malaria, Plasmodium falciparum, and its mosquito vector, Anopheles gambiae, to human and primate histories, thereby discerning long-term parallels and anthropogenic population explosions. The growth history of Trichinella spiralis, a zoonotic parasite disseminated by swine, proved regionally specific, paralleling distinctive growth histories for wild boar in Asia and Europe. Parallel histories were inferred for an anemone and its algal symbiont (Exaiptasia pallida and Symbiodinium minutum). Concerted growth in potatoes and the agent of potato blight (Solanum tuberosum and Phytophthora infestans) did not commence until the age of potato domestication. Through these examples, we illustrate the utility of comparative historical demography as a new exploratory tool by which to interrogate the origins and durability of myriad ecological relationships. To facilitate future use of this approach, we introduce a tool called C-PSMC to align and evaluate the similarity of demographic history curves.

Figures referred to in Methods: Figure S1. Comparison of minimum average slope difference for all considered host-symbiont pairs. Bars representing the fit of biologically plausible hosts are highlighted, while apocryphal matches are faded, but shown for context. In all cases, the true host is preferred from the arena of plausible hosts (e.g. European boar > Asian boar for European T. spiralis), but in some cases an apocryphal match scores better (e.g. European boar > potato for P. infestans).
Supplemental Table 3: Range of plausible generation times (in years) assumed for each symbiont genus.

Effect of variant-calling methodology on PSMC
We tested the sensitivity of PSMC to alternative variant-calling methods, which may influence the number and distribution of heterozygotic base calls in the consensus sequence. We observed little difference in shape or scaling between PSMC plots of T. spiralis isolates based on heterozygote-calling performed, alternatively, by Geneious or the mpileup function of Samtools (samtools mpileup -uf Genome.fasta test.sorted.bam | bcftools view -cg -| vcfutils.pl vcf2fq > consensus.fq) (Supplemental Figure S2).
The heterozygote frequencies found in the PSMC input files (.psmcfa) generated from these consensus fasta files were also similar irrespective of the variant-calling program, with the Geneious consensuses generally having a greater proportion of missing data, suggesting the filtering thresholds we applied (75% consensus, >5x coverage) were slightly more stringent than the Samtools defaults (Supplemental Table 4). Nadachowska-Brzyska et al. (2016) [1] reported that increasing their per-site filtering thresholds shifted their PSMC curve toward the present, and this effect is also apparent to a minor degree in our plots (Supplemental Fig. S2). The same authors found that more stringent filtering helped to resolve population fluctuations in PSMC plots, though the significance of this faded for average read depths greater than ~18x.
Supplementary  Figure S2. PSMC plots based on consensus sequences generated by Geneious (orange) and Samtools (grey), for two T. spiralis isolates from America (left) and two from China (right) In general, joint histories of species pairs could be traced from about 500,000 years ago (where years result from joint estimates of mutation rate and generation time) until a few thousand years ago. In considering the evolutionary longevity of particular parasitic or symbiotic relationships, our methods appear well-suited to that interval (and no longer). It was often the case that one, but not both, genomes in a given pair provided temporally deeper inferences on demographic history.

Malaria system
The origin and expansion of the agent of malaria, P. falciparum, is tied to human demographic history [2,3]. The most recent evidence indicates that P. falciparum shifted from gorillas to humans between 60,000 and 130,000 years ago [4], though chimpanzees were previously suspected as the proximate source [5,6]. When humans began clearing land for agriculture, anopheline mosquito numbers increased accordingly, resulting in increased malaria transmission and presumably increased P. falciparum population sizes. Therefore, we sought to examine the demographic changes written in the genomes of hosts, vectors, and parasites to confirm our understanding of falciparum malaria in humans.
Plasmodium falciparum is only briefly diploid in its anopheline mosquito definitive hosts, and to our knowledge has not been sequenced in that life-history stage. Therefore, we mapped the distribution of potentially heterozygous positions by creating "pseudodiploid" genomes in silico by mapping short reads from two haploid sequencing projects derived from P. falciparum isolated in the same geographic region to the P. falciparum 3D7 chromosomes (the most complete representation of the P. falciparum genome to date). We took these assemblies to represent the joining of haploid gametes into transient diploid ookinetes. We selected sequencing projects with approximately equivalent numbers of reads in order that heterozygous bases could be determined based on read depths as above, assuming comparable coverage of chromosomes by each deeply sequenced genome. The resulting consensus chromosomes were used for downstream analysis, including estimates of relative generation time (Supplemental Figure S3) that well match empirical estimates based on clinical and entomological observation. Producing pseudodiploids in this way, we note, could render any haploid organisms amenable to PSMC analysis.
Plasmodium falciparum notably parallels the bottleneck evident for human beings and their subsequent ascendency in the last ten thousand years; by comparison, other great apes have experienced stable or declining effective population sizes since the Neolithic (Supplemental Figure S4). Similar demographic patterns have been inferred for primates circa 100 kya-1 Mya, and would be expected to eventually converge on the shared demography of our common ancestors; subsequent to divergence, it would not be surprising to find that various primate species were subjected to common beneficial and injurious environmental conditions. Figure S3. PSMC curve-fit to French Human as a function of generation time for A.gambiae (green) and P. falciparum (purple), assuming a generation time of 25 years for humans and identical per-generation mutation rates of 2.5e -8 . The lowest average slope difference was found at 1.43 years/generation for A.gambiae and 0.48 years/generation (~23 weeks/generation) for P. falciparum. However, the slower Drosophila mutation rate (1.1e -9 ) may be more appropriate for mosquitoes [7], which would imply an A. gambiae generation time of 0.057 years, or approximately one generation every 3 weeks, in line with empirical estimates based on knowledge of mosquito ecology.
Chimpanzees had previously been considered the penultimate host for the ancestors of falciparum malaria until more closely related parasites were discovered in gorillas.  Figure S4). The demographic history of P. falciparum tracks human demography best, especially during the most recent historical interval, when human populations grew exponentially while other great apes did not. However, of the considered chimpanzee populations, P.t. schweinfurthii was the most consistent with human and malaria demography in the range of 80K -1M years ago. Figure S4. Estimated demographic histories for P. falciparum falciparum (violet) and A. gambiae gambiae (green) in relation to that for human (black) and gorilla (brown) and multiple subspecies of chimpanzee (blue). Gorilla populations grew and shrank in parallel with humans, Anophelines, and P. falciparum, with steep declines from 200,000 ybp to 70,000 ybp, followed by a small stable population between 70,000 and 11,000 ybp. From ~100-11 kya, gorilla populations declined in concert; prior to then, their populations were notably larger and since then notably smaller. P. falciparum and mosquito closely parallel human demography for a considerable period, including growth in the most recent past unique, among primates, to human populations. Recent growth in the mosquito population, while present, is much less dramatic than that observed in the human and parasite populations. Note that human and gorilla effective population sizes were uniformly scaled up by a factor of 10 to be more easily visible in this figure.

Trichinella spiralis system
Trichinella spiralis is contracted by the consumption of meat in which larval parasites have encysted. Adult pairs mate within weeks, and disseminate new larvae to the tissues whereupon the next generation of larvae establish chronic infections that endure for the life of the host. Thus, infections are acquired after weaning and are transmitted at death. The mean generation time of the parasite would be shortened to the extent that infection reduced swine longevity and to the extent that the parasite also exploits shorter-lived hosts, such as rats. The curve-fitting procedures described above suggested ~ 4 parasite generations per wild boar generation (or 1.1 years, assuming 5-year generations for wild boar).
The inferred relative generation time of T. spiralis was inconsistent in Asia and in Europe (Supplemental Figure S5), probably owing to limited complexity and depth in the demographic history reconstructions of European T. spiralis. Since domestication, pigs have undergone a population explosion, and our analysis illustrate commensurate exponential growth in the size of the parasite population (Supplemental Figure S6). Wild boar have interbred with domesticated pigs, leaving tracts of homozygosity in wild boar genomes [9] that may inflate estimates of recent population growth in wild boar. Figure S5. Fit of European (blue) and Chinese (red) T. spiralis to sympatric S. scrofa as a function of their generation time. Given 5 years per boar generation, the lowest average slope difference between the average of the two regional fits (grey line) stabilizes at approximately 1.1 years per T. spiralis generation, which is the regional optimum for the Asian T. spiralis. Because the Asian genomes record a deeper and more complex demographic history than the European ones, we elected to proceed with 1.1 years/generation as the species optimum generation time for T. spiralis. Note that although the Asian T. spiralis slope difference is minimized at <6 weeks, this biologically implausible result derives from few overlapping time intervals (between parasite and host).
As stated, the demographic plots for western T. spiralis were generally lacking in complexity, depriving curve-fitting procedures of information needed to achieve precise estimates (Supplemental Figure S5). Nonetheless, their optimal fit is consistent with a decline in pig populations in both Europe and Asia that has previously been interpreted as having been caused by the LGM [10]. The inferred histories of Asian T. spiralis reach deeper in time and show greater demographic variability over time. Most notably, in contrast to European T. spiralis and to both European and Asian S. scrofa, the effective population size of Asian T. spiralis does not decline over the last ~50,000 years, but stabilizes and grows. All pigs and T. spiralis eventually gave way to dramatic growth over the Holocene.

Symbiodinium system
A long-term demographic association between the anemone Aiptasia pallida and its symbiont Symbiodinium minutum might not have been expected owing to expulsion and die-offs of symbionts during periods of environmental stress [11]. Cnidarians, including anemones and corals, sometimes acquire new species after heat-induced bleaching [11]. Although such associations might therefore prove ephemeral, the population of S. minimum tracked that of A. pallida remarkably closely (Supplemental Figure S8). The demographic histories of all three species were similar until approximately 10-20 kya, following the LGM, when the coral population began a period of extreme growth while the anemone and symbiont populations began to decline.
We employed equivalent curve-fitting procedures so as to determine the best-possible demographic match over time, and derived a relative generation time of ~0.5 against the actual host, Aiptasia pallida, and 0.97 against the apocryphal host, Acropora digitifera (Supplemental Figure S8). Notably, the estimated slope difference between host and symbiont demographic curves was less for the natural pair at almost every considered value of relative generation time (Supplemental Figure S7); the average slope difference between S. minutum and that of its actual host (~.09) was notably less than to that of its apocryphal host (~0.17) when optimally scaled to each. Supplemental Figure S8 depicts the histories of all three species, assuming a relative generation time of 0.52 for S. minutum (fit to Aiptasia pallida).
Their shared history in the Red Sea may partially explain the evidently shared demographic histories of the anemone and its symbiont (in contrast to that of the coral, obtained from Okinawa in the East China Sea). Although oceanic coral systems generally prosper when the global climate cools, glacier formation during the last glacial maximum resulted in a 120m drop in sea level, halving the surface area of the Red Sea, significantly reducing exchange with the Indian Ocean through the Strait of Bab el Mandab, and increasing the sea's salinity by about 50% [12,13]. These events would likely have reduced the census size of the anemone population, as well as fragmenting it within the Red Sea and between the sea and ocean, further reducing the effective population size of the anemone. Smaller, fragmented host populations would likely have constrained growth of the symbiont. Subsequent rising seas may have engendered population growth of the anemone and its algal symbiont owing to increased habitat availability and increased gene flow among once-separated populations (Supplemental Figure  S8). Figure S7. The slope difference between S. minutum and its natural host (Aiptasia pallida; blue line) was less than to another marine invertebrate, the coral Acropora digitifera (green line). The fit to Aiptasia pallida was optimized at a relative generation time of approximately half (0.52x) that of the host. Optimal fitting to Acropora digitifera would require almost identical generation times for host and symbiont (0.97x). Figure S8. Demographic histories of coral (Acropora digitifera), anemone (Aiptasia pallida), and the anemone's symbiont (S. minutum), assuming annual generations for the coral and anemone and biannual generations for the symbiont. Rises and falls of S. minutum (dotted blue line) approximate those for its extant host, Aiptasia pallida (solid blue line) for approximately 1 million years. A similar pattern of growth and decline is observed, albeit out of sync, for Acropora digitifera, but the demographic trajectory of Acropora digitifera separates in the last 20,000 years, experiencing rapid population growth while the anemone and its symbiont decline and stabilize.

Potato system
Infections of potatoes (Solanum tuberosum) with the potato blight Phytophthora infestans caused one of the greatest human disasters in agricultural history, resulting in the death or emigration of 20% of the Irish in the mid-1800's. Potatoes originated in the Andes and may have been domesticated twice. Biogeographic data suggest the central Mexican highlands as the ancestral home of P. infestans, achieving global reach only after the potato domestication in the Neolithic and the globalization of agriculture beginning 500 years ago. P. infestans appears restricted to tuber-forming species of Solanum and does not easily cross with other species of Phytophthora. The closest known relatives to P. infestans (P. andina, P. ipomoea, P. mirabilis, and P. phaseoli) are also pathogens of plants native to the Neotropical highlands. However, other closely related congeners (P. iranica, P. clandestina, P. tentacultate) afflict diverse plant types (eggplant, clover, chrysanthemum) in diverse regions (Asia, Europe); an even greater diversity of host types and biogeographic regions characterize more inclusive groups of Phytophthora, including species infecting orchids native to Indonesia, lilacs in the Balkans, rhododendron native to the Himalayas, raspberry native to Northern Europe and Asia, and Douglas fir in Western North America. This broader evidence for host switching provides valuable context to our results, suggesting parallel demography only during the most recent past.
Prior to potato domestication and cultivation, potatoes and the agent of potato blight shared no obvious episodes of concerted population growth or contraction. Lacking such compelling information, and with a relatively short and simple reconstructed history of P. infestans to overlap, our curve-fitting algorithm easily identified a clear preferred relative generation time for the two at ~0.14 P. infestans generations per potato generation (which was assumed to be one year) (Supplemental Figure S9). Adopting this temporal scalar resulted in an average slope difference (0.38 on a scale of 0-2) that was as poor as the best fit between apocryphal species pairs (such as wild boar), even with relatively tight constraints on plausible generation times for the pathogen (<3 months [14]). Among poor alternatives, the best-fit curves register by reference only to the recent (post-domestication) interval of population growth in potatoes. This is preceded by ~20,000 years of decline in potatoes and relative stasis in P. infestans (Supplemental Figure S10). These findings are consistent with a close affinity of P. infestans with potatoes only since their domestication. Figure S9. Fit of P. infestans demographic history to that of potatoes as a function of relative generation time. The optimum is reached when a single potato generation equals ~10 generations of P. infestans. This short generation time is based on concerted growth in the host and pathogen only since the domestication of potatoes in the last ~8,000 years ( Figure S10), suggesting that shared demography may have been limited to this recent temporal interval. Figure S10. Contrasting demographic histories of potato and P. infestans the agent of potato blight excepting since domestication. The only demographic feature indicating shared history in this plant/pathogen system relates to recent growth in each. Here, the demographic history of potato (solid line) is depicted against the history of the fungal pathogen assuming the optimal relative generation time, which does not convincingly match the history of potatoes prior to potato domestication.
Noting an ancient tripling of chromosomal content in dicot plants [15], we were careful to assess whether, how, and to what extent demographic plots might be biased if reads derived from paralogs in the potato genome were mismapped to the same position, thereby inflating estimates of heterozygosity. We examined this potential source of bias by filtering regions of the potato genome sequenced to especially high coverage (arguably over-represented owing to their occurrence as multiple copies). Prior to filtering, the potato genome dataset had an average sequencing depth of 31.8 reads per site with a large standard deviation (148.5). After filtering, the mean sequencing depth and standard deviation was reduced to 12.4 and 13.3, respectively. Filtering resulted in removal of approximately 1/5 of each scaffold for consideration by PSMC. After filtering, about 40 million bases of contiguous sequence in 10 scaffolds were available for analysis. Filtering high coverage areas reduced the overall heterozygosity in the potato genome from 1.19% to 0.895%. It also reduced the length of the scaffolds used for secondary PSMC analysis from 43,353,824 bp to 36,083,464 bp (17% removed from consideration). Confirming simulations entertained by Li and Durbin (2011) [16], removing suspected paralogous regions did not affect the overall shape of curve, but did reduce the estimate of the most ancient populations (Supplemental Figure S11). The essential demographic signal to remain stable, further substantiating the conclusion that potatoes and the fungal agent of potato blight did not experience parallel demographic histories until the age of potato domestication, beginning ~9,000 years ago. Figure S11. Effect of removing suspected paralogs from the potato genome prior to estimating its demographic history using PSMC. The PSMC plot derived from the potato genome assembly (dotted line) was compared to an equivalent analysis performed on a dataset from which anomalously deeply sequenced regions were removed (solid line), owing to the suspicion that paralogous gene copies derived from genome duplications might be inflating estimates of heterozygosity. As predicted (Li and Durbin, 2011) [16], removing these sequences did reduce estimates of ancient population size. However, this quality control step did not fundamentally alter the inference as to when ancestral potato populations grew or contracted, and specifically had no appreciable effect on demographic estimates within the last 30,000 years.