RNA Extraction followed a modified protocols from RNeasy plant mini kit (Qiagen), below we share only the modifications used to extract high qualilty RNA from difficult plant species.

http://www.qiagen.com/us/products/discovery-and-translational-research/dna-rna-purification/rna-purification/total-rna/rneasy-plant-mini-kit

Note: Amount of Starting Material: the amount of starting material to use can vary between species, we found that for more troublesome species containing a lot of secondary metabolites, such as Plantago , using less starting (~30-50mg) material resulted in higher quality RNA output.

Note: Grind with mortar and pestle: while instruments to disrupt tissue can work, we found that grinding for an extended period by hand in liquid nitrogen produced higher quality results for RNA extractions on non model plants

Note: When adding 450 ul of RLT Buffer perform a 3 min incubation before vortexing, proper mixing here is key to a successful extraction

Modification: The wash step with RW1 was performed twice on difficult samples, this helped remove excess metabolites that were still present in plant samples

Modification: The wash step with RPE was performed a minimum of three times but more washes were performed if color of the sample was still green or brown. The final sample should be a clear liquid.

After sequencing the first step is to set up a clean compute space to analyze the sequencing data to avoid issues with mixing versions of software.

There are a few different ways to manage a compute environment, here conda/mamba was used to create a clean environment as it works with all of the downstream programs and has an extensive manual. Mamba is a newer alternative to conda that can install programs with greater speed.

conda doc: http://conda.io/projects/conda/en/latest/user-guide/getting-started.html

mamba doc: http://mamba.readthedocs.io/en/latest/

mamba and conda can both be used in the same way, just replace mamba with conda in the codes below.

To set up an environment with conda or mamba you can use the following code to create an environment and install the required programs:

#First create environment with mamba
mamba create -n rnaseq_env
#Next activate our new environment 
mamba activate rnaseq_env
#Install tools needed in future analysis
mamba install fastqc
mamba install trinity

When obtaining data from a sequencing facility, it is essential to check the quality of data received for each sample and remove any unwanted (see adaptors) or low-quality reads. FastQC is the standard QC tool for sequencing experiments and other options are mainly command-line-only options that require

multiple steps or tools (Lo and Chain 2014).

Fast qc can be ran through command line or by downloading their GUI interface here:

http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

Read quality tends to decrease as it progresses

to tail ends of the reads, and FastQC should be used to check the quality at

the read tail. If the read quality is low or adapters are still present, then

the reads should be removed through one of multiple available trimming tools

Low quality reads should be removed by one of the many available read tools (Table 2), below is an example of how this is done with Trimmomatic (A. M. Bolger, Lohse, and Usadel 2014).

A	B
trimmomatic	http://www.usadellab.org/cms/?page=trimmomatic
cutadapt	http://cutadapt.readthedocs.io/en/stable/
fastp	http://github.com/OpenGene/fastp
AdapterRemoval	http://adapterremoval.readthedocs.io/en/stable/

Table 2: : List of possible tools for trimming data.

When no reference genome is available, or the

available genome is of low quality, RNA-Seq reads can be assembled into a de

novo transcriptome using a variety of tools. A few of the available tools are listed below.

A	B
Trinity	http://github.com/trinityrnaseq/trinityrnaseq/wiki
SOAP	http://github.com/aquaskyline/SOAPdenovo-Trans
Oases	http://github.com/dzerbino/oases
Trans-ABySS	http://github.com/bcgsc/transabyss
rnaSPAdes	http://github.com/ablab/spades

No tool is perfect for every organism, with specific tools performing better or worse depending on the ploidy level and repetitiveness of the reads and so a comprehensive comparison should be done before starting assembly (Madritsch, Burg, and Sehr 2021; Hölzer and Marz 2019).

For this protocol we will be using Trinity due to its high performance and its use in

previous publications for similar species to be used in this study (Ponnuchamy et al. 2020; Boehm et al., 2017).

Trinity \
--seqType fq \
--SS_lib_type RF  \
--left left.fq \
--right right.fq \
--CPU n \ #change n to the number of cpus in use
--max_memory 220G \
--trimmomatic #use this option to run trimmomatic with trinity, can be removed if trimming was completed already

