A draft reference assembly of the Psilocybe cubensis genome

DNA isolation

Dried stems from Psilocybe cubensis strain P.envy. The strain name is anecdotal reported to have been grown axenically (unknown conditions) and obtained in Somerville, MA, US. These samples were used for isolation of high molecular weight (HMW) DNA using a modified CTAB/Chloroform and SPRI protocol. Briefly, 300mg of stem sample were ground to a fine powder using a -80C frozen mortar and pestle. 150 mg of powder was then aliquoted into 2 mL conical tubes (USA Scientific) with 1.5 mL cetrimonium bromide. These tubes were then incubated at room temperature on a tube rotator for 10 minutes. 6 uL of RNase A (Promega 4 mg/mL) was then added and both tubes were incubated at 37°C for one hour, vortexing every 15 minutes. Following this incubation, 7.5 uL Proteinase K (New England Biolabs 20 mg/mL) was added and the tubes were incubated at 60°C for 30 minutes, vortexing every 10 minutes. At the conclusion of the Proteinase K incubation, both tubes were incubated on ice for 10 minutes. The samples were then centrifuged for 5 minutes at 14000 rpm. 600 uL of supernatant was removed from each tube and added to 600 uL chloroform. The tubes were then vortexed until opaque and spun for 5 minutes at 14000 rpm. 400 uL of the aqueous layer was removed using a wide bore tip and added to a 1.5 mL Eppendorf tube. 400 uL MIP (marijuana infused products) Solution B and 400 uL DNA Binding Beads (Medicinal Genomics PN 420004) were added to the Eppendorf tube and inverted to homogenize. The tubes were then incubated at room temperature on the tube rotator for 15 minutes. The tubes were then removed from the rotator and placed on a magnetic tube rack for 3 minutes. The supernatant was removed, the beads were washed twice with 1 mL of 70% ethanol and allowed to dry for 5 minutes. The beads were then eluted in 100 uL of 56°C Elution Buffer (Medicinal Genomics PN 420004) using a wide bore tip and incubated at 56°C for 5 minutes. Following this incubation, the tubes were returned to the magnetic rack, the supernatant of both tubes were removed using a wide bore tip and pooled in a fresh Eppendorf tube. HMW DNA length was quantified on an Agilent TapeStation and produced a DIN of 8.1. Qubit Fluorometer (Thermo Fisher Scientific) quantified 55ng/ul. Nanodrop Spectrophotometer (Thermo Fisher Scientific) quantified 104ng/ul with 260/280nm ratio of 1.85 and 260/230nm of 0.95.

Sequencing

Sequencing libraries were constructed according to the manufacturer’s instructions for Pacific Biosciences Sequel II HiFi sequencing. 773,735 CCS reads were generated. Quast (Gurevich et al. 2013) was used to assess the quality of the input fasta sequence file (N50 = 13.9Kb) and the output assembly fasta file (3.33Mb N50).

Assembly and annotation

The unfiltered CCS data was assembled using the Peregrine assembler (pg_asm 0.3.5,arm_config5e69f3d+) (Chin 2019). Reads were assembled into 32 contigs with lengths ranging from 32 kilobases to 4.6 megabases (Figure 1). The Peregrine assembler requires at least 2 HiFi reads to substantially overlap to contribute to a contig and as a result we did not observe any bacterial contamination in the assembly BUSCO v3.0.2 completeness scores (97.6%) were measured using agaricales_odb10.2020-08-05 BUSCO lineage database (Table 1) (Simao et al. 2015; Waterhouse et al. 2018). FunAnnotate v1.8.4 was used to annotate the genome (Li and Wang 2021) resulting in 13,478 genes.

Figure 1. Psilocybe cubensis P.envy contig length distribution (n = 32).

