DOI: 10.5061/dryad.sj3tx968b
title: Comparative transcriptomics reveals divergence in pathogen response gene families amongst 20 forest tree species

description: Forest trees provide critical ecosystem services for humanity that are
 under threat due to ongoing global change. Measuring and characterizing
 genetic diversity is key to understanding adaptive potential and
 developing strategies to mitigate negative consequences arising from
 climate change. In the area of forest genetic diversity, genetic
 divergence caused by large-scale changes at the chromosomal level has been
 largely understudied. In this study, we used the RNA-seq data of twenty
 co-occurring forest tree species from genera including Acer, Alnus,
 Amelanchier, Betula, Cornus, Corylus, Dirca, Fraxinus, Ostrya, Populus,
 Prunus, Quercus, Ribes, Tilia, and Ulmus sampled from Upper
 Peninsula of Michigan. These data were used to infer the origin
 and maintenance of gene family variation, species divergence time, as well
 as gene family expansion and contraction. We identified a signal of common
 whole genome duplication events shared by core eudicots. We also found
 rapid evolution, namely fast expansion or fast contraction of gene
 families, in plant-pathogen interaction genes amongst the studied diploid
 species. Finally, the results lay the foundation for further research on
 the genetic diversity and adaptive capacity of forest trees, which will
 inform forest management and conservation policies.
description: Leaves of each of twenty forest
 species were sampled in August 2021, at the University of Notre Dame
 Environmental Research Center (UNDERC) near Land O' Lakes, Wisconsin,
 U.S.A. (Fig. 1). The UNDERC site is located at latitude: 46.23391,
 longitude: -89.537254, and encompasses approximately ~8,000 acres of mixed
 deciduous forest on both sides of the state line between Wisconsin and
 Michigan’s upper peninsula. The twenty species used in this study are
 dominant components of the woody vegetation at the site and are:
 <em>Acer rubrum</em> (red maple), <em>Acer
 saccharum</em> (sugar maple), <em>Acer spicatum</em>
 (mountain maple), <em>Alnus incana</em> (Gary alder),
 <em>Amelanchier laevis</em> (smooth serviceberry),
 <em>Betula alleghaniensis</em> (yellow birch),
 <em>Betula papyrifera</em> (paper birch), <em>Cornus
 alternifolia</em> (pagoda dogwood), <em>Corylus
 cornuta</em> (beaked hazelnut), <em>Dirca palustris</em>
 (Eastern leatherwood), <em>Fraxinus nigra</em> (black ash),
 <em>Ostrya virginiana</em> (American hophornbeam),
 <em>Populus grandidentata</em> (bigtooth aspen),
 <em>Populus tremuloides</em> (quaking aspen), <em>Prunus
 serotina</em> (black cherry), <em>Prunus virginiana</em>
 (chokecherry), <em>Quercus rubra</em> (northern red oak),
 <em>Ribes cynosbati</em> (prickly gooseberry), <em>Tilia
 americana</em> (American basswood), <em>Ulmus
 americana</em> (American elm). They are all angiosperms (flowering
 plants) and eudicots, belonging to the clades of asterids, saxifragales,
 and rosids.  These species were chosen as they are the most dominant at
 the study site and in large national forests in the region (i.e. Ottawa,
 Nicolet, and Hiawatha National Forests). Specifically, these species
 compose &gt;95% of the tree individuals and above-ground biomass in
 these forests (Swenson et al. 2017). We focus on the most dominant species
 in order to provide insights into the ecological and evolutionary
 processes driving transcriptomic variation within and across populations
 and communities (Swenson and Jones 2017). Upon collection, the samples
 were flash-frozen in liquid nitrogen, then transported and stored in an
 ultracold freezer at the field station. The samples were sent to Novogene
 Corporation Inc. (Chula Vista, CA, USA) for mRNA extraction, library
 preparation, mRNA sequencing, and quality control. The samples were
 sequenced on an Illumina NovaSeq 6000 platform to produce 150 bp
 paired-end reads. Raw reads filtering was performed by Novogene to get the
 clean reads including removing reads containing adapters, removing reads
 containing N &gt;10%, and removing reads containing low quality base
 (Phred quality score ≤5) which is over 50% of the total base. An average
 of 6 Gb paired-end clean reads were acquired for each sample and deposited
 in NCBI Sequence Read Archive (SRA accession number: PRJNA896247).
 <em>De novo</em> transcriptome assembly was performed for each
 sample using the software TransLiG v1.3 with default parameters (Liu et
 al. 2019). The clean reads were then mapped to the assembled transcriptome
 using the software RSEM v1.2.25 with default parameters (Li and Dewey
 2011). The longest transcript (isoform) with a transcript per million
 (TPM) value ≥1.5 for each gene was retained. Redundancy was removed from
 the retained transcripts using the software CD-HIT v4.8.1 with parameters
 of -c 0.95 -n 8 -T 32 -M 1540 (Fu et al. 2012). To remove potential
 contamination, the retained transcripts were aligned against the
 non-redundant (nr) protein database (downloaded on January 22, 2022) using
 the blastx function implemented by the software DIAMOND v2.0.13.151
 (Buchfink et al. 2021) with parameter of --outfmt 102 to assign the
 best-known taxon ID to each query transcript. The file that was used to
 convert accession numbers to taxonomy was downloaded from https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/, and the taxon ID file was downloaded from https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/new_taxdump/. The query transcripts that could be assigned to green plant (Viridiplantae) taxon IDs were retained. The filtered transcripts were aligned to the nr database that only includes green plant subjects using the blastx function implemented by the software DIAMOND v2.0.13.151 with default parameters. The aligned transcripts were combined with those transcripts assigned with green plant taxon IDs to form the final transcriptomes. The completeness of the final transcriptomes was examined using the 1,614 orthologues in the BUSCO (Benchmarking Universal Single Copy Orthologs) v5.2.2 set of embryophyta_odb10 (Manni et al. 2021). The transcriptomes were annotated with Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms. The transcripts were aligned to the nr database that only includes green plant subjects using the blastx function implemented by the software DIAMOND v2.0.13.151 with parameters of -taxonlist 33090 --very-sensitive -f5. The .xml outputs were used as input for mapping and annotation using the software Blast2GO (Götz et al. 2008). The transcripts were translated to protein sequences using the software TransDecoder v5.5.0 (https://github.com/TransDecoder/TransDecoder/wiki) with default parameters to get the longest open reading frame (ORF) per transcript and then were annotated by the software InterProScan v5.54-87.0 (Jones et al. 2014) with the parameters of -dra -f XML -goterms -pa. The final GO terms were acquired by merging Blast2GO and InterProScan annotation outputs using the software Blast2GO. To get the KEGG terms, the longest ORF per transcript were also annotated using EggNOG-mapper v2.1.7 (Cantalapiedra et al. 2021) with parameters of -m diamond -sensmode more-sensitive -tax_scope 33090 -tax_scope_mode 33090. The assembled transcriptome and annotation were used in the downstream analyses as shown in the flowchart (Supplementary Fig. S1).