Herbaria provide a valuable resource for obtaining informative mRNA

RNA recovery, quality, and sequencing results

We extracted RNA from three herbarium specimens (“-Herbarium” samples) whose vouchers are stored at MICH herbarium, with collection dates ranging from 1950 to 1995. We also investigated the quality differences between samples sequenced during the traditionally recommended time frame (“-Fresh” samples) and their respective 4-year-old silica-dried vouchers (“-Silica” samples), which were newly sequenced for this study (Supplemental Table 1). The RNA integrity (RIN) scores of the silica-dried and herbarium samples were consistently below the recommended level for sequencing (Fig. 1A; Supplemental Table 1), unlike those of the fresh samples (Supplemental Table 2). Despite low RIN values and in two cases (Antistrophe-Herbarium and Clavija-Herbarium) RNA concentrations that were undetectable using a Qubit (see Methods), RNA-seq libraries were successfully constructed, and the expected number of raw sequence reads was recovered.

Figure 1.

Extracted RNA from fresh, silica-dried, and herbarium specimens provide usable transcriptomes based on bioinformatic metrics. (A) The RNA integrity (RIN) scores for samples used in this study. (B) The proportion of contigs inferred to be “good contigs” compared with the proportion inferred as “bad contigs” based upon the TransRate metrics. (C) The inferred completeness of the transcriptomes based upon sequenced single-copy orthologs. (D) The total number of inferred orthologs from each sample used to infer the phylogeny. (E) The coalescent-based inferred phylogenetic relationships and sample source for transcriptomes in this study. Branch lengths converted to molecular rates (subs/bp).

Open in new tabDownload PowerPointLink to figure

To ensure uniformity in the transcriptome assembly procedure, the raw reads obtained from Carruthers et al. (2024) and the newly sequenced raw reads were all assembled using the same read cleaning, assembly, and postprocessing procedure outlined in the Methods. One assessment of RNA quality we used was the ratio of paired to unpaired reads in these assemblies. Theoretically, fewer paired reads should be present as sequences degrade. When comparing results from samples collected from the same plant, the silica-dried samples had a higher proportion of unpaired sequence reads after trimming, ranging from 26.49% (Navarretia-Silica) to 43.69% (Collomia-Silica) compared with their respective fresh sample: 15.52% (Navarretia-Fresh) to 15.79% (Collomia-Fresh) (Supplemental Table 3). Herbarium samples had an even greater proportion of unpaired reads after trimming, with values ranging from 70.08% (Clavija-Herbarium) to 79.95% (Antistrophe-Herbarium). The loss of read pairs likely reflects a decrease in RNA quality over time, as evidenced by the shifted read-length distribution of the remaining reads that lost their pairs (Supplemental Table 4).

To further investigate signatures of RNA degradation, we examined the length of the insert size for the fragments before sequencing. We found that the insert length for the herbarium samples was shorter than that of the fresh and silica-dried samples (Fig. 2A; Supplemental Fig. 1). Finally, to test for age-related damage, we analyzed the raw sequencing data for signs of deamination. None of the samples, regardless of preservation approach, exhibited typical damage patterns (Fig. 2B; Supplemental Fig. 2). The maximum level of adenosine to inosine (A > I) degradation was 0.12% in Primula-Silica, and the maximum for cytidine to uridine (C > U) was 0.13% in Ipomopsis-Silica.

Figure 2.

Quality assessment of sequenced and assembled RNA from representative samples for each preservation method. The quality analysis results for one representative from each preservation treatment are shown here; results for all samples are available in Supplemental Figures 1 through 3. (A) Deamination levels for adenosine to inosine (A > I) and cytidine to uridine (C > U). The plots show position-specific nucleotide misincorporations from the 5′ (left) and 3′ (right). (B) The length of the insert for the sequenced fragments. (C) The scores for TransRate2 broken into different predicted error types: reliability of the base pair calls (sCnuc), the coverage level of the transcript (sCov), the pairing of the mapped reads (sCord), and evidence of chimerism within the transcripts (sCseg). (D) The BUSCO scores for each sample.

Open in new tabDownload PowerPointLink to figure

Transcripts assembled from sequenced archival RNA are inferred to be high quality

