Here, we present the findings of a phylogenomic investigation of hallucinogenic mushroom genomes. We sequenced three known PS + mushrooms ( Psilocybe cyanescens , Gymnopilus dilepis , and Panaeolus cyanescens ) representing diverse lineages, searched for clustered genes that were simultaneously associated with functions and phylogenetics of PS production, and converged on a single gene cluster that shows signatures of HGT among the sequenced fungi (Fig. 1 B and C). We then confirmed the chemical function and specificity of the presumed first step in the PS pathway (conversion of tryptophan to tryptamine). We assessed further evidence of HGT within a dung environment, and investigated the ecological trends in Agaricales genome content, to elucidate the ecological context in which the PS cluster transfer may have taken place. Together, our work suggests that shared environmental selection pressures may have favored the transfer of the PS gene cluster among hallucinogenic fungi, and provides a methodological roadmap for the future discovery of novel fungal pharmaceuticals.

A common feature of fungal secondary metabolite biosynthesis is the organization of most or all of the required anabolic, transport, and regulatory genes in gene clusters. Gene clusters are often discontinuously distributed among fungal taxa, partly due to horizontal gene transfer (HGT) among species with overlapping ecological niches (Gluck‐Thaler and Slot 2015 ). The limited phylogenetic distribution of psilocybin (Fig. 1 A), coupled with the requirement for multiple enzymatic steps for its biosynthesis (in order: tryptophan‐decarboxylation, tryptamine‐4‐hydroxylation, 4‐hydroxytryptamine O‐phosphorylation, and N ‐methylation; Fricke et al. 2017 ) suggested that the psilocybin pathway might have dispersed via HGT of a gene cluster. We therefore predicted that the genetic mechanism for psilocybin biosynthesis would be identified in searches for gene clusters with a common phylogenetic history and distribution restricted to psilocybin producing (PS + ) mushrooms.

Rationale and approach to search for psilocybin gene clusters. (A) The patchy distribution of psilocybin in Agaricales from different fungal lifestyles . Psilocybin (PS) production has a limited phylogenetic distribution among Agaricales fungi, and is associated with ecological lifestyles (dung decay [DG], wood decay [WD], and mycorrhizal [MR]) with similarly spotty distributions. (B) Pipeline for psilocybin gene cluster discovery . Protein models from PS + (blue shading) and matched PS − (gray shading) genomes (Psicy2 = Psilocybe cyanescens ; Pancy2 = Panaeolus cyanescens; Gymlu2 = Gymnopilus dilepis ; Gymch1 = Gymnopilus chrysopellus ; Galma1 = Galerina marginata ; Hypsu1 = Hypholoma sublateritium ) were sorted into psilocybin discovery homolog groups (PDHs), which were filtered by taxon representation. Gymch1 was retroactively inferred PS + (dotted gray outline). Four taxon filters allow for 0, 1, 2, or 3 PS − genomes per PDH. Taxonomically filtered PDHs contained few candidates for decarboxylation, methylation, oxygenation, phosphorylation, and transport‐related functions. A single cluster of taxonomically filtered PDHs was identified in each of four genomes (two partial clusters in Pa. cyanescens ), and corresponded to expected functions identified in the PS + +1 PDH set.

Secondary metabolites are small molecules that are widely employed in defense, competition, and signaling among organisms (Raguso et al. 2015 ). Due to their physiological activities, secondary metabolites have been adopted by both ancient and modern human societies as medical, spiritual, or recreational drugs. Psilocin is a psychoactive agonist of the serotonin (5‐hydroxytryptamine, 5‐HT) ‐2A receptor (Halberstadt and Geyer 2011 ) and is produced as the phosphorylated prodrug psilocybin by a restricted number of phylogenetically disjunct mushroom forming families of the Agaricales (Bolbitiaceae, Inocybaceae, Hymenogastraceae, Pluteaceae, Fig. 1 A) (Allen 2010 ; Dinis‐Oliveira 2017 ). Hallucinogenic mushrooms have a long history of religious use, particularly in Mesoamerica, and were a catalyst of cultural revolution in the West in the mid‐20th century (Nyberg 1992 ; Letcher 2006 ). Psilocybin was structurally described and synthesized in 1958 by Albert Hoffman (Hofmann et al. 1958 ), and a biosynthetic pathway was later proposed based on the transformation of labeled precursor molecules by Psilocybe cubensis (Agurell and Nilsson 1968 ). However, prohibition since the 1970s (21 U.S. Code § 812—schedules of controlled substances) has limited advances in psilocybin genetics, ecology, and evolution. There has been a recent resurgence of research on hallucinogens in the clinical setting; brain state imaging studies of psilocin exposure have identified changes in neural activity and interconnectivity that underlie subjective experiences, and therapeutic trials have investigated psilocybin's potential for treating major depression and addictive disorders (Griffiths et al. 2011 ; Carhart‐Harris et al. 2012 ; Petri et al. 2014 ; Carhart‐Harris et al. 2016 ; Johnson et al. 2017 ). Although the ecological roles of psilocybin, like most secondary metabolites, remain unknown, psilocin's mechanism of action suggests metazoans may be its principal targets.

