Regulatory dynamics distinguishing desiccation tolerance strategies within resurrection grasses

Abstract Desiccation tolerance has evolved recurrently in grasses using two unique strategies of either protecting or dismantling the photosynthetic apparatus to minimize photooxidative damage under life without water (anhydrobiosis). Here, we surveyed chromatin architecture and gene expression during desiccation in two closely related grasses with distinguishing desiccation tolerance strategies to identify regulatory dynamics underlying these unique adaptations. In both grasses, we observed a strong association between nearby chromatin accessibility and gene expression in desiccated tissues compared to well‐watered, reflecting an unusual chromatin stability under anhydrobiosis. Integration of chromatin accessibility (ATACseq) and expression data (RNAseq) revealed a core desiccation response across these two grasses. This includes many genes with binding sites for the core seed development transcription factor ABI5, supporting the long‐standing hypothesis that vegetative desiccation tolerance evolved from rewiring seed pathways. Oropetium thomaeum has a unique set of desiccation induced genes and regulatory elements associated with photoprotection, pigment biosynthesis, and response to high light, reflecting its adaptation of protecting the photosynthetic apparatus under desiccation (homoiochlorophyly). By contrast, Eragrostis nindensis has unique accessible and expressed genes related to chlorophyll catabolism, scavenging of amino acids, and hypoxia, highlighting its poikilochlorophyllous adaptations of dismantling the photosynthetic apparatus and degrading chlorophyll under desiccation. Together, our results highlight the complex regulatory and expression dynamics underlying desiccation tolerance in grasses.


| INTRODUCTION
Water deficit was the most pervasive challenge faced by charophyte green algae during the colonization of land, and it has been a continual force shaping the evolution and diversification of plants. Plants have evolved a versatile arsenal of strategies to avoid or overcome water limitations. At the extreme, a small group of plants can survive prolonged desiccation for months to years until the return of water.
Desiccation tolerance, or the ability to survive atmospheric drying, has been investigated for more than a century (Bewley, 1979;Bristol, 1919). Anhydrobiosis, or "life without water," causes cells to enter a glassy or solid state and puts tremendous stress on all of the macromolecules and organelles. Polysomes are lost during drying, effectively halting protein synthesis, but they are quickly restored during rewetting (Bewley, 1979). Mitochondria and chloroplasts of most species swell and become ill defined during desiccation, and only species able to restore these organelles are associated with recovery from dry conditions (Wellburn & Wellburn, 1976;Sherwin & Farrant, 1998). Biochemical, physiological, and molecular mechanisms underlying desiccation tolerance have been thoroughly described, but comparatively little is known about the regulatory dynamics controlling this trait and what distinguishes desiccation tolerance from more typical drought responses (Gechev et al., 2021).
Desiccation tolerance requires tight coordination of numerous cellular processes and pathways to protect plants from the effects of water deprivation and photooxidative damage under excess light.
Various osmoprotectants, heat shock proteins, reactive oxygen species scavengers, changes in membrane lipid composition, and late embryogenesis abundant (LEA) proteins have well-documented roles in desiccation tolerance (Hoekstra et al., 2001). Along with these conserved responses, desiccation tolerant plants utilize two distinct strategies to mitigate photo-oxidative damage. Homoiochlorophyllous species retain and protect their chlorophyll and thylakoids during desiccation, whereas poikilochlorophyllous species break down and rebuild chlorophyll, thylakoids, and components of the photosynthetic apparatus under anhydrobiosis using specialized membrane free plastids (desiccoplasts) (Tuba et al., 1994). These divergent strategies are associated with distinct rehydration and photosynthetic kinetics and have utility in different environments.
Vegetative desiccation tolerance is thought to have evolved through rewiring preexisting pathways, perhaps convergently across independent lineages. Desiccation tolerance mechanisms are broadly conserved between vegetative tissues and seeds, prompting the long-standing hypothesis that desiccation tolerance evolved through rewiring existing seed pathways (Costa et al., 2017;Illing et al., 2005;VanBuren, 2017;VanBuren, Man Wai, et al., 2018). However, recent comparative experiments suggest this connection is more nuanced, as seed pathways also play an important role in typical drought responses of desiccation sensitive species . Water deficit responses are broadly regulated by the phytohormone abscisic acid (ABA), and the ABA regulon has a major role in desiccation tolerance (Gaff & Oliver, 2013;Manfre et al., 2009;Shinozaki and Yamaguchi-Shinozaki & Yamaguchi-Shinozaki, 2007) and drought responsive pathways (Daszkowska-Golec, 2016). Several transcription factor families including dehydration-responsive element-binding factor (DREB), basic leucine zipper (bZIP), and NAM-ATAF-CUC2 (NAC) are regulated by ABA and play integral roles in drought and desiccation responses (Nakashima et al., 2014;Takasaki et al., 2015;Wang et al., 2019;Yoshida et al., 2015). The unique and overlapping pathways between drought and desiccation and the regulatory machinery underlying these stress responses remain poorly understood.
The regulation of gene activity based on proximity to cis-sequences or entire regions has been a long standing interest of the biological community (Baker, 1968) and has been the source of significant advancement in our understanding of eukaryotic gene regulation (Eissenberg, 1989). High quality analysis of plant chromatin structure and accessibility has been under investigation for several decades (Bowler et al., 2004). Recently, the use of assay for transposase accessible chromatin sequencing (ATAC-seq) has become a powerful tool to reveal binding sites of transcription factors and correlates well with changes in gene expression in various tissues and cell types (Lu et al., 2017). ATAC-seq has even been used to draw correlations between species to identify evolutionarily conserved chromatin regions (Lu et al., 2019). Chromatin dynamics have been profiled under several abiotic stresses including cold, heat, and flooding Liang et al., 2021;Reynoso et al., 2019;Wang et al., 2021), but little is known about chromatin changes during desiccation. One small-scale investigation identified a desiccation induced promoter, but this study did not identify underlying cis-regularoty elements or provide a (Smith-Espinoza et al., 2007) clear picture of how chromatin is affected by the loss of water or the link between chromatin changes and gene expression.
Here, we surveyed changes in chromatin architecture, gene expression, and regulatory dynamics under desiccation in two related resurrection grasses with contrasting photoprotective strategies. We collected parallel datasets from the homoiochlorophyllous (chlorophyll retaining) model resurrection grass Oropetium thomaeum and the poikilochlorophyllous (chlorophyll degrading) grass Eragrostis nindensis.
These species are found in the Chloridoideae subfamily of grasses, a group of diverse, stress tolerant species that includes the important underutilized crops teff (Eragrostis tef ) and finger millet (Eleusine coracana). O. thomaeum also has the smallest diploid genome among grasses ($240 Mbp) , and the compact monoploid genomes within this subfamily enable association between accessible chromatin regions (ACRs) and nearby genes. Previous studies in O. thomaeum showed an induction of pathways and gene families related to seed development such as ELIPs and LEAs during desiccation . Here, we sought to identify regulatory elements that enabled these unique pathways compared with typical drought responses. With the use of this comparative system, we show that significant shifts in chromatin architecture are associated with desiccation and photoprotective responses in resurrection grasses. We identified a core set of genes, accesable chromatin peaks, and cis-regulatory elements that underly conserved desiccation tolernace responses across these two grasses. We also found distinct regulatory dynamics that are associated with the different strategies of retaining or degrading the photosynthetic apparatus under anhydrobiosis.

