# Gene and Transposable Element Expression Evolution Following Recent and Past Polyploidy Events in Spartina (Poaceae)

Delphine Giraud, Oscar Lima, Mathieu Rousseau-Gueutin, Armel Salmon, Malika Ainouche

<figure><img src="image_1.png" /></figure>

## To cite this version:

Delphine Giraud, Oscar Lima, Mathieu Rousseau-Gueutin, Armel Salmon, Malika Ainouche. Gene and Transposable Element Expression Evolution Following Recent and Past Polyploidy Events in Spartina (Poaceae). Frontiers in Genetics, 2021, 12, pp.589160. 10.3389/fgene.2021.589160. hal-03216905

HAL Id: hal-03216905

https://hal.science/hal-03216905v1

Submitted on 4 May 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

Check for
updates

# Gene and Transposable Element Expression Evolution Following Recent and Past Polyploidy Events in *Spartina* (Poaceae)

Delphine Giraud¹, Oscar Lima¹, Mathieu Rousseau-Gueutin², Armel Salmon¹ and Malika Aïnouche¹*

¹ UMR CNRS 6553 Ecosystèmes, Biodiversité, Evolution (ECOBIO), Université de Rennes 1, Rennes, France, ² IGEPP, INRAE, Institut Agro, Univ Rennes, Le Rheu, France

Gene expression dynamics is a key component of polyploid evolution, varying in nature, intensity, and temporal scales, most particularly in allopolyploids, where two or more sub-genomes from differentiated parental species and different repeat contents are merged. Here, we investigated transcriptome evolution at different evolutionary time scales among tetraploid, hexaploid, and neododecaploid *Spartina* species (Poaceae, Chloridoideae) that successively diverged in the last 6–10 my, at the origin of differential phenotypic and ecological traits. Of particular interest are the recent (19th century) hybridizations between the two hexaploids *Spartina alterniflora* (2n = 6x = 62) and S. maritima (2n = 6x = 60) that resulted in two sterile F1 hybrids: *Spartina × townsendii* (2n = 6x = 62) in England and *Spartina × neyrautii* (2n = 6x = 62) in France. Whole genome duplication of S. × townsendii gave rise to the invasive neo-allododecaploid species *Spartina anglica* (2n = 12x = 124). New transcriptome assemblies and annotations for tetraploids and the enrichment of previously published reference transcriptomes for hexaploids and the allododecaploid allowed identifying 42,423 clusters of orthologs and distinguishing 21 transcribed transposable element (TE) lineages across the seven investigated *Spartina* species. In 4x and 6x mesopolyploids, gene and TE expression changes were consistent with phylogenetic relationships and divergence, revealing weak expression differences in the tetraploid sister species *Spartina bakeri* and *Spartina versicolor* (<2 my divergence time) compared to marked transcriptome divergence between the hexaploids S. alterniflora and S. maritima that diverged 2–4 mya. Differentially expressed genes were involved in glycolysis, post-transcriptional protein modifications, epidermis development, biosynthesis of carotenoids. Most detected TE lineages (except *SINE* elements) were found more expressed in hexaploids than in tetraploids, in line with their abundance in the corresponding genomes. Comparatively, an astonishing (52%) expression repatterning and deviation from parental additivity were observed following recent reticulate evolution (involving the F1 hybrids and the neo-allododecaploid S. anglica), with various patterns of biased homoeologous gene expression, including genes involved in epigenetic regulation. Downregulation of TEs was observed in both

**OPEN ACCESS**

**Edited by:**
Yves Van de Peer,
Ghent University, Belgium

**Reviewed by:**
Pamela Soltis,
University of Florida, United States
Clayton J. Visger,
California State University,
Sacramento, United States

*Correspondence:
Malika Aïnouche
malika.ainouche@univ-rennes1.fr

**Specialty section:**
This article was submitted to
Plant Genomics,
a section of the journal
Frontiers in Genetics

**Received:** 30 July 2020
**Accepted:** 23 February 2021
**Published:** 25 March 2021

**Citation:**
Giraud D, Lima O, Rousseau-Gueutin M, Salmon A and Aïnouche M (2021) *Gene and Transposable Element Expression Evolution Following Recent and Past Polyploidy Events in *Spartina* (Poaceae).* Front. Genet. 12:589160. doi: 10.3389/fgene.2021.589160
hybrids and accentuated in the neo-allopolyploid. Our results reinforce the view that allopolyploidy represents springboards to new regulatory patterns, offering to worldwide invasive species, such as *S. anglica*, the opportunity to colonize stressful and fluctuating environments on saltmarshes.

**Keywords: allopolyploidy, hybridization, transcriptome evolution, transposable elements (TE), *Spartina***

INTRODUCTION

Gene expression dynamics is a key component of polyploid evolution, and has received considerable interest in the last decades, where studies on various polyploid systems revealed an important role of whole genome duplication (WGD) in modulating diverse and novel gene expression patterns (Chen, 2007; Flagel et al., 2008; Yoo et al., 2014). Comparative analyses indicate that the evolution of duplicated gene expression has a temporal dimension, and that immediate responses to polyploidy may not only last over long periods, but also contribute to the long-term processes of diploidization and fractionation (Wendel et al., 2018). One of the key parameters of this dynamics is the divergence level between the genomes being merged and duplicated during the polyploid speciation process (Buggs et al., 2012; Tayalé and Parisod, 2013). Neopolyploids are usually classified as autopolyploids, where homologous genomes (i.e., within species) are duplicated, or as allopolyploids involving duplication of more or less divergent (homoelogous) genomes reunited in the same nucleus following interspecific hybridization (Stebbins, 1971; reviewed in Doyle et al., 2008). Most autopolyploids so far explored exhibit moderate transcriptome or proteome alteration compared to their diploid progenitors (Albertin et al., 2005; Parisod et al., 2010b; del Pozo and Ramirez-Parra, 2014; Visger et al., 2019; Song et al., 2020). In contrast, allopolyploidy seems to be accompanied by profound parental expression repatterning in naturally formed neopolyploids (Hegarty et al., 2006; Chelaifa et al., 2010b; Buggs et al., 2011), experimentally resynthesized allopolyploids and/or their naturally established counterparts (Wendel, 2000; Chen and Ni, 2006; Ha et al., 2009; Soltis et al., 2015), such as in oilseed rape, wheat, cotton, or coffee (Akhunova et al., 2010; Higgins et al., 2012; Combes et al., 2013; Yoo et al., 2013; Wu et al., 2018). These dynamics reflect both the effects of the reunion of divergent regulatory networks (resulting from hybridization) and the effects of genetic redundancy (resulting from WGD), all these contributing to a "transcriptomic shock" (Hegarty et al., 2006), as a functional extension of what was earlier termed as "genomic shock" (McClintock, 1984), including its associated genetic and epigenetic regulatory processes (Hu and Wendel, 2019).

Duplicated genes may undergo partitioning of ancestral functions (subfunctionalization, Force et al., 1999) or benefit from mutations conferring new functionality (neofunctionalization, Ohno, 1970), which in turn affects the long-term retention of the duplicated copies (duplication-degenerescence-complementation model, Lynch and Force, 2000). The enhanced functional plasticity resulting from gene copy redundancy and allelic diversity facilitates the

emergence of novel variation, traits and phenotypes (Finigan et al., 2012), thus contributing to species adaptation and long-term diversification (Comai, 2005; Conant and Wolfe, 2008; Jackson and Chen, 2010; Madlung, 2013; Tank et al., 2015; Van de Peer et al., 2017).

Expression evolution in polyploids has been extensively explored by testing their deviation from an expected transcriptomic "parental additivity." This non-additive parental expression may be perceived by considering either the overall gene expression level (usually measured by comparisons with the average expression of both parental species, i.e., Mid-Parent Value; MPV) or by considering the relative contribution of each homoeologous copy to the total expression level (Grover et al., 2012). The high frequency of biases in the respective contributions of homeologs reported in various allopolyploids lead to the "genomic dominance" concept whereby one parental (homoeologous) subgenome expression state is exhibited in strong preference over the other parental expression state (Flagel and Wendel, 2009).

Recent studies provide accumulating evidence that the "dominant" subgenome is early established following hybridization and may be maintained over generations (Edger et al., 2017; Bird et al., 2018). This, very interestingly, permits to connect the short-term genome dominance phenomenon to the long-term fractionation process affecting polyploid genomes (Wendel et al., 2018), where biases in gene loss between duplicated homoeologous genomes result from the selection against loss of the most expressed copy (Schnable et al., 2011). The corollary is that following allopolyploidy, the "dominant" subgenome (the highest expressed) is more likely to be retained than the "recessive" subgenome (the lowest expressed) that becomes the most fractionated (Cheng et al., 2018). As expected, subgenome dominance is not observed in neo-autopolyploids (Garsmeur et al., 2014) suggesting a major role of the composition of the two subgenomes being merged in this process, notably in their content of transposable or regulatory elements.

Allopolyploidy merges and duplicates more or less differentiated genomes, including repetitive sequences that represent a dynamic component of plant genomes (Bennetzen and Wang, 2014). Transposable elements (TEs) affect in various ways the structure and expression of polyploid genomes. Several studies (reviewed in Vicient and Casacuberta, 2017) have reported transcriptional reactivation of some TE lineages in hybrids and allopolyploids or new insertions as found in tobacco, sunflower, wheat, or *Brachiaria* sp. (Ungerer et al., 2006; Petit et al., 2010; Yaakov and Kashkush, 2012; Santos et al., 2015). As TE expression is controlled by epigenetic regulation (e.g., DNA hypermethylation and siRNAs), new TE regulation
following polyploid speciation may result in altered neighbor gene expression including gene silencing and novel gene expression patterns (Kashkush et al., 2003; Hollister and Gaut, 2009; Parisod et al., 2010a; Zhang et al., 2015; Lerat et al., 2019). In allopolyploids, divergence between the parental genomes in terms of TE abundance and distribution is expected to increase subgenome dominance, as the TE-rich subgenome is expected to be targeted by epigenetic regulation and then less expressed than the TE-sparse subgenome (Bird et al., 2018; Bottani et al., 2018).

This complex interplay between parental genome legacy, the challenges faced by the needed coordination of divergent regulatory networks in the same nucleus and the connected processes taking place at the early and latest stages of polyploid species evolution raise critical questions about the rules and mechanisms shaping modern extant genomes and their predictability among polyploid lineages.

The *Spartina* genus (Poaceae, Chloridoideae) provides excellent opportunities to explore such processes at different evolutionary time scales, providing a phylogenetically well-understood framework of polyploidization that occurred in the last 10–12 my (Ainouche et al., 2012; Rousseau-Gueutin et al., 2015). Hybridization and polyploidy have recurrently shaped the history of this genus (also considered as a subsection in the *Sporobolus* genus; Peterson et al., 2014 but see Bortolus et al., 2019) and have resulted in modern species with various ploidy, ranging from tetraploids to dodecaploids, with a basic chromosome number of *x* = 10 (Marchant, 1968a; Ainouche et al., 2012). No diploid species was reported to date in this genus which most likely evolved from a tetraploid ancestral lineage (*2n* = 4*x* = 40) that diversified mainly in North and South America (and more recently in Europe) where they colonize coastal and inland marshes (Mobberley, 1956; Baumel et al., 2002a). In this study, we focused on seven species representing different ploidy, selected on the basis of their known phylogenetic history (Baumel et al., 2002a; Fortune et al., 2007; Rousseau-Gueutin et al., 2015). In tetraploids, two sister species were selected: *S. bakeri* and *Spartina patens*. *S. bakeri* (sand cordgrass) is native to the southeastern United States, where it grows on sandy beaches along the coast and in inland freshwater habitat in Florida and southern Georgia. It is morphologically highly similar to *S. patens* (saltmeadow grass), from which it can be distinguished by its unique vegetative habit (wanting rhizomes) and tolerance to freshwater (Mobberley, 1956). *S. patens* is distributed along the eastern United States to the Caribbean and northeast Mexico and has been introduced to the Mediterranean coast since the mid-19th century (Fabre, 1849) where it has been first considered as an endemic species and named *S. versicolor* Fabre. Recent molecular studies revealed that this latter actually represents an introduced, invasive genotype of *S. patens* (Prieto et al., 2011; Baumel et al., 2016). In the hexaploid clade, which diverged from the tetraploids 6–10 mya (Rousseau-Gueutin et al., 2015), we analyzed two species that diverged c.a. 2–4 mya: *S. alterniflora* (*2n* = 6*x* = 62) native to the Atlantic American coasts and *S. maritima* (*2n* = 6*x* = 60) native to the European Atlantic coast (Baumel et al., 2002a). During the 19th century, introductions

of *S. alterniflora* to Europe led to hybridization with the native species and had important ecological and evolutionary consequences (reviewed in Ainouche et al., 2009; Strong and Ayres, 2013). Hybridization with *S. maritima* (as male parent) resulted in two natural F1 hybrids: *S. × townsendii* (*2n* = 6*x* = 62; Groves and Groves, 1880) in Southern England and *S. × neyrautii* (*2n* = 6*x* = 62; Foucaud, 1897) in south-west France (Marchant, 1963; Baumel et al., 2003; Ainouche and Wendel, 2014). Chromosome doubling of *S. × townsendii* at the end of the 19th century resulted in the formation of the new fertile allododecaploid (*2n* = 12*x* = 124) *S. anglica* (Hubbard, 1965; Marchant, 1968a). This system is now considered as one of the textbook examples of recent allopolyploid speciation, which allows examining the immediate effects of hybridization and genome duplication in natural populations (Ainouche et al., 2004).

Polyploidy and hybridization have important adaptive consequences in *Spartina*. Hexaploid species evolved anatomical and physiological traits [including salt and xenobiotic tolerance increased by specific cells such as salt glands or trichomes, production of DMSP (Dimethylsulfoniopropionate), an osmoprotectant molecule; Céccoli et al., 2015; Paredes-Páliz et al., 2017; Robertson et al., 2017; Rousseau et al., 2017; Cavé-Radet et al., 2019b] providing the ability to colonize low marsh zones and tolerate several hours of tidal immersion (Maricle et al., 2009; Bedre et al., 2016). Invasive abilities appear as an immediate consequence of hybridization and/or genome duplication in *Spartina* (Strong and Ayres, 2013; Ainouche and Gray, 2016). The neopolyploid *S. anglica* has rapidly expanded in range in its native region (Europe) and has been deliberately introduced in several continents for land reclamation and estuary stabilization (e.g., China, Australia, and New Zealand), where several attempts of eradication are being conducted to control the species (Gray and Benham, 1990; Hammond and Cooper, 2002; Ainouche et al., 2009). This recently formed allopolyploid species exhibits higher ecological amplitude and higher stress tolerance than its parental species (Cavé-Radet et al., 2019b). The highly duplicated nature of its allododecaploid genome is thought to play an important role in the plasticity and the observed adaptive features. First investigations of global gene expression using microarrays on this system revealed various patterns of non-additive parental expression affecting genes involved in stress tolerance (Chelaifa et al., 2010b; Ferreira de Carvalho et al., 2017). Significant DNA methylation changes were also observed notably in regions flanking transposable elements, as evaluated from methylation-sensitive AFLP (Salmon et al., 2005) and methylation sensitive transposon display (Parisod et al., 2009). More recently, reference transcriptomes were built using Next Generation Sequencing for the hexaploid and allododecaploid species (Ferreira de Carvalho et al., 2013b; Boutte et al., 2016) and small-RNAs (Cavé-Radet et al., 2019a) and repetitive sequences (Giraud et al., 2020) from tetraploid and hexaploid *Spartina* genomes were annotated.

In this study, gene and TE expressions were analyzed at different evolutionary time scales (Figure 1). The two
<figure><img src="image_4.png" /><figcaption>**FIGURE 1** | Phylogenetic relationships among *Spartina* species (redrawn from Baumel et al., 2002a), and position of the analyzed tetraploid (*Spartina bakeri* and *Spartina versicolor*), hexaploid (*Spartina maritima*, *Spartina alterniflora*, *Spartina × townsendii*, *Spartina × neyrautii*) and dodecaploid (*Spartina anglica*) taxa analyzed in this study. Divergence times (my) represent upper limits estimated so far (Rousseau-Gueutin et al., 2015). Stars indicate whole genome duplication (WGD) events during *Spartina* evolution.</figcaption></figure>

sister tetraploid species and two hexaploid species were first compared to understand how gene and TE expressions evolved since their divergence 6–10 mya. Then, the evolution of gene and TE expression following recent hybridization and neopolyploidy was analyzed by comparing the species involved in the recent allopolyploid speciation, including the hexaploid parental species, their reciprocal F1 hybrids and the allododecaploid formed during the end of the 19th century. More specifically, we addressed the following questions: (i) how similar are the tetraploid and hexaploid transcriptomes? (ii) how additive are F1 hybrids/allopolyploid transcriptomes compared to parental transcriptomes? (iii) how has the hybridization vs. genome duplication affected genes and the expression of TEs? (iv) are the transcriptomes of two independently formed hybrids affected in the same way? (v) what are the functions affected by interspecific hybridizations and/or genome duplication? To answer these questions, reference transcriptomes were assembled for each investigated species (providing new reference transcriptomes for the tetraploid species and enriching previously assembled transcriptomes of hexaploids and hybrid/allopolyploid species Ferreira de Carvalho et al., 2013b; Boutte et al., 2016). Homologous (putative orthologous) transcripts were identified among species and the relative contribution of homeologs was assessed for a