The rate of horizontal gene transfer (HGT) between species of microorganisms is thought to be higher for genes located in gene clusters, which often encode all of the enzymatic, regulatory, and transport‐related steps required for a metabolic pathway to function in a single genomic locus. Such clusters may enhance the evolvability of fungi by facilitating the rapid loss or gain of multigene traits such as the production of bioactive molecules. Although developmentally complex mushroom‐forming fungi are thought to experience little HGT compared with morphologically simpler fungi, a scattered distribution of the hallucinogenic molecule psilocybin among diverse “magic” mushrooms led us to hypothesize that its biosynthetic pathway has been dispersed by HGT of a gene cluster. To test our hypothesis, we sequenced the genomes of three distantly related hallucinogenic mushroom species for comparison with closely related, nonhallucinogenic species. We identified a homologous multigene cluster in each hallucinogenic species by searching for clustering among all genes with a psilocybin‐like distribution among mushroom species. The enzymatic functions of genes within this cluster were confirmed here and in another concurrent study, and phylogenetic analyses support HGT of the cluster between divergent dung decomposers in the genera Psilocybe and Panaeolus , a first for mushroom‐forming fungi. Bioactive molecules like psilocybin are often presumed to have niche‐specific roles, but the ecological contexts in which they evolved are rarely known. We found that distantly related dung‐ and wood‐decay fungi have less variation in their genome content compared to close relatives in alternative niches, suggesting that this content is shaped in part by shared ecological pressures. Coupled with the inheritance patterns of the psilocybin cluster, these data support the hypothesis that psilocybin production is part of a larger adaptive strategy to dung and late wood‐decay niches, which harbor abundant invertebrates that eat or compete with fungi. We speculate that neuroactive compounds like psilocybin that target broadly conserved neurotransmitter receptors may have evolved as a strategy to influence arthropod activity in these niches, and that fungi within these niches could be further sources of neuroactive molecules.

Results and Discussion

IDENTIFICATION OF PSILOCYBIN GENE CLUSTERS IN THREE GENERA BY WHOLE GENOME SEQUENCING We identified candidate psilocybin genes by sequencing three diverse PS+ mushroom homokaryon genomes—Ps. cyanescens, Pa. ( = Copelandia) cyanescens, and Gy. dilepis (Table 1, GenBank MG548652‐MG548659), then comparing them to three related mushrooms not known to produce psilocybin (PS−): Galerina marginata, Gymnopilus chrysopellus, and Hypholoma sublateritium. Of the 37 PDHs that we identified to be consistent with a PS+ distribution among these taxa, only five genes were clustered, all in PS+ genomes. We retroactively designated Gy. chrysopellus, potentially PS+ because it possesses a cluster identical to Gy. dilepis, which is not surprising given inconsistent identifications, and geographical variation among Gymnopilus spp. phenotypes. Predicted functions of the five clustered genes were also consistent with psilocybin biosynthesis and metabolite transport, and were putatively designated tryptophan decarboxylase (PsiD), psilocybin‐related N‐methyltransferase (PsiM), psilocybin‐related hydroxylase (PsiH), psilocybin‐related phosphotransferase (PsiK), and psilocybin‐related transporter (PsiT). The orthologs of these enzymes shared 75–95% sequence similarity with those from a concurrently discovered psilocybin gene cluster in Ps. cubensis (Fricke et al. 2017), so we have adopted the same naming conventions here. Table 1. Genome assembly and annotation of psilocybin‐producing mushrooms Accession Length (nt) Scaffolds Contigs Average depth of coverage (x) N50 (nt) Complete BUSCOs (%) Total proteins Decarboxylases/PSD1 P450's1 Methyltransferases/DUF890 domain—proteins1 Kinases/phosphotransferases/PDH term 0PNAW1 MFS1 Gymnopilus dilepis2 SAMN07169108 47,177,497 8,423 10,681 16.5 33,540 73.43 16,257 28/9 151 89/4 275/29/1 37 Panaeolus cyanescens SAMN07166494 44,965,162 9,521 11,850 25.7 32,751 75.66 13,420 28/8 148 91/2 267/16/2 34 Psilocybe cyanescens SAMN07169033 53,483,841 18,721 38,006 44.7 46,250 72.18 15,973 38/17 178 102/2 298/23/1 44