| Dynamic patterns of chromatin accessibility under desiccation
We surveyed the changes in chromatin architecture and gene expression underlying desiccation tolerance using a combination of ATACseq and RNAseq in O. thomaeum and E. nindensis. Mature plants were slowly dried over a period of 10 days without watering to trigger desiccation, and samples for RNAseq and ATAC-seq were collected in parallel from leaves of desiccated and well-watered plants (Figure 1).
RNAseq reads were quality trimmed, pseudo-aligned to the O. thomeaum and E. nindensis gene models using Kallisto (Bray et al., 2016), and differentially expressed genes were identified using Sleuth (Pimentel et al., 2017). In O. thomaeum, 5,187 genes had higher expression under desiccation, and 5,354 genes had lower expression compared with well-watered (q < 0.05). Similar timepoints in E. nindensis had 8,071 genes with higher expression under desiccation and 8,971 genes with lower expression compared with well-watered.
This represents roughly one third of the genes in the O. thomaeum genome and almost half of the genes with detectable expression. In E. nindensis, the differentially expressed genes represent about 15% of all genes and 30% of the genes with detectable expression. These expression dynamics highlight the complex transcriptional reprogramming required for the successful deployment of desiccation tolerance  (Table S1). Genes with the highest differential expression in desiccated E. nindensis plants included two late embryogenesis-related proteins (LEAs; En_0084954, En_0002557), an aldehyde dehydrogenase (En_0060658), and a glucosyl transferase (En_0017075 ; Table S2).
We utilized ATAC-seq to profile chromatin openness and link expression dynamics with putative cis-regulatory elements. Nuclei were isolated from flash frozen samples, and $50,000 nuclei were used to construct ATACseq libraries for desiccated and well-watered samples in each species. Naked DNA controls were also collected to mask regions showing TN5 insertional bias. Trimmed reads were aligned to the O. thomeaum V2 or E. nindensis V2 genomes using bow-tie2 (Langmead & Salzberg, 2012), and peaks of accessible chromatin were called using Genrich (Gaspar, 2018).  (Tables S3 and S4).
In total, we identified 26,225 ACRs in O. thomaeum and 83,455 ACRs in E. nindensis occupying 6.4% and 4.9% of the genome, respectively ( Figure 1a,c). Roughly 85% and 82% of all ACRs were near expressed genes in O. thomaeum and E. nindensis, respectively. For both species, the 0-500 bp bin upstream of the transcription start site (TSS) has an over-representation of ACRs based on the relative proportion of the genome occupied by this range (Figure 1b,d). We also observed that ACR density decreases with distance from the TSS. The O. thomaeum and E. nindensis genomes are highly compact compared to other grassses VanBuren et al., 2015;VanBuren et al., 2020), and expressed genes and their flanking regions (i.e., 10 kb upstream to 1 kb downstream) cover 65.8% and 59.4% of their genomes, respectively. Not all ACRs were found near canonical genes. Of the open regions, 7.2% and 14.2% are near non-expressed genes, and 7.5% and 21.5% of the regions were not near any genes.