subset of non-additively expressed genes in the hybrids and/or the allopolyploid.

# MATERIALS AND METHODS

## Plant Material, RNA Extraction, and Sequencing

Among the seven *Spartina* species analyzed in this study (Figure 1), two individuals per species were collected from their natural habitat along the French coasts for *S. maritima*, *S. alterniflora*, *S. × neyrautii*, *S. anglica*, and *S. versicolor*, in Hythe (England) for *S. × townsendii* and in Florida (United States) for *S. bakeri* (Table 1). Previous transcriptome analyses (Ferreira de Carvalho et al., 2017) indicated that low intraspecific expression variation was encountered among plants collected along a tidal gradient in *Spartina* populations, but in order to limit potential plastic responses among samples, all individuals were transplanted to the same natural site (Morbihan, France on June, 19, 2013) in order to maintain all plants in the same environmental conditions. After 1 month of acclimation (2013/07/18), newly formed leaves from each individual were harvested between 10–12 am and flash frozen in liquid nitrogen. Thereafter, their RNAs
<table><caption>**TABLE 1** | Studied polyploid *Spartina* species and summary of RNA-Seq Illumina datasets used.</caption><thead><tr><th>Species</th><th>Synonyms</th><th>Ploidy, number of chromosomes</th><th>Sample origin</th><th>Number of replicates</th><th>Total library depth (cumulative number of paired reads)</th></tr></thead><tbody><tr><td><b>S. maritima</b> Curtis</td><td><b>Sporobolus maritimus</b> (Curtis)<br>P.M.Peterson & Saarela, comb. nov.</td><td>6x = 60</td><td>Le Hezo (Morbihan, France)</td><td>4</td><td>146,438,783</td></tr><tr><td><b>S. alterniflora</b> Loisel.</td><td><b>Sporobolus alterniflorus</b> (Loisel.)<br>P.M.Peterson & Saarela, comb. nov.</td><td>6x = 62</td><td>Le Faou (Finistère, France)</td><td>4</td><td>179,792,458</td></tr><tr><td><b>S. x townsendii</b> H.<br>Groves & J. Groves</td><td><b>Sporobolus x townsendii</b> (H.<br>Groves & J. Groves) P.M.Peterson<br>& Saarela, comb. nov.</td><td>6x = 62</td><td>Hythe (Hampshire, United Kingdom)</td><td>4</td><td>171,773,818</td></tr><tr><td><b>S. x neyrautii</b><br>Foucaud</td><td>–</td><td>6x = 62</td><td>Hendaye (Pyrénées Atlantiques, France)</td><td>3</td><td>139,586,516</td></tr><tr><td><b>S. anglica</b> C. E. Hubb.</td><td><b>Sporobolus anglicus</b> (C.E.Hubb.)<br>P.M.Peterson & Saarela, comb. nov.</td><td>12x = 124</td><td>Le Hezo (Morbihan, France)</td><td>8</td><td>296,374,022</td></tr><tr><td><b>S. versicolor</b> Fabre</td><td><b>Sporobolus versicolor</b> (Fabre)<br>P.M.Peterson & Saarela, comb. nov.</td><td>4x = 40</td><td>Aresquiers (Hérault, France)</td><td>3</td><td>76,965,982</td></tr><tr><td><b>S. bakeri</b> Merr.</td><td><b>Sporobolus bakeri</b> (Merr.)<br>P.M.Peterson & Saarela, comb. nov.</td><td>4x = 40</td><td>Florida (United States)</td><td>2</td><td>85,839,501</td></tr></tbody></table>

were directly extracted according to Chelaifa et al. (2010a). Illumina libraries were prepared according to the Tru-Seq PCR-Free Protocol provided by Illumina. One Illumina library was prepared per individual (two individuals per species) and was sequenced on two independent lanes. The *S. anglica* libraries were sequenced twice more than the hexaploid species because of its larger genome (each individual was sequenced on 4 independent lanes; see **Table 1** and **Supplementary Table 1**). Sequencing was performed at the BGI (Beijing Genomics Institute) using the Hi-Seq 2000 technology (100 bp paired-ends reads; insert-size of 500 bp). Raw sequence data have been deposited at GenBank under the accession number PRJNA657699.

## Enrichment of Hexaploid Reference Transcriptomes and New Tetraploid Transcriptomes

Prior to differential gene expression analyses, the reference transcriptomes previously obtained by co-assembling pyrosequencing (Ferreira de Carvalho et al., 2013b) and Illumina data (Boutte et al., 2016) were enriched with the set of newly sequenced reads to build updated references for hexaploid and allododecaploid species (data used for reference transcriptome assemblies are available on GenBank, BioProject PRJNA338100). Considering the complex assembly process of transcriptomes in polyploid species, bioinformatics tools used were configured in order to obtain contigs representing consensus sequences of all expressed homeologous copies (or homeologs) from one given gene (Ferreira de Carvalho et al., 2013b). The expression level estimates for each contig was thus the sum of expressions of homeologs (and possibly, slightly divergent paralogs, resulting from recent individual gene duplication). Individual copies of a given gene were detected after remapping reads to these contigs, detecting SNPs (Single Nucleotide Polymorphisms) and phasing them to build

haplotypes (see section Parental Origin of Transgressively Expressed Orthologs).

First, mappings of new Illumina reads on reference transcriptomes were performed using Bowtie2 (parameters: local alignment; authorized mismatch; G, 52, 8; Langmead and Salzberg, 2012). Approximately, 60% of reads per library were aligned on respective transcriptoms. Non-mapped reads, corresponding to the expression of non-assembled transcripts, were then pooled by species and assembled following a method previously described (Boutte et al., 2016). Briefly, reads were assembled using the Trinity algorithm (v. 2.8.4; Grabherr et al., 2011) with a k-mer size of 25 bp and a minimum assembled contig length of 48 bp. To improve the assembly, all new contigs with a length higher than 100 bp were co-assembled with Newbler software (parameters: ml = 40 bp; mi = 90%; v.2.8 Margulies et al., 2005). Finally, a selfBLAST+ of all new contigs was performed using BLAST+ (v. 2.5.0; Altschul et al., 1990) to identify redundant sequences. Contigs included in other contigs were removed and contigs with a minimum overlapping of 50 bp, with an identity higher or equal to 90% were assembled using custom python scripts (Boutte et al., 2016). The newly assembled contigs were then added to previous reference transcriptomes (Boutte et al., 2016) for each hexaploid species (*S. maritima*, *S. alterniflora*, *S. × townsendii*, and *S. × neyrautii*) and the allopolyploid *S. anglica*. As no reference transcriptomes were available for the analyzed tetraploid species, Illumina libraries sequenced for expression analyses were used to assemble *de novo* reference transcriptomes for each tetraploid species. Reads were assembled using the method described above (Trinity assembly followed by Newbler co-assembly and custom scripts to remove redundancy). Quality of assemblies was tested using BUSCO (v. 4.1.4; Seppey et al., 2019). Transcriptome completeness was measured by comparing new assembled transcriptomes with pre-selected BUSCO genes specific to Poales species (lineage dataset: poales_odb10; creation date: 2020-08-05; number of species: 12; number of BUSCO genes: 4,896; e-value < 1e-06).
## Functional Annotations

Prior to contig annotations, potential contamination during RNA extraction and/or library sequencing was verified using BLASTn (Altschul et al., 1990) against the NCBI nucleotide database (accessed 08/12/2018). Such contigs displaying similarities with Prokaryote, Fungi, Protist or Metazoan sequences were removed from the analyses (*e*-value < 1e-6 and 90% of identity).

New assembled and reference transcriptomes were annotated as detailed in Ferreira de Carvalho et al. (2013b) and Boutte et al. (2016). Contigs were aligned via BLASTx analyses (one-way Best Blast Hits with an *e*-value < 1e-6) against protein databases from six species: *Arabidopsis thaliana* (genome version TAIR10), *Oryza sativa* (v. IRGSP-1.0), *Brachypodium distachyon* (v. 3.0), *Panicum hallii* (v. 3.1), *Sorghum bicolor* (v. 3), *Zea mays* (v. B73 RefGen v4). Gene Ontologies (GO terms) of each gene identified in our transcriptomes were retrieved from plant protein databases and included in our annotations. Similarities with known protein domains were also found using hmmer [a Hidden Markov Model (HMM) algorithm hmmer.org, v. 3.2] against the Pfam database (v. 32.0 Finn et al., 2014). Figures representing results from annotation analyses were drawn using the VennDiagram (v. 1.6.20; Chen and Boutros, 2011) and the UpSetR v. 1.4 R packages (Conway et al., 2017).

To complete annotations, contigs were also aligned against the *Spartina* transposable element databases built by Giraud et al. (2020) as well as against ribosomal DNA and chloroplast genomes previously assembled (Boutte et al., 2015; Rousseau-Gueutin et al., 2015) using BLASTn (*e*-value < 1e-6; Altschul et al., 1990) and Bowtie2 (parameters: local alignment; authorized mismatch; G, 52, 2; Langmead and Salzberg, 2012).

## Identification of Orthologous Transcripts/Regions Among Species

To compare gene expression among *Spartina* species, the OrthoVenn2 program (Xu et al., 2019) was used to identify orthologous contigs or (orthologs) in the different transcriptomes. This program provides clusters of orthologous sequences by sequence comparison analyses (BLASTp) and graph-based clustering among up to eight plant species simultaneously. In order to reduce the computational time, only annotated contigs were conserved for analyses. Prior to OrthoVenn2 analysis, contigs were translated into protein sequences (in the 6 reading frames) using the EMBOSS transeq tool (v. 6.6.0; Madeira et al., 2019). OrthoVenn2 was run with default parameters (*e*-value 1e-5; inflation value 1.5) to identify orthologs among the different transcriptomes. OrthoVenn2 allowed grouping orthologs found in the studied species into the same cluster. Within each cluster, regions or windows where orthologs showed sequence similarities were delimited in order to measure expression levels in each species from the same sequence length. To achieve this, all contigs or orthologs within the same cluster were aligned against each other using BLASTn (*e*-value < 1e-6; Altschul et al., 1990) and windows displaying a maximal number of aligned contigs and a maximal size

were selected using custom python scripts for differential expression analyses.

## Differential Expression Analyses

Gene expression levels were measured by mapping reads of each sequenced library on newly assembled reference transcriptomes using the Bowtie2 software with the following parameters: local alignment; authorized mismatch; G, 52, 8 corresponding to 90% of minimum identity (Langmead and Salzberg, 2012). For each cluster of orthologs, gene expression levels were given by the number of reads mapped on contigs and included in windows or regions previously selected.

Expression of TEs was also estimated according to the same method as used for gene expression. RNA-seq libraries were mapped using Bowtie2 (same settings as mentioned above) on their respective *Spartina* TE database previously built by Giraud et al. (2020). The expression of each TE family was measured by the total number of reads aligned on the full-length sequence of the elements.

Read counts for genes and TEs were calculated separately for each sequenced library representing biological or technical replicates. Data counts were extracted with samtools (v. 1.3.1; Li et al., 2009) and were normalized using the EDASeq Bioconductor R-package (v. 2.18.0; Risso et al., 2011) in order to correct sequencing bias within each library (GC-content and gene length) and between libraries (sequencing depth) (Bullard et al., 2010; Li et al., 2015). These biases were corrected according to the full-quantile normalization implemented in EDASeq. Normalization factor with DESeq2 was skipped because differences in sequencing depth were already taken into account with EDASeq. After normalization, data quality assessment and sample homogeneity were analyzed by clustering (Euclidean distances between samples) and principal component analysis (PCA).

Statistical analyses of gene and TE expression were performed with the DESeq2 package (Love et al., 2014). This package calculates log2 fold-change between species (from normalized count), estimates the dispersion and uses Wald tests (negative binomial linear model) to determine significant differential expression. Before DESeq2 analyses, clusters of orthologs weakly expressed in all studied species (mean expression of all replicates lower than 5 reads after normalization step) were removed because (i) they potentially represented noise from RNA sequencing and (ii) their low counts could affect the statistical tests (biased dispersions and variance in log2 fold change). Regarding the limited number of biological replicates per species we considered for each species both biological and technical replicates without distinction, in order to maximize the intraspecific variance. Considering the large number of tests, adjusted p-values were estimated according to Bonferroni correction to minimize false-positive rates. In our analyses, genes and TEs were considered as differentially expressed (DE) between two species when the adjusted p-value was lower than 0.01.

Several comparisons were performed: (i) expression levels between species presenting the same ploidy or of the same hybrid origin (*S. versicolor* vs. *S. bakeri*; *S. alterniflora* vs.
S. maritima; S. × townsendii vs. S. × neyrautii), (ii) expression levels between the two tetraploid species (S. versicolor vs. S. bakeri) and the two hexaploid species (S. alterniflora vs. S. maritima), (iii) expression levels between each F1 hybrid (S. × townsendii vs. S. × neyrautii) and the MPV (Mid Parent Value; expected under additive parental expression) and (iv) expression levels between S. × townsendii and S. anglica. The MPV was calculated by averaging the expression levels found in both parental species (mean of expression levels found in all parental replicates after normalization). Genes that deviate from parental additivity in hybrids and allopolyploid (statistically DE compared to the estimated MPV) were classified in different patterns: under parental dominance when the gene expression level found in the hybrids or allopolyploid was not statistically DE compared to the expression level found in the maternal or paternal species (maternal or paternal dominance, respectively); transgressively expressed when the gene expression level found in the hybrids or allopolyploid was statistically DE compared to the expression levels found in maternal and in paternal species (transgressively up-regulated or down-regulated).

## GO Enrichment Analysis

Gene ontology analyses were conducted to identify metabolic pathways or biological processes that were found differentially expressed following recent hybridization and genome doubling or between tetraploid and hexaploid species. GO term enrichment analyses were performed with the online platform agriGO (v. 2.0; Tian et al., 2017), considering the GO terms associated with A. thaliana (GO annotations available on TAIR10). The set of reference genes included only a subset of Arabidopsis genes for which homologous genes were found in *Spartina* transcriptomes. Biological processes over-represented in DE genes were determined using hypergeometric statistical tests using Singular Enrichment Analysis tool (Bonferroni multi-test adjustment, *p*-value < 0.01).

## Parental Origin of Transgressively Expressed Orthologs

We selected a set of non-additive, transgressively expressed orthologs in the F1 hybrid S. × townsendii and/or the allopolyploid S. anglica compared to their parental species (S. alterniflora and S. maritima) and identified the parental origin of the transcribed hybrid/allopolyploid haplotypes (homeologous copies) involved in such an expression pattern. The detection of haplotypes was performed as previously described for high ploidy-level species without diploid genome reference (Boutte et al., 2016) by detecting polymorphisms (in regions with minimum read-depth of 30; with minimum proportion of reads displaying the alternative character state of 10%) and phasing reads sharing at least two character states in sliding windows of 120 bp.

## RESULTS

### New Reference Transcriptomes and Functional Annotations

The *de novo* transcriptome assembly process allowed building of new reference transcriptomes for the two tetraploid species containing 103,101 contigs for S. versicolor and 107,319 contigs for S. bakeri. Non-aligned reads to previously assembled transcriptomes of the hexaploid and allododecaploid species (40% of reads) were separately assembled, which resulted in the enrichment of reference transcriptomes from 158,825 to 240,710 contigs according to species (Table 2). Consequently, the new reference transcriptomes used for expression analyses contained 217,800 and 206,348 contigs for the parental hexaploid species S. maritima and S. alterniflora, respectively, 281,153 and 266,906 contigs in the F1 hybrids S. × townsendii and S. × neyrautii, and 297,327 contigs in the allododecaploid S. anglica. Their GC-contents range from 45 to 49% and their N50 varies from 417 bp in the F1 hybrids and the allopolyploid to 724 bp in the tetraploid species, with an intermediate value in the hexaploid parental species (471–494 bp). Quality of the new