CONFIRMATION OF PSILOCYBIN GENE CLUSTER ENZYME ACTIVITY To confirm gene cluster function, we profiled the enzymology of heterologously expressed full‐length cDNAs of Ps. cyanescens PsiD and PsiK in bacterial expression systems, and assayed their activities by LC‐MS/MS analyses. We determined that PsiD, the first committed step in the reaction and the only one not producing a drug‐scheduled compound, has specific decarboxylase activity on tryptophan. PsiD reactions produced tryptamine, identified at the characteristic m/z 144.1 [M + H]+, (Fig. S1; Supporting Information Data 1). PsiD did not decarboxylate phenylalanine, tyrosine, or 5‐hydroxy‐l‐tryptophan (5‐HTP) under the same conditions. We note that PsiD is similar to type II phosphatidylserine decarboxylases (PSDs), but has no significant sequence similarity with a pyridoxal‐5’‐phosphate‐dependent decarboxylase recently characterized in Ceriporiopsis subvermispora as specific for l‐tryptophan and 5‐HTP (Kalb et al. 2016). A unique GGSS sequence in a conserved C‐terminal motif (Fig. 2), suggests tryptophan decarboxylation is a previously unknown derived function among PSDs (Wriessnegger et al. 2009; Choi et al. 2015). We detected no activity of PsiK on 5‐HT or 4‐hydroxyindole (4‐HI) as alternatives to the psilocin substrate, possibly due to requirements for the 4‐hydroxyl and the methylated amine groups of psilocin. Our further characterization of PsiK and other enzymes was prevented by regulatory restrictions on the possession of substrates and products. However, a recent study has characterized complete biosynthesis of psilocybin by the homologous Ps. cubensis cluster (Fricke et al. 2017), serendipitously confirming the function of the gene clusters presented here. The pathway inferred by Fricke et al. (2017) proceeds as follows: tryptophan decarboxylation, tryptamine 4‐hydroxylation, 4‐hydroxytryptamine O‐phosphorylation, and sequential N‐methylations to produce psilocybin without a psilocin intermediate. This pathway was unexpected, because previous isotopic studies suggested tryptamine N‐methylation precedes hydroxylation of dimethyltryptamine and O‐phosphorylation of psilocin (Agurell and Nilsson 1968). Gene duplications among the clusters we have identified could suggest alternate or reticulated pathways also exist. Figure 2 Open in figure viewer PowerPoint Substrate specificity and C‐terminal sequence signatures of tryptophan decarboxylases (TDCs). (A) Table and sequence alignment of PsiD and other TDC homologs, previously characterized phosphatidylserine decarboxylase (PSD) type I and II enzymes, and closely related fungal TDC‐like proteins indicating conserved C‐terminal residues (386–409). The unique GGSS sequence in the conserved C‐terminal motif is outlined. (B) MEME‐derived sequence logo of conserved C‐terminal residues in fungal homologs of PsiD.