The quality of the assembled transcripts was assessed using the TransRate algorithm (Smith-Unna et al. 2016) as implemented in TransRate2 (see Methods). TransRate2 is a tool designed to evaluate the quality of de novo transcriptome assemblies, expanding on TransRate by implementing the splice-aware HISAT2 mapping approach (Kim et al. 2019). In this approach, contig quality is inferred by analyzing common assembly errors (Smith-Unna et al. 2016). All transcripts that did not exhibit errors in assembly were considered “good quality,” whereas those with issues that exceed the allowable threshold (see methods) were considered “bad quality” (Fig. 1B). Between 65.8% (Leucothoe-Fresh) and 79.9% (Pyrola-Fresh) of the transcripts from the fresh tissue were regarded as good quality by this metric. For the silica-preserved samples, these numbers varied from 46.4% (Ipomopsis-Silica) to 93.7% (Collomia-Silica), and for the herbarium samples, they ranged from 66.3% (Antistrophe-Herbarium) to 94.1% (Clavija-Herbarium) (Supplemental Table 5). Five of the seven fresh samples produced a higher proportion of good-quality transcripts than their silica-dried counterpart. However, overall quality determinations and individual quality metrics were similar across preservation methods (Fig. 2C; Supplemental Fig. 3). TransRate2 was also used to calculate the proportion of bases covered across each transcriptome (Supplemental Table 5). This value for fresh tissue ranged from 97.70% to 98.27%; for silica-dried tissue, from 93.85%–96.58%; and for herbaria tissue, from 95.58%–97.00%. We investigated the proportion of bases covered on a contig-by-contig basis and found 100% of the bases covered ranged from 69.99%–76.23% for fresh, 48.27%–66.92% for silica-dried, and 59.05%–64.59% for herbarium samples. Contigs that did not have 100% coverage in general had >95% coverage (Supplemental Fig. 4).

We examined transcriptome completeness using Benchmarking Universal Single-Copy Orthologs (BUSCOs) (Simão et al. 2015). In all comparisons between silica-dried and fresh samples, the silica-dried samples had fewer missing BUSCOs (Figs. 1C, 2D; Supplemental Table 6). The sequencing depth for the silica-dried samples was, on average, much greater than that for fresh-collected samples (Supplemental Table 3). Thus, sequencing depth likely explains the greater number of complete BUSCOs recovered from silica-dried samples. These results indicate that single-copy orthologs were still identifiable in silica-dried samples, despite being stored at room temperature for 5 years. For transcriptomes derived from herbarium specimens, we found that the 1950s samples had fewer complete BUSCOs despite having similar sequencing depth to the silica-dried samples and higher sequencing depth than the fresh samples. Missing BUSCOs ranged from 198 to 227 for fresh-derived, 26–102 for silica-derived, and 131–340 for herbarium-derived transcriptomes. However, two of the three counts of missing BUSCOs for herbarium samples were still within the range of those for fresh-collected samples (counts of 131 and 211 missing BUSCOs), highlighting the importance of sequencing depth when assembling transcriptomes from herbarium specimens. We fitted a multinomial logistic regression of BUSCO recovery as a function of preservation method and sequencing depth (Supplemental Fig. 5). As expected, the proportion of complete loci increased with sequencing depth. In contrast, proportions of fragmented and missing loci decreased. Interestingly, the fresh samples in this analysis showed a marginally lower proportion of complete BUSCOs compared with the silica-dried samples and showed a lower proportion of fragmented loci.

Archival RNA may be used for phylotranscriptomic analyses

Using all samples in the study, we inferred gene trees and a coalescent-based species tree using a phylotranscriptomic approach with the branch lengths converted to molecular rates (see Methods). In total, 3341 orthologs were inferred, which contained at least 10 taxa (Supplemental Table 7). The lowest matrix occupancy was for Antistrophe-Herbarium with 842 orthologs recovered. The highest occupancy was for Darlingtonia1-Silica with 2992 orthologs recovered (Fig. 1D). The overall topology of the inferred species tree matched the current consensus relationships among species (Larson et al. 2020; Zuntini et al. 2024), and samples collected from the same plant were always sisters in the tree (Fig. 1E).

We found no evidence of longer branches subtending silica-dried or herbarium specimens, as is sometimes observed for degraded museum samples owing to damage to nucleic acids not being completely accounted for during analysis (Derkarabetian et al. 2019). This suggests that damaged nucleic acids did not erroneously increase branch lengths in our study. Although assemblies from the same individual were not identical, this is not unexpected as stochastic differences introduced during the library prep, sequencing, and assembly processes can result in differences in the final assembled transcript, especially for transcripts containing heterozygous sites. Importantly, this demonstrates that mRNA sequencing of herbarium samples can be used to infer evolutionary relationships among species.

Archival RNA uncovers genes associated with plant immune systems