<table><caption>TABLE 2 | *Spartina* reference transcriptomes and functional annotations.</caption><thead><tr><th>Species</th><th>S. maritima</th><th>S. alterniflora</th><th>S. × townsendii</th><th>S. × neyrautii</th><th>S. anglica</th><th>S. versicolor</th><th>S. bakeri</th></tr></thead><tbody><tr><td>Number of contigs in reference transcriptomes (Boutte et al., 2016)</td><td>58,975</td><td>43,521</td><td>58,120</td><td>62,101</td><td>56,617</td><td>-</td><td>-</td></tr><tr><td>Number of new assembled contigs</td><td>158,825</td><td>162,827</td><td>223,033</td><td>204,805</td><td>240,710</td><td>103,101</td><td>107,319</td></tr><tr><td>Total number of contigs in reference transcriptome</td><td>217,800</td><td>206,348</td><td>281,153</td><td>266,906</td><td>297,327</td><td>103,101</td><td>107,319</td></tr><tr><td>Number of annotated contigs</td><td>76,916 (35%)</td><td>69,980 (34%)</td><td>93,645 (33%)</td><td>91,798 (34%)</td><td>94,130 (32%)</td><td>66,924 (65%)</td><td>62,827 (59%)</td></tr><tr><td>Number of unigenes (detected using the O. sativa genome)</td><td>22,186</td><td>22,133</td><td>23,756</td><td>23,139</td><td>23,544</td><td>19,328</td><td>19,779</td></tr><tr><td>GC-content</td><td>45.18%</td><td>45.25%</td><td>46.10%</td><td>45.43%</td><td>45.33%</td><td>49.30%</td><td>47.23%</td></tr><tr><td>N50 (bp)</td><td>471</td><td>494</td><td>418</td><td>417</td><td>417</td><td>700</td><td>724</td></tr><tr><td>Mean contig length (bp)</td><td>345</td><td>351</td><td>313</td><td>322</td><td>315</td><td>559</td><td>565</td></tr><tr><td>Median contig length (bp)</td><td>230</td><td>226</td><td>211</td><td>223</td><td>209</td><td>358</td><td>356</td></tr><tr><td>Minimum contig length (bp)</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td></tr><tr><td>Maximum contig length (bp)</td><td>10,313</td><td>10,839</td><td>9,423</td><td>14,705</td><td>8,833</td><td>17,112</td><td>16,705</td></tr></tbody></table>
reference transcriptomes was assessed using BUSCO with a panel of genes representing Poales species. Results are given in **Supplementary Figure 1**.

Functional annotations were assigned to 35 and 34% of contigs in *S. maritima* and *S. alterniflora*, respectively (Table 2). The vast majority of annotations was performed with the Blastx method (e.g., in *S. maritima*, Blastx allowed us to assign a function to 86% of annotated contigs and 48% with the HMM-Pfam method; **Supplementary Figure 2**). In comparison with parental species, the same proportions of contigs were also annotated in the F1 hybrids and the allopolyploid with 93,645 contigs (33%) in *S. × townsendii*, 91,798 contigs (34%) in *S. × neyrautii* and 94,130 contigs (32%) in *S. anglica*. Finally, within tetraploid transcriptomes, 65% of contigs in *S. versicolor* and 59% of contigs in *S. bakeri* were annotated (representing, respectively, 66,924 and 62,827 contigs). In order to reduce the computational time of analyses, we explored gene expression evolution for this set of annotated contigs, which represent a large part of the global gene expression in *Spartina* leaves (83% of mapped reads for each library **Supplementary Table 2**).

## Orthologous Contig Identification

Orthologous contigs among reference transcriptomes were identified and 62,830 clusters of orthologs were formed using OrthoVenn2. It is important to note that a few clusters may sometimes belong to a single gene if it could not be fully assembled. Among these clusters, 12,520 contained at least one ortholog of each of the seven *Spartina* species, 1,725 contained only orthologs of hexaploid species and their derived taxa (orthologs found in transcriptomes of the parental species, the

F1 hybrids and the allododecaploid) and 2,640 contained only orthologs of tetraploid species (Figure 2). In addition, 4,856 clusters contained exclusively orthologs identified in both F1 hybrids, and 2,399 orthologs of only *S. × townsendii* and *S. anglica*. Overall, the identification of orthologs with the OrthoVenn2 method allowed finding orthologs for more than 62% of annotated contigs [from 62% in *S. versicolor* (41,832 contigs) to 66% in *S. alterniflora* (50,245 contigs)]; **Supplementary Table 2**]. The composition of each cluster of orthologs (contig name, annotation) is detailed in **Supplementary Table 3**.

## Gene Expression

After data normalization and elimination of low counts, 42,423 clusters of orthologs were retained for statistical analyses (representing 19,417 unigenes detected using the *O. sativa* genome). Data quality was assessed via hierarchical clustering (Euclidean distances) and principal component analysis of counts from each library (**Supplementary Figure 3**). Given the absence of aberrant results, all libraries from each species were conserved as replicates for analyses.

## Expression Evolution in Tetraploid and Hexaploid Contexts

Comparisons between tetraploids showed that in leaves, 6,111 contigs (or orthologs) were significantly differentially expressed (DE) between *S. versicolor* and *S. bakeri*, representing 14.4% of the studied contigs (**Figure 3A**). Almost half of them were over-expressed in *S. bakeri* (3,178; 7.5%) the other half being over-expressed in *S. versicolor* (2,933; 6.9%). Orthologs more expressed in *S. bakeri* appear to be involved in cell development, modification and protein transport, whereas orthologs more expressed in *S. versicolor* are involved in epidermis development, response to biotic and abiotic stresses (salt, light, cold, and bacterium), gene silencing and post-transcriptional protein modifications (**Supplementary Table 4**).

More orthologs (15,278; 36.0%) were found DE between the hexaploid species *S. maritima* and *S. alterniflora* than between the tetraploids. Among them, 7,482 orthologs (48.9% of DE contigs) were over-expressed in *S. alterniflora* compared to *S. maritima* and were mainly involved in cell development, post-translational protein modifications (protein desumoylation, protein amino acid myristoylation), fatty acid and starch metabolic processes, response to salt and cold stresses and gene silencing (**Supplementary Table 4**). Conversely, 7,805 orthologs (51.1% of DE contigs) were over-expressed in *S. maritima* and were involved in chloroplast and epidermis development, protein catabolic process, post-translational protein modifications (amino acid phosphorylation and methylation), glycolysis process. When comparing tetraploids with hexaploids, the number of DE orthologs varies from 11,573 (27.3%) between *S. bakeri* and *S. alterniflora* to 14,911 (35.1%) between *S. versicolor* and *S. maritima* (**Figure 3A**). GO-term enrichment analyses showed that these DE orthologs have notably a key role in glycolysis process, epidermis development, pigment biosynthetic process and gene silencing (**Figure 3B and Supplementary Table 5**).

<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /></figure>

## Expression Evolution Following Recent Hybridization Events

In order to examine the transcriptomic effects of interspecific hybridization in the two independently formed F1 hybrids (with *S. alterniflora* as maternal parent and *S. maritima* as paternal parent), gene expression levels were compared

first with the *in silico* MPV (average expression between parental species) (Figure 4A). Then orthologs that were found differentially expressed between the two hybrids were examined to (i) identify the potentially affected functions and (ii) in which way these DE genes deviate from parental additivity.

<figure><img src="image_7.png" /></figure>
Comparisons of each hybrid with the MPV revealed that among 42,423 studied clusters of orthologs, 48.1% of clusters (20,390 orthologs) in S. × townsendii and 42.8% of clusters (18,170 orthologs) in S. × neyrautii deviated from parental additivity. For each hybrid, an equivalent number of DE orthologs was under-expressed compared to the MPV (10,218 and 8,871 orthologs, respectively in S. × townsendii and S. × neyrautii) and over-expressed compared to the MPV (10,172 and 9,299 orthologs, respectively in S. × townsendii and S. × neyrautii). Expression patterns in hybrids (Table 3) revealed that 70% of non-additive genes were expressed similarly to one of the parental species (parental expression dominance). In S. × townsendii, 6,908 orthologs (16.3%) mimicked the maternal parent (S. alterniflora) whereas 7,151 orthologs (16.9%) mimicked the paternal parent S. maritima. In S. × neyrautii, 6,421 orthologs (15.1%) mimicked the maternal parent S. alterniflora whereas 6,466 orthologs (15.2%) mimicked the paternal parent S. maritima. Other DE genes were found transgressively expressed in the hybrids compared to the MPV. Interestingly, in both hybrids a higher number of orthologs was more transgressively up-regulated (3,884 and 3,452 in S. × townsendii and S. × neyrautii, respectively) than transgressively down-regulated (2,447 and 1,831, respectively).

We found that among the genes deviating from parental additivity in hybrids (20,390 orthologs in S. × townsendii and 18,170 orthologs S. × neyrautii as indicated above), 13,971 of them were DE compared to MPV in both hybrids (Figure 4B). These genes were mainly involved in post-transcriptional protein modification processes, nucleotide metabolic processes and glycolysis (Figure 4B and Supplementary Table 6). Enriched GO-terms also belongs to gene silencing mechanisms (e.g., GO:0006342 chromatin silencing; GO:0016571 histone methylation; GO:0035196 production of miRNAs involved in gene silencing by miRNA). More specifically, 63% of DE

orthologs (8,848/13,971) shared a similar pattern in both hybrids: 6,322 orthologs were under the same parental dominance and 2,526 orthologs were transgressively expressed in S. × townsendii and S. × neyrautii (Table 4). Enriched GO-terms for genes under paternal dominance in S. × townsendii and S. × neyrautii were mainly involved in gene silencing processes whereas genes under maternal dominance were involved in glycolysis process, post-transcriptional protein modification and protein transport (Figure 5). Other genes involved in the glycolysis process also appeared transgressively expressed in both F1 hybrids.

When comparing the two hybrids, we identified 10,703 DE orthologs (25.2% of clusters) between S. × townsendii and S. × neyrautii (Figure 4A). A large majority of these 25.2% orthologs (10,037; 93.8%) exhibits different expression patterns in each hybrid with regard to parental additivity (Supplementary Table 7). Specifically, 3,321 orthologs (31% of DE orthologs) were under parental dominance in both hybrids but did not mimic the same parent: 1,639 orthologs appeared under maternal dominance in S. × townsendii but under paternal dominance in S. × neyrautii and conversely, 1,682 orthologs appeared under paternal dominance in S. × townsendii and under maternal dominance in S. × neyrautii. These orthologs were involved in trichome development and response to stresses (salt, light, and cold) (Supplementary Table 7). Other DE orthologs between hybrids were mainly additive in one hybrid and expressed transgressively in the other hybrid (4,034 orthologs corresponding to 38% of DE orthologs). In total, 2,711 orthologs (25.3% of DE contigs) were transgressively up-regulated in one hybrid and additive in the other and 1,323 orthologs (12.4% of DE orthologs) were transgressively down-regulated in one hybrid and additive in the other (Supplementary Table 7). GO-term enrichment analyses reveal that orthologs transgressively up-regulated in S. × townsendii and additive in S. × neyrautii (1,378 orthologs) were involved in biotic/abiotic stress responses (virus,

<table><caption>TABLE 3 | Gene expression patterns in the two F1 hybrids and the allopolyploid compared to the parental species S. alterniflora (the maternal parent) and S. maritima (the paternal parent).</caption><thead><tr><th rowspan="3"></th><th colspan="6">No DE from the MPV</th><th colspan="6">DE from MPV</th></tr><tr><th rowspan="2">Additivity</th><th colspan="2">Maternal expression dominance</th><th colspan="2">Paternal expression dominance</th><th colspan="2">Transgressive down-regulation</th><th colspan="2">Transgressive up-regulation</th></tr><tr><th>H-H-H</th><th>H-H-O</th><th>H-H-O</th><th>H-H-O</th><th>H-H-O</th><th>H-H-O</th><th>H-H-O</th><th>H-H-O</th><th>H-H-O</th></tr></thead><tbody><tr><td>S. × townsendii</td><td>21,129</td><td>422</td><td>482</td><td>3,028</td><td>3,880</td><td>3,260</td><td>3,891</td><td>1,842</td><td>317</td><td>288</td><td>3,404</td><td>256</td><td>224</td></tr><tr><td>F1 hybrid</td><td></td><td>22,033<br>(51.9%)</td><td></td><td>6,908<br>(16.3%)</td><td></td><td>7,151<br>(16.9%)</td><td></td><td>2,447<br>(5.8%)</td><td></td><td></td><td>3,884<br>(9.2%)</td><td></td><td></td></tr><tr><td>S. × neyrautii</td><td>23,488</td><td>340</td><td>425</td><td>2,860</td><td>3,561</td><td>2,987</td><td>3,479</td><td>1,323</td><td>273</td><td>235</td><td>3,113</td><td>199</td><td>140</td></tr><tr><td>F1 hybrid</td><td></td><td>24,253<br>(57.2%)</td><td></td><td>6,421<br>(15.1%)</td><td></td><td>6,466<br>(15.2%)</td><td></td><td>1,831<br>(4.3%)</td><td></td><td></td><td>3,452<br>(8.1%)</td><td></td><td></td></tr><tr><td>S. anglica</td><td>18,621</td><td>667</td><td>782</td><td>3,300</td><td>4,540</td><td>3,022</td><td>4,057</td><td>2,186</td><td>335</td><td>337</td><td>3,946</td><td>361</td><td>269</td></tr><tr><td>allopolyploid</td><td></td><td>20,070<br>(47.3%)</td><td></td><td>7,840<br>(18.5%)</td><td></td><td>7,079<br>(16.7%)</td><td></td><td>2,858<br>(6.7%)</td><td></td><td></td><td>4,576<br>(10.8%)</td><td></td><td></td></tr></tbody></table>

For each taxon, the first line indicates the number of contigs from each pattern, the second line the total number of contigs from each category and the third line indicates the percentage of contigs from each category.
bacterium, salt, and cadmium) and in gene silencing (production of small RNA involved in gene silencing by RNA). In comparison, orthologs transgressively up-regulated in S. × *neyrautii* and additive in S. × *townsendii* (1,333 orthologs) were involved in gravitropism, cellular development and also in responses to abiotic stresses (salt and cadmium).

## Expression Evolution Following Recent Allododecaploid Speciation

Gene expression in S. anglica was compared (i) with the MPV and (ii) following the genome duplication *per se* by comparing gene expression of S. anglica with S. × *townsendii*. Analyses showed that 22,353 orthologs (52.7% of clusters) deviated from parental additivity in the allopolyploid, which is more than what was found in S. × *townsendii* (20,390 DE orthologs, 48.1% from MPV) (Figure 4A). The comparison against MPV revealed that 68% of DE orthologs in S. anglica were also DE in S. × *townsendii* (Figure 4B). These genes were mainly involved in glycolysis and fatty acid process, post-transcriptional protein modification process, nucleotide metabolic process

and gene silencing (Figure 4B and Supplementary Table 6). As reported for the F1 hybrids, the majority of identified DE orthologs (67%) were under parental expression dominance (70% in S. × *townsendii*; Table 3). However, in S. anglica, the number of genes under maternal dominance was superior to genes under paternal dominance (7,840 orthologs expressed similarly as S. alterniflora against 7,079 expressed similarly as S. maritima). This is in accordance with species-species comparisons (Figure 6A), which showed a smaller number of DE orthologs between S. anglica and S. alterniflora (15,409; 36.3%) than between S. anglica and S. maritima (16,143; 38.1%). In addition, gene expression patterns revealed an increase of transgressively expressed orthologs in S. anglica (7, 434) compared to both F1 hybrids (6,331 in S. × *townsendii* and 5,283 in S. × *neyrautii*), with a large majority of them (62%) being transgressively up-regulated.

We found that 56.7% of the studied orthologs in S. anglica conserved the same pattern as the F1 hybrid after genome doubling (Table 5). Among these, 35.02% of the analyzed orthologs were additive in S. × *townsendii* and conserved the additivity in the allopolyploid S. anglica. Interestingly, in the F1 hybrid and in the allopolyploid, genes mimicking the expression found in the maternal parent S. alterniflora (3,074 orthologs; 7.2%) were mainly involved in trichome development (Table 5 and Figure 6B). Conversely, genes expressed as the paternal parent S. maritima (3,050 orthologs; 7.2%) were involved in gene silencing and down-regulation processes. Among genes expressed transgressively in S. × *townsendii* and S. anglica, those involved in epidermis development appear downregulated whereas those involved in starch metabolism appear upregulated.

When examining gene expression evolution following genome duplication *per se*, 14,219 orthologs (33.5%) were found DE between the F1 hybrid S. × *townsendii* and the allopolyploid S. anglica (Figure 6A). An equivalent number of genes were up-regulated in each species: 6,989 orthologs were over-expressed in S. × *townsendii* and 7,230 orthologs were over-expressed in S. anglica. Genes more expressed in S. × *townsendii* appeared involved in chloroplast development, sulfur biosynthetic process, glycolysis process and response to stress (water deprivation, bacterium, and fungus infection) whereas genes more expressed in S. anglica appeared involved in epidermis development, negative regulation/gene silencing, response to salt stress (Supplementary Table 8). Most of the expression changes affected (i) genes that were additively expressed in S. × *townsendii* and found under maternal dominance or transgressively up-regulated in S. anglica (5.0% and 4.9% of orthologs, respectively) or, (ii) genes that were under paternal dominance in S. × *townsendii* and found under maternal dominance in S. anglica (5.0%) (Table 5). Finally, more genes shifted toward maternal dominance in S. anglica (7,840 orthologs) than to paternal dominance (7,079 orthologs).

Detection of haplotypes was performed for a set of 58 orthologous contigs showing an expression level clearly different following hybridization and duplication (transgressively expressed in S. × *townsendii* and S. anglica) and with enriched

