Ascomycetes

RESCIPt supports retrieval of SSU and LSU marker-gene data from SILVA via an automated method, “get-silva-data”, or manual import of the necessary sequence and taxonomy files (Fig 1). The “get-silva-data” pipeline allows selection of (a) which version of the database to download, (b) whether to download LSU, SSU sequences, or the SSU NR99 sequences, and (c) which taxonomic ranks to use and other options for taxonomy parsing (see software documentation for more details). These options are all stored in the data provenance of the output files, for later retrieval and reproducibility. RESCRIPt parses the SILVA taxonomy, using three files as input:
The “parse-silva-taxonomy” method utilizes the taxrank, taxmap, and taxtree files to generate a consistent user-defined rank-associated taxonomy. Although the set of ranks can be configured by the user, the following ranks are extracted by default: domain (d_), phylum (p_), class (c_), order (o_), family (f_), and genus (g_). Any ranks not associated with taxonomy have their upper-level taxonomic lineage propagated downward (i.e. the values are forward filled with the last observed taxonomic value) towards lower-level ranks. This ensures general compatibility with downstream taxonomy classification tools, many of which may require non-empty fields at each rank. Rank propagation can be optionally disabled.
Finally, the user can choose to append the organism name (from the taxmap file) for use as the species (s_) rank taxonomy. We generally warn against this due to the myriad of inconsistent information found within the organism name field (based on our benchmarking results described herein), but it can occasionally be useful. If the user does decide to leverage the organism name, we currently only return the first two words, to remove subspecies-level information that is often included in the given organism name and which can degrade classification accuracy (e.g., because the extra information causes that species to be interpreted as a unique label).
Rank propagation is provided to allow users to extract more taxonomic information, rather than explicitly pulling down only the ranks of interest. For example, if a user opted to download sequence data along with only the six standard taxonomic ranks (see above), they may obtain the following taxonomic output when rank propagation is not used:
Z27393.1.1722 d__Eukaryota; k__Fungi; p__Ascomycota; c__; o__; f__; g__
AB671439.1.2071 d__Eukaryota; k__Fungi; p__Ascomycota; c__; o__; f__; g__
The user might assume that query sequences that “hit” either of these reference sequences would be unable to classify beyond the phylum level. However, applying rank propagation will yield the following for these same accessions:
Z27393.1.1722 d__Eukaryota; k__Fungi; p__Ascomycota; c__Taphrinomycotina; o__Taphrinomycotina; f__Taphrinomycotina; g__Taphrinomycotina
AB671439.1.2071 d__Eukaryota; k__Fungi; p__Ascomycota; c__Pezizomycotina; o__Pezizomycotina; f__Pezizomycotina; g__Pezizomycotina
This is because intermediate ranks not selected by the user (e.g., sub-phyla Taphrinomycotina and Pezizomycotina) were propagated downward and used to fill in the unannotated ranks. Hence, forward filling allows users to disambiguate incompletely annotated reference sequences. The drawback is the conflation of taxonomy by mixing ranks from other levels.
The RESCRIPt project page (https://github.com/bokulich-lab/RESCRIPt) lists several tutorials describing how to use various RESCRIPt functions, including methods to import and parse SILVA data.

Optimal annealing temperatures for the selected fungal ITS primer pairs were explored in vitro. Three ascomycete and four basidiomycete fungal species collected on Yakushima Island, Kagoshima Prefecture, Japan, were subjected to the in vitro examination (A1–3 and B1–4: Table 2). DNA was extracted from fruiting body tissues of the fungal specimens using the CTAB method, as described elsewhere [32] . PCR was conducted using the buffer system of Ampdirect Plus (Shimadzu) with BIOTAQ HS DNA Polymerase (Bioline) under a temperature profile of 95°C for 10 min, followed by 35 cycles at 94°C for 20 s, 47°C, 50°C, 53°C or 56°C for 30 s, and 72°C for 20 s (40 s for the entire ITS region), followed by 72°C for 7 min. The concentration of MgCl₂, dNTPs, PCR primers and the template DNA in the reaction buffer were 1.5 mM, 200 µM, 0.5 µM and 1 ng/µl, respectively. Amplification of the DNA fragments was confirmed using the Flash Gel System for DNA (Lonza).
To evaluate the coverage of the selected fungal ITS primer pairs in vitro, an additional PCR assay was conducted under the optimal annealing temperature tested above using seven ascomycete and seven basidiomycete specimens (Table 2). All necessary permits of the sample collection were issued by Kyushu Regional Forest Office, Ministry of Agriculture, Forestry and Fisheries, Japan.

We tested Twisst on two published genomic datasets from Neurospora spp. (ascomycete fungi) and Heliconius spp. (butterflies), selected to represent different sampling strategies (four and five taxa, respectively), as well as different levels of evolutionary complexity. The Neurospora dataset (Corcoran et al. 2016 (link)) consisted of 22 aligned haploid genome sequences from Neurospora tetrasperma samples (10 of mating type A and 12 of mating type a), along with single genomes representing two related species: Neurospora crassa and Neurospora hispaniola. Whole genome alignments were obtained from http://datadryad.org/resource/doi:10.5061/dryad.162mh. We used Lineage-10 (UK) samples of N. tetrasperma, as these had been shown to carry a strong signal of introgression from N. hispaniola (Corcoran et al. 2016 (link)). Trees were constructed for sliding windows of 50 SNPs using BIONJ as described above, with the requirement that each sample had to be genotyped at ≥40 of the 50 SNPs per window. Topology weightings were computed using Twisst, with four defined taxa: N. tetrasperma mat a (12 sequences), N. tetrasperma mat A (10 sequences), N. crassa (one sequence), and N. hispaniola (one sequence).
The Heliconius dataset consisted of 18 resequenced genomes (or 36 haploid genomes) from Martin et al. (2013 (link)). These samples comprised five populations: two geographically isolated races of Heliconius melpomene, from Panama (H. m. rosina, n = 4) and Peru (H. m. amaryllis, n = 4), and their respective sympatric relatives Heliconius cydno chioneus from Panama (n = 4) and Heliconius timareta thelxinoe from Peru (n = 4), with which they are known to hybridize; along with two additional samples of the more distant silvanifrom clade to serve as outgroups. We limited our analysis to two chromosomes: 18, which carries the gene optix, known to be associated with red wing pattern variation; and 21, the Z sex chromosome, which has been shown to experience reduced gene flow between these species, probably due to genetic incompatibilities (Martin et al. 2013 (link)). Fastq reads were downloaded from the European Nucleotide Archive, study accession no. ERP002440. Reads were mapped to the H. melpomene reference genome version 2 (Davey et al. 2016 (link)) using BWA-mem (Li and Durbin 2009 (link); Li 2013 ), with default parameters. Genotyping was performed using the Genome Analysis Toolkit (DePristo et al. 2011 (link)) version 3 HaplotypeCaller and GenotypeGVCFs, with default parameters except that heterozygosity was set to 0.02. Phasing and imputation was performed using Beagle version 4 (Browning and Browning 2007 (link)). Trees were inferred as described above, and weightings were computed using Twisst, with the five taxa described above.

Flavonoid biosynthesis proteins from A. thaliana were used to search for homologs in S. cerevisiae using DELTA-BLAST with an e-value < 1e-15. To search for CHI homologs in yeasts, a HMMER profile was constructed using S. cerevisiae Aim18p and Aim46p and used to conduct a HMMER search 332 annotated budding yeast genomes (15 (link)) to identify hits with score >50 and e-value <0.001. To search for CHI homologs in more divergent fungi, tBLASTn was used to identify hits in over 700 unannotated Ascomycota genomes (16 (link)) with e-value <1e-6. The full plant and fungal CHI homolog phylogeny was built by adding the fungal sequences to a structural alignment of plant CHI homolog sequences (7 (link)) using the —add function of MAFFT (55 (link)) using default settings. The resulting alignment was truncated to the regions matching the initial structural alignment that included the CHI domain and filtered to sites with less than 50% gaps using trimAl (56 (link)). A phylogenetic tree was constructed using fasttree (57 ) with the settings (d -lg -gamma -spr 4 -slownni -mlacc 2), midpoint rooted, and edited to collapse clades using iTOL (58 (link)). The absence of any CHI homologs in metazoan lineages was determined by repeated search attempts using both blastp and DELTA-BLAST at permissive e-value thresholds restricted to relevant taxonomic groups. The determination that the single CHI homolog in the published genome of S. arboricola was due to a fusion of paralogs was made by aligning the coding sequence of all hits from Saccharomyces species via MAFFT using default settings (55 (link)) and looking for signatures of recombination between AIM18 and AIM46 homologs in the S. arboricola gene using RDP4 (59 (link)). A single recombination event was found in this sequence (all tests of recombination were significant at p<1e-13), which supports a fusion with AIM46 providing roughly the 5′ 60% and AIM18 the final 3′ 40% of the sequence.
A subset of protein sequences of plant and fungal CHIs were aligned via MAFFT using default settings (55 (link)). Protein sequence alignments were visualized using Jalview 2.11.2.5 and colored by sequence conservation (60 (link)). Mid-point–rooted phylogenetic trees of plant and fungal CHIs were generated through the Phylogeny.fr “one click” mode workflow (http://phylogeny.lirmm.fr/) using default settings (61 , 62 (link), 63 (link), 64 (link), 65 ). Alignments and phylogenetic trees were exported as svg files and annotated in Adobe Illustrator 26.5.

The first 25 amino acids of Mco10 or Atp19 of S. cerevisiae were used as templates to identify homologs of these proteins from the published fungal genomes using BlastP. The phylogenetic tree was constructed using homologous sequences from model Ascomycota and Basidiomycota. Maximum Likelihood phylogenetic inference was performed in MEGA X software^{52 (link)} with a Poisson correction model^{53 (link)}. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the Poisson model, and then selecting the topology with superior log likelihood value.