After we confirmed the quality of the assembled transcriptomes, we examined the nucleotide-binding and leucine-rich repeat (NLR) immune receptor repertoires in transcriptomes from each of the three collection types (i.e., fresh, silica-dried, and herbarium samples). We selected this gene family because (1) NLRs are fast evolving and complex in architecture; (2) their expression often correlates with function; and (3) herbarium samples represent an untapped resource to interrogate NLR diversity. Using the NLRtracker algorithm (Kourelis et al. 2021), we first defined the total set of expressed NLRs in transcriptomes from fresh, silica-dried, and herbarium samples, which revealed an overall increase in detectable NLRs in silica-dried samples relative to freshly collected controls for each species (Fig. 3A; Supplemental Table 8), potentially owing to differences in sequencing depth (Supplemental Table 3). By comparison, fewer NLR-related genes were detected for archival RNA from the 1950s, although a full-length receptor was recovered (Fig. 3A). Next, we generated a maximum likelihood phylogeny using the conserved NB-ARC domain to examine NLR relationships among all samples relative to a reference set of functionally validated receptors (Kourelis et al. 2021). Consistent with known NLR phylogenetic relationships (Chia et al. 2024), we observed three distinct clades corresponding to NLR subtype as defined by their N-terminal resistance to powdery mildew 8 (RPW8), toll/Interleukin-1 receptor (TIR), or coiled coil (CC) domains (Fig. 3B). Importantly, these analyses revealed overlap between the NLRs from each of the three collection types, highlighting the diversity of NLRs present within archival RNA and providing new candidates for functional investigation.

Figure 3.

An NLR immune receptor sourced from a 68-year-old herbarium specimen of Clavija costaricana is functionally transferable to Nicotiana benthamiana. (A) Total number of NLRs identified in each transcriptome sample. (B) Maximum likelihood phylogeny of NLR immune receptors based on the central NB-ARC regulatory domain. The outer ring indicates whether each NLR was sourced from herbaria samples, silica samples, or fresh samples or belongs to a reference set of functional NLRs from flowering plants. Branch color represents NLR lineage: RPW8-NLRs (RPW8), TIR-NLRs (TIR), or CC-NLRs (CC). The three tested CC-NLRs (CNLs) are labeled on the tree. Tree scale = substitutions/site. (C) HR cell death phenotypes induced by wild-type (WT) or autoactive (DV) variants of DcCNL-3xFLAG (derived from silica-dried Darlingtonia), CcCNL-3xFLAG (derived from herbarium-pressed Clavija), or the NbZAR1-3xFLAG control in N. benthamiana. Leaves were imaged 7 days post agroinfiltration under bright light as well as UV light, in which cell death responses appear as green fluorescence. (D) Quantification of HR cell death performed 7 d post agroinfiltration. Data from three independent experimental replicates are shown (n ≥ 5 infiltrations per replicate). Statistically significant differences between the means of WT and DV treatments are indicated with an asterisk (Kruskal–Wallis test, P < 0.01; ns = not significant). (E) Immunoblots of all tested NLRs 2 d post agroinfiltration. Total protein levels were visualized with Ponceau S staining. Immunoblotting was repeated twice with similar results. Red asterisks indicate proteins of interest.

Open in new tabDownload PowerPointLink to figure

NLR genes recovered from archival RNA retain function

To determine whether NLRs from archival RNA retain functions in immunity, we synthesized two full-length CC-NLRs tagged with 3xFLAG epitopes for comparison against the well-established NbZAR1 CC-NLR from Nicotiana benthamiana. The CC-NLRs were DcCNL from a Darlingtonia californica sample dried with silica in 2018 and CcCNL from a Clavija costaricana herbarium specimen pressed in 1956. To assess their function, we synthesized and compared wild-type (WT) versions of each receptor against “autoactivated” mutants carrying aspartic acid (D) to valine (V) mutations in the conserved “MHD” motif known to sensitize CC-NLRs and produce immune-related cell death (Adachi et al. 2019). When transiently expressed in N. benthamiana leaves, autoactivated NbZAR1^DV and CcCNL^DV produced appreciable immune-related cell death phenotypes observable under typical light and UV illumination, whereas DcCNL^DV and all WT controls were inactive (Fig. 3C,D). Importantly, immunoblotting indicated sufficient protein accumulation for each construct in N. benthamiana (Fig. 3E). Our results demonstrate that the CcCNL receptor sequence, collected in 1956 and archived at the University of Michigan Herbarium (MICH), is indeed viable and supports the idea that herbaria can be used to explore functionally relevant NLR immune receptors.

Herbarium preserved tissue as a source for coding sequence data

In this study, we present the first demonstration that plant mRNA can be obtained from herbarium specimens. Importantly, this work builds on the small but growing body of knowledge on mRNA stability in eukaryotic tissues (Fordyce et al. 2013; Smith et al. 2014, 2017, 2019; Mármol-Sánchez et al. 2023, 2025). These studies collectively overturn the notion that archival mRNA degrades too fast to be usable for transcriptomics.