<table><caption>**TABLE 4** | Gene expression pattern evolution in the two independently formed hybrids S. × *townsendii* and S. × *neyrautii*.</caption><thead><tr><th>Expression pattern in S. × *townsendii*</th><th>Expression pattern in S. × *neyrautii*</th><th>Number of genes</th><th>% of genes</th></tr></thead><tbody><tr><td rowspan="5">Additive</td><td>Additive</td><td>17,834</td><td><strong>42.0%</strong></td></tr><tr><td>Maternal expression dominance</td><td>1,069</td><td>2.5%</td></tr><tr><td>Paternal expression dominance</td><td>984</td><td>3.3%</td></tr><tr><td>Transgressive down-regulation</td><td>638</td><td>1.5%</td></tr><tr><td>Transgressive up-regulation</td><td>1,508</td><td>3.6%</td></tr><tr><td rowspan="5">Maternal expression dominance</td><td>Additive</td><td>1,730</td><td>4.1%</td></tr><tr><td>Maternal expression dominance</td><td>3,045</td><td><strong>7.2%</strong></td></tr><tr><td>Paternal expression dominance</td><td>1,826</td><td>4.3%</td></tr><tr><td>Transgressive down-regulation</td><td>155</td><td>0.4%</td></tr><tr><td>Transgressive up-regulation</td><td>152</td><td>0.4%</td></tr><tr><td rowspan="5">Paternal expression dominance</td><td>Additive</td><td>1,733</td><td>4.1%</td></tr><tr><td>Maternal expression dominance</td><td>1,868</td><td>4.4%</td></tr><tr><td>Paternal expression dominance</td><td>3,277</td><td><strong>7.7%</strong></td></tr><tr><td>Transgressive down-regulation</td><td>163</td><td>0.4%</td></tr><tr><td>Transgressive up-regulation</td><td>110</td><td>0.3%</td></tr><tr><td rowspan="5">Transgressive down-regulation</td><td>Additive</td><td>1,150</td><td>2.7%</td></tr><tr><td>Maternal expression dominance</td><td>216</td><td>0.5%</td></tr><tr><td>Paternal expression dominance</td><td>209</td><td>0.5%</td></tr><tr><td>Transgressive down-regulation</td><td>858</td><td><strong>2.0%</strong></td></tr><tr><td>Transgressive up-regulation</td><td>14</td><td>0.1%</td></tr><tr><td rowspan="5">Transgressive up-regulation</td><td>Additive</td><td>1,806</td><td>4.3%</td></tr><tr><td>Maternal expression dominance</td><td>223</td><td>0.5%</td></tr><tr><td>Paternal expression dominance</td><td>170</td><td>0.4%</td></tr><tr><td>Transgressive down-regulation</td><td>17</td><td>0.1%</td></tr><tr><td>Transgressive up-regulation</td><td>1,668</td><td><strong>3.9%</strong></td></tr></tbody></table>

The percentage of genes that conserved the same expression pattern in both hybrids are in bold.
functions. Among these 58 contigs, eight were too weakly expressed in the parents, preventing the identification of parental haplotypes in hybrids/allopolyploids; five were displaying no polymorphism within and between species; 28 were displaying an additive pattern (in terms of presence/absence) of parental haplotypes in the hybrids and the allopolyploid (but expressed in a transgressive manner); five were non-additive in S. × *neyrautii*, six in S. × *townsendii*; one in both hybrids; and five were non-additive in both F1 hybrids and

the allopolyploid (three of them are outlined in **Table 6**). The gene encoding the subunit SSRP1 of the FACT complex (Facilitates Chromatin Transcription) was transgressively upregulated in S. × *townsendii* and S. *anglica*. For this gene, only the haplotype originating from S. *alterniflora* was detected in the F1 hybrid and the allopolyploid, together with haplotypes shared by the two parents (i.e., of “undetermined origin”). The two other genes involved in salt stress response, showed different contributions of the

**Maternal dominance (mimicked S. *alterniflora*)**

<figure><img src="image_8.png" /></figure>

**Paternal dominance (mimicked S. *maritima*)**

<figure><img src="image_9.png" /></figure>

**Transgressive up-regulation**

<figure><img src="image_10.png" /></figure>

**Transgressive down-regulation**

<figure><img src="image_11.png" /></figure>

**FIGURE 5** | Bar graph of the most enriched GO-terms (biological process) corresponding to genes found DE compared to the MPV and with the same expression patterns in both F1 hybrids S. × *townsendii* and S. × *neyrautii* (e.g., bar graph of maternal dominance represents the enriched functions of the 3,045 genes found under maternal dominance in S. × *townsendii* and in S. *neyrautii*). GO-terms were grouped according to their metabolic pathway affiliation. X-axis corresponds to the negative logarithm of adjusted *p*-value.

**A**
**S. alterniflora**
2n = 6x = 62
Maternal parent

**B**
**S. maritima**
2n = 6x = 60
Paternal parent

**S. × townsendii**
2n = 6x = 62
F1 hybrid

**Maternal dominance (mimicked S. *alterniflora*)**

<figure><img src="image_12.png" /></figure>

**Transgressive up-regulation**

<figure><img src="image_13.png" /></figure>

**Paternal dominance (mimicked S. *maritima*)**

<figure><img src="image_14.png" /></figure>

**Transgressive down-regulation**

<figure><img src="image_15.png" /></figure>

**FIGURE 6** | Transcriptomic changes following hybridization and whole genome duplication in *Spartina*. (A) Number of differentially expressed genes between species (pairwise comparisons). The numbers in grey circles represent the number of orthologs differentially expressed among the 42,423 analyzed clusters. The numbers of upregulated orthologs are represented near species names. (B) Bar graph of the most enriched GO-terms (biological process) for non-additively expressed orthologs in the F1 hybrid and in the allopolyploid (expressed under parental dominance or transgressive). GO-terms were grouped according to their metabolic pathway affiliation. The X-axis corresponds to the negative logarithm of adjusted *p*-value. The complete list of enriched GO terms is available in **Supplementary Tables 6, 7**.
<table><caption>**TABLE 5** | Evolution of gene expression patterns between the F1 hybrid S. × townsendii and the allopolyploid S. anglica.</caption><thead><tr><th>Expression pattern in S. × townsendii</th><th>Expression pattern in S. anglica</th><th>Number of genes</th><th>% of genes</th></tr></thead><tbody><tr><td rowspan="5">Additive</td><td>Additive</td><td>14,931</td><td>35.2%</td></tr><tr><td>Maternal expression dominance</td><td>2,140</td><td>5.0%</td></tr><tr><td>Paternal expression dominance</td><td>1,693</td><td>4.0%</td></tr><tr><td>Transgressive down-regulation</td><td>1,208</td><td>2.8%</td></tr><tr><td>Transgressive up-regulation</td><td>2,061</td><td>4.9%</td></tr><tr><td rowspan="5">Maternal expression dominance</td><td>Additive</td><td>1,384</td><td>3.3%</td></tr><tr><td>Maternal expression dominance</td><td>3,074</td><td>7.2%</td></tr><tr><td>Paternal expression dominance</td><td>1,896</td><td>4.5%</td></tr><tr><td>Transgressive down-regulation</td><td>268</td><td>0.6%</td></tr><tr><td>Transgressive up-regulation</td><td>286</td><td>0.7%</td></tr><tr><td rowspan="5">Paternal expression dominance</td><td>Additive</td><td>1,455</td><td>3.4%</td></tr><tr><td>Maternal expression dominance</td><td>2,109</td><td>5.0%</td></tr><tr><td>Paternal expression dominance</td><td>3,050</td><td>7.2%</td></tr><tr><td>Transgressive down-regulation</td><td>268</td><td>0.6%</td></tr><tr><td>Transgressive up-regulation</td><td>269</td><td>0.6%</td></tr><tr><td rowspan="5">Transgressive down-regulation</td><td>Additive</td><td>795</td><td>1.9%</td></tr><tr><td>Maternal expression dominance</td><td>284</td><td>0.7%</td></tr><tr><td>Paternal expression dominance</td><td>256</td><td>0.6%</td></tr><tr><td>Transgressive down-regulation</td><td>1,081</td><td>2.5%</td></tr><tr><td>Transgressive up-regulation</td><td>31</td><td>0.1%</td></tr><tr><td rowspan="5">Transgressive up-regulation</td><td>Additive</td><td>1,505</td><td>3.5%</td></tr><tr><td>Maternal expression dominance</td><td>233</td><td>0.5%</td></tr><tr><td>Paternal expression dominance</td><td>184</td><td>0.4%</td></tr><tr><td>Transgressive down-regulation</td><td>33</td><td>0.1%</td></tr><tr><td>Transgressive up-regulation</td><td>1,929</td><td>4.5%</td></tr></tbody></table>

The percentage of genes that conserved the same pattern of expression in both species are in bold.

parental copies in their transgressive over-expression in the F1 hybrid and in the allopolyploid. For the CLC-d gene encoding a chloride channel protein, all haplotypes retrieved in the parental species were detected in S. × townsendii and S. anglica. However, for the gene encoding the salt-stress inducible tonoplast aquaporin 2, only haplotypes inherited from S. maritima were detected in S. × townsendii, whereas in S. anglica, haplotypes from both parental species were detected.

## TE Expression

A strong transcriptional activity was detected for nine TE families in the analyzed *Spartina* species (Figure 7A). Most of them were Class I elements such as LTR-retrotransposons [2 Gypsy lineages (Tekay and Ogre) and four Copia lineages (SIRE, Ivana, Ikeros, and Ale), LINE and SINE]. Concerning Class II elements, only CACTA lineages were notably expressed. LINE elements appear as the most highly expressed TE (Figure 7A) with an average expression level two times higher than Gypsy Tekay (the second highest expressed TE) and four times higher than Copia lineages. Expression levels of other

<table><caption>**TABLE 6** | Origin of haplotypes identified for three genes transgressively up-regulated following hybridization and genome doubling in S. × townsendii and S. anglica.</caption><thead><tr><th>Haplotype origin</th><th>Haplotypes inherited from S. alterniflora</th><th>Haplotypes inherited from S. maritima</th><th>Undetermined parental origin (with no polymorphism)</th></tr></thead><tbody><tr><td colspan="4"><strong>Subunit SSRP1 of FACT complex</strong></td></tr><tr><td>Haplotypes identified in parental species</td><td>✓ (2)</td><td></td><td>✓ (1)</td></tr><tr><td>Origin of haplotypes identified in S. × townsendii (F1 hybrid)</td><td>✓ (2)</td><td></td><td></td></tr><tr><td>Origin of haplotypes identified in S. anglica (Allopolyploid)</td><td>✓ (2)</td><td></td><td></td></tr><tr><td colspan="4"><strong>Chloride channel protein CLC-d</strong></td></tr><tr><td>Haplotypes identified in parental species</td><td>✓ (6)</td><td></td><td>✓ (5)</td></tr><tr><td>Origin of haplotypes identified in S. × townsendii (F1 hybrid)</td><td>✓ (3)</td><td></td><td>✓ (3)</td></tr><tr><td>Origin of haplotypes identified in S. anglica (Allopolyploid)</td><td>✓ (4)</td><td></td><td>✓ (4)</td></tr><tr><td colspan="4"><strong>Salt-stress inducible tonoplast intrinsic protein 2</strong></td></tr><tr><td>Haplotypes identified in parental species</td><td>✓ (1)</td><td>✓ (2)</td><td>✓ (1)</td></tr><tr><td>Origin of haplotypes identified in S. × townsendii (F1 hybrid)</td><td>✓ (1)</td><td></td><td></td></tr><tr><td>Origin of haplotypes identified in S. anglica (Allopolyploid)</td><td>✓ (1)</td><td>✓ (2)</td><td></td></tr></tbody></table>

The number of haplotypes detected is indicated in brackets.

less expressed TEs were shown in Supplementary Figure 4. Among these, mutator elements are significantly more expressed in the tetraploid S. bakeri than in the other tetraploid or hexaploid species.

These nine TE families were transcriptionally active in all *Spartina* species but their expression levels varied according to species (Figure 7B). The two tetraploid species showed similar and moderate TE expression levels compared to the hexaploid species. Indeed, all *Copia* lineages, as well as the *Gypsy Tekay*, *LINE* and *CACTA*, were significantly less expressed in S. versicolor and S. bakeri than in S. maritima and/or in S. alterniflora. Only *SINE* elements were more expressed in both tetraploids than in hexaploids (twofold more expressed in tetraploids).

Expression comparisons between the parental hexaploid species S. maritima and S. alterniflora revealed the greatest number of differentially expressed TEs (six of the nine highest expressed TEs; Figure 7B) and the greatest variation in terms of expression levels. The *Gypsy Tekay*, *Ogre*, the *LINEs*
and the *CACTAs* elements were clearly more expressed in *S. maritima* than in *S. alterniflora* (twofold more expressed for *Gypsy Tekay*, *LINEs* and ninefold more for *CACTAs*) whereas the *Copia Ikeros* and *Ale* were more expressed in *S. alterniflora* than in *S. maritima* (eightfold and threefold more expressed, respectively).

In contrast, both F1 hybrids $S. \times townsendii$ and $S. \times neyrautii$ shared similar expression levels for all studied TEs. No significant differences were found except for the *Gypsy Tekay*, the *SINE* and *CACTA* elements (slightly more expressed in $S. \times townsendii$ than in $S. \times neyrautii$). When comparing TE expression levels of the F1 hybrids with those of their parental species, statistical analyses demonstrated that most TEs were non-additively expressed in the F1 hybrids. Four of the nine studied TEs (Figure 7B) exhibited parental dominance in both or at least one F1 hybrid. The *Gypsy Tekay*, *Copia SIRE*, and *LINE* elements mimicked the expression of the maternal parent *S. alterniflora* whereas only *Gypsy Ogre* mimicked the expression of the paternal parent *S. maritima*. TEs under parental dominance in F1 hybrids displayed expressions similar to the parent with the lowest expression, especially when the parental species exhibited large expression differences. The same pattern was also observed for *SINE* elements that showed a transgressive downregulated expression in both F1 hybrids (slightly less expressed in $S. \times townsendii$ and $S. \times neyrautii$ than in the parental species).

Finally, analysis of TE expression of the allopolyploid *S. anglica* revealed that the down-regulation patterns observed in $S. \times townsendii$ were reinforced in *S. anglica*. Indeed, *Gypsy Tekay*, *Ogre*, *Copia Ivana*, *Ikeros* as well as *SINE* and *CACTA* elements were significantly less expressed in *S. anglica* than in $S. \times townsendii$. Comparisons with expression levels detected in the parental species revealed that the expression dominance, biased toward the parent with the lowest expression level, was also strengthened in *S. anglica*: *Gypsy Tekay*, *Ogre* and *Copia Ivana* elements were found under maternal dominance in *S. anglica*, in addition to *Copia SIRE* and *LINEs* already found under the same dominance in $S. \times townsendii$.

## DISCUSSION

Transcriptomic changes affecting genes and transposable elements are major responses to hybridization and polyploidy. These changes that may occur immediately after (allo)polyploid speciation and persist over long-term of evolutionary time, as well as their consequences in terms of species adaptation and ecology, are timely but challenging questions. Taking advantage from the *Spartina* system where recurrent hybridization and polyploidy are well-documented in a well-understood phylogenetic context (Marchant, 1968b; Baumel et al., 2002a; Ainouche et al., 2009), we evaluated transcriptomic changes following past and recent polyploid speciation events by comparing gene and TE expressions among seven *Spartina* species with ploidy ranging from $4 \times$ to $12 \times$. New reference transcriptomes were

assembled for the tetraploid *S. bakeri* and *S. versicolor* (syn. *S. patens*), and the previously annotated transcriptomes were enriched for the hexaploid species (*S. alterniflora*, *S. maritima*, *S. × townsendii*, and *S. × neyrautii*) and the allododecaploid *S. anglica* (Ferreira de Carvalho et al., 2013b; Boutte et al., 2016). Ortholog detection and statistical analyses of gene and TE expression allowed comparing the expression levels of orthologous contigs included in 42,423 clusters and 21 TE lineages.

## Contrasted Gene Expression and TE Dynamics in Tetraploid and Hexaploid *Spartina* Subclades

The tetraploid lineage ($2n = 4x = 40$) is composed of eight *Spartina* species native to North or South America that so far have been poorly investigated regarding transcriptomic analyses (except the prairie cordgrass *Spartina pectinata*; Gedye et al., 2010; Nah et al., 2016). Here we analyzed two additional tetraploid species, *S. bakeri* and *S. versicolor*. These species show different ecological preferences (freshwater habitat and high marsh, respectively) but are weakly divergent morphologically (Mobberley, 1956) and genetically (Baumel et al., 2002a, 2016; Fortune et al., 2007; Peterson et al., 2014). Our results also show a weak transcriptome divergence as only a limited number of genes ($6,111$ orthologs; 14.4% of clusters) were found differentially expressed between these species (compared with other studied *Spartina*, see below). Genes differentially expressed were mainly involved in conserved metabolic processes such as cell development, transport and post-transcriptional protein modifications. However, in *S. versicolor* we observed a higher expression of genes involved in epidermis development. This result is consistent with epidermis differences identified by Maricle et al. (2009) who showed the presence of silica cells and a higher cuticle thickness in *S. versicolor* compared to *S. bakeri*. Specific analyses of genes involved in cuticle and silica cell development are necessary to confirm this causal relationship. In addition, *S. versicolor* displayed a higher expression of genes involved in abiotic stress including salt, light and cold stress response. *S. versicolor* is adapted to salt marsh conditions (Casolo et al., 2015) whereas *S. bakeri* occupies freshwater habitats. Higher tolerance to salt stress, as well as the presence of rhizomes in *S. versicolor* (absent in *S. bakeri*, Mobberley, 1956), may explain its invasiveness and its greater worldwide distribution compared to *S. bakeri*.