HORIZONTAL TRANSFER OF A PSILOCYBIN GENE CLUSTER Phylogenetic analyses of PS homologs from a local database of 618 fungal proteomes yielded congruent gene tree topologies with respect to PS+ taxa, and clades of clustered PS genes from all gene trees excluded the PS− taxa in the database, suggesting the clustered genes are coordinately inherited (Figs. 3 and S2A–E). The gene trees also suggest HGT of the cluster from Psilocybe to Panaeolus and HGT of most PS genes between Atheliaceae and Agaricaceae when compared to a phylogenomic tree of related Agaricales (Fig. 3). The direction of the latter HGT is ambiguous, and not strongly supported by all five genes. Analyses with additional PsiD and PsiK amplicon sequences retrieved by degenerate PCR of unsequenced Psilocybe and Conocybe genomes (Supporting Information Data 1, GenBank Accessions MG548652‐MG548659) suggest the dung fungus Ps. cubensis vertically inherited the cluster, and Pa. cyanescens acquired the cluster from Psilocybe sp., possibly from a dung‐associated lineage. Alternative hypotheses of vertical inheritance in these lineages were rejected; constrained topologies that exclude Pa. cyanescens and Conocybe smithii (AU test, P = 0.004) or Pa. cyanescens alone (P = 0.036) from putative donor clades were rejected (Supporting Information Data 1). Furthermore, a PsiD gene tree‐species tree reconciliation model allowing duplication, HGT, and loss (six events: D = 1, HT = 3, L = 2) is more parsimonious than a model that only allows duplication and loss (28 events: D = 3, L = 25) (Fig. S3). PS gene orthologs were not detected in Ps. fuscofulva; Ps. fuscofulva is sister to the rest of the genus, a pattern consistent with the ancestor of Psilocybe being PS− (Borovička et al. 2015). Conservation of synteny flanking the Ps. cyanescens PS cluster (Fig. 3B) suggests it may have been recently acquired in Psilocybe as well, or alternatively lost as a unit in close relatives. A genome wide scan did not identify any additional HGTs of genes or clusters between Psilocybe and Panaeolus (Fig. S4); however, additional supported HGTs from Agaricales to distant fungal lineages (Fig. 3C) suggest the constituent gene families may have experienced similar ecological distribution prior to the origin of the psilocybin gene cluster. Figure 3 Open in figure viewer PowerPoint Psilocybe cyanescens version 2 scaffold 5617 and Galerina marginata version 1 scaffold 9. PS clusters consist of a tryptophan decarboxylase (PsiD), one to two P450 monooxygenases (PsiH), a methyltransferase (PsiM), a phosphotransferase (PsiK), and one to two MFS transporters (PSiT). (C) RAxML phylogeny of TDC indicating putative HGT branches; Eutypa lata is in Xylariales (Ascomycota, violet shading), an order correlated with absence of termites in coarse woody debris (Kirker et al. 2012 Basidiobolus meristosporus (Zoopagomycota, pink shading) is commonly associated with amphibian dung and arthropods; members of the Phallomycetidae (green shading) commonly associated with dung, decayed wood, and/or insect spore dispersal. Gray taxon names = PCR sequences, black = whole genome. Support is percent of 100 ML bootstraps. **PsiH exists in high copy number in Fibulorhizoctonia sp. for which 54 similar PsiH gene homologs are not shown here. Phylogenies of all psilocybin genes are in Figure The evolution of psilocybin. (A) Phylogenomic tree of Agaricales (tan shading) with Atheliales (blue shading) outgroup. Support values = (internode certainty, tree certainty). Psilocybin‐producing taxa are indicated by a blue outline box. (B) PS locus synteny relative toversion 2 scaffold 5617 andversion 1 scaffold 9. PS clusters consist of a tryptophan decarboxylase (PsiD), one to two P450 monooxygenases (PsiH), a methyltransferase (PsiM), a phosphotransferase (PsiK), and one to two MFS transporters (PSiT). (C) RAxML phylogeny of TDC indicating putative HGT branches;is in Xylariales (Ascomycota, violet shading), an order correlated with absence of termites in coarse woody debris (Kirker et al.), with members that produce a white rot of wood;(Zoopagomycota, pink shading) is commonly associated with amphibian dung and arthropods; members of the Phallomycetidae (green shading) commonly associated with dung, decayed wood, and/or insect spore dispersal. Gray taxon names = PCR sequences, black = whole genome. Support is percent of 100 ML bootstraps. **PsiH exists in high copy number insp. for which 54 similar PsiH gene homologs are not shown here. Phylogenies of all psilocybin genes are in Figure S2 A–E. Recent studies have suggested HGT is pervasive in the fungi, especially among lifestyle‐associated genes (Wisecaver et al. 2014; Gluck‐Thaler and Slot 2015), and may occur along “highways” (frequent partners in gene exchange) that could correspond to shared environments (Szöllősi et al. 2015; Qiu et al. 2016). However, HGT has been found to be rare in Basidiomycota compared to Ascomycota (Wisecaver et al. 2014), suggesting the transfer of the PS cluster may have provided a significant fitness benefit to the recipient, and is to our knowledge, the first report of HGT of a secondary metabolite gene cluster between lineages of mushroom‐forming fungi (Agaricomycetes). A number of secondary metabolism gene cluster HGT events have been previously reported in Ascomycota, where gene clusters have been much more frequently identified than in Basidiomycota; however, any causal associations among rates of gene clustering, rates of HGT, and the strength of selection among fungal lineages remain largely uninvestigated (Slot 2017).