| Chromatin architecture is tightly linked to the regulation of desiccation responsive genes
Because chromatin architecture is believed to impact the expression level of nearby genes, we investigated the relationship between chromatin openness and gene expression using two different methods.
Expressed genes typically have ACRs near the TSS (Figure 1b,d), and we used deepTools (Ramírez et al., 2016) to aggregate the genomic regions near differentially expressed genes. We found that ACRs near genes are significantly more open when the gene has higher expression ( Figure 2). This correlation was strongest for genes with higher expression under desiccation in both O. thomaeum and E. nindensis (Figure 2a,b). Genes with upregulated expression in well-watered tissues had a distinct peak in openness near the TSS, but this was less pronounced than desiccation induced genes. We then looked at differentially expressed genes that had a differentially accessible region in each bin and plotted the best-fit line for each combination. The R 2 value for that line-fit drops off inside the 5000 to 2000 bp region and is highest for the 0 to 500 bp region upstream of the TSS ( Figure S5).
To quantify the impact of chromatin accessibility on expression dynamics, we split the differentially expressed genes into four categories depending on the direction of differential gene expression and To enable detailed comparisons of expression patterns and chromatin dynamics across species, we utilized a set of syntenic orthologs between the two grasses . Roughly half of the dif- and S8). GO terms associated with syntelogs that are similarly downregulated in desiccation are related to photosynthesis, core and secondary metabolism, cell wall processes, and hormone metabolism (Table S8). Broadly, this suggests these two grasses utilize fundamentally similar mechanisms to survive desiccation.

| Cis-regulatory elements associated with desiccation tolerance
With the use of these datasets of ACRs and gene expression, we identified putative cis-regulatory elements that are involved in the regulation of desiccation responses. We extracted sequences within differential  (Bailey, 2020). This motif may play a central role in anhydrobiosis related processes in grasses.
F I G U R E 3 Contrasting enrichment patterns of desiccation-associated genes in Oropetium thomaeum and Eragrostis nindensis. The proportion of syntelogs with similar, unique, or opposing expression patterns is plotted for O. thomaeum (a) and E. nindensis (b). Enriched gene ontology terms (GO terms) for genes that have more accessible chromatin and higher expression under desiccation are plotted for O. thomaeum (c) and E. nindensis (d). GO terms are transformed using multidimensional scaling to reduce dimensionality, and terms are grouped by semantic similarities. GO terms with previously characterized roles in desiccation and photoprotection responses are highlighted. The color of the circles represent significance and size represents the number of genes in that group.

| Regulatory dynamics distinguishing photoprotective strategies under desiccation
O. thomaeum and E. nindensis utilize different strategies to mitigate photooxidative damage under anhydrobiosis. O. thomaeum retains and protects, whereas E. nindensis degrades and resynthesizes chlorophyll, thylakoid membranes, and components of the photosynthetic apparatus during desiccation and rehydration cycles (Bartels & Mattar, 2002). We searched for unique ACRs, expression dynamics, and putative cis-regulatory elements that distinguish the distinct desiccation tolerance strategeis in these two grasses. The sets of syntenic orthologs and ACRs between O. thomaeum and E. nindensis were used for comparative analyses.
Syntenic orthologous or species-specific genes that are uniquely upregulated in O. thomaeum compared with E. nindensis have functions related to cellular response to heat and high light, photoprotection, protein localization to chloroplast, regulation of seed germination, and regulation of cell cycle (Table S8)