The majority of studied TE lineages exhibited similar expression levels in the two tetraploids. Previous investigations of repetitive sequences in the *S. bakeri* and *S. versicolor* genomes (Giraud et al., 2020) indicated that these species share the same quantity and diversity of TEs. Similar TE expression patterns thus indicate that no new significant transcriptional and/or transposition activity occurred following their divergence <2 mya, which occurred relatively recently in the *Spartina* clade (Rousseau-Gueutin et al., 2015). One exception can be made for Class II Mutator elements that appear to be ninefold more expressed in *S. bakeri* than in *S. versicolor* (Supplementary Figure 3), suggesting their recent transcriptional activation and
<figure><img src="image_16.png" /></figure>

insertion in the *S. bakeri* genome. Quantitative estimations of *Mutator* elements indicate that they are notably more abundant in *S. bakeri* (9.4 Mb) than in *S. versicolor* (1.5 Mb) genomes (Giraud *et al.*, 2020.), which would be consistent with their reported high rates of transposition (Dupeyron *et al.*, 2019).

In contrast to the two tetraploids investigated, the hexaploids *S. maritima* (2n = 6x = 60) and *S. alterniflora* (2n = 6x = 62) are morphologically (Mobberley, 1956) and genetically (Baumel et al., 2001, 2002a; Salmon et al., 2005; Parisod et al., 2009) well-differentiated species. They diverged about 2-4 mya along the European and American Atlantic coasts, respectively (Rousseau-Gueutin et al., 2015). About 1-5% nuclear nucleotide divergence was reported between these species (Fortune et al., 2007; Chelaifa et al., 2010a; Ferreira de Carvalho et al., 2013a,b; Boutte et al., 2016). Genome wide expression patterns are consistent with this divergence, as 15,278 differentially expressed orthologs (36.0% of studied clusters) were detected between *S. alterniflora* and *S. maritima*. Genes DE were involved in several conserved metabolic processes (cell development, protein catabolic process, glycolysis process, transport, and post-transcriptional protein modifications), but interestingly also in processes related to the physiology and ecology of *Spartina*. As previously observed using quantitative PCR analyses on target genes (Ferreira de Carvalho et al., 2017), our study confirmed that genes involved in response to salt stress were up-regulated in *S. alterniflora* compared to *S. maritima*. This result was directly in line with better tolerance to abiotic stress and invasiveness of

*S. alterniflora* highlighted in several studies that investigated its leaf morphology and anatomy (Maricle et al., 2009) and its tolerance to salt and hydrocarbon stress conditions (Watts et al., 2006; Bedre et al., 2016; Cavé-Radet et al., 2019b). These results somehow differ from previous transcriptomic comparisons between *S. maritima* and *S. alterniflora* performed using rice microarrays (Chelaifa et al., 2010a), which identified only 13.3% of DE genes (1,247 among the 9,353 examined genes in both species), most of them (957, belonging to developmental and cellular growth genes) being upregulated in *S. alterniflora* and downregulated in *S. maritima*. Various causes may explain the observed differences between our results and Chelaifa et al. (2010b) findings, such as: (i) rice microarray specificity, which tends to reduce the number of analyzed genes (70% of rice genes were hybridized with *Spartina* RNA on microarray), (ii) the potentially variable expression patterns in different conditions (i.e., plants maintained in the Greenhouse in Chelaifa's study and plants in natural conditions in this study) and different tissues (i.e., leaves and roots in Chelaifa's study and leaves in this study).

Divergent evolution between *S. maritima* and *S. alterniflora* appears to have also resulted in significant TE expression divergence. *Gypsy*, *Ogre*, *LINEs*, and *CACTAs* elements were significantly more expressed in *S. maritima* than in *S. alterniflora*, whereas *Copia Ikeros* and *Ale* elements were more expressed in *S. alterniflora* than in *S. maritima*. The expression levels seem positively correlated with their relative
abundance in both genomes. For example, *Copia Ikeros* elements, which represent 25 Mb of the *S. alterniflora* genome and 8.9 Mb of the *S. maritima* genome (Giraud et al., 2020), were eightfold more expressed in *S. alterniflora* than in *S. maritima*. These results clearly show that TE activity in *S. alterniflora* and *S. maritima* has differently evolved since their divergence 2–4 mya and probably led to some TE insertions in one or the other species according to TE lineage. Cavé-Radet et al. (2019a) identified 3,730 ra-siRNAs involved in the TE regulation in *S. maritima* and *S. alterniflora* preferentially targeting *Copia Ivana* and *SIRE*, *Gypsy Tekay* and *LINE* elements. These findings as well as our expression analyses indicate that in both hexaploids, these four active TEs are post-transcriptionally regulated (via small RNA synthesis) preventing their accumulation in the genome. However, highly expressed TEs such as *Copia Ikeros* and *Ale* in *S. alterniflora* and *CACTAs* in *S. maritima*, were not clearly under smallRNA control (Cavé-Radet et al., 2019a). This would suggest that both *Copia* elements benefit from *trans*-regulation of other expressed *Copia* elements (such as *Ivana* and *SIRE*). The limited post-transcriptional regulation of *CACTAs* despite their high transcriptional level may thus explain their accumulation in *S. maritima* (50.7 Mb vs. 26.0 Mb in *S. alterniflora*; Giraud et al., 2020).

All-to-all comparisons between tetraploid and hexaploid species revealed that the number of DE genes ranged from 27.3% (between *S. bakeri* and *S. alterniflora*) to 35.1% (between *S. versicolor* and *S. maritima*). GO-terms enrichment analyses showed that differences between tetraploid and hexaploid mainly concern genes involved in glycolysis, post-transcriptional protein modification, epidermis development. Interestingly, genes involved in the biosynthesis of carotenoids were less expressed in tetraploids than in hexaploids. Studies on plant responses to zinc and phenanthrene stresses in *Spartina densiflora* or to chromium stress in *Spartina argentinensis* (syn. *S. spartinae*; $2n = 4x = 40$) exhibit a decrease of carotenoids in stress condition (Mateos-Naranjo et al., 2008; Redondo-Gomez et al., 2011; Redondo-Gómez et al., 2011). Observed variations of genes involved in the carotenoid metabolic pathway between tetraploid and hexaploid species may be due to stress conditions of salt marsh (salt stress, flooding, and pollution) or linked to functions of carotenoids in plants (antioxidant during photosynthesis, precursors for the abscisic acid synthesis). To date, it remains unclear whether transcriptomic changes observed between tetraploid and hexaploid species appeared (i) following ploidy increase, (ii) more progressively during *Spartina* evolution or (iii) both. The auto- or allo-polyploid origin of the tetraploid and hexaploid *Spartina* lineages is not yet fully elucidated, and divergence of the duplicated ancestral genomes (including regulatory elements) and their subsequent evolution must have affected the transcriptome fates of the studied species. Nuclear gene phylogenies or haplotype detection from RNA-Seq data performed so far have revealed the presence of differentiated homeologs in both tetraploid and hexaploids (Fortune et al., 2007; Boutte et al., 2016; Ferreira de Carvalho et al., 2017) which would suggest a reticulate (i.e., allopolyploid) origin of these lineages, although

differentiating autopolyploidy (followed by duplicate gene divergence) from allopolyploidy in the meso-tetraploid species is a challenging task with no known related diploid species. Moreover, which ancestral tetraploid genome(s) contributed to the hexaploid ancestor of *S. maritima* and *S. alterniflora* remains an open question.

## Rapid Transcriptome Evolution Following Interspecific Hybridization: Alteration of Gene Expression and TE Silencing

The *Spartina* system offers unique opportunities to explore two components of the allopolyploid speciation process in natural conditions (hybridization versus genome duplication), a situation rarely met among recent and natural allopolyploid models (Ainouche and Wendel, 2014). Moreover, the first steps of the allopolyploid speciation process in natural conditions, i.e., consequences of divergent genome merger, can be explored in two independently formed hybrids $S. \times townsendii$ and $S. \times neyrautii$, which share the same maternal (*S. alterniflora*) and paternal (*S. maritima*) species (Baumel et al., 2003; Ainouche et al., 2004). These two hybrids exhibit different morphologies, $S. \times townsendii$ bearing intermediate traits between the parental species, and $S. \times neyrautii$ being highly similar to *S. alterniflora* (Marchant, 1977; Baumel et al., 2003). In spite of high pollen sterility, $S. \times townsendii$ still forms vigorous populations at the hybridizing site in England (Renny-Byfield et al., 2010; Huska et al., 2016), whereas only remnant sterile $S. \times neyrautii$ individuals are surviving in south-west France as a result of site disturbance and urbanization (Hubbard, 1968; Baumel et al., 2003). The different traits of these two hybrids sharing the same genetic origin have always been puzzling. Our analyses reveal that hybridization resulted in consistent transcriptomic changes, with slightly more genes deviating from parental expression additivity in $S. \times townsendii$ (48.1%) than in $S. \times neyrautii$ (42.8%). These results confirmed the gene expression alteration following hybridization previously found in $S. \times townsendii$ and $S. \times neyrautii$ by microarray analyses (Chelaifa, 2010; Chelaifa et al., 2010b). The majority of genes considered DE compared to MPV in $S. \times townsendii$ and $S. \times neyrautii$ were found DE in both hybrids indicating consistent effects of two independent hybridization events on gene regulation evolution. More specifically, in both $S. \times townsendii$ and $S. \times neyrautii$ genes under paternal dominance were mainly involved in gene silencing processes whereas genes under maternal dominance were involved in glycolysis process, post-transcriptional protein modification and protein transport. Our results also revealed that genes involved in negative regulation and chromatin silencing were overexpressed in *S. alterniflora* compared to *S. maritima*. This negative regulation mimicking the paternal parent in hybrids, not previously reported, is consistent with DNA methylation repatterning identified following hybridization using methylation sensitive amplified polymorphism analyses (Salmon et al., 2005; Parisod et al., 2009). Results thus suggest that epigenetic modifications, which appear rapidly in the newly formed F1 hybrids, led to gene repression or silencing in $S. \times townsendii$ and $S. \times neyrautii$. Similar transcriptomic
differences were observed in rice hybrids and were shown to result from the differential production of small interfering RNA (siRNA) in the parental lines. Thus, the F1 hybrids produced siRNA at an intermediate level compared to the two parents and if the level of a siRNA produced was sufficient to methylate both homeologs, this led to the hypermethylation of the locus from the unmethylated parent and its lower transcription (Chodavarapu et al., 2012), as potentially observed in our *Spartina* hybrids.

In contrast with previous microarray analyses (Chelaifa, 2010), no significant parental expression dominance was observed for most genes in both F1 hybrids, with an equivalent number of genes under paternal and maternal dominance (16.9 vs. 16.3% in S. × townsendii and 15.2% vs. 15.1% in S. × neyrautii). However, our results agree with Chelaifa (2010) regarding genes differentially expressed between the two hybrids where over-expression of genes involved in development and cellular growth in S. × neyrautii compared to S. × townsendii was reported. Our study also found that genes transgressively up-regulated in S. × neyrautii and additive in S. × townsendii compared to MPV were mainly involved in gravitropism, development and cellular organization. These results are consistent with the known morphological differences between these hybrids: S. × neyrautii plants exhibit high stems, long and fleshy leaves and long rhizomes whereas S. × townsendii usually display small height with reduced leaves (Marchant, 1977). In agreement with the microarray study, we also found that genes involved in salt stress response were up-regulated in S. × townsendii compared to S. × neyrautii. Comparison of salt stress tolerance between these hybrids was not performed to date but could help to better understand their physiology, ecology and contrasted distribution. Finally, genes involved in gene silencing also appear more expressed in S. × townsendii than in S. × neyrautii suggesting that hybridization induced a gene expression repression more important in S. × townsendii. Interestingly, when examining the non-additive patterns (compared to the MPV) of the genes DE between the two hybrids, we surprisingly found contrasted patterns in S. × townsendii and S. × neyrauti. Thus, our results provide new insights regarding the consequences of divergent genome merger, where non-additive gene expression in independent hybridization events (involving similar parental genotypes, Baumel et al., 2003) may entail differential parental gene expression repatterning, which very well illustrates the myriads of possible outcomes resulting from the “genomic shock” of hybridization (Nieto Feliner et al., 2020). The observed expression changes could also reflect post-hybridization or post-genome duplication evolution, although this might be limited regarding the generation time in these perennial young (c.a. 150 years old) plants. Hybridization is widespread and recurrent in natural populations, which increases the hybrid population genetic diversity. Recurrent and reciprocal hybridization between diploid *Tragopogon dubius* and *T. pratensis* resulted in the formation of morphologically diverse neo-allotetraploid *T. miscellus* individuals, which exhibit consistent transcriptomic differences, departure from parental additivity and differential homeologous expression bias (Shan et al., 2020). In our case, S. × neyrautii and S. × townsendii

share the same maternal parent (*S. alterniflora*, Baumel et al., 2003), in spite of morphological and transcriptomic differences, indicating different outcomes from the “replay of the evolution tape” (Gould, 1989).

Contrasted patterns are observed for repetitive sequences where TE transcriptional activity was similar in both F1 hybrids. Indeed, *Gypsy*, *Copia*, *LINEs*, *CACTAs* lineages expressed in the parental species were also highly expressed in S. × townsendii and S. × neyrautii. Both F1 hybrids exhibited similar levels and patterns of expression suggesting that both hybridization events induced similar consequences on TE activity. Comparisons with parental expression allowed classifying TEs into two main categories. On one hand, *Gypsy Ogre*, *Copia Ikeros* and *Ale* as well as *CACTA* elements were additively expressed in the F1 hybrids indicating no expression evolution following hybridization. On the other hand, *Gypsy Tekay*, *Copia SIRE*, *Ivana*, and *LINE* elements were repressed in both hybrids (mimicking the parental species displaying the lowest expression level). Small RNA analysis conducted by Cavé-Radet et al. (2019a) revealed that this second category of TEs was specially targeted by rasiRNAs in hexaploid species. Decrease of their transcriptional activity thus seems directly assignable to small RNA regulation inherited from parental species. In addition, Parisod et al. (2009) showed, using methyl-sensitive transposon display, that following hybridization in *Spartina*, DNA methylation increased drastically near TEs. Small RNAs can also act upstream to induce DNA methylation (RdDM pathways; Axtell, 2013; Borges and Martienssen, 2015; Wendte and Schmitz, 2018) and thus reinforce TE silencing. Consequently, all these results suggest that activity of several TEs in *Spartina* hybrids was constrained by various epigenetic regulations established rapidly following hybridization, causing TE silencing and preventing putative TE burst. This is consistent with the “genomic quiescence” with no transposition burst reported for a few targeted transposable elements in *Spartina* hybrids by Baumel et al. (2002b) and Parisod et al. (2009).

## Maternal Dominance and Strengthening of TE Silencing in the Neo-Allopolyploid

In the recently formed allododecaploid *S. anglica*, we analyzed the superimposed effects of hybridization and genome duplication on one hand, and the effect of genome duplication *per se* on the other hand, by comparing expression levels of the allododecaploid *S. anglica*, with those of the hexaploid parental species *S. alterniflora* and *S. maritima* and with those of *S. × townsendii*. About 33.5% contigs were DE between *S. × townsendii* and *S. anglica*. This result demonstrates that additional gene regulation changes affected gene expression after genome doubling *per se* as previously shown in *Spartina* (Chelaifa et al., 2010b), and in other allopolyploids such as *Senecio* (Hegarty et al., 2006), *Arabidopsis* (Wang et al., 2004; Madlung et al., 2005), *Triticum* (Shaked et al., 2001) or *Tragopogon* (Buggs et al., 2011; Boatwright et al., 2018). Among the main transcriptomic changes following genome duplication, it appeared that the number of genes mimicking the maternal expression pattern (*S. alterniflora*) increased in *S. anglica*
with 7,840 genes (18.5%) under maternal dominance. This is novel finding compared to previous analyses using microarrays (Chelaifa et al., 2010b) where the maternal dominance in S. × townsendii was found attenuated in S. anglica. Interestingly, genes involved in trichome development were under maternal dominance in S. × townsendii (mimicking S. alterniflora) and stayed under the same parental dominance in S. anglica. Several studies (e.g., in *Arabidopsis* and *Medicago sativa*; Alkio et al., 2005; Alves et al., 2017; Cavé-Radet et al., 2020) have shown the key role of trichomes in organic xenobiotics detoxification. In *Spartina*, tolerance to phenanthrene (PAH, polycyclic aromatic hydrocarbons) was increased following allopolyploidization (Cavé-Radet et al., 2019b). S. alterniflora was found more tolerant to phenanthrene than S. maritima but less tolerant than S. × townsendii and S. anglica. Then the trichome development regulation in the hybrid and allopolyploid, inherited from the maternal parent S. alterniflora may explain in part their better tolerance to PAH.