Once the de novo assembly is completed, the next step is to check assembly quality. The quality of an assembly will affect multiple downstream analysis such as differential expression analysis, gene structure prediction, and isoform identification (Raghavan et al. 2022). Multiple tools and methods were researched to find the best practices for quality assessment, but no universal protocol is currently used to allow for comparison across studies (Haas et al. 2013; Xie et al. 2014).

A large percentage (70-90%) of a high-quality transcriptome will also map back to the reads, and this is an important step before proceeding with downstream analysis. For this protocol multiple ways to characterize the quality assembly were tested, starting with aligning reads to the assembly to count the number of proper and improper pairs.

A	B
STAR	http://github.com/alexdobin/STAR
Bowtie2	http://github.com/BenLangmead/bowtie2
BWA	http://github.com/lh3/bwa
HISAT2	http://daehwankimlab.github.io/hisat2/
TopHat2	http://github.com/infphilo/tophat

Table 3: List of tools that can be used to map transcripts back to sequencing reads and links to their repositories.

The code for QA can be found below

#Check BUSCO scores 
#install BUSCO 
mamba install -c conda-forge -c bioconda busco=5.7.1 #set lastest version here 
#run BUSCO ; set linage for your species, threads, and name of input and output files 
busco -l <viridiplantae_odb10> -m transcript -c <n> -i <your_genome.fasta> -o <your_output> 

#map genome back to input reads
#install bowtie2
mamba install bioconda::bowtie2
# first build index based on assembled transcriptome 
bowtie2 -build <input.fasta> <output.fasta>
#perform read alignment  
bowtie2 -p 10 -q --no-unal -k 20 -x Trinity.fasta -1 reads_1.fq -2 reads_2.fq  \
     2>align_stats.txt| samtools view -@10 -Sb -o bowtie2.bam
#combine stats 
cat 2>&1 align_stats.txt