The final genome assembly was aligned to three other public Psilocybe cubensis datasets (Fricke et al. 2017; Torrens-Spence et al. 2018; Reynolds et al. 2018) and one different Psilocybe species (Psilocybe cyanescens) to verify taxonomic identification (Table 2). In total, 96-98.75% of these Psilocybe cubensis sequences align to the new HiFi generated Psilocybe cubensis P.envy reference using minimap2 and bwa-mem (Li and Durbin 2010; Li 2018). Mapping rates were determined using samtools flagstat on bam files (Li et al. 2009). Alignments were visualized with MUMmer V4.0.0beta2 and Integrative Genomics Viewer v2.4.16 (Delcher et al. 2003; Robinson et al. 2011; Thorvaldsdottir et al. 2013).

Three Illumina genome assemblies (Reynolds et al., McKernan et al., Fricke et al.) were additionally aligned using MUMmer for whole genome alignment plots (Figure 2).

Figure 2. Whole genome alignments of short read Illumina assemblies to Psilocybe cubensis strain P. envy.

Left is Psilocybe cyanescens from Reynolds et al. Middle is McKernan et al. (MGC) Illumina assembly. Right is Fricke et al. or JGI assembly.

Polymorphisms

Illumina whole-genome shotgun data (McKernan et al. NCBI Project: PRJNA687911) was mapped to the P. envy HiFi reference assembly using bwa-mem (version0.7.17-r1188), samtools (version 1.8), sorted with sambamba (version 0.6.7) and variants were identified using GATK HaplotypeCaller (version 4.1.6.0) with default arguments.  The annotation from the funannotate pipeline was converted from gff3 format into SnpEff (version 4.3t 2017-11-24) database as described here (https://pcingola.github.io/SnpEff/se_buildingdb/) and the variants that came out of HaplotypeCaller were annotated. 553,716 variants (471,443 SNPs and 82,273 small insertions and deletions) were called and annotated equating to aSNP every 99bp and a variant every 83bp including indels. Of these, 375,896 (67.9%) are heterozygous and 177,820 (32.1%) are homozygous with a ratio of just over 2 to 1 heterozygous:homozygous variants. Lastly, as a quality check, the original Pacific Biosciences CCS corrected shotgun reads were mapped back to the reference with minimap2 (version 2.17-r941) and variants were called again using GATK HaplotypeCaller. A total of 15,963 variants are identified and 15,674 (98.2%) are heterozygous with only 289 homozygous variants called. Whole genome shotgun reads mapped back to their consensus reference should produce predominantly heterozygous calls in a diploid organism. Scripts utilized to for variant calling are in github and described in the Data availability section.

Structural variation

The N-methyltransferase gene responsible for Psilocybin production in P.envy contains a structural variation not seen in previous P. cubensis surveys (Figure 3). Illumina read mapping of the McKernan et al. P. cubensis assembly in NCBI (NCBI Project: PRJNA687911) demonstrates multiple read pairs spanning a 4.6kb insertion in the HiFi P. cubensis strain P.envy (SRA submission SRP299420). This insertion extends the 3’ end of the P.envy N-methyltransferase gene. The 4.6kb insertion is also observed as a deletion in Psilocybe cyanescens and as a deletion in RNA-Seq data from Torrens-Spence et al. (NCBI Project: PRJNA450675) (Reynolds et al. 2018). Other SNPs also exist in these genes and should be considered in context of this deletion. Further work is required to understand the biological significance of this variation.

Figure 3. IGV display of Illumina reads mapped to HiFi Psilocybe cubensis P.envy assembly.

Top track is Medicinal Genomics Illumina whole genome shotgun data of a different P. cubensis (strain name unknown: NCBI Project: PRJNA687911) mapped to the HiFi P. cubensis strain P.envy. Second track contains RNA-Seq data from a third P. cubensis genome (strain name also unknown: NCBI Project: PRJNA450675) hosted at JGI. Third track is Psilocybe cyanescens genome mapped to HiFi P. cubensis P.envy reference genome. Fourth track is FunAnnotate GFF3 annotation of the HiFi P. cubensis P.envy genome.