Conversely, in S. × townsendii and S. anglica, genes involved in silencing (chromatin silencing, histone methylation, and production of smallRNA) mimicked expression levels found in the paternal parent S. maritima. This result confirms that the transcriptomic repression observed following hybridization was inherited after genome doubling. This is consistent with a previous MSAP study that showed in S. anglica the inheritance of epigenetic marks, appearing following hybridization (Salmon et al., 2005; Parisod et al., 2009). Further comparisons between genes downregulated in S. × townsendii and S. anglica and miRNA-target genes identified in *Spartina* (Cavé-Radet et al., 2019a) will be interesting, to explore the link between miRNA production and decrease of such gene expression in hybrid and allopolyploid compared to parental species.

In addition to maternal dominance, an increase of transgressive genes was identified in S. anglica (7,434 genes in S. anglica vs. 6,331 in S. × townsendii), a phenomenon already observed by Chelaifa et al. (2010b). Genes transgressively expressed in the F1 hybrids and S. anglica (up or down-regulated compared to the parental lines) were involved in different biological processes, including epidermis development, starch metabolic process, post-transcriptional protein modifications. Among a selected set of 58 transgressively expressed genes, 28 of them showed that copies from both parental species contribute to the transgressive pattern. But in contrast, for 30 of them only one parental subgenome was involved in transgressive pattern.

Detecting homeologs in an allododecaploid species such as *Spartina anglica* is particularly challenging (Boutte et al., 2015, 2016), as the parental species are hexaploids, most likely of hybrid origins (Fortune et al., 2007) and as there is no known diploid *Spartina* that could be used as reference. The origins of the hexaploid lineage as well as the number of differentiated genomes that have been merged in the hexaploid ancestor are not fully elucidated yet. Notwithstanding this complexity and the absence of *Spartina* reference genome, we were able to take advantage from the divergence that occurred between the two hexaploid parents S. maritima and S. alterniflora in the

last 2–4 my to identify polymorphic orthologous regions and detect *maritima* versus *alterniflora* haplotypes, using parental reference transcriptomes, a procedure successfully employed in recent allopolyploids [e.g., *Capsella bursa-pastoris*, Kryvokhyzha et al. (2019); *Mimulus peregrinus*, Edger et al. (2017); *Tragopogon mirus* and *T. miscellus*, Boatwright et al. (2018)]. Further analyses on the *Spartina* genomes, aiming at exploring the nature and history of monoploid genomes in tetraploids and hexaploids, will allow distinguishing the meso-homeologs in the modern polyploids.

Non-additive patterns of parental expression contribute to enhance plasticity and adaptive responses to fluctuating environments in natural allopolyploid populations (Ferreira de Carvalho et al., 2017; Shimizu-Inatsugi et al., 2017). Our results also show an over-expression of genes involved in salt stress response in S. anglica compared to the F1 hybrid and the parental species. Several studies reported the intercontinental invasiveness of S. anglica populations that cope with severe chemical or physical constraints on salt marshes and high salinity levels (An et al., 2007; Strong and Ayres, 2013; Ainouche and Gray, 2016; Wong et al., 2018). This increased tolerance in S. anglica could be explained by the up-regulation of genes involved in salt stress response after genome doubling. For example, haplotype detection on two up-regulated genes in S. anglica (compared to S. × townsendii) and specially involved in salt stress response showed two different ways of haplotype expression evolution. For the CLC-d gene encoding the chloride channel protein (Wang et al., 2015), increase of expression was induced by the up-regulation of all or only a part of parental haplotypes retrieved in S. anglica. However, for the gene encoding the salt-stress inducible tonoplast aquaporin 2 (Wang et al., 2014), the restoration of “*alterniflora*-type” haplotype expression may account for its up-regulation.

In addition, we observed an over-expression of genes involved in silencing in S. anglica compared to S. × townsendii. For example, the FACT gene is known to encode a histone chaperone that can mediate nucleosome disassembly and reassembly (Grasser, 2020). Modification of the chromatin states via such histone chaperone was shown to mediate gene expression programs and help plants to more efficiently cope with stressful conditions (Probst and Mittelsten Scheid, 2015). Thus, it is tempting to speculate that such a gene may have played a key role during the allopolyploidization process, which merged two different genomes with divergent regulatory elements in the same nucleus. The identification of the parental origin of haplotypes suggests that the up-regulation of this gene after genome doubling was due to the increase of the expression of copies from S. alterniflora (maternal subgenome) and silencing of “*maritima*-type” haplotypes (already observed in S. × townsendii). The upregulation of genes involved in silencing in S. anglica suggested not only a parental legacy of repression as indicated above but also a strengthening of repression after genome duplication.

In this paper, we compared the relative expression per transcriptomes in the hexaploid parents and F1 hybrids and
the allododecaploid, by normalizing the number of reads relative to the total transcripts in each species. This assumes that mRNA transcriptome size (total number of transcripts per cell) is constant, which holds true when cells produce similar levels of RNA/cell (Lovén et al., 2012). Although widely used in comparative RNA-Seq studies, this approach does not take into account the potential absolute mRNA transcriptome size variation (reviewed Coate and Doyle, 2015) that may result from increased ploidy, which affects both gene copy number and cell size for some tissues as reported when comparing the allotetraploid *Glycine dolichocarpa* to its diploid progenitors (Coate and Doyle, 2010). These authors developed a quantitative PCR (qPCR) assay that normalizes individual gene expression to the genomic copy number when comparing the allotetraploid *G. dolichocarpa* to its diploid progenitors. Relative mRNA transcriptome sizes could be estimated by coupling this essay with transcriptome-normalized expression data (RNA-Seq) and revealed that the allotetraploid leaf transcriptome was approximately 1.4-fold larger than either diploid progenitor transcriptome. In *Spartina anglica*, gene dosage response is more complex as the parental genomes are hexaploid, and their genome composition as well as the extent of the fractionation process has to be elucidated as mentioned above. Future studies paralleling co-extracted genomic DNA and RNA Illumina sequencing should allow exploring this question. A larger sampling of additional plant organs and tissues, including root and inflorescences would be of particular interest for exploring previously reported differential tolerance to stresses (i.e., HAP and salt-stress) of investigated species and with regard to their fertility. As we are aware that transcriptome changes do not obviously correlate to the translatome (e.g., Coate et al., 2014; Soltis et al., 2016), another interesting question would be the impact of these transcriptional changes on the proteome of *S. anglica*, and its correlation to transcriptome responses.

Regarding repetitive sequences, TE repression, already observed in F1 hybrids, seems reinforced following genome doubling in *S. anglica*. Among the nine highest expressed TEs, six of them (*Gypsy Tekay*, *Ogre*, *Copia Ivana*, *Ikeros*, *SINES*, and *CACTAs*) were less expressed in *S. anglica* than in *S. × townsendii*. This indicates that additional regulation takes place following genome doubling to repress TE transcriptional activity. Studies on *Arabidopsis* or *Triticum* polyploids showed that allopolyploidization induced rapid methylation changes near TEs, avoiding TE burst (Madlung et al., 2005; Parisod et al., 2010a; Ben-David et al., 2013). Analyses in *Spartina* also revealed that DNA methylation that appeared near TEs following hybridization were conserved after genome doubling (Parisod et al., 2009). Therefore, in *S. anglica*, TE transcriptional activity was partly under the control of repressive epigenetic marks already implemented following hybridization and inherited from *S. × townsendii*.

In conclusion, our comparative transcriptomic analyses among *Spartina* species allowed us to understand the evolution of gene and TE expression following recent and

past polyploidization events. Gene expression changes were consistent with phylogenetic relationships and divergence time between species. Comparisons of tetraploid and hexaploid species showed that the TE dynamics was clearly different, reflecting a complex evolutionary history in both lineages since their divergence 6–10 mya. Particularly remarkable is the significant transcriptome repatterning following reticulate evolution, where expression changes (consistent with epigenetic and regulatory mechanisms alterations) that took place in 150–170 years old hybrids and neo-allopolyploid *S. anglica* far exceeded long term divergent transcriptome evolution in the meso-tetraploid and meso-hexaploid lineages. The superimposed polyploidization events which took place in the *Spartina* clade during the last 10–12 mya offered increased opportunity to partition parental expression. Recent allopolyploidy provided springboards for new regulatory and expression patterns that played a central role in the species traits and ecology, including abilities to colonize stressful and fluctuating environments on saltmarshes, as particularly well illustrated in the worldwide invasive allododecaploid *S. anglica*. The genomic and transcriptomic resources being developed on this system now open new perspectives to explore the deepest history of the parental species, the extent of fractionation affecting the ancestral tetraploid and hexaploid genomes, and the way this dynamic affects adaptation and invasiveness of the modern species.

## DATA AVAILABILITY STATEMENT

The datasets generated for this study are available at the NCBI SRA database on Bioproject PRJNA657699. Data used for reference transcriptome assemblies by Boutte et al. (2016) are available on BioProject PRJNA338100 (www.ncbi.nlm.nih.gov/sra/?term=PRJNA338100).

## AUTHOR CONTRIBUTIONS

MA, AS, and MR-G designed the experiments. OL contributed to molecular experiments for RNA seq. DG, AS, MA, and MR-G analyzed the results and wrote the manuscript. All authors contributed to the article and approved the submitted version.

## FUNDING

This study was supported by University of Rennes 1 (grant allocation to DG), European Union Seventh Framework Program: project 'Genomerge' [FP7-CIG-2013-2017; Grant No. 333709 to MR-G (transcriptome sequencing)].

## ACKNOWLEDGMENTS

The following institutions are thanked for their support: University of Rennes 1, UMR 6553 Ecosystèmes, Biodiversité, Evolution (Ecobio), Observatoire des Sciences de l'Univers (Rennes), INEE CNRS. This work benefited from the technical
platforms PEM and Ecolex (UMR Ecobio), Genouest. T. Fontaine and R. Bodiguel (UMR Ecobio) are thanked for their assistance in maintaining the *Spartina* plant collection. The reviewers are thanked for helpful comments and suggestions on the earliest version of this manuscript.

## REFERENCES

Ainouche, M., Chelaifa, H., Ferreira, J., Bellot, S., Ainouche, A., and Salmon, A. (2012). "Polyploid evolution in *Spartina*: dealing with highly redundant hybrid genomes," in *Polyploidy and Genome Evolution*, eds P. S. Soltis and D. E. Soltis (Berlin: Springer), 225-243.

Ainouche, M., and Gray, A. (2016). Invasive *Spartina*: lessons and challenges. *Biol. Invasions* 18, 2119-2122. doi: 10.1007/s10530-016-1201-7

Ainouche, M. L., Baumel, A., and Salmon, A. (2004). *Spartina anglica* C. E. Hubbard: a natural model system for analysing early evolutionary changes that affect allopolyploid genomes. *Biol. J. Linn. Soc.* 82, 475-484. doi: 10.1111/j.1095-8312.2004.00334.x

Ainouche, M. L., Fortune, P. M., Salmon, A., Parisod, C., Grandbastien, M.-A., Fukunaga, K., et al. (2009). Hybridization, polyploidy and invasion: lessons from *Spartina* (Poaceae). *Biol. Invasions* 11, 1159-1173. doi: 10.1007/s10530-008-9383-2

Ainouche, M. L., and Wendel, J. F. (2014). "Polyploid speciation and genome evolution: lessons from recent allopolyploids," in *Evolutionary Biology: Genome Evolution, Speciation, Coevolution and Origin of Life*, ed. P. Pontarotti (Cham: Springer International Publishing), 87-113.

Akhunova, A. R., Matniyazov, R. T., Liang, H., and Akhunov, E. D. (2010). Homoeolog-specific transcriptional bias in allopolyploid wheat. *BMC Genomics* 11:505. doi: 10.1186/1471-2164-11-505

Albertin, W., Brabant, P., Catrice, O., Eber, F., Jenczewski, E., Chèvre, A.-M., et al. (2005). Autopolyploidy in cabbage (*Brassica oleracea* L.) does not alter significantly the proteomes of green tissues. *Proteomics* 5, 2131-2139. doi: 10.1002/pmic.200401092

Alkio, M., Tabuchi, T. M., Wang, X., and Colón-Carmona, A. (2005). Stress responses to polycyclic aromatic hydrocarbons in *Arabidopsis* include growth inhibition and hypersensitive response-like symptoms. *J. Exp. Bot.* 56, 2983-2994. doi: 10.1093/jxb/eri295

Altschul, S. F., Gish, W., Miller, W., Myers, E. W., and Lipman, D. J. (1990). Basic local alignment search tool. *J. Mol. Biol.* 215, 403-410. doi: 10.1016/S0022-2836(05)80360-2

Alves, W. S., Manoel, E. A., Santos, N. S., Nunes, R. O., Domiciano, G. C., and Soares, M. R. (2017). Detection of polycyclic aromatic hydrocarbons (PAHs) in *Medicago sativa* L. by fluorescence microscopy. *Micron* 95, 23-30. doi: 10.1016/j.micron.2017.01.004

An, S. Q., Gu, B. H., Zhou, C. F., Wang, Z. S., Deng, Z. F., Zhi, Y. B., et al. (2007). *Spartina* invasion in China: implications for invasive species management and future research. *Weed Res.* 47, 183-191. doi: 10.1111/j.1365-3180.2007.00559.x

Axtell, M. J. (2013). Classification and comparison of small RNAs from plants. *Annu. Rev. Plant Biol.* 64, 137-159. doi: 10.1146/annurev-arplant-050312-120043

Baumel, A., Ainouche, M. L., Bayer, R. J., Ainouche, A. K., and Misset, M. T. (2002a). Molecular phylogeny of hybridizing species from the genus *Spartina* Schreb. (Poaceae). *Mol. Phylogenet. Evol.* 22, 303-314. doi: 10.1006/mpev.2001.1064

Baumel, A., Ainouche, M. L., Kalendar, R., and Schulman, A. H. (2002b). Retrotransposons and genomic stability in populations of the young allopolyploid species *Spartina anglica* C.E. Hubbard (Poaceae). *Mol. Biol. Evol.* 19, 1218-1227. doi: 10.1093/oxfordjournals.molbev.a004182

Baumel, A., Ainouche, M. L., and Levasseur, J. E. (2001). Molecular investigations in populations of *Spartina anglica* C.E. Hubbard (Poaceae) invading coastal Brittany (France). *Mol. Ecol.* 10, 1689-1701. doi: 10.1046/j.1365-294X.2001.01299.x

Baumel, A., Ainouche, M. L., Misset, M. T., Gourret, J.-P., and Bayer, R. J. (2003). Genetic evidence for hybridization between the native *Spartina maritima* and the introduced *Spartina alterniflora* (Poaceae) in South-West France: *Spartina*

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.589160/full#supplementary-material

x neyrautii re-examined. *Plant Syst. Evol.* 237, 87-97. doi: 10.1007/s00606-002-0251-8

Baumel, A., Rousseau-Gueutin, M., Sapienza-Bianchi, C., Gareil, A., Duong, N., Rousseau, H., et al. (2016). *Spartina versicolor* fabre: another case of *Spartina* trans-Atlantic introduction? *Biol. Invasions* 18, 2123-2135. doi: 10.1007/s10530-016-1128-z

Bedre, R., Mangu, V. R., Srivastava, S., Sanchez, L. E., and Baisakh, N. (2016). Transcriptome analysis of smooth cordgrass (*Spartina alterniflora* Loisel), a monocot halophyte, reveals candidate genes involved in its adaptation to salinity. *BMC Genomics* 17:657. doi: 10.1186/s12864-016-3017-3

Ben-David, S., Yaakov, B., and Kashkush, K. (2013). Genome-wide analysis of short interspersed nuclear elements SINES revealed high sequence conservation, gene association and retrotranspositional activity in wheat. *Plant J.* 76, 201-210. doi: 10.1111/tpj.12285

Bennetzen, J. L., and Wang, H. (2014). The contributions of transposable elements to the structure, function, and evolution of plant genomes. *Annu. Rev. Plant Biol.* 65, 505-530. doi: 10.1146/annurev-arplant-050213-035811

Bird, K. A., VanBuren, R., Puzey, J. R., and Edger, P. P. (2018). The causes and consequences of subgenome dominance in hybrids and recent polyploids. *New Phytol.* 220, 87-93. doi: 10.1111/nph.15256

Boatwright, J. L., McIntyre, L. M., Morse, A. M., Chen, S., Yoo, M.-J., Koh, J., et al. (2018). A robust methodology for assessing differential homeolog contributions to the transcriptomes of allopolyploids. *Genetics* 210, 883-894. doi: 10.1534/genetics.118.301564

Borges, F., and Martienssen, R. A. (2015). The expanding world of small RNAs in plants. *Nat. Rev. Mol. Cell Biol.* 16, 727-741. doi: 10.1038/nrm4085

Bortolus, A., Adam, P., Adams, J. B., Ainouche, M. L., Ayres, D., Bertness, M. D., et al. (2019). Supporting *Spartina*: interdisciplinary perspective shows *Spartina* as a distinct solid genus. *Ecology* 100:e02863. doi: 10.1002/ecy.2863

Bottani, S., Zabet, N. R., Wendel, J. F., and Veitia, R. A. (2018). Gene expression dominance in allopolyploids: hypotheses and models. *Trends Plant Sci.* 23, 393-402. doi: 10.1016/j.tplants.2018.01.002

Boutte, J., Aliaga, B., Lima, O., Ferreira, de Carvalho, J., Ainouche, A., et al. (2015). Haplotype detection from next-generation sequencing in high-ploidy-level species: 45S rDNA gene copies in the hexaploid *Spartina maritima*. *G3* 6, 29-40. doi: 10.1534/g3.115.023242

Boutte, J., Ferreira, de Carvalho, J., Rousseau-Gueutin, M., Poulain, J., Da Silva, C., et al. (2016). Reference transcriptomes and detection of duplicated copies in hexaploid and allododecaploid *Spartina* species (Poaceae). *Genome Biol. Evol.* 8, 3030-3044. doi: 10.1093/gbe/evw209

Buggs, R. J. A., Renny-Byfield, S., Chester, M., Jordon-Thaden, I. E., Viccini, L. F., Chamala, S., et al. (2012). Next-generation sequencing and genome evolution in allopolyploids. *Am. J. Bot.* 99, 372-382. doi: 10.3732/ajb.1100395

Buggs, R. J. A., Zhang, L., Miles, N., Tate, J. A., Gao, L., Wei, W., et al. (2011). Transcriptomic shock generates evolutionary novelty in a newly formed, natural allopolyploid plant. *Curr. Biol.* 21, 551-556. doi: 10.1016/j.cub.2011.02.016

Bullard, J. H., Purdom, E., Hansen, K. D., and Dudoit, S. (2010). Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments. *BMC Bioinform.* 11:94. doi: 10.1186/1471-2105-11-94

Casolo, V., Tomasella, M., De Col, V., Braidot, E., Savi, T., and Nardini, A. (2015). Water relations of an invasive halophyte (*Spartina patens*): osmoregulation and ionic effects on xylem hydraulics. * Funct. Plant Biol.* 42, 264-273. doi: 10.1071/FP14172

Cavé-Radet, A., Giraud, D., Lima, O., El Amrani, A., Ainouche, M., and Salmon, A. (2019a). Evolution of small RNA expression following hybridization and allopolyploidization: insights from *Spartina* species (Poaceae, Chloridoideae). *Plant Mol. Biol.* 102, 55-72. doi: 10.1007/s11103-019-00931-w
Cavé-Radet, A., Rabhi, M., Gouttefangeas, F., and El Amrani, A. (2020). Do specialized cells play a major role in organic xenobiotic detoxification in higher plants? *Front. Plant Sci.* 11:1037. doi: 10.3389/fpls.2020.01037

Cavé-Radet, A., Salmon, A., Lima, O., Ainouche, M. L., and El Amrani, A. (2019b). Increased tolerance to organic xenobiotics following recent allopolyploidy in *Spartina* (Poaceae). *Plant Sci.* 280, 143–154. doi: 10.1016/j.plantsci.2018.11.005

Céccoli, G., Ramos, J., Pilatti, V., Dellaferrera, I., Tivano, J. C., Taleisnik, E., et al. (2015). Salt glands in the Poaceae family and their relationship to salinity tolerance. *Bot. Rev.* 81, 162–178. doi: 10.1007/s12229-015-9153-7

Chelaifa, H. (2010). *Spéciation Allopolyploïde et Dynamique Fonctionnelle du Génome Chez les Spartines*. Dissertation thesis. Rennes: Université de Rennes 1.

Chelaifa, H., Mahé, F., and Ainouche, M. (2010a). Transcriptome divergence between the hexaploid salt-marsh sister species *Spartina maritima* and *Spartina alterniflora* (Poaceae). *Mol. Ecol.* 19, 2050–2063. doi: 10.1111/j.1365-294X.2010.04637.x

Chelaifa, H., Monnier, A., and Ainouche, M. (2010b). Transcriptomic changes following recent natural hybridization and allopolyploidy in the salt marsh species *Spartina x townsendii* and *Spartina anglica* (Poaceae). *New Phytol.* 186, 161–174. doi: 10.1111/j.1469-8137.2010.03179.x

Chen, H., and Boutros, P. C. (2011). VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. *BMC Bioinform.* 12:35. doi: 10.1186/1471-2105-12-35

Chen, Z. J. (2007). Genetic and epigenetic mechanisms for gene expression and phenotypic variation in plant polyploids. *Annu. Rev. Plant Biol.* 58, 377–406. doi: 10.1146/annurev.arplant.58.032806.103835

Chen, Z. J., and Ni, Z. (2006). Mechanisms of genomic rearrangements and gene expression changes in plant polyploids. *BioEssays* 28, 240–252. doi: 10.1002/bies.20374

Cheng, F., Wu, J., Cai, X., Liang, J., Freeling, M., and Wang, X. (2018). Gene retention, fractionation and subgenome differences in polyploid plants. *Nat. Plants* 4, 258–268. doi: 10.1038/s41477-018-0136-7

Chodavarapu, R. K., Feng, S., Ding, B., Simon, S. A., Lopez, D., Jia, Y., et al. (2012). Transcriptome and methylome interactions in rice hybrids. *Proc. Natl. Acad. Sci. U.S.A.* 109, 12040–12045. doi: 10.1073/pnas.1209297109

Coate, J. E., Bar, H., and Doyle, J. J. (2014). Extensive translational regulation of gene expression in an allopolyploid (*Glycine dolichocarpa*). *Plant Cell* 26, 136–150. doi: 10.1105/tpc.113.119966

Coate, J. E., and Doyle, J. J. (2010). Quantifying whole transcriptome size, a prerequisite for understanding transcriptome evolution across species: an example from a plant allopolyploid. *Genome Biol. Evol.* 2, 534–546. doi: 10.1093/gbe/evq038

Coate, J. E., and Doyle, J. J. (2015). Variation in transcriptome size: are we getting the message? *Chromosoma* 124, 27–43. doi: 10.1007/s00412-014-0496-3

Comai, L. (2005). The advantages and disadvantages of being polyploid. *Nat. Rev. Genet.* 6, 836–846. doi: 10.1038/nrg1711

Combes, M.-C., Dereeper, A., Severac, D., Bertrand, B., and Lashermes, P. (2013). Contribution of subgenomes to the transcriptome and their intertwined regulation in the allopolyploid *Coffea arabica* grown at contrasted temperatures. *New Phytol.* 200, 251–260. doi: 10.1111/nph.12371

Conant, G. C., and Wolfe, K. H. (2008). Turning a hobby into a job: how duplicated genes find new functions. *Nat. Rev. Genet.* 9, 938–950. doi: 10.1038/nrg2482

Conway, J. R., Lex, A., and Gehlenborg, N. (2017). UpSetR: an R package for the visualization of intersecting sets and their properties. *Bioinformatics* 33, 2938–2940. doi: 10.1093/bioinformatics/btx364

del Pozo, J. C., and Ramirez-Parra, E. (2014). Deciphering the molecular bases for drought tolerance in *Arabidopsis* autotetraploids. *Plant Cell Environ.* 37, 2722–2737. doi: 10.1111/pce.12344

Doyle, J. J., Flagel, L. E., Paterson, A. H., Rapp, R. A., Soltis, D. E., Soltis, P. S., et al. (2008). Evolutionary genetics of genome merger and doubling in plants. *Annu. Rev. Genet.* 42, 443–461. doi: 10.1146/annurev.genet.42.110807.091524

Dupeyron, M., Singh, K. S., Bass, C., and Hayward, A. (2019). Evolution of *Mutator* transposable elements across eukaryotic diversity. *Mob. DNA* 10:12. doi: 10.1186/s13100-019-0153-8

Edger, P. P., Smith, R., McKain, M. R., Cooley, A. M., Vallejo-Marin, M., Yuan, Y., et al. (2017). Subgenome dominance in an interspecific hybrid, synthetic allopolyploid, and a 140-year-old naturally established neo-allopolyploid monkeyflower. *Plant Cell* 29, 2150–2167. doi: 10.1105/tpc.17.00010

Fabre, M. E. (1849). Description d'une nouvelle espèce de *Spartina*, abondante sur une portion du littoral mééditerranéen. *Ann. Sci. Nat. Bot. Biol.* 3, 122–125.

Ferreira de Carvalho, J., Boutte, J., Bourdaud, P., Chelaifa, H., Ainouche, K., et al. (2017). Gene expression variation in natural populations of hexaploid and alldodecaploid *Spartina species* (Poaceae). *Plant Syst. Evol.* 303, 1061–1079. doi: 10.1007/s00606-017-1446-3

Ferreira de Carvalho, J., Chelaifa, H., Boutte, J., Poulain, J., Couloux, A., et al. (2013a). Exploring the genome of the salt-marsh *Spartina maritima* (Poaceae, Chloridoideae) through BAC end sequence analysis. *Plant Mol. Biol.* 83, 591–606. doi: 10.1007/s11103-013-0111-7

Ferreira de Carvalho, J., Poulain, J., Da Silva, C., Wincker, P., Michon-Coudouel, S., et al. (2013b). Transcriptome de novo assembly from next-generation sequencing and comparative analyses in the hexaploid salt marsh species *Spartina maritima* and *Spartina alterniflora* (Poaceae). *Heredity* 110, 181–193. doi: 10.1038/hdy.2012.76

Finigan, P., Tanurdzic, M., and Martienssen, R. A. (2012). "Origins of novel phenotypic variation in polyploids," in *Polyploidy and Genome Evolution*, eds P. S. Soltis and D. E. Soltis (Berlin: Springer), 57–76.

Finn, R. D., Bateman, A., Clements, J., Coggill, P., Eberhardt, R. Y., Eddy, S. R., et al. (2014). Pfam: the protein families database. *Nucleic Acids Res.* 42, 222–230. doi: 10.1093/nar/gkt1223

Flagel, L., Udall, J., Nettleton, D., and Wendel, J. (2008). Duplicate gene expression in allopolyploid *Gossypium* reveals two temporally distinct phases of expression evolution. *BMC Biol.* 6:16. doi: 10.1186/1741-7007-6-16

Flagel, L. E., and Wendel, J. F. (2009). Gene duplication and evolutionary novelty in plants. *New Phytol.* 183, 557–564. doi: 10.1111/j.1469-8137.2009.02923.x

Force, A., Lynch, M., Pickett, F. B., Amores, A., Yan, Y., and Postlethwait, J. (1999). Preservation of duplicate genes by complementary, degenerative mutations. *Genetics* 151, 1531–1545.

Fortune, P. M., Schierenbeck, K. A., Ainouche, A. K., Jacquemin, J., Wendel, J. F., and Ainouche, M. L. (2007). Evolutionary dynamics of Waxy and the origin of hexaploid *Spartina species* (Poaceae). *Mol. Phylogenet. Evol.* 43, 1040–1055. doi: 10.1016/j.mycopv.2006.11.018

Foucaud, J. (1897). *Un *Spartina* inédit*. *Ann. Soc. Sci. Nat. Char. Inf.* 32, 220–222.

Garsmeur, O., Schnable, J. C., Almeida, A., Jourda, C., D'Hont, A., and Freeling, M. (2014). Two evolutionarily distinct classes of paleopolyploidy. *Mol. Biol. Evol.* 31, 448–454. doi: 10.1093/molbev/mst230

Gedye, K., Gonzalez-Hernandez, J., Ban, Y., Ge, X., Thimmapuram, J., Sun, F., et al. (2010). Investigation of the transcriptome of prairie cord grass, a new cellulosic biomass crop. *Plant Genome J.* 3:69. doi: 10.3835/plantgenome2010.06.0012

Giraud, D., Lima, O., Huteau, V., Coriton, O., Boutte, J., Kovarik, A., et al. (2020). Evolutionary dynamics of transposable elements and satellite DNAs in polyploid *Spartina species*. *Plant Sci.* 302:110671. doi: 10.1016/j.plantsci.2020.110671

Gould, S. J. (1989). *Wonderful Life: The Burgess Shale and The Nature of History*. New York, NY: W. W. Norton & Company.

Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D. A., Amit, I., et al. (2011). Full-length transcriptome assembly from RNA-Seq data without a reference genome. *Nat. Biotechnol.* 29, 644–652. doi: 10.1038/nbt.1883

Gray, A. J., and Benham, P. E. M. (eds) (1990). *Spartina anglica: A Research Review*. London: HMSO.

Grasser, K. D. (2020). The FACT histone chaperone: tuning gene transcription in the chromatin context to modulate plant growth and development. *Front. Plant Sci.* 11:85. doi: 10.3389/fpls.2020.00085

Grover, C. E., Gallagher, J. P., Szadkowski, E. P., Yoo, M. J., Flagel, L. E., and Wendel, J. F. (2012). Homoeolog expression bias and expression level dominance in allopolyploids. *New Phytol.* 196, 966–971. doi: 10.1111/j.1469-8137.2012.04365.x

Groves, H., and Groves, J. (1880). *Spartina townsendii*. *Reports of the Botanical Society and Exchange Club of the British Isles* 1:37.

Ha, M., Kim, E.-D., and Chen, Z. J. (2009). Duplicate genes increase expression diversity in closely related species and allopolyploids. *Proc. Natl. Acad. Sci. U.S.A.* 106, 2295–2300. doi: 10.1073/pnas.0807350106

Hammond, M. E. R., and Cooper, A. (eds) (2002). *Turning the Tide: The Eradication of Invasive Species?: Proceedings of the International Conference on Eradication of Island Invasives*. Gland: IUCN.

Hegarty, M. J., Barker, G. L., Wilson, I. D., Abbott, R. J., Edwards, K. J., and Hiscock, S. J. (2006). Transcriptome shock after interspecific hybridization
in *Senecio* is ameliorated by genome duplication. *Curr. Biol.* 16, 1652-1659. doi: 10.1016/j.cub.2006.06.071

Higgins, J., Magusin, A., Trick, M., Fraser, F., and Bancroft, I. (2012). Use of mRNA-seq to discriminate contributions to the transcriptome from the constituent genomes of the polyploid crop species *Brassica napus*. *BMC Genomics* 13:247. doi: 10.1186/1471-2164-13-247

Hollister, J. D., and Gaut, B. S. (2009). Epigenetic silencing of transposable elements: a trade-off between reduced transposition and deleterious effects on neighboring gene expression. *Genome Res.* 19, 1419-1428. doi: 10.1101/gr.091678.109

Hu, G., and Wendel, J. F. (2019). *Cis* – *trans* controls and regulatory novelty accompanying allopolyploidization. *New Phytol.* 221, 1691-1700. doi: 10.1111/nph.15515

Hubbard, C. E. (1968). *Grasses*. London: Penguin Books.

Hubbard, J. C. E. (1965). *Spartina* marshes in southern England: VI. Pattern of invasion in Poole Harbour. *J. Ecol.* 53, 799-813. doi: 10.2307/2257637

Huska, D., Leitch, I. J., de Carvalho, J. F., Leitch, A. R., Salmon, A., Ainouche, M., et al. (2016). Persistence, dispersal and genetic evolution of recently formed *Spartina* homoploid hybrids and allopolyploids in Southern England. *Biol. Invasions* 18, 2137-2151. doi: 10.1007/s10530-015-0956-6

Jackson, S., and Chen, Z. J. (2010). Genomic and expression plasticity of polyploidy. *Curr. Opin. Plant Biol.* 13, 153-159. doi: 10.1016/j.pbi.2009.11.004

Kashkush, K., Feldman, M., and Levy, A. A. (2003). Transcriptional activation of retrotransposons alters the expression of adjacent genes in wheat. *Nat. Genet.* 33, 102-106. doi: 10.1038/ng1063

Kryvokhyzha, D., Milesi, P., Duan, T., Orsucci, M., Wright, S. I., Glémin, S., et al. (2019). Towards the new normal: transcriptomic convergence and genomic legacy of the two subgenomes of an allopolyploid weed (*Capsella bursa-pastoris*). *PLoS Genet.* 15:e1008131. doi: 10.1371/journal.pgen.1008131

Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. *Nat. Methods* 9, 357-359. doi: 10.1038/nmeth.1923

Lerat, E., Casacuberta, J., Chaparro, C., and Vieira, C. (2019). On the importance to acknowledge transposable elements in epigenomic analyses. *Genes* 10:258. doi: 10.3390/genes10040258

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and SAMtools. *Bioinformatics* 25, 2078-2079. doi: 10.1093/bioinformatics/btp352

Li, P., Piao, Y., Shon, H. S., and Ryu, K. H. (2015). Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data. *BMC Bioinform.* 16:347. doi: 10.1186/s12859-015-0778-7

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. *Genome Biol.* 15:550. doi: 10.1186/s13059-014-0550-8

Lovén, J., Orlando, D. A., Sigova, A. A., Lin, C. Y., Rahl, P. B., Burge, C. B., et al. (2012). Revisiting global gene expression analysis. *Cell* 151, 476-482. doi: 10.1016/j.cell.2012.10.012

Lynch, M., and Force, A. (2000). The probability of duplicate gene preservation by subfunctionalization. *Genetics* 154, 459-473.

Madeira, F., Park, Y. M., Lee, J., Buso, N., Gur, T., Madhusoodanan, N., et al. (2019). The EMBL-EBI search and sequence analysis tools APIs in 2019. *Nucleic Acids Res.* 47, W636-W641. doi: 10.1093/nar/gkz268