#run rnaQUAST 
#install 
mamba install -c bioconda rnaquast
#run 
rnaQUAST \
--transcripts /PATH/TO/transcripts1.fasta /PATH/TO/ANOTHER/transcripts2.fasta /PATH/TO/MULTIPLE/*.fasta [...] \
--reference /PATH/TO/reference_genome.fasta --gtf /PATH/TO/gene_coordinates.gtf
#gtf file can be obtained from code that is available in downstream analysis. 

After transcriptome assembly and quality assessment, the putative functionality of the contained sequences can be elucidated. Most tools follow a similar procedure comprised of three steps: homology transfer and assignment via a search algorithm, sequence feature annotation, and gene ontology assignment.

A	B
Trinotate	http://github.com/Trinotate/Trinotate/releases
Annocript	http://github.com/frankMusacchia/Annocript
EnTAP	http://entap.readthedocs.io/en/v0.8.1-beta/index.html
OmixBox	http://www.biobam.com/omicsbox/

Table 4: List of tools that can be used for functional annotation and links to their repositories.

This protocol will use trinotate as it works directly from trinity, used above for de novo transcriptome assembly, the code to run trinotate is below.

#First we run transdecoder to determine likely coding regions 
#install 
mamba install transdecoder

#run transdecoder
/Pathtotransdecoder/TransDecoder.LongOrfs -t trinity.fasta
/Pathtotransdecoder/TransDecoder.Predict -t trinity.fasta 

# the output that we want from this is the transdecoder.pep file
# we can also create a gtf/gff3 file with transdecoder if needed with gff3_gene_to_gtf_format.pl 

# Next we will start the annotation process, this will be done in trinotate and can be run in parallel with the DE analysis

#make sure environement is activate and install 
mamba activate trinityenv
mamba install trinotate

#compile the sequence database, this only needs to be done once
/Path_to_trinotate/Build_Trinotate_Boilerplate_SQLite_db.pl Trinotate

#once completed you will have the following files
Trinotate.sqlite 
uniprot_sprot.pep
Pfam-A-hmm.gz

#prepare the protein database for blast searches, set working directory to file 
makeblastdb -in uniprot_sprot.pep -dbtype prot

#uncompress and prepare the pfam database
gunzip Pfam-A.hmm.gz
hmmpress Pfam-A.hmm

#next we will blast our fasta to the database to find annotations
#this is best to run on an HPC as it will take up to a few days for the blastx 
blastx -query your_transcriptome.fasta -db uniprot_sprot.pep -num_threads <n> -max_target_seqs 1 -outfmt 6 -evalue 1e-3 > blastx.outfmt6 	
blastp -query transdecoder.pep -db uniprot_sprot.pep -num_threads <n>-max_target_seqs 1 -outfmt 6 -evalue 1e-3 > blastp.outfmt6

#identify protein domains (optional)
#run this on HPC as well if possible 
hmmscan --cpu 12 --domtblout TrinotatePFAM.out Pfam-A.hmm trinity.fasta.transdecoder.pep > pfam.log

#load trinotate sqlite database 
Trinotate Trinotate.sqlite init --gene_trans_map Trinity.fasta.gene_trans_map --transcript_fasta Trinity.fasta --transdecoder_pep trinity.fasta.transdecoder.pep

# load blast homologies
Trinotate Trinotate.sqlite LOAD_swissprot_blastp blastp.outfmt6

Trinotate Trinotate.sqlite LOAD_swissprot_blastx blastx.outfmt6

# load pfam doamain entires
Trinotate Trinotate.sqlite LOAD_pfam TrinotatePFAM.out

#generate output report
Trinotate Trinotate.sqlite report [opts] > trinotate_annotation_report.xls

#extract go assignments per gene and make into annotations file 
/path_to_trinity_install/bin/extract_GO_assignments_from_Trinotate_xls.pl \
                         --Trinotate_xls trinotate.xls \
                         -G --include_ancestral_terms \
                         > go_annotations.txt

#get a gene lengths file, this can be used for downstream analysis 
/path_to_trinity_install/opt/trinity-2.13.2/util/misc/fasta_seq_length.pl  Trinity.fasta > Trinity.fasta.seq_lens

/path_to_trinity_install/opt/trinity-2.13.2/util/misc/TPM_weighted_gene_length.py  \
         --gene_trans_map /media/chase/Samsung_T5/Trinotate_Lotus/Trinity.fasta.gene_trans_map \
         --trans_lengths /media/chase/Samsung_T5/Trinotate_Lotus/Trinity.fasta.seq_lens \
         --TPM_matrix kallisto.isoform.TMM.EXPR.matrix > Trinity.gene_lengths.txt

EnTAP was also performed to compare annotation results across platforms, the code for this is below.

#install EnTAP based on the site instructions (http://entap.readthedocs.io/en/latest/)

#Run database configuration, this can be run for all databases that your will search against eg: nr, swissprot, etc.. 
EnTAP --config -d  /Path_to_database_fasta/nr_database.fasta --out-dir output_folder -t <n>--ini entap_config.ini
#Run EnTAP 
EnTAP --runP -i protiens.pep -d /Path_to_database/nr.dmnd -t <n>--ini Path_to_Entap/EnTAP-v0.10.8-beta/entap_config.ini

As with the protocols described above, multiple tools are currently available to estimate transcript abundance in the absence of a genome. Recent comparisons of these tools can be found by (Wu et al. 2018; Corchete et al. 2020). These studies found multiple tools perform well and the choice of tool comes down to multiple factors that will differ between experiments.

A	B
STAR	http://github.com/alexdobin/STAR
TopHat	http://github.com/infphilo/tophat
Kallisto	http://github.com/pachterlab/kallisto
Salmon	http://github.com/COMBINE-lab/salmon

Table 5: List of tools that can be used for transcript quantification and links to their repositories.

Kallisto was used for this protocol due to speed and low computational requirments with high accuracy.

#install 
mamba install kallisto
#create index file 
kallisto index -i <index.idx> <transcriptome_file.fasta>

#Can be run with trinity
path_to_trinity/util/align_and_estimate_abundance.pl \
--transcripts transcriptome_file.fasta
--seqType fq \
--left left.fq.gz \
--right right.fq.gz \
--est_method kallisto \
--thread_count <n>\
--gene_trans_map Trinity.fasta.gene_trans_map \
--SS_lib_type RF \
--output_dir kal_out

#create count file for multiple samples with kallisto 
path_to_trinity/utils/abundance_estimates_to_matrix.pl \
--est_method kallisto \
--gene_trans_map Trinity.fasta.gene_trans_map \
--out_prefix kallisto \
--name_sample_by_basedir \
--basedir_index -3 \
Sample1/kal_out/abundance.tsv \
Sample2/kal_out/abundance.tsv 

#can add all samples, two are given as example 

Some plant species already have a reference genome available, if this is the case, mapping can be done directly to the reference. This can be done with multiple mapping tools (Table 3). A quick example of mapping with STAR is given below.

#First create index file 
STAR \
--runThreadN <n> \
--runMode genomeGenerate \
--genomeDir Path_to_genomefile/Genome \
--genomeFastaFiles Genome/genome_fasta \
--sjdbGTFfile Genome/genome_gtf \
--sjdbOverhang 99 \
--genomeChrBinNbits 7

#next run mapping 
STAR \
--readFilesIn Path_to_forwardread.fastq Path_to_reverseread.fastq 
--runThreadN <n> \
--genomeDir Path_to_genomefile/Genome \
--readFilesCommand zcat \
--quantMode GeneCounts \
--outFileNamePrefix ${i/"_Clipped_R1.fastq.gz"/""} \
--sjdbGTFfile $genomedir$genome_gtf  \
--sjdbOverhang 99 \
--outSAMtype BAM SortedByCoordinate
done

Once the transcripts have been quantified, the next step is to investigate the changes in gene expression caused by the experimental stimuli. There are a variety of tools that can be used for differential gene expression analysis with DeSeq2, EdgeR, and Limma currently being the most used (Costa-Silva, Domingues, and Lopes 2017). This protocol uses DeSeq2, which can be run seperately in R or following the trinity analysis protocol.

A	B
DeSeq2	http://bioconductor.org/packages/release/bioc/html/DESeq2.html
edgeR	http://bioconductor.org/packages/release/bioc/html/edgeR.html
Limma	http://bioconductor.org/packages/release/bioc/html/limma.html
ROTS	http://www.btk.fi/research/research-groups/elo/software/rots/

Table 6: List of tools that can be used for differential expression analysis and links to their repositories.

After differential gene analysis, the next step is to place the change in genes to a biological context. GO enrichment analysis is a technique that is widely used to highlight important biological processes in complex outputs of RNA-Seq experiments (Ashburner et al. 2000).

A	B
GSEA	http://www.gsea-msigdb.org/gsea/index.jsp
Panther	http://geneontology.org/
G: Profiler	http://biit.cs.ut.ee/gprofiler/gost
GOseq	http://bioconductor.org/packages/release/bioc/html/goseq.html

Table 7: List of tools that can be used for GO enrichment analysis and links to their repositories.

Example code for DEG and GO enrichment analyisis is included below:

# This code is ran following the trinity analysis pipeline, however all analysis can also be run in R/Rstudio seperately
# setup is based on the sample file and contrasts file, see trinity wiki for setup 
/Path_to_trinityfile/Analysis/DifferentialExpression/analyze_diff_expr.pl \
--matrix /media/chase/Samsung_T5/GOA_RNA/Holcus/kallisto.gene.TMM.EXPR.matrix \
-P 0.05 \
-C 1 \
--max_genes_clust 20000 \
--samples /media/chase/Samsung_T5/GOA_RNA/Holcus/sample_file.txt \
--examine_GO_enrichment --GO_annots /Path_to_annotationfile/go_annotations.txt --gene_lengths /Path_to_genelengthfile/Trinity.gene_lengths.txt

#In R
library("DESeq2")
#Creating a DESeq2 dataset, here only considering the treatment (control vs. treated).
dds_treatment <- DESeqDataSetFromMatrix(countData = cts,
                                        colData = metadata_file,
                                        design = ~ Treatment(variable))
#run DeSeq2
dds_treatment <- DESeq(dds_treatment)
# get results 
res_trt <- results(dds_treatment)
# Export Results 
write.csv(x = as.data.frame(res_trt[order(res_trt$pvalue),]), 
          file = "treatment_results.csv")