https://doi.org/10.5061/dryad.qz612jmng

Scripts and code to reproduce RNAseq analysis for looking at changes in expression due to captivity effects in shrews; specifically in the cortex, hippocampus, and olfactory bulb. These regions will help to understand region specific changes to the brain, as well as impacts on the whole brain.

Description of the data and file structure

First, we assess RNA-seq data quality, filtering low quality reads and trimming adapters in fastp. RNA-seq sample lists (samples_hippoc.txt, samples_olf_bulb.txt, samples_cortex.txt) are described below.

• Rows =Sample_ids. Note: Samples named HFS* correspond to wild shrews, and those named Eto* correspond to captive shrews as shown in the file.
• Columns = Organ, experimental group (wild or captive), sex, toothrow length (cm)

Then, we map reads pruned reads and mapped to the reference transcriptome, quantifying transcript abundance in kallisto. Quality of these processes can be found in Captivity_BioinforTable_Full.xlsx. NAs represent data that is not applicable to the row.

• Rows = Samples
• Columns = Region, Condition, Capture (Date), Euthanized (Date), Weight Caught (g), Weight Euthanized, Sex, Toothrow (cm), RIN, Reads pre filtered, Reads post filtered, Reads pre filtered (pairwise), Reads post filtered (pairwise), Reads removed, Reads removed (pairwise), K processed, K aligned, and mapping percent (%).

Gene count files for each tissue type (Olf_bulb_counts.txt, Hippocampus_counts.txt, Cortex_counts.txt), are described below.
Gene Counts (counts)

• Rows = Genes
• Columns = Samples

After this, we normalize data and analyze differential expression between conditions using DESeq2. Below are the results for the comparisons between wild and captive individuals (Olf_bulb_results_all.csv, Hippocampus_results_all.csv, Cortex_results_all.csv). NAs in p-values represent genes were differential expression was not tested, as they were filtered out if within a row, all samples have zero counts, if a row contains a sample with an extreme count outlier, or if a row is filtered by DESeq2 automatic independent filtering.
DESeq2 Results

• Rows = Genes
• Columns = means , log-fold changes, p-values

Additionally, here are the significant results (p<0.05) for each tissue (Olfactory_bulb__significant_gene_list.txt, Hippocampus_significant_gene_list.txt, Cortex_significant_gene_list.txt), as well as those genes sorted into unregulated and downregulated lists without loci (Downregulated_genes_cortex.txt, Downregulated_genes_hippoc.txt, Downregulated_genes_olf_bulb.txt, Upregulated_genes_cortex.txt, Upregulated_genes_hippoc.txt and Upregulated_genes_olf_bulb.txt) used for pathway analyses.

• Rows = List of significant genes

With the resultant genes we also perform a pathway enrichment analysis using DAVID Functional Enrichment Tools, to test whether they enrich KEGG pathways (Cortex_pathways_down.txt, Cortex_pathways_up.txt, Olf_bulb_pathways_down.txt, Olf_bulb_pathways_up.txt, Hippocampus_pathways_down.txt, Hippocampus_pathways_up.txt).

• Rows = Pathway
• Columns = genes tested, total hits, percent hits, enrichment value, p-values, adjusted p-values

Finally, an excel file where percentage of complex I NADH:ubiquinone oxidoreductase subunit genes in 12 of the enriched pathways was calculated.
NADH genes 2024.xlsx. NAs represent blank cells, as no more complex I NADH:ubiquinone oxidoreductase subunit genes were found in the pathway.

• Rows = Gene symbols in pathways, Number of NADH:ubiquinone oxidoreductase subunit genes (total), NADH:ubiquinone oxidoreductase subunit genes (%)
• Columns = Pathways

Code/Software

Required packages
Kallisto version 0.46.1
R version 4.2.0
DESeq2 version 1.36.0
Database for Annotation, Visualization and Integrated Discovery
(DAVID) (2022)

RNA-seq analyses first align to a reference and then quantify of reads. Genome and original unfiltered reads can be downloaded as described below. However, these steps could be skipped if reproducing, as count data has been provided above.
The reference (sorAra2; GCF_000181275.1) can be download from straight from NCBI, or using code provided.

mkdir ./data/ref/
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/181/275/GCF_000181275.2_SorAra2.0/GCF_000181275.2_SorAra2.0_genomic.gff.gz ./data/ref/
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/181/275/GCF_000181275.2_SorAra2.0/GCF_000181275.2_SorAra2.0_rna.fna.gz ./data/ref/
gunzip ../ref/GCF_000181275.2_SorAra2.0_rna.fna.gz

RNA-seq data from this project can be found on NCBI SRA. These can be downloaded manually, or using the getter.sh script with the help of sratoolkit (https://github.com/ncbi/sra-tools). Note, all scripts use indirect paths and may need to be updated depending on file structure.

bash get_rawseq.sh

Quality control, filtering, trimming

These scripts can be skipped if reproducing from counts. Trim adapters from downloaded reads and remove low quality reads using default settings and fastp. You will need to download fastp to your environment (https://github.com/OpenGene/fastp).

bash fastp.sh

Mapping and quantification

Reads that have went through quality control are then mapped to the reference transcriptome and quantified . This method does not directly map reads to the genome. (https://pachterlab.github.io/kallisto/about).

bash kallisto.sh

Note: This will create new transcript abundances separate ffrom the ones used in this analysis, but should be the same just stored in a different location. Change names in code below if required.

Analyses

Each analysis was conducted using the R code below for each tissue type. For best results, run in RStudio, as each matrix and figure is not set to print out in a best attempt to not overwrite results. Additionally, some intersection tables require data from the other codes, so 2024_2Cortex_captivity.R and 2024_2Olf_bulb_captivity.R are supposed to be run before 2024_2Hippocampus_captivity.R code, since hippocampus script includes codes for plots that need outputs from the two previous scripts.

2024_2Cortex_captivity.R
2024_2Olf_Bulb_captivity.R
2024_2Hippocampus_captivity.R

DAVID Geneset Enrichment

Both the above programs were done online at the below links. In a perfect world this should be scripted, however, due to conflicts in packages and Rversions they were not. https://david.ncifcrf.gov/summary.jsp