Madlung, A. (2013). Polyploidy and its effect on evolutionary success: old questions revisited with new tools. *Heredity* 110, 99-104. doi: 10.1038/hdy.2012.79

Madlung, A., Tyagi, A. P., Watson, B., Jiang, H., Kagochi, T., Doerge, R. W., et al. (2005). Genomic changes in synthetic *Arabidopsis* polyploids. *Plant J.* 41, 221-230. doi: 10.1111/j.1365-313X.2004.02297.x

Marchant, C. J. (1963). Corrected chromosome numbers for *Spartina x townsendii* and its parent species. *Nature* 199:929. doi: 10.1038/199929a0

Marchant, C. J. (1968a). Evolution in *Spartina* (Gramineae): II. Chromosomes, basic relationships and the problem of *S. x townsendii* agg. *J. Linn. Soc. Lond. Bot.* 60, 381-409. doi: 10.1111/j.1095-8339.1968.tb00096.x

Marchant, C. J. (1968b). Evolution in *Spartina* (Gramineae): III. Species chromosome numbers and their taxonomic significance. *J. Linn. Soc. Lond. Bot.* 60, 411-417. doi: 10.1111/j.1095-8339.1968.tb00097.x

Marchant, C. J. (1977). Hybrid characteristics in *Spartina x neyrautii* Fouc., a taxon rediscovered in northern Spain. *Bot. J. Linn. Soc.* 74, 289-296. doi: 10.1111/j.1095-8339.1977.tb01182.x

Margulies, M., Egholm, M., Altman, W. E., Attiya, S., Bader, J. S., Bemben, L. A., et al. (2005). Genome sequencing in microfabricated high-density picolitre reactors. *Nature* 437, 376-380. doi: 10.1038/nature03959

Maricle, B. R., Koteyeva, N. K., Voznesenskaya, E. V., Thomasson, J. R., and Edwards, G. E. (2009). Diversity in leaf anatomy, and stomatal distribution and conductance, between salt marsh and freshwater species in the C4 genus *Spartina* (Poaceae). *New Phytol.* 184, 216-233. doi: 10.1111/j.1469-8137.2009.02903.x

Mateos-Naranjo, E., Redondo-Gómez, S., Cambrollé, J., Luque, T., and Figueroa, M. E. (2008). Growth and photosynthetic responses to zinc stress of an invasive cordgrass. *Spartina densiflora. Plant Biol.* 10, 754-762. doi: 10.1111/j.1438-8677.2008.00098.x

McClintock, B. (1984). The significance of responses of the genome to challenge. *Science* 226, 792-801. doi: 10.1126/science.15739260

Mobberley, D. G. (1956). *Taxonomy and Distribution of the Genus *Spartina*. Retrospective Theses and Dissertations. Ames. IL: IOWA State University.

Nah, G., Lee, M., Kim, D.-S., Rayburn, A. L., Voigt, T., and Lee, D. K. (2016). Transcriptome analysis of *Spartina pectinata* in response to freezing stress. *PLoS One* 11:e0152294. doi: 10.1371/journal.pone.0152294

Nieto Feliner, G., Casacuberta, J., and Wendel, J. F. (2020). Genomics of evolutionary novelty in hybrids and polyploids. *Front. Genet.* 11:792. doi: 10.3389/fgene.2020.00792

Ohno, S. (1970). *Evolution by Gene Duplication*. Berlin: Springer.

Paredes-Páliz, K. I., Mateos-Naranjo, E., Doukkali, B., Caviedes, M. A., Redondo-Gómez, S., Rodríguez-Llorente, I. D., et al. (2017). Modulation of *Spartina densiflora* plant growth and metal accumulation upon selective inoculation treatments: a comparison of gram negative and gram positive rhizobacteria. *Mar. Pollut. Bull.* 125, 77-85. doi: 10.1016/j.marpolbul.2017.07.072

Parisod, C., Alix, K., Just, J., Petit, M., Sarilar, V., Mhiri, C., et al. (2010a). Impact of transposable elements on the organization and function of allopolyploid genomes. *New Phytol.* 186, 37-45. doi: 10.1111/j.1469-8137.2009.03096.x

Parisod, C., Holderegger, R., and Brochmann, C. (2010b). Evolutionary consequences of autoploidy. *New Phytol.* 186, 5-17. doi: 10.1111/j.1469-8137.2009.03142.x

Parisod, C., Salmon, A., Zerjal, T., Tenaillon, M., Grandbastien, M.-A., and Ainouche, M. (2009). Rapid structural and epigenetic reorganization near transposable elements in hybrid and allopolyploid genomes in *Spartina*. *New Phytol.* 184, 1003-1015. doi: 10.1111/j.1469-8137.2009.03029.x

Petersson, P. M., Romaschenko, K., Arrieta, Y. H., and Saarela, J. M. (2014). A molecular phylogeny and new subgeneric classification of *Sporobolus* (Poaceae: Chloridoideae: Sporobolinae). *Taxon* 63, 1212-1243. doi: 10.12705/636.19

Petit, M., Guidat, C., Daniel, J., Denis, E., Montoriol, E., Bui, Q. T., et al. (2010). Mobilization of retrotransposons in synthetic allotetraploid tobacco. *New Phytol.* 186, 135-147. doi: 10.1111/j.1469-8137.2009.03140.x

Prieto, J., Cires, E., Sánchez Corominas, T., and Vázquez, V. (2011). Systematics and management of natural resources: the case of *Spartina* species on European shores. *Biologia* 66:1011. doi: 10.2478/s11756-011-0109-z

Probst, A. V., and Mittelsten Scheid, O. (2015). Stress-induced structural changes in plant chromatin. *Curr. Opin. Plant Biol.* 27, 8-16. doi: 10.1016/j.pbi.2015.05.011

Redondo-Gomez, S., Andrades-Moreno, L., Parra, R., Valera-Burgos, J., Real, M., Mateos-Naranjo, E., et al. (2011). *Spartina densiflora* demonstrates high tolerance to phenanthrene in soil and reduces its concentration. *Mar. Pollut. Bull.* 62:9. doi: 10.1016/j.marpolbul.2011.05.018

Redondo-Gómez, S., Mateos-Naranjo, E., Vecino-Bueno, I., and Feldman, S. R. (2011). Accumulation and tolerance characteristics of chromium in a cordgrass Cr-hyperaccumulator, *Spartina argentinensis*. *J. Hazard. Mater.* 185, 862-869. doi: 10.1016/j.jhazmat.2010.09.101

Renny-Byfield, S., Ainouche, M., Leitch, I. J., Lim, K. Y., Le Comber, S. C., and Leitch, A. R. (2010). Flow cytometry and GISH reveal mixed ploidy populations and *Spartina* nonaploids with genomes of *S. alterniflora* and *S. maritima* origin. *Ann. Bot.* 105, 527-533. doi: 10.1093/aob/mcq008

Risso, D., Schwartz, K., Sherlock, G., and Dudoit, S. (2011). GC-content normalization for RNA-Seq data. *BMC Bioinform.* 12:480. doi: 10.1186/1471-2105-12-480

Robertson, M., Schrey, A., Shayter, A., Moss, C. J., and Richards, C. (2017). Genetic and epigenetic variation in *Spartina alterniflora* following the
deepwater horizon oil spill. *Evol. Appl.* 10, 792-801. doi: 10.1111/eva.
12482

Rousseau, H., Rousseau-Gueutin, M., Dauvergne, X., Boutte, J., Simon, G., Marnet, N., et al. (2017). Evolution of DMSP (dimethylsulfoniopropionate) biosynthesis pathway: origin and phylogenetic distribution in polyploid *Spartina* (Poaceae, Chloridoideae). *Mol. Phylogenet. Evol.* 114, 401-414. doi: 10.1016/j.ympev. 2017.07.003

Rousseau-Gueutin, M., Bellot, S., Martin, G. E., Boutte, J., Chelaifa, H., Lima, O., et al. (2015). The chloroplast genome of the hexaploid *Spartina maritima* (Poaceae, Chloridoideae): comparative analyses and molecular dating. *Mol. Phylogenet. Evol.* 93, 5-16. doi: 10.1016/j.ympev.2015.06.013

Salmon, A., Ainouche, M. L., and Wendel, J. F. (2005). Genetic and epigenetic consequences of recent hybridization and polyploidy in *Spartina* (Poaceae). *Mol. Ecol.* 14, 1163-1175. doi: 10.1111/j.1365-294X.2005. 02488.x

Santos, F. C., Guyot, R., do Valle, C. B., Chiari, L., Techio, V. H., Heslop-Harrison, P., et al. (2015). Chromosomal distribution and evolution of abundant retrotransposons in plants: gypsy elements in diploid and polyploid *Brachiaria* forage grasses. *Chromosome Res.* 23, 571-582. doi: 10.1007/s10577-015- 9492-6

Schnable, J. C., Springer, N. M., and Freeling, M. (2011). Differentiation of the maize subgenomes by genome dominance and both ancient and ongoing gene loss. *Proc. Natl. Acad. Sci. U.S.A.* 108, 4069-4074. doi: 10.1073/pnas. 1101368108

Seppey, M., Manni, M., and Zdobnov, E. M. (2019). "BUSCO: assessing genome assembly and annotation completeness," in *Gene Prediction Methods in Molecular Biology*, ed. M. Kollmar (New York, NY: Springer), 227-245.

Shaked, H., Kashkush, K., Ozkan, H., Feldman, M., and Levy, A. A. (2001). Sequence elimination and cytosine methylation are rapid and reproducible responses of the genome to wide hybridization and allopolyploidy in wheat. *Plant Cell* 13, 1749-1759.

Shan, S., Boatwright, J. L., Liu, X., Chanderbali, A. S., Fu, C., Soltis, P. S., et al. (2020). Transcriptome dynamics of the inflorescence in reciprocally formed allopolyploid *Tragopogon miscellus* (Asteraceae). *Front. Genet.* 11:888. doi: 10. 3389/fgenet.2020.00888

Shimizu-Inatsugi, R., Terada, A., Hirose, K., Kudoh, H., Sese, J., and Shimizu, K. K. (2017). Plant adaptive radiation mediated by polyploid plasticity in transcriptomes. *Mol. Ecol.* 26, 193-207. doi: 10.1111/mec.13738

Soltis, D. E., Misra, B. B., Shan, S., Chen, S., and Soltis, P. S. (2016). Polyploidy and the proteome. *Biochim. Biophys. Acta Proteins Proteomics* 1864, 896-907. doi: 10.1016/j.bkapap.2016.03.010

Soltis, P. S., Marchant, D. B., Van de Peer, Y., and Soltis, D. E. (2015). Polyploidy and genome evolution in plants. *Curr. Opin. Genet. Dev.* 35, 119-125. doi: 10.1016/j.gde.2015.11.003

Song, M. J., Potter, B. I., Doyle, J. J., and Coate, J. E. (2020). Gene balance predicts transcriptional responses immediately following ploidy change in *Arabidopsis thaliana*. *Plant Cell* 32, 1434-1448. doi: 10.1105/tpc.19.00832

Stebbins, G. L. (1971). *Chromosomal Evolution in Higher Plants*. London: Edward Arnold.

Strong, D. R., and Ayres, D. R. (2013). Ecological and evolutionary misadventures of *Spartina*. *Annu. Rev. Ecol. Evol. Syst.* 44, 389-410. doi: 10.1146/annurev- ecolsys-110512-135803

Tank, D. C., Eastman, J. M., Pennell, M. W., Soltis, P. S., Soltis, D. E., Hinchliff, C. E., et al. (2015). Nested radiations and the pulse of angiosperm diversification: increased diversification rates often follow whole genome duplications. *New Phytol.* 207, 454-467. doi: 10.1111/nph.13491

Tayalé, A., and Parisod, C. (2013). Natural pathways to polyploidy in plants and consequences for genome reorganization. *Cytogenet. Genome Res.* 140, 79-96. doi: 10.1159/000351318

Tian, T., Liu, Y., Yan, H., You, Q., Yi, X., Du, Z., et al. (2017). agriGO v2.0: a GO analysis toolkit for the agricultural community, 2017 update. *Nucleic Acids Res.* 45, 122-129. doi: 10.1093/nar/gkx382

Ungerer, M. C., Strakosh, S. C., and Zhen, Y. (2006). Genome expansion in three hybrid sunflower species is associated with retrotransposon proliferation. *Curr. Biol.* 16, 872-873. doi: 10.1016/j.cub.2006.09.020

Van de Peer, Y., Mizrachi, E., and Marchal, K. (2017). The evolutionary significance of polyploidy. *Nat. Rev. Genet.* 18, 411-424. doi: 10.1038/nrg.2017.26

Vicient, C. M., and Casacuberta, J. M. (2017). Impact of transposable elements on polyploid plant genomes. *Ann. Bot.* 120, 195-207. doi: 10.1093/aob/mcx078

Visger, C. J., Wong, G. K.-S., Zhang, Y., Soltis, P. S., and Soltis, D. E. (2019). Divergent gene expression levels between diploid and autotetraploid *Tolmiea* relative to the total transcriptome, the cell, and biomass. *Am. J. Bot.* 106, 280-291. doi: 10.1002/ajb2.1239

Wang, J., Tian, L., Madlung, A., Lee, H.-S., Chen, M., Lee, J. J., et al. (2004). Stochastic and epigenetic changes of gene expression in *Arabidopsis* polyploids. *Genetics* 167, 1961-1973. doi: 10.1534/genetics.104.027896

Wang, L.-L., Chen, A.-P., Zhong, N.-Q., Liu, N., Wu, X.-M., Wang, F., et al. (2014). The *Thellungiella salsuginea* tonoplast aquaporin TsTIP1:2 functions in protection against multiple abiotic stresses. *Plant Cell Physiol.* 55, 148-161. doi: 10.1093/pcp/pct166

Wang, S., Su, S. Z., Wu, Y., Li, S. P., Shan, X. H., Liu, H. K., et al. (2015). Overexpression of maize chloride channel gene ZmCLC-d in *Arabidopsis thaliana* improved its stress resistance. *Biol. Plant* 59, 55-64. doi: 10.1007/ s10535-014-0468-8

Watts, A. W., Ballestero, T. P., and Gardner, K. H. (2006). Uptake of polycyclic aromatic hydrocarbons (PAHs) in salt marsh plants *Spartina alterniflora* grown in contaminated sediments. *Chemosphere* 62, 1253-1260. doi: 10.1016/j. chemosphere.2005.07.006

Wendel, J. F. (2000). Genome evolution in polyploids. *Plant Mol. Biol.* 42, 225-249. doi: 10.1023/A:1006392424384

Wendel, J. F., Lisch, D., Hu, G., and Mason, A. S. (2018). The long and short of doubling down: polyploidy, epigenetics, and the temporal dynamics of genome fractionation. *Curr. Opin. Genet. Dev.* 49, 1-7. doi: 10.1016/j.gde.2018. 01.004

Wendte, J. M., and Schmitz, R. J. (2018). Specifications of targeting heterochromatin modifications in plants. *Mol. Plant* 11, 381-387. doi: 10.1016/j.molp.2017.10.002

Wong, J. X. W., Costantini, F., Merloni, N., Savelli, L., Geelen, D., and Airoldi, L. (2018). The widespread and overlooked replacement of *Spartina maritima* by non-indigenous *S. anglica* and *S. townsendii* in north-western Adriatic saltmarshes. *Biol. Invasions* 20, 1687-1702. doi: 10.1007/s10530-017-1654-3

Wu, J., Lin, L., Xu, M., Chen, P., Liu, D., Sun, Q., et al. (2018). Homoeolog expression bias and expression level dominance in resynthesized allopolyploid *Brassica napus*. *BMC Genomics* 19:586. doi: 10.1186/s12864-018-4966-5

Xu, L., Dong, Z., Fang, L., Luo, Y., Wei, Z., Guo, H., et al. (2019). OrthoVenn2: a web server for whole-genome comparison and annotation of orthologous clusters across multiple species. *Nucleic Acids Res.* 47, 52-58. doi: 10.1093/nar/ gkz333

Yaakov, B., and Kashkush, K. (2012). Mobilization of Stowaway-like MITEs in newly formed allohexaploid wheat species. *Plant Mol. Biol.* 80, 419-427. doi: 10.1007/s11103-012-9957-3

Yoo, M.-J., Liu, X., Pires, J. C., Soltis, P. S., and Soltis, D. E. (2014). Nonadditive gene expression in polyploids. *Annu. Rev. Genet.* 48, 485-517. doi: 10.1146/ annurev-genet-120213-092159

Yoo, M.-J., Szadkowski, E., and Wendel, J. F. (2013). Homoeolog expression bias and expression level dominance in allopolyploid cotton. *Heredity* 110, 171-180. doi: 10.1038/hdy.2012.94

Zhang, J., Liu, Y., Xia, E.-H., Yao, Q.-Y., Liu, X.-D., and Gao, L.-Z. (2015). Autotetraploid rice methylome analysis reveals methylation variation of transposable elements and their effects on gene expression. *Proc. Natl. Acad. Sci. U.S.A.* 112, 7022-7029. doi: 10.1073/pnas.15151 70112

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2021 Giraud, Lima, Rousseau-Gueutin. *Salmon and Ainouche*. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.