| Chromatin architecture and anhydrobiosis
Chromatin dynamics have not been described for desiccated vegetative tissues, but chromatin is highly condensed and compacted in desiccated seeds (van Zanten et al., 2011). Differential chromatin accessibility is established in developing seeds, and gene level modifications enable rapid transcription for germination related processes upon hydration (Fransz & de Jong, 2011). We observed a high signalto-noise ratio (FRiP) of the ACRs in desiccated leaf tissue stemming from a significantly lower background level of reads compared with well-watered leaf tissue. This enrichment resulted in a higher abundance of reads in open peaks, potentially indicating a major difference in chromatin state in the nucleus of desiccated cells.
We hypothesize several factors that may contribute to the low background signal and robust peaks observed under desiccation. The first is that all nuclei in the desiccated plants are more condensed, and this compaction and associated factors prevent Tn5 accessibility in background regions. The more parsimonious explanation is that desiccated cells and nuclei are highly uniform with fewer cells actively undergoing transient processes that open chromatin, such as transcription or replication. Under desiccation, all cell types must arrest their normal developmental or metabolic functions and enact a series of highly coordinated responses to successfully prepare for anhydrobiosis and subsequent rehydration. In turn, this would "synchronize" cells and produce the observed tight coordination between expression dynamics and chromatin architecture that is typically only found in single cell data (Farmer et al., 2021). Drying of this magnitude is experienced by all cells, and risky processes such as growth and photosynthesis during the stresses of anhydrobiosis would be detrimental and preferentially avoided. Thus, a clear, tightly regulated signal from all cells in the desiccated samples is achieved. This is supported by the massive transcriptional reprogramming we observed, and the clear enrichment of GO terms related to photoprotective and anhydrobiosis processes with comparatively few background pathways. This pattern of clearer peaks and more accessible chromatin under desiccation contrasts what has been observed under drought and salinity (Raxwal et al., 2020) but is similar to cold stress (Zeng et al., 2019).

| Evolutionary dynamics of photoprotective strategies under desiccation
O. thomaeum and E. nindensis are found within the Chloridoideae subfamily of grasses, a group of stress tolerant C4 species with superior drought, heat, and salinity tolerance (Marcum, 1999;Peterson et al., 2001). Evolutionary access to existing resilience traits within this subfamily likely enabled the independent evolution of desiccation tolerance in O. thomaeum and E. nindensis (Pardo & VanBuren, 2021).
With the use of detailed comparative genomics approaches with syntenic orthologs, we identified a core set of genes and cis-regulatory regions that are induced during desiccation in O. thomaeum and E. nindensis. Integration of ATAC-seq and RNAseq data allowed us to filter out spurious associations, and we identified numerous enriched genes with previously described roles in desiccation. Among the stress responsive cis-elements, we identified numerous binding motifs for the seed desiccation transcription factor ABI5, supporting the hypothesis that vegetative desiccation tolerance evolved from rewiring existing seed pathways (Oliver et al., 2000;Oliver et al., 2005). This hypothesis is contentious, and previous work in the monocot Xerophyta humilis failed to find evidence linking desiccation processes to the canonical seed development transcription factor network of LEC1, ABI3, ABI5, and others (Lyall et al., 2020). Our results provide a direct link between expression dynamics and regions of accessible chromatin with seed related regulatory motifs, but this association needs to be further tested using transcription factor binding data.
Species-specific changes in chromatin architecture and expression dynamics during desiccation reflect the distinct strategies that E. nindensis and O. thomaeum utilize to mitigate photooxidative damage. We identified a cluster of desiccation associated GO terms related to chlorophyll processes, UV-A responses, photoprotection, and pigment metabolism that was unique to O. thomaeum. These orchestrated responses in O. thomaeum reflect homoiochlorophyly, or the strategy to protect chlorophyll and photosynthesis related macromolecules during desiccation. Consistent with this, O. thomaeum has more ELIPs than E. nindensis (when accounting for ploidy), and all of the ELIPs have massive shifts in chromatin accessibility and high expression under desiccation. By comparison, only two ELIPs in E. nindensis have associated ACRs, and few are differentially expressed under desiccation. Instead, ELIPs in E. nindensis have high expression during rehydration and likely function in protecting leaves as they resynthesize and repair their photosynthetic apparatus .
We identified a cis-regulatory motif associated with the drought transcripton factor DREB upstream of each ELIP in O. thomaeum. This motif is missing from the Arabidopsis ELIP orthologs and may repre- Two grams of leaf tissue was collected from each plant and frozen in liquid nitrogen immediately. Tissues were ground into coarse powder and aliquoted into multiple 1.5 ml tubes with $100 mg per tube. The samples were then subjected to nuclei isolation for ATAC-seq or RNA extraction for RNAseq.

| ATAC-seq library construction
Crude nuclei extracts were prepared using the protocol reported in Lu et al. (Lu et al., 2017). Briefly, approximately 0.5 g of liquid nitrogen ground tissue was resuspended in 10 ml of isolation buffer (

| ATAC-seq analysis
Paired-end raw reads were trimmed using Trim_galore (v0.6.6) to remove Nextera adapters and low-quality bases (Krueger, 2018). Filtered reads were then aligned to the O. thomaeum v2.1, or E. nindensis Peaks associated with ACR were called using Genrich (v0.6.1) (Gaspar, 2018) with options specific for ATAC-seq. Additionally, reads mapped from naked DNA were used as control for O. thomaeum. The resulting narrowPeak files and the read files (bam format) were used in the R library DiffBind (v3.4.11) (Ross-Innes et al., 2012;Stark et al., 2011) to determine differentially open regions. DESeq2 (Love et al., 2014) was used to quantify differential ACRs, with normalization using the background library size.

| Integrating ATACseq and RNAseq datasets
A genome-wide assessment of ACR distribution was performed using deepTools (v3.4.3) (Ramírez et al., 2016). ACRs were classified into different groups based on their proximity to genes, and comparisons of peak openness and differential gene expression were performed.
These correlations were used as input for gene ontology and ciselement enrichment analyses.
We identified genes near ACRs using bedTools (v2.29.2) (Quinlan, 2014). ACRs were grouped into different categories based on their distance to genic elements including: 10,000 to 5,001, 5,000 to 2,001, 2,000 to 1,001, 1,000 to 501, 500 to 0 bp upstream (5 0 ) of the TSS; overlapping with the gene; and 0 to 1,000 bp downstream (3 0 ) of the transcription termination site (TTS). On the basis of the distribution and correlation of chromatin openness and gene expression patterns, we chose the range of 0-3,000 bp upstream of gene TSSs as the cutoff of associating ACRs with specific genes. The O. thomaeum and E. nindensis genomes are relatively compact compared to other grass species, and putative cis-regulatory regions are in close proximity to genes.
Genes were grouped by their differential expression patterns and ACR, and GO term enrichment analyses were performed using the R library topGO (v2.46) (Alexa & Rahnenführer, 2009  Motifs associated with ACRs within 3 Kbp upstream of differentially expressed genes were identified, and the Motif Discovery and Enrichment Analysis (XSTREME) algorithm from Multiple Em for Motif Elicitation (MEME) (v5.4.1) was used to identify enriched motifs in our data and associate these motifs with known transcription factor binding sites .

| Comparative genomic analyses
Comparative genomic analyses were conducted between O. thomaeum and E. nindensis to identify shared and species-specific sets of differentially expressed genes and ACRs dynamics. Syntenic orthologs were identified between the two closely related grasses using the python version of MCScan (https://github.com/tanghaibao/jcvi/wiki/MCscan-( Python-version)) (Wang et al., 2012). The chromosome scale O. thomaeum genome was used as an anchor, and the two sets of homeologous genes in the E. nindensis genome were mapped to the single corresponding syntenic ortholog in O. thomaeum. Genes in each species were then classified as syntenic or non-syntenic, and these two designations were incorporated with differential expression analyses, ACR dynamics, and enriched motif analyses. Enriched GO terms and putative cis-element motifs were identified for syntenic orthologs with conserved and species-specific expression patterns as described above.

ACKNOWLEDGMENTS
This work is supported by the US National Science Foundation (NSF) Grant MCB-1817347 (to R.V.). We thank members of the Water and Life Interface Institute (WALII), supported by NSF DBI grant 2213983, for helpful discussions.

CONFLICT OF INTEREST
The Authors did not report any conflict of interest.

DATA AVAILABILITY STATEMENT
The raw RNAseq and ATACseq data are available from the National Center for Biotechnology Information (NCBI) Short Read Archive.
Raw data for this project can be found under BioProject accession no. PRJNA807505.