The expanding field of historical transcriptomics can begin to provide a source of gene annotation to complement paleogenomics and supply empirical evidence about when and where particular genes were/are expressed. Unlocking the transcriptome of extinct, rare, and/or endangered plants from the convenience of an herbarium can facilitate the discovery of novel genetic variation not evident in the genome sequence alone, such as isoform production. Although this study is limited to data from a single tissue, previous work in mammals has shown that some tissue-specific expression is maintained postmortem (Smith et al. 2019; Mármol-Sánchez et al. 2023, 2025), and future work with herbaria tissue will benefit from similar investigations.

Sequence data show decreased insert length but not the characteristic deamination seen in previous studies

To help avoid cross-contamination and procedural errors that can arise from laboratory work, the RNA extractions were performed on three different days. The low RIN scores recovered for many samples (Supplemental Table 1; Supplemental Document 1) suggest that extracted RNA of low predicted quality can provide reliable transcriptomes. The lack of detectable RNA in two samples that were successfully sequenced has also been recently demonstrated for studies on the woolly mammoth (Mármol-Sánchez et al. 2025). Regardless, libraries were successfully prepared following the manufacturer's protocols. In the libraries, small amounts of DNA were detected in the RNA samples despite DNase treatment. DNA is almost always present in plant samples, even when DNase treatments are performed (e.g., Morales-Briones et al. [2021] retained enough DNA to assemble a complete plastome despite a DNase treatment). Furthermore, the recent work on the mammoth transcriptome showed RNA was extracted despite not using a DNase treatment (Mármol-Sánchez et al. 2025).

We examined the sequences for the characteristic patterns of damage often seen in studies of unpreserved mRNA (Smith et al. 2019; Mármol-Sánchez et al. 2023). Read-pairing loss was more prevalent in the herbarium samples than in the silica-dried or fresh ones. Loss of read pairs could occur when a shorter insert size causes the sequencer to read into adapters. Information would then be lost as adapters are trimmed during bioinformatic processing. An examination of insertion length, calculated as the distance between paired-end reads, revealed the expected pattern: The insert length shortened with increasing sample age, with fresh tissue exhibiting the longest insert size and herbarium tissue the shortest. Read-pairing loss and insert length provide complementary evidence that the sequences shortened, and future work that examines the rate at which these decrease and whether they influence expression patterns in a biologically meaningful way will help researchers better understand the potential insights historical transcriptomics may provide. As researchers examine older samples, it may be important to use modified protocols specifically designed for shorter reads, such as the approach taken by Smith et al. (2019).

To further assess the sequence data, we used a standard approach to examining deamination (Jónsson et al. 2013). The most common sources of deamination are cytidine to uridine (C > U), which is converted to C > T for sequencing, and adenosine to inosine (A > I), converted to A > G for sequencing, which have been detected in previous ancient RNA work (Smith et al. 2019; Mármol-Sánchez et al. 2023). In the present study, we did not find the same characteristic pattern of damage (Fig. 2A; Supplemental Fig. 2). Although our current understanding of archival RNA does not allow for a detailed synthesis of what may cause this difference in deamination patterns, the sampling in this study was ∼70 years old compared with the ∼130-year-old thylacine tissue obtained from museums. Another possible explanation for the observed damage patterns could be our use of the Rcorrector (Song and Florea 2015) approach to correcting base pair errors. Rcorrector is commonly used for transcriptomic assemblies but is not typically used for the study of ancient biomolecules; the program is designed to correct random sequencing errors, and its use may have corrected errors arising from deamination. Furthermore, this tissue is derived from plant material, which may, for a yet-to-be-determined reason, help protect against deamination. Further work on ancient RNA should help uncover the reasons for this pattern of deamination.

Reliable transcriptome assemblies are obtainable from archival RNA

The individual transcriptomes used here were independently and uniformly assembled using many field-standard programs and approaches for sequence read cleaning, de novo transcriptomic assembly, and postprocessing of assembled transcripts (Woodcock-Girard et al. 2025). Despite the loss of reads during trimming, each sample yielded more than 10,000 predicted coding sequences, demonstrating that herbaria are a powerful source of data. The reliability of the inferred coding sequences was evaluated for accuracy using the TransRate procedure (Smith-Unna et al. 2016), which identifies errors common with misassembled contigs and therefore assesses the quality of assembled contigs. The proportion of contigs considered “good” based on the overall quality assessment and scores for the individual quality metrics was similar across the samples, regardless of preservation method. The proportion of bases that were covered by sequence reads showed a 1%–2% increase for the fresh samples compared with the silica-dried and herbarium samples, indicating slightly better coverage of the assemblies from fresh tissue. On a per-contig basis, the samples showed that in both cases the majority of contigs had 100% of their bases covered, indicating that the assembled contigs were concordant with their sequencing reads. The main difference among preservation types was the number of BUSCOs recovered, with the herbarium samples from the 1950s providing the fewest complete BUSCOs. These comparisons, when taken together, demonstrate the reliability of the assembled transcripts assembled from all methods, but herbarium samples provided fewer contigs for downstream analysis. It remains to be seen whether sequencing to greater depth or deploying programs specifically designed for archival RNA will produce assemblies similar to those of fresh tissue.

Silica preservation provides a functional alternative for projects seeking to obtain sequence data

Silica-preserved tissue stored at room temperature for up to 6 months has been shown to provide reliable coding sequence data for population genetic and phylotranscriptomic pursuits (Ruiz-Vargas et al. 2024). This study expanded upon that investigation and examined transcriptome quality after 5 years using a pairwise comparison. The seven specimen pairs consisted of one sample collected and sequenced using standard practices for field-collected RNA, as well as comparable tissue collected at the same time in silica gel and processed 5 years later. Phylogenetically, pairs always clustered together and showed almost no difference in subs/bp (Fig. 1E). Based on our results, we found no advantage to storing the tissue in RNAlater for phylotranscriptomics as opposed to collecting in silica gel. Although this is a small sample size, transcriptomics from silica-collected samples deserve further investigation as they can facilitate collaborations with researchers who lack funding for specialized storage (Tyszka et al. 2025). Furthermore, the herbarium-derived samples also appeared in the expected phylogenetic position (Larson et al. 2020; Carruthers et al. 2024), providing further evidence that they were not the result of cross-contamination. They also did not exhibit abnormally long branches and, thus, appear to be suitable for phylotranscriptomic pursuits.

Expressing an immune system gene from the 1950s

Although bioinformatic methods are powerful tools for data quality assessment, experimental validation remains an important standard. To provide a nonbioinformatic demonstration of the reliability of the obtained coding sequences, we validated them by synthesizing two CC-NLRs and quantifying their ability to produce cell death phenotypes when expressed in N. benthamiana. We demonstrate that an NLR immune receptor from a sample collected in the 1950s can be introduced and expressed in a living plant. This highlights that the predicted contigs are of reliable quality for molecular analyses and that RNA extracted from herbaria is of sufficient quality to yield reliable transcriptomes for researchers interested in studying older specimens. The placement of historic NLRs in a phylogenetic context, combined with functional validation in living plants, serves as a confirmation that the transcriptomes obtained and assembled are accurate representations of genes that existed within the plants sampled. This finding also demonstrates the possibility of synthesizing genes from herbarium specimens, opening a new source of tissue for RNA extraction for researchers seeking coding sequences from older samples.

Historical transcriptomics as a data source for inferring the tree of life

A body of literature has slowly developed, highlighting that RNA in seeds maintains its integrity for germination under extreme conditions (Sallon et al. 2008; Gros-Balthazard et al. 2021), and this can even be leveraged for RNA extraction (Rollo 1985; Fordyce et al. 2013). Ancient mummies have had their mRNA extracted (Venanzi and Rollo 1990), and canids found in permafrost have had their transcriptomes profiled (Smith et al. 2019). The expression profiles of museum tissue show that some signals of tissue-specific expression are maintained (Mármol-Sánchez et al. 2023). Furthermore, microRNAs have been extracted from a variety of eukaryotic species (Keller et al. 2017; Shaw et al. 2019; Fromm et al. 2021), and with these data, scientists are gaining new insights into the gene regulation of extinct organisms (Mármol-Sánchez et al. 2023, 2025).

Here, we show that transcriptomes from herbarium specimens can be used to infer the tree of life. Phylogenetics is a field constantly adapting to incorporate novel sources of data. The One Thousand Transcriptomes project (One Thousand Plant Transcriptomes Initiative 2019) generated a wealth of data that researchers have been able to build on and complement with their own data sets. Furthermore, transcriptomes generated for differential expression experiments are often made publicly available and can be repurposed for phylotranscriptomics (Wong and Peakall 2022; Arreguin et al. 2025). When pursuing new phylotranscriptomic analyses, these existing data sources may be complemented with herbarium data to include phylogenetically contentious and difficult-to-acquire lineages.

The predominant data source for phylogenomic analyses using herbaria has been probe data sets, which rely on DNA extraction (Johnson et al. 2019; Zuntini et al. 2024). Although probe-based DNA sequencing provides a powerful method of targeting specific genes, these methods are largely limited to the set of genes targeted a priori, and unless the same probe set is used, the data sets are not compatible. As probe data sets are often built using transcriptomes as a guide, the coding regions sequenced from transcriptomes can be used to complement existing probe data sets and transcriptomic data. The inclusion of transcriptomic data alongside probe data sets has allowed researchers to refine not only phylogenetic hypotheses but also patterns of gene and genome duplication (Lagou et al. 2024). The ability to obtain transcriptomes from herbaria has the potential not only to expand our knowledge of the tree of life but also to provide insights into the processes that have shaped and continue to shape it.

Limitations of this study and future directions for understanding herbaria as a source of transcriptomic data

Despite the insights provided by the present study, several questions about historical RNA remain unanswered. Our focus on mRNA raises the question of how reliably microRNAs may also be obtained from samples like ours. Furthermore, in this study, we lacked replication of our samples; therefore, we cannot reliably analyze biologically meaningful differences in expression patterns. Expression patterns can help answer questions such as whether the immune response to past exposure to pathogens is still identifiable in herbarium specimens. It could also be assessed whether the stress of desiccation is reflected in expression patterns. Another limitation we have is the lack of reference genomes for the samples in question. We are unable to analyze exon–exon junctions and exonic versus intronic reads. Without replicates and complete genomes, we were unable to determine whether there is biased retention of genes and exonic regions. Further research may help answer this and other questions, such as whether different plant tissue types, such as succulent versus papery tissues, vary in how well they preserve mRNA. How does the length and type of preservation influence degradation, and at what age (and to what extent) will deamination play a role in herbarium-based transcriptomic studies? The ability to access mRNA from herbaria opens the door for these and many other questions. Here, we provide the first step by establishing that usable plant mRNA may be obtained from herbarium-sourced tissue.

In summary, our work (1) demonstrates that mRNA may be obtained from pressed herbarium specimens and years-old silica-dried samples; (2) this mRNA may be sequenced and assembled into reliable transcriptomes that can be used for a variety of pursuits including phylotranscriptomics; and (3) using transcriptomes from herbarium specimens, we can predict NLR genes, indicating the continued presence of genetic information including essential variation underlying disease resistance. To date, work on historical RNA remains almost nonexistent across the tree of eukaryotes, with great potential for insights across broad disciplines in biology (Friedländer and Gilbert 2024; Tyszka et al. 2025). This work further highlights the value of the world's increasingly threatened herbaria (Davis 2024) by explicitly combining the historical collections with insights from high-throughput sequencing. Our results emphasize the importance of protecting these invaluable resources as novel uses for them continue to arise.

Sampling scheme

For this study, we chose to sample members of the flowering plant order Ericales, a clade of more than 13,000 species that contains important crops such as tea, kiwi, and blueberries; ornamentals; parasitic lineages; and carnivorous plants (Rose et al. 2018; Larson et al. 2020; Carruthers et al. 2024). We used samples stored under three preservation approaches: RNAlater followed by freezing (i.e., “-Fresh” samples), silica-dried tissue collected in abundant silica gel in the field (i.e., “-Silica” samples), and herbarium tissue (i.e., “-Herbarium” samples). Our sampling included seven paired silica-dried leaves and corresponding fresh leaves, in six cases collected from the same individual plant and in all cases at the same time and from the same tissue type. Our sampling also included three herbarium specimens collected between 1950 and 1995 and three additional silica-dried samples (Supplemental Table 9).

RNA extraction

We obtained mRNA from three herbarium samples and 10 silica-dried samples, which were stored at room temperature for ∼5 years (Supplemental Table 1). Tissue grinding, RNA extractions with the DNase treatment, and quality assessment were performed using the same procedure of Ruiz-Vargas et al. (2024). To ensure the extraction procedure did not introduce cross-contamination, the extractions were performed at three different times. The silica-dried samples were extracted on June 4, 2023, and June 11, 2023, and the herbarium samples were extracted on June 12, 2023. For two of the samples, the RNA concentrations measured with a Qubit fluorometer (Thermo Fisher Scientific) were below the minimum concentration required for analysis with a TapeStation (Supplemental Table 1; Supplemental Document 1). Ribosomal RNAs were removed using the QIAseq FastSelect kit (Qiagen), and the RNA-seq libraries were prepared with the Kapa hyperstranded mRNA library kit (Roche). Paired-end 150 bp was carried out on a single lane of an Illumina NovaSeq X.

Transcriptome assembly and quality assessment

Transcriptomes for all samples, including the fresh samples from Carruthers et al. (2024), were assembled using the same procedure described below. The quality of the raw sequence reads was initially assessed using the program FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Single base pair errors were corrected using the program Rcorrector v1.0.4 (Song and Florea 2015), with reads deemed unfixable removed using the unfixable_filter.py script from Morales-Briones et al. (2021). Trimmomatic v0.39 (Bolger et al. 2014) was used to remove adapter sequences and trim low-quality reads.

Kraken2 v2.1.2 (Wood et al. 2019) was used to classify any material matching the standard protist and fungal database, and any classified material was removed as a contaminant. The same procedure was then used to classify reads that matched the chloroplast and then the mitochondrion database from the program Kakapo (Ramanauskas and Igić 2023). The remaining reads were classified as nuclear reads.

Trinity v2.1.1 (Haas et al. 2013) was used to assemble the nuclear sequence reads into transcripts. Chimeric transcripts were detected using the procedure of Yang and Smith (2013) with a database composed of the peptide sequences for Arabidopsis thaliana (TAIR10) and Actinidia chinensis (Red5_PS1_1.69.0). The abundances of the nuclear reads were then quantified using Salmon v1.3.0 (Patro et al. 2017) and clustered into the highest supported transcript using Corset v1.09 (Davidson and Oshlack 2014). The open reading frames were predicted using TransDecoder v5.5.0 (https://github.com/TransDecoder/TransDecoder), and the best supported ORF was retained. In total, we generated 13 new transcriptomes, three of which originated from herbarium material, and reassembled seven previously sequenced transcriptomes.

The completeness of the assemblies was assessed using BUSCO v.5.8.2 (Manni et al. 2021) with the “‐‐mode transcriptomes” option and the viridiplantae_odb10 database. The quality of the transcripts was assessed using TransRate2 (https://github.com/ericbretz/transrate2), an implementation of the TransRate procedure (Smith-Unna et al. 2016). TransRate2 is a tool designed to evaluate the quality of de novo transcriptome assemblies, expanding on TransRate by implementing the splice-aware HISAT2 mapping approach (Kim et al. 2019). Specifically, the score a transcript is penalized based on the following criteria: (1) the transcriptome shows evidence of chimeric sequences; (2) the transcript has low coverage from the raw reads; (3) the reads mapping to the transcript are not properly paired; and (4) a base in the transcript does not correspond to the base of the sequences that are mapped to it. Given sufficient penalties, the transcript will be considered “bad quality”; otherwise, transcripts will be considered “good quality.” The mapping rates obtained from TransRate2 were used to calculate the total coverage across all bases and the per-contig base coverage.

RNA damage analysis

To determine the length of the insert size associated with each read pair, we used HISAT2 to map the cleaned raw reads back to the assembled transcriptome. The output BAM file was then used to calculate the distance between the reads using a novel Python script (see Data access). The program mapDamage2 v.2.3.0a0-10bc442 (Jónsson et al. 2013) with “‐‐merge-libraries” and all other settings set to default was used to test for deamination, with the same BAM file generated for the insert size analysis used as input.

Phylogenetic analysis

We conducted a phylogenomic analysis by identifying homologous gene clusters from the assembled transcriptomes of fresh, silica-dried, and herbarium samples. The phylogenomic analysis also included a previously assembled transcriptome of Cornus florida (NCBI RefSeq assembly [https://www.ncbi.nlm.nih.gov/refseq/]: GCF_030987335.1). CD-Hit v4.8.1 (Fu et al. 2012) with a sequence identity threshold of 0.99 was used to reduce redundancy in the transcriptomes. An all-by-all BLASTN (Altschul et al. 1990) search on the transcriptomes was conducted for initial homology detection, with an e-value of 10 and the maximum target sequences set to 1000. The homolog clusters, as detected by BLASTN, were then further refined using mcl v14-137 (Van Dongen 2008), with an inflation value of 1.4. All remaining homolog clusters with at least 10 taxa were retained for phylogenomic analysis.

The putative homologs were aligned using MAFFT v7.490 (Katoh and Standley 2013) with the “‐‐auto” and “‐‐maxiterate 1000” settings. The phyx program pxclsq v1.2 (Brown et al. 2017) was used to remove all alignment columns with occupancy <10%. Maximum likelihood, as implemented in IQ-TREE v1.6.12 (Nguyen et al. 2015), with the GTR model of evolution and gamma rate variation, was used to infer a phylogenetic tree for each homolog. We chose the GTR+G model because it is the most parameter-rich within the GTR family, and other commonly used time-reversible models are nested within it. From the inferred trees, all tips that had branch lengths of 0.2 subs/bp and/or were 10× longer than their sister or 0.3 subs/bp absolute value were removed, and all clades in trees that consisted of gene sequences from a single sample were collapsed to retain only the one with the most information in the alignment. The refined homologs were then run through the same procedure three more times to refine the phylogenetic data set further.

Orthologous groups were extracted from the homolog trees using the maximum inclusion method (Yang and Smith 2014), requiring a minimum of 10 taxa and removing any branch lengths with a relative value of 0.2 subs/bp or an absolute value of 0.3 subs/bp. The sequences corresponding to the extracted ortholog trees were then aligned using MAFFT, cleaned using pxclsq as before, and used to estimate gene trees with IQ-TREE using the same settings as the homology procedure. A species tree was then inferred using the maximum quartet support species tree method implemented in ASTRAL v5.7.8 (Zhang et al. 2018). The coalescent units were converted to subs/bp using the median branch length value of all concordant branches in the gene trees as implemented in the program BES (Walker et al. 2022).

NLR identification and phylogenetic analysis

NLR identification was performed using the default settings of the NLRtracker pipeline, which is available online (https://github.com/slt666666/NLRtracker). Phylogenetic analysis was performed by comparing the central NB-ARC domains of fresh, silica-dried, and herbarium samples against a reference set of flowering plant NLRs (Kourelis et al. 2021). Isolated domains were aligned using MAFFT and trimmed using trimAl (Capella-Gutiérrez et al. 2009). The resulting alignments were used for maximum likelihood phylogenetic analysis with IQ-TREE v2.0.3 using the “JTT+F+G4” model selected by ModelFinder (Kalyaanamoorthy et al. 2017) and 1000 ultrafast bootstrap replicates to assess support (Hoang et al. 2018). The phylogenetic tree was subsequently rendered using TVbot (Xie et al. 2023).

Plant growth conditions

N. benthamiana plants were grown in soil under controlled conditions, with a long-day photoperiod (16-h light; 160–200 µE, using Sylvania F58W/GRO fluorescent lighting) at 22°C.

Cloning, immunoblotting, and cell death assays

Constructs for cell death assays were assembled using the Modular Cloning (MoClo) system with the following components: pICH85281 (mannopine synthase promoter + W [MasWpro], Addgene 50272), pICSL50007 (AddGene 50308), pICSL60008 (Arabidopsis heat shock protein terminator [HSPter], TSL SynBio), binary vector pICH47742 (Addgene 48001), and the appropriate synthesized NLR gene or variant. The CC-NLRs synthesized were 3xFLAG epitope tagged WT and autoactive (DV) variants of DcCNL (identified from the Darlingtonia-Silica transcriptome) and CcCNL (identified from the Clavija-Herbarium transcriptome). NbZAR1 constructs were generated using parts from an established Golden Gate toolkit (Harant et al. 2022). All resulting constructs were transformed into Agrobacterium tumefaciens GV3101 (pMP90) by electroporation. Cell death assays were performed by Agrobacterium-mediated transient expression, in which fully expanded leaves of 4-week-old N. benthamiana plants were infiltrated with A. tumefaciens strains containing a binary expression plasmid mixed at a 1:1 ratio with a suspension of A. tumefaciens containing the p19 silencing suppressor. Bacterial cells were resuspended in a fresh infiltration buffer (10 mM MES-KOH, 10 mM MgC₂, and 150 mM acetosyringone at pH 5.6) and adjusted to an OD600 of 0.3. Immune-related cell death phenotypes were assessed using an HR index ranging from zero (no visible symptoms) to seven (fully confluent cell death). For immunoblotting, protein samples were prepared from six N. benthamiana leaf discs (8 mm diameters) at 2 days after agroinfiltration. Leaf samples were then ground to a fine powder and homogenized in 150 µL of 2× LDS sample buffer (NuPAGE LDS sample buffer). Immunoblotting was performed with antiflag antibody coupled with HRP (Sigma-Aldrich A8592, 1:5000 dilution), and total protein loading was visualized with Ponceau S solution (Sigma-Aldrich P7170).

Statistical analysis

We estimated differences in the recovery rates of BUSCO loci from transcriptomes derived from fresh, silica-dried, and herbarium tissue, accounting for differences in sequencing depth. We fitted a multinomial logistic regression using nnet v7.3-20 (Venables and Ripley 2013), with counts of complete (summing single-copy and duplicated loci), fragmented, and missing BUSCOs as responses and with storage condition and read depth as predictors. To assess the inferences of this model, we visualized the predicted proportions of recovery from each category per sequencing depth. All analyses were conducted in R v4.5.1 (R Core team 2021).