Robustness of Transposable Element Regulation but No Genomic Shock Observed in Interspecific Arabidopsis Hybrids

Robustness of Transposable Element Regulation but No Genomic Shock Observed in Interspecific... The merging of two divergent genomes in a hybrid is believed to trigger a “genomic shock”, disrupting gene regulation and transposable element (TE) silencing. Here, we tested this expectation by comparing the pattern of expression of transposable elements in their native and hybrid genomic context. For this, we sequenced the transcriptome of the Arabidopsis thaliana genotype Col-0, the A. lyrata genotype MN47 and their F1 hybrid. Contrary to expectations, we observe that the level of TE expression in the hybrid is strongly correlated to levels in the parental species. We detect that at most 1.1% of expressed transposable elements belonging to two specific subfamilies change their expression level upon hybridization. Most of these changes, however, are of small magnitude. We observe that the few hybrid-specific modifications in TE expression are more likely to occur when TE insertions are close to genes. In addition, changes in epigenetic histone marks H3K9me2 and H3K27me3 following hybridization do not coincide with TEs with changed expression. Finally, we further examined TE expression in parents and hybrids exposed to severe dehydration stress. Despite the major reorganization of gene and TE expression by stress, we observe that hybridization does not lead to increased disorganization of TE expression in the hybrid. Although our study did not examine TE transposition activity in hybrids, the examination of the transcriptome shows that TE expression is globally robust to hybridization. The term “genomic shock” is perhaps not appropriate to describe transcrip- tional modification in a viable hybrid merging divergent genomes. Key words: hybridization, genomic shock, transposable element, stress, Arabidopsis. Introduction Yet, a recent study of the homoploid hybrid of Arabidopsis Interspecific hybridization is an important factor in plant thaliana and Arabidopsis lyrata did not report a major disrup- evolution: While fertile allopolyploid hybrids may become tion of the transcriptome and epigenome. Moderate differ- the founders of new species (Soltis and Soltis 2009), the ences in expression of protein-coding genes were reported merging of two divergent genomes in a hybrid can also and the parental pattern of the histone mark H3K27me3 cause mild outbreeding depression or even result in com- appeared to be maintained (Zhu et al. 2017). These species, plete incompatibility (Bomblies and Weigel 2007; Todesco however, differ in both genomic structure and TE content. et al. 2016). At the molecular level, hybridization can lead A recent increase in TE number has been reported for the to a genome-wide misregulation of the transcriptome and A. lyrata genome (Hu et al. 2011). Arabidopsis thaliana TEs epigenome (Lafon-Placette and Ko¨hler 2015). The mobili- are concentrated in pericentromeric regions, rarely venturing zation of up-regulated transposable elements (TE) can fur- in gene rich chromosome arms. By contrast, in A. lyrata TEs ther result in chromosomal rearrangements. Such account for most of the 40% larger genome and are broadly phenomena, first described by Barbara McClintock in distributed, occurring in closer vicinity to expressed genes (Hu maize (McClintock 1984) and termed “genomic shock”, et al. 2011). The relatively low frequency of insertion poly- are thought to provide a postzygotic barrier against gene morphisms within species revealed evolutionary tensions on flow between species. insertions in gene rich regions, where permanent TE silencing The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non- commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com Genome Biol. Evol. 10(6):1403–1415. doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 1403 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Go ¨ bel et al. GBE can have deleterious consequences on the expression of 1-week intervals. Two micrograms of total RNA were used for neighboring genes (Hollister and Gaut 2009; Lockton and library preparation following the TruSeq Illumina RNA Sample Gaut 2010; Hollister et al. 2011). Finally, several studies Preparation v2 Guide. This includes poly-Aþ RNA selection have observed that A. lyrata TEs tended to be more highly and the use of random primers for reverse-transcription. expressed than A. thaliana alleles in hybrids (He et al. 2012). In Sequencing was performed on Illumina HiSeq2000 following view of these differences in number, distribution, and regula- the manufacturer’s protocols and paired-end 100-bp-long tion of transposable elements, we hypothesized that for reads were obtained. In total, 15–20 million paired reads for hybrids between these species, a genomic shock is more likely the parents and 30–40 million reads for the hybrid were pro- to occur at the level of TE expression. To test this hypothesis, duced (supplementary table 1, Supplementary Material on- we quantified the impact of parental versus hybrid genomic line). RNA-seq reads were filtered and trimmed for quality background on TE regulation. controlasin He et al. (2016) and mapped on the hybrid ge- We studied transposable elements in the A. thaliana (Col- nome, a concatenation of the A. thaliana Col-0 reference 0)  A. lyrata (MN47) hybrid and its respective parental lines, genome (TAIR10, www.arabidopsis.org; last accessed may using RNA-Seq expression data and ChIP-seq data of the his- 25, 2018) and the A. lyrata MN47 reference genome tone marks H3K27me3 and H3K9me2. It is difficult to estab- (Araly1, Nordberg et al. [2014]), using STAR with default lish orthology relationships for TEs, because many of them parameters (Dobin et al. 2013). duplicated after the species separation (Hu et al. 2011; de la Chaux et al. 2012). Thus, our analysis focuses on trans-acting effects experienced by TEs of each genome after hybridiza- tion. Contrary to our primary expectations, we find TE silenc- ing to be largely unaffected in the hybrid. No systematic TE and Gene Read Counting relationship could be found between a change of TE expres- sion in the hybrid and the gain or loss of histone marks. We The gene annotations TAIR10 (Col-0, Berardini et al. 2015), further exposed interspecific hybrids and their parents to se- Alyrata_107_v1.0 (MN47, Hu et al. 2011) and TE annotations vere dehydration stress and confirm that TE regulation (e.g., position and TE family memberships) for A. lyrata and A. remains robust to hybridization in stress conditions. thaliana (Pietzenuk et al. 2016) were merged to form the genome annotation file used in this analysis. Aligned reads were filtered as follows: for the MN47 genome, only the main Materials and Methods scaffolds 1–8 were considered. Read or fragment alignments Plant Materials and Growth Conditions shorter than 20 or longer than 700 base pairs (bp) were dis- carded. A minimum alignment score of 33 was required in Seeds of A. thaliana accession Col-0 were obtained from the order to accept an RNA-seq match. We allowed multiple Arabidopsis Biological Resource Center (ABRC, USA). RNA-seq read matches on TEs if they were entirely within a Arabidopsis lyrata ssp. lyrata genotype MN47 was obtained single (super)family, following Pietzenuk et al. (2016), and Jin from the lab of D. Weigel (Max Planck Institute for et al. (2015). A read was assigned to a (super)family F on Development, Tu ¨ bingen, Germany). Parental lines were genome G if its primary match was in F on G, and any sec- crossed by pollinating emasculated A. thaliana flowers with ondary matches on G were in F, too. In order to assign a A. lyrata pollen, as described in De Meaux et al. (2006).We match to a (super)family F, a minimum overlap of 50 bp to obtained approximately 40 seeds per silique, 90% of which a member of F was required. Secondary matches on the al- aborted. We did not use embryo rescue to generate more ternative genome were allowed but not counted. Each read hybrids as in Zhu et al. (2017). Reciprocal crosses using thus either contributed a single count to a (super)family or A. thaliana as pollen donor were not successful. was discarded as multiply matching. If a read matched to TE .. TE of family F, the read was defined to contribute 1 n RNA Sampling and Sequencing in Standard and Stress 1/n counts to each of these TEs. Reads from a region of over- Conditions lap a TE and a gene were not discarded, but counted with the Seeds were stratified for 3 days at 4 C, germinated on soil and TE. These summed fractional per-element counts were the grown for 4 weeks in a growth chamber at 20 C under 14 h final output from the counting procedure. They were either 1 2 light/16 C 10 h dark under dim light (100 mmol s m ). A used directly as counts per individual TE or the counts of the dehydration treatment was applied following Seki et al. members of (super)families were summed up to yield aggre- (2002). Plants were removed from the soil and dehydrated gated per-(super)family counts. We implemented the count- on paper for 0 and 3h, in the same growth chamber. The ing procedure in the R programming language. After the aerial part of the plant was sampled, flash-frozen in liquid counting step, annotated TEs without any read count were nitrogen, and RNA was extracted as previously described excluded from the analysis. Genes were counted with the (He et al. 2016). Four biological replicates were collected in same procedure as the TEs. 1404 Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Robustness of Transposable Element Regulation GBE Read Count Normalization each bin of the reference sample, the same number of TEs as contained in the bin was randomly selected from the corre- Read counts from genes and TEs were simultaneously nor- sponding distance bin of all TEs. The union of the resampled malized. Each count was first multiplied by a factor of 1,000/l, bins then formed one resampled set. The width of the bins where l is the length in base pairs of the respective annotated increased with increasing distance: up to a distance of 1 kb, element (gene, individual TE, or entire TE (super)family). This the bin width was 200 base pairs (bp), between 1 and 5 kb it removes the dependence of the read count on element length was 500 bp, and beyond 5 kb it was 1,000 bp, to avoid having or (super)family size. The total length of a TE (super)family was bins with sparse numbers of TEs. A median, 5% and 95% defined to equal the summed length of all (super)family mem- quantile of each superfamily frequency was extracted from bers considered, discounting any regions of overlap with the 10,000 random samples and compared with the observed other annotated elements. frequency of TEs with hybrid-dependent expression. Differential Expression Analysis The Generalized Linear Models (GLM) support of the R pack- Chromatin Immune-Precipitation Experiments age edgeR (Robinson et al. 2010; McCarthy et al. 2012)was used to test the significance of differences in length- Plants were germinated and grown on germination medium normalized counts between the four biological replicates of containing = Murashige and Skoog salts, 3% sucrose, and hybrid and parental transcriptomes. P values were adjusted 0.8% agar. Seeds were stratified for 5 days at 4 C, and then with the FDR correction (Benjamini and Hochberg 1995). For transferred to a growth chamber under long-day conditions visualization in scatterplots and for the stringent definition of (14 h cool white light supplemented with red light at 20 C, differential TE expression, length-normalized counts were and 10 h dark at 18 C). One week after germination, plants convertedtoCounts per Million(CPM) andscaled were transferred to soil and grown in the same chamber. (Robinson and Oshlack 2010) using the calcNormFactors Finally, leaves larger than 0.5 cm were sampled from function of the R package edgeR. This transformation mirrors 3-week-old Col-0, and 6-week-old MN47 and F1 hybrids what edgeR does internally when computing differential ex- for chromatin immunoprecipitation (ChIP). Leaves from 10 pression. It removes global count biases due to library size to 15 plants were pooled. Sampling was shifted in the hybrid differences, easing the interpretation of the scatterplots and and the A. lyrata individual, in order to sample material of allowing the set-up of a sample-independent cutoff on the similar development (3 cm rosette diameter). absolute expression value at 0 and 10 CPM. ChIP experiments were performed with the MAGnify Chromatin Immunoprecipitation System (49-2024, Thermo Fisher Scientific) according to the manufacturer’s instructions, Analysis of Distance Distribution of TEs to Neighboring with the following modifications. Plant material was fixed in Genes 20 ml Crosslinking Buffer (0.4 M Sucrose, 10 mM Tris–HCl, pH Distances of differential and nondifferential TEs to the nearest 8.0, 10 mM MgCl2, 1% formaldehyde) by applying vacuum gene or gene 5 -end were binned into classes 0–1, 1–2, 2–3, twice for 15 min. Nuclei and chromatin were purified as in He and 3 kb. (TEs overlapping a gene were not considered.) A et al. (2012). Chromatin was sonicated with a BioRuptor de- pseudocount of 1 was added to each bin count. A saturated vice (Diagenode, Lie ` ge, Belgium) for 30 times 30 s at high log-linear model of the counts was computed using the R setting with 30 s intermittent cooling in ice-water. A DNA frag- glm() function (family ¼ “poisson”), using sum-to-zero con- ment size of 300–600 bp was controlled by running an aliquot trasts. The P value for distance d in supplementary table 1, of decrosslinked and purified DNA on 1.5% agarose gels. The Supplementary Material online, is the Pr(>jzj) value of the following antibodies were used in immunoprecipitation: anti- interaction coefficient between the factors distance category rat IgG (R9255, Sigma; St. Louis, MO), as well as anti- and significance of hybrid-dependent TE expression. H3K27me3, anti-H3K9me2 and anti-H3 (07-449, 07-441, and 07-690, respectively, Millipore; Temecula, CA). To reach Random Expectations of Superfamily Frequencies enough precipitated DNA, the ChIP-DNA precipitation proto- To evaluate the enrichment of TEs changing their expression col was performed independently for eight aliquots of the in hybrid background within each superfamily, we computed collected leaf material, pooled and purified by Qiagen a random expectation. For this, we used the set of TEs show- MinElute Reaction Cleanup Kit. Finally, we obtained ChIP- ing differential expression in the hybrid background as a ref- DNA samples from Col-0, MN47, and Col-0xMN47 with anti- erence set and performed 10,000 random resamplings of the bodies targeting H3K27me3 and H3K9me2 marks. In addi- same number of TEs from all TEs with evidence of expression. tion, we generated a ChIP-DNA sample targeting Histone 3 To account for the skewed distribution of these TEs in the in a Col-0xMN47 sample to verify that both genomes could be vicinity of non-TE genes, both the reference set and the set of efficiently immuno-precipitated and assess ChIP-seq quality all TEs were binned by distance to the closest gene 5 -end. For (see supplementary Methods, Supplementary Material online). Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 1405 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Go ¨ bel et al. GBE Table 1 Numbers and Fractions of Differentially Expressed TEs at Two Different False-Discovery Rates (FDR 0.01 or0.05) and Two Different Thresholds on the Normalized Counts Per Million (CPM) (all: 0CPM; 10 CPM) FDR £ 0.01 FDR £ 0.05 All 10 CPM All 10 CPM Arabidopsis thaliana 3 þ 2(0.3%) 1 þ 1(0.1%) 5 þ 5(0.6%) 1 þ 2(0.2%) Arabidopsis lyrata 57 þ 2(0.8%) 14 þ 0(0.2%) 73 þ 6(1.1%) 18 þ 1(0.3%) NOTE.—Numbers are given as “number of up-regulations” þ “number of down-regulations.” Fractions were computed relative to all TEs with evidence of expression in standard conditions, discounting members of genome-specific families. reads were filtered with samtools (Li et al. 2009)for a map- ChIP-Seq Library Preparation and Sequencing ping quality of 3 (-q 3), which discards multimapping reads DNA was sheared with Ion Shear reagents and the Ion Xpress with two equally good mappings but retains those with a Fragment Library Kit (Thermo Fisher Scientific, Waltham, MA) single best mapping location. Duplicate reads from the am- was used for subsequent library preparation. [AQ10] After plification step during library preparation were removed with ligation of the respective Ion Xpress barcode adapters, sam- the MarkDuplicates program from Picard Tools (http://broad- ples were amplified by 17 or 23 PCR cycles. Amplified libraries institute.github.io/picard; last accessed May 25, 2018). ChIP- were quantified on Bioanalyzer DNA High-Sensitivity Chips seq read counts for TEs and non-TE genes were counted fol- (Agilent Technologies, Santa Clara, CA) and diluted to lowing the same pipeline as for RNAseq reads. A first analysis 9 pM. Emulsion PCR was performed on the Ion OneTouch of ChIP-seq count distribution confirmed that the distribution System, followed by enrichment for template positive Ion of marks recovered known chromatin domains but revealed Sphere Particles. Sequencing sample was loaded on an Ion 64 genomic regions giving aberrant ChIP signals (see supple- Proton chip and sequenced on the Ion proton sequencer mentary Methods, Supplementary Material online). These 64 (Thermo Fisher Scientific) according to the manufacturer’s regions were removed from all analyses presented in the instructions. manuscript. ChIP-Seq Read Processing, Mapping and Peak Detection Results Raw data were preprocessed using the TORRENT SUITE ver- Only a Minor Fraction of All Expressed TEs Is Differentially sion 4.0.1 (Thermo Fisher Scientific) to trim adapter sequen- Expressedinthe Hybrid ces. Using the FASTX-Toolkit (http://hannonlab.cshl.edu/ fastx_toolkit/; last accessed May 25, 2018), reads were then Our annotation contains 10863 TEs on the A. thaliana Col- trimmed from the 3 end to remove low quality bases (phred 0 genome and 30,868 elements on the A. lyrata MN47 <15). Reads shorter than 30 bases and reads of poor overall genome. After excluding genome-specific TE families, quality (more than half of the bases with phred< 20) were 10,805 and 29,150 elements remained on the Col-0 and discarded. The mapping of the reads was performed with MN47 genomes, respectively. Of these, 1,592 (15%) of the bowtie2 (Langmead and Salzberg 2012), with the end-to- A. thaliana TEs and 7,045 (24%) of the A. lyrata TEs end parameter, without allowing mismatches in the seed (- showed evidence of expression in at least one replicate. N 0), and allowing for up to two alignments per read (-k 2); We quantified TE expression in normalized Counts Per although then only the best alignment was kept. All samples Million (CPM, see Materials and Methods). TEs contributed were mapped to a hybrid reference genome obtained by but a minor fraction of all reads (on average, only 0.07% in concatenating the genomes of the parent as described above A. thaliana samples and 0.24% in A. lyrata samples). This for RNAseq read mapping. Cross-hybridization between fraction remained comparable in the hybrid transcriptome genomes was negligible (0.5–1% of mapped reads). A sum- (on average 0.06% and 0.30% of A. thaliana and A. lyrata mary of ChIP seq read data are given in supplementary table reads, respectively in the hybrid). Only a small fraction 2, Supplementary Material online. (maximally 2%) of the expressed TEs of both genomes Quality assessment of the assays was performed using were significantly differentially expressed between paren- phantompeaktools (Kharchenko et al. 2008; Landt et al. tal and hybrid backgrounds, using different cutoffs to de- 2012) to calculate NSC and RSC metrics (normalized and rel- fine significance (table 1). Moreover, this percentage ative strand coefficients, respectively), and the bioconductor dropped below 0.5% when the median expression was ChIPQC package (Carroll et al. 2014) to calculate the square required to be above ten CPM in at least one of the two sum of deviation. NSC was low but RSC metrics were well backgrounds (table 1). These observations show that TE within the recommended range (Landt et al. 2012). Mapped transcription is not severely disturbed in the hybrid. 1406 Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Robustness of Transposable Element Regulation GBE FIG.1.—Correlation of TE expression levels, estimated as the median of normalized counts per million over the replicate plants grown in standard conditions, in parent and hybrid for individual TEs, TE families, and superfamilies. TEs or (super)families with significant expression change are marked in red. In A. lyrata, the 59 (0.8%) TEs that respond significantly to Specific Rather than Global Effects Act on TE Expression in hybridization at FDR 0.01 were generally induced. In A. the Hybrid thaliana, the 5 (0.3%) affected elements were either induced In order to investigate the effect of hybridization at the levels or repressed (supplementary fig. S1, Supplementary Material of individual elements, families and superfamilies of TEs, we online). SINEs dominated by Sadhu elements (Rangwala and correlated expression of individual TEs and combined expres- Richards 2010) and Harbingers dominated by the ATIS112A sion levels of TE families between parents and hybrids (fig. 1). family (Kapitonov and Jurka 2004) stood out as enriched Although TE transcription in the two genetic backgrounds among the responding A. lyrata TEs (fig. 2). These superfami- correlated well at all levels, the correlation increased from lies were significantly overrepresented relative to the set of the individual TEs to the family and superfamily levels. expressed TEs (1 RNA-seq count, permutation test, Spearman’s rank correlation of the median expression of rep- FDR< 0.05). The fact that only a small number of TE families licates in the two backgrounds under control conditions was showed enrichment for differential expressed elements indi- 0.73, 0.89, 0.97 on the three levels of individual TE, family, cates that hybridization does not have a global effect on TE and superfamily, respectively in A. thaliana, and 0.65, 0.86, expression. 0.91, respectively, in A. lyrata. Including four replicates in the RNA-seq analysis was necessary to establish these relation- ships with increased confidence. Using subsets of only two Rare TE Expression Changes Coincide with Alterations in replicates, for example, reduced the parent–hybrid correlation Nearby Gene Expression coefficient of TE expression to 0.53, 0.77, 0.87 on the three levels for A. thaliana, and to 0.51, 0.84, 0.89, respectively, for TEs can change expression of nearby genes by influencing A. lyrata. local chromatin accessibility. We therefore asked whether The small subset of TEs with significantly different expres- hybrid-specific expression changes of TEs predict a corre- sion level between parent and hybrid was restricted to lower sponding change in neighboring genes. level units (individual TEs and families), indicating that the Compared with TEs with expression not affected by hybrid- differences are case-specific rather than caused by a global ization, A. lyrata TEs expressed differentially in the hybrid were disruption of TE silencing. depleted in genome regions distant more than 3kb from Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 1407 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Go ¨ bel et al. GBE FIG.2.—Observed and expected superfamily distribution of differentially expressed TEs (parent vs. hybrid). Colored bars show the observed numbers of TEs showing hybrid-dependent expression in each TE superfamily. Each color corresponds to one superfamily as indicated in the legend. Superfamily names are as given in the annotation of (Pietzenuk et al. 2016). Names with trailing question marks were annotated as deviant/uncertain members of the respective superfamily. Transparent bars indicate expected median count for random sampling within each superfamily. Error bars indicate the 5% and 95% quantile of the 10,000 random draws. Single black stars mark superfamilies whose observed count is significantly higher than expected from random expectations. Gray stars indicate depletion (maximally 5% of the universe have counts the observed count of the respective superfamily). Double stars indicate significant enrichment or depletion that is robust to FDR correction. Table 2 Distance Distributions to Nearby Genes of Arabidopsis lyrata TEs with Differential Expression in Parent and Hybrid Under Control Conditions Distance to Closest Gene Closest 5 End of a Gene <1 kb 1–2 kb 2–3 kb 3kb <1 kb 1–2 kb 2–3 kb 3kb Number of TEs Diff. exp. PvsH 15 18 13 13 7 24 12 16 Not Diff. Expr. 1,238 1,675 913 3,160 654 1,455 1,078 3,799 P (frequency in distance bin) 0.239 0.539 0.061 4.0e-4 0.554 0.01 0.517 1.5e-4 NOTE.—The table combines distance information of TEs with a significant change in expression in the hybrid to either the nearest gene (any orientation) or the nearest 5 of any gene. Numbers of significant and nonsignificant TEs are listed for each distance bin (FDR 0.01). In the lower part, P values are given for the contribution of each distance bin to the difference between the distributions of significant and nonsignificant TEs (see Materials and Methods). Bins with significant depletion or enrichment in significant TEs are highlighted in blue and orange, respectively. 0 4 the closest gene 5 ends (log-linear model, P 5  10 , preferential location of hybridization-sensitive TEs in euchro- table 2). In addition, counts of hybridization-sensitive TEs matic, gene-rich regions in A. lyrata.In A. thaliana, there were were enriched in the bins 1–2 kb upstream of gene 5 -ends too few differentially expressed TEs in the hybrid background (log-linear model, P 0.01, table 2). This indicates a to reliably compare the distance distributions. 1408 Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Robustness of Transposable Element Regulation GBE We further analyzed if changed expression of TEs and their neighboring genes in a hybrid background were correlated. Hybrid-specific changes in transcription levels of non-TE genes have already been reported (Zhu et al. 2017). We also quan- tified non-TE gene expression levels in parents and hybrids. Note that in our experiment, the timing of sampling accounted for the developmental differences between the fast developer A. thaliana and both A. lyrata and the hybrid, which are late flowering. Markedly fewer genes were af- fected by hybridization in our experiment (234 in A. thaliana and 395 in A. lyrata), yet, we observed a significant overlap in genes reported to change expression in the hybrid back- ground in Zhu et al. (2017) (Hypergeometric test of overlap in genes changing their expression level between hybrid and parents, P< 10 for all comparisons). In A. lyrata 153 TEs with detectable expression were within 2 kb of a gene with hybridization-sensitive expression (FDR 0.01). Among these, 8 (5%) TEs displayed a change in expression in the hybrid background. This percentage dropped to 1% (71/6,892) for TEs with a gene neighbor showing no significant change in expression in hybrid versus parent (P< 5  10 ,hypergeo- metric test). In the A. thaliana genome, the number of expressed TEs located in the vicinity of genes was too small to calculate correlations. We conclude that hybrid-specific TE expression can coincide with a modification of the expression of nearby protein-coding genes, but this only affects a very small number of expressed loci. Histone Marks Are Not Modified in the Hybrid Background Two epigenetic marks are associated with transcriptional silencing in plants: the plastic mark H3K27me3 and the permanent mark H3K9me2 (Kim and Zilberman 2014; Sequeira-Mendes et al. 2014). We profiled these two histone modifications by Chromatin Immuno-Precipitation followed by sequencing (ChIP-Seq) in hybrids and A. thaliana and A. lyrata parents. Although the ChIP-seq H3K27me3 yielded a comparatively lower number of reads, especially for Col-0 (supplementary tables 2, 5,and 6, Supplementary Material online), the distribution of marks on TEs was in agreement with previous reports (fig. 3A, Seymour et al. 2014) and cap- tured well the chromatin domains that were defined in A. thaliana (Sequeira-Mendes et al. 2014, supplementary Methods, Supplementary Material online). In our data, H3K9me2 marks showed a high correlation of normalized read counts between parental and hybrid genomes, both genome-wide (calculated across 10-kb genomic windows) and for individual TEs, TE families, and superfamilies. The FIG.3—Correlation of histone mark levels on TEs in parent and hybrid. genome-wide Spearman correlation coefficients of (A) Genome-wide distribution of H3K9me2 and H3K27me2 on the pa- H3K9me2 read count along the genomes were 0.66 and rental and hybrid genomes. The circular plots represent the hybrid genome used as reference, with the five A. thaliana chromosomes in green and the eight A. lyrata chromosomes in red. Inner tracks represent gene, TE, FIG.3—Continued H3K9me2 and H3K27me3 densities in 500 kb bins. Top left: Col-0, Top counts. The color code indicates superfamily membership, using the right: MN47, Bottom: hybrid. (B and C) Scatterplots of ChIP-seq read palette of figure 2. (B) H3K9me2; (C) H3K27me3. Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 1409 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Go ¨ bel et al. GBE FIG.4—No change in histone mark occupancy is detected for differentially expressed TEs. Shown are scatterplots of each histone marks on each of the two genomes for TEs with significantly different expression level in parent versus hybrid under standard conditions. Superfamily membership of the TEs is indicated using the palette of figure 2. Symbol shapes indicate the direction of gene expression change in hybrid versus parent (upward triangle: up, downward: down). 0.92 for A. thaliana or A. lyrata, respectively. Spearman’s rank Nevertheless, the correlation of H3K27me3 occupancy in correlation coefficients of H3K9me2 for individual TEs, TE- parent and hybrid remained high even for TEs, in A. lyrata families, and superfamilies were 0.76, 0.94, and 0.96 for A. (Spearman’s rank correlation 0.58, 0.78, 0.92 for individual thaliana and 0.87, 0.97, and 0.99 for A. lyrata (fig. 3B and C). TEs, families, and superfamilies). For the A. thaliana/hybrid Although the high correlation values pointed to a largely ad- comparison, the correlation dropped to 0.2, 0.27, 0.21 for ditive inheritance of the mark in the hybrid on all levels of the individual TEs, families, and superfamilies, respectively annotation, the fact that correlation was lowest on the level (fig. 3C). Note however that the Col-0 H3K27me3 ChIP- of individual elements again underscores the role of element- seq experiments yielded a comparatively lower number of specific responses in the hybrid. reads (supplementary table 2, Supplementary Material While the function of H3K9me2 is to permanently silence online). TEs, H3K27me3 mainly serves to reversibly silence genes We tested then whether changes in histone mark occu- whose expression is restricted to specific conditions or de- pancy could explain TE expression differences between velopmental time windows (Gan et al. 2015). The distribu- parents and hybrids. TEs with a modified expression level in tion of H3K27me3 marks was strongly correlated genome- hybrid background generally showed similar H3K27me3 level wide between parental and hybrid backgrounds (Spearman in parent and hybrid (fig. 4). When a departure in epigenetic correlation coefficients: 0.63 and 0.86 for A. thaliana and A. levels was detected between parents and hybrids, the direc- lyrata, respectively, fig. 3A). H3K27me3 is not a typical si- tion of the change in expression of the TE was not necessarily lencing mark of TEs (Sequeira-Mendes et al. 2014). This was in the expected direction (fig. 4). A. lyrata TEs with a modified most apparent in our data for the LTR superfamily, which expression level in hybrid background also generally showed tended to be associated with the highest level of H3K9me2 similar H3K9me2 mark levels in parents and hybrids. We and the lowest level of H3K27me3 (fig. 3B and C). detected a concordant loss of H3K9me2 marks and increased 1410 Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Robustness of Transposable Element Regulation GBE Table 3 Numbers and Fractions of Differentially Expressed TEs in Stress Conditions (3h) and in Response to Stress (3h specific), at Two Different False-Discovery Rates (FDR 0.01 or 0.05) and Two Different Thresholds on the Normalized Counts Per Million (CPM) (all: 0CPM; 10 CPM) FDR £ 0.01 FDR £ 0.05 All 10 CPM All 10 CPM Arabidopsis thaliana 3h 4 þ 14 (1.1%) 2 þ 4(0.4%) 7 þ 26 (2.1%) 2 þ 5(0.4%) 3h specific 3 þ 12 (0.9%) 2 þ 3(0.3%) 4 þ 23 (1.7%) 2 þ 4(0.4%) Arabidopsis lyrata 3h 33 þ 18 (0.7%) 4 þ 4(0.1%) 54 þ 37 (1.3%) 10 þ 15 (0.4%) 3h specific 9 þ 17 (0.4%) 3 þ 4(0.1%) 12 þ 34 (0.7%) 7 þ 14 (0.3%) NOTE.—Fractions were computed relative to all TEs with evidence of expression under stress conditions, discounting members of genome-specific families. expression in the hybrid for 6 A. lyrata TEs (fig. 4). But there conditions, the correlation of TE expression levels between again, we observed for 2 TEs that the change in H3K9me2 parents and hybrids increased at the family and superfamily occupancy did not reflect the change in expression. levels. The Spearman’s rank correlation of the replicates’ me- Altogether, in our data, a differential deposition of dian expression in the two backgrounds was (0.7, 0.81, 0.83) H3K27me3 and H3K9me2 histone methylation marks in the for A. thaliana, and (0.67, 0.91, 0.94) for A. lyrata, for indi- hybrid is rare and generally does not coincide with hybrid- vidual, family and superfamily levels, respectively. dependent changes in TE expression. Interestingly, we observed that more TEs were down- regulated in the hybrids compared with the parents under stress conditions than under standard conditions (fig. 5A and table 3). Only the Harbinger and LINE_L1 superfamiles Robustness of TE Expression to Hybridization Is Maintained were enriched among expressed A. lyrata TEsinstresscon- under Stress Conditions ditions (fig. 5B; permutation test using TEs with 1 RNA-seq It is assumed that transposition activity is induced by stress count as reference, FDR< 0.05). This further corroborates the conditions (McClintock 1984; Cavrak et al. 2014). conclusion that in our two species, TE regulation is robust to Dehydration stress is associated with a vast transcriptome hybridization even under stress conditions. reprogramming (Matsui et al. 2008). We analyzed the depen- dence of transcript abundance upon interspecific hybridiza- tion for TEs and non-TE genes under strong drought stress Discussion conditions in the same way as described for control condi- No evidence for a genome shock upon hybridization of tions. Having found that TEs with a hybrid-specific expression A. lyrata and A. thaliana change tend to be close to genes, and given the reprogram- The term “genomic shock,” was initially used to describe the ming of the gene transcriptome by the stress, we asked breakage and large scale rearrangement of chromosomes in whether TE regulation was robust to hybridization even under maize (McClintock 1946, 1984). This term is thus appropriate stress conditions. In our data, 30–40% of the genes of both to describe genome instability following interspecific hybridi- species changed expression after 3h of dehydration, with sim- zation, that is, unbalanced segregation of homoelogous chro- ilar proportions in the parental and hybrid backgrounds. The mosomes, chromosome loss or translocations (Xiong et al. fraction of stress responsive TEs was much smaller: 3–4% on 2011; Chester et al. 2012; Henry et al. 2014; Spoelhof the A. thaliana genome and 2–3% for A. lyrata (not shown). et al. 2017). Over the years, the term “genome shock” has In addition, we observed no increase in the variance of TE also been used to describe the anticipated alteration in ge- expression in hybrid context (supplementary fig. 3, nome activity immediately following interspecific hybridiza- Supplementary Material online). Although the fraction of tion or allopolyploidization, although transposon A. thaliana TEs with a significant expression difference be- reactivation has seldom been rigorously documented imme- tween parent and hybrid is higher under dehydration stress diately after genome merging (Parisod et al. 2010). In plants, than under control conditions, both the absolute numbers (4 DNA cytosine methylation (mC) and the histone mark up and 14 down) and the fraction (1.1% of all TEs with read H3K9me2 embed TEs in a chromatin context unfavorable count1) remained small (table 3). The fraction of respond- for transcription (Kim and Zilberman 2014). The establishment ing A. lyrata TEs decreased slightly compared with control and maintenance of these marks is not only inherited through conditions, and again the numbers remained small (33 up, DNA replication, it is also guided by small interfering RNAs 18 down, fraction 0.7%). The percentage dropped below (siRNAs) (see Kim and Zilberman 2014 for a review of plant TE 0.5% if an expression level of10 CPM in at least one of silencing mechanisms). Because of their potential to exert the backgrounds was required (table 3). Like under control Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 1411 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Go ¨ bel et al. GBE genome-wide trans effects via base pair homology of the short 24 nt siRNAs, the population of siRNAs expressed by a genome must be adapted to its gene content to ensure proper gene expression. Interspecific hybridization, therefore, brings together trans-regulators of TE expression that are not mutually adapted and could deeply disorganize gene regula- tion. Gene expression modifications following interspecific hy- bridization or allopolyploidization have been monitored across plant genera as diverse as Gossypium, Senecio, Brassica, Tragopogon, or Oryza (e.g., Hegarty et al. 2006; Buggs et al. 2009;Yoo et al. 2013). Many studies have concluded that epigenome incompatibilities are important contributors to hybrid sterility (Parisod et al. 2009; Ghani et al. 2014; Wang et al. 2014; Xu et al. 2014; Senerchia et al. 2015; Wu et al. 2016). In Aegilops hybrids, for example, parental divergence in TE content was associated with the emergence of new TE copies in hybrids (Senerchia et al. 2015). Allopolyploidy is also believed to be associated with drastic epigenetic reorganiza- tion although evidence for a significant reactivation of TE transposition remains scant (Parisod et al. 2010). Our analysis contradicts the conclusions of these previous reports. Comparing the TE expression and epigenetic profile of the A. lyrata and A. thaliana genomes in parental and hy- brid context, we find no evidence for a major transcriptomic or epigenomic “shock” driven by inappropriate TE silencing and epigenetic mechanisms. Both TE transcription and the decoration of TEs by the silencing histone mark H3K9me2 are globally unaffected in the hybrid. We observe that only a handful of specific TEs show a change in expression in hybrids. In addition, the changes in TE expression are generally of small magnitude. We further show that the distribution of silencing epigenetic marks is strongly correlated in parents and hybrids, and its variation does not match the changes in TE expression. Finally, under severe dehydration stress, a FIG.5.—(A) Correlation of TE expression levels in stress conditions, large number of genes change their transcriptional activity estimated as Median normalized counts per million over the replicates, in (Matsui et al. 2008). This broad genome-wide remodeling parent and hybrid for individual TEs, TE families, and superfamilies. TEs or of transcriptional activity, however, far from inducing a sev- (super)families with significant expression change are marked in red (sig- ered genome shock, tends to increase the robustness of TE nificant only under control conditions), orange (significant both under expression patterns to hybridization. These findings thus sug- standard and stress conditions), or yellow (significant only in stress con- gest that stress does not increase the impact of hybridization, ditions). (B) Observed and expected superfamily distribution of TEs whose further supporting the idea that the control of TE elements is expression is hybrid-dependent specifically in stress conditions. Colored robust to genome merging. Although these observations par- bars show the observed numbers of TEs showing hybrid-dependent ex- tially contradict previous reports, they are not isolated. In re- pression in each TE superfamily. Each color corresponds to one superfamily cent years, an increasing number of studies, mostly as indicated in the legend. Superfamily names are as given in the annota- conducted in the Brassicaceae family, has begun to report tion of Pietzenuk et al. (2016). Names with trailing question marks were annotated as deviant/uncertain members of the respective superfamily. cases where hybridization does not massively alter gene reg- Transparent bars indicate expected median count for random sampling ulation (Akama et al. 2014; Shi et al. 2015; He et al. 2016; within each superfamily. Error bars indicate the 5% and 95% quantile of Zhang et al. 2016; Zhu et al. 2017). In addition, in tomato, the 10,000 random draws. Single black stars mark superfamilies whose transgressive expression of small RNAs was associated with observed count is significantly higher than expected from random expect- the dysregulation of stress responsive genes, but not with that ations. Gray stars indicate depletion (maximally 5% of the universe have of TEs (Shivaprasad et al. 2012). counts the observed count of the respective superfamily). Double stars The contradiction between these studies and previous indicate significant enrichment or depletion that is robust to FDR reports may have various explanations. First, it is striking to correction. note that many of the studies that conclude on the absence of 1412 Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Robustness of Transposable Element Regulation GBE a global gene expression dysregulation were conducted in the A. thaliana allele of FLC, which is up-regulated by species with well-known genomes, where gene expression FRIGIDA, is one of the alleles most strongly impacted by hy- levels can be quantified with the best accuracy. They were bridization (supplementary table 3, Supplementary Material not conducted in polyploids with large and imperfectly assem- online; Zhu et al. 2017). Interspecific differences in trans-act- bled genomes. Second, biological and experimental variance ing factor activity will therefore determine the extent to which between replicates affects the strength of the correlation of gene expression changes upon hybridization. If F1 hybrids are expression in parental versus hybrid backgrounds, and thus fertile, these changes can segregate in subsequent F2 gener- the apparent global impact of hybridization. Finally, it is nec- ations as in, for example, Shivaprasad et al. (2012). essary to contrast diverse sources of information to test Here, focusing on the trans-acting effect of the passage whether hybridization has a broad impact on gene regulation. from a native parental background to a hybrid background on In the allotetraploid A. suecica, modification in siRNA levels TE expression, we show that at most 1–2% of expressed TEs was not associated with nonadditive gene expression in newly have expression that is differentially controlled in trans. formed allopolyploids, suggesting that changes in siRNA ex- Several elements indicate that these trans-acting effects may pression were not the result of epigenetic disorganization in be associated with interspecific differences in mechanisms of the hybrid (Ha et al. 2009). Similarly, in Arabidopsis homo- gene regulation in euchromatic regions. When combined in a ploid hybrids, we observed that hybrid-specific TE expression hybrid, these trans-acting differences modify the regulatory changes do not coincide with obvious epigenetic changes, a environment of a few TE subfamilies. Indeed, the small pro- pattern also reported for protein-coding genes (Zhu et al. portion of TEs with differential expression in hybrids tend to 2017). Taken together, the different molecular data sets be located in the vicinity of genes. Interspecific differences in used in this study show that TE expression and the distribution the regulation of euchromatic regions are however likely to be of epigenetic marks, which control gene expression, are not minor, given the small number of elements affected. massively disorganized in hybrids. Interspecific differences in gene regulation controlled in cis can also create a pattern of genome dominance in hybrids, with a majority of orthologous transcripts showing a transcrip- trans-and cis-Acting Differences between Species Can tion bias towards one parental genome (Yoo et al. 2013; Manifest in F1 Hybrids Cheng et al. 2016; Zhang et al. 2016). For example, in The comparative analysis of gene expression in hybrids and A. thaliana  A. lyrata hybrids, a bias towards preferential their parents reflects the variability of trans-and cis-regulatory expression of the A. lyrata allele has been reported by inde- mechanisms controlling gene expression in the two species pendent studies (He et al. 2012, 2016; Zhu et al. 2017). (Wittkopp et al. 2004). After hybridization, changes in total Cis-effects have not been examined in this study, because expression will simply reflect the number and magnitude of orthologous relationships cannot be established for most trans-regulatory differences between species, whereas the TEs, but for the subset of TEs that are located in orthologous specific cis-regulatory variant of each transcript will result in positions in both species, the A. lyrata TE copy tended to be allelic expression differences. These differences will manifest more highly expressed than the A. thaliana copy in the hybrid in hybrids on a restricted number of genes or elements, which (He et al. 2012). does not imply that there is a major “genomic shock”. Only trans-acting differences will have nonadditive conse- Incompatibilities Can Also Have an Oligogenic Basis quences on expression. In fungi and cotton allopolyploid, >50% of homeologous genes inherit the expression pattern It may seem at first surprising that TE silencing shows so little of the parents (Yoo et al. 2013; Cox et al. 2014). In rice, allelic disruption in interspecific A. thaliana  A. lyrata hybrids, given expression bias in hybrids between Oryza sativa and O. japon- the estimated 6 My of divergence and their vastly different TE ica also reflects mostly interspecific differences. In hybrids be- content (Hu et al. 2011; Hohmann et al. 2015). The proba- tween the diploid O. sativa and the tetraploid O. punctata, bility to observe specific disruptions in gene expression in a only 16% of the genes showed an expression pattern differ- hybrid background, however, does not necessarily scale with entfrom thatofthe parents (Wu et al. 2016). The analysis of evolutionary distance. The proportion of genes with hybrid- protein-coding gene expression in A. thaliana  A.lyrata specific changes in expression is not greater for hybrids be- hybrids reported a larger proportion of nonadditive gene ex- tween the closely related subspecies O. sativa and O. japonica pression changes for the A. thaliana genome compared with than between O. sativa and the more distant species O. punc- the A. lyrata genome (Zhu et al. 2017). This pattern, like in our tata (Wu et al. 2016). A simple oligogenic combination of study, could not be associated with chromatin modifications alleles, in fact, can cause what is known as a Dobzhansky– induced by hybridization or any other genome-wide change Muller incompatibility, even in the absence of large in gene regulatory mechanisms (Zhu et al. 2017). However, evolutionary distances between parents. Dobzhansky–Muller the A. thaliana parental genotype Col-0 is missing an active incompatibilities, which can cause major disruption at the allele of the major developmental regulator FRIGIDA.In fact, phenotypic and, sometimes at the transcriptome level, are Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 1413 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Go ¨ bel et al. GBE Berardini TZ, et al. 2015. The Arabidopsis Information Resource: making not seldom within species, and were even detected within and mining the “gold standard” annotated reference plant genome. populations (Alcazar et al. 2010; Bomblies and Weigel Genesis 53(8):474–485. 2010; Durand et al. 2012; Wright et al. 2013; Chae et al. Bomblies K, Weigel D. 2007. Hybrid necrosis: autoimmunity as a potential 2014). We conclude that the term “shock” is not appropriate gene-flow barrier in plant species. Nat. Rev. Genet. 8(5):382. to describe the transcriptional consequences of the merging Bomblies K, Weigel D. 2010. Arabidopsis and relatives as models for the study of genetic and genomic incompatibilities. Philos Trans R Soc B of divergent genomes in a viable F1 hybrid. Our data rather Biol Sci. 365(1547):1815–1823. supports a less dramatic scenario, where the merging of two Buggs RJ, et al. 2009. Gene loss and silencing in Tragopogon miscellus genomes creates a new trans-acting background, in which (Asteraceae): comparison of natural and synthetic allotetraploids. oligogenic Dobzhansky–Muller incompatibilities can manifest. Heredity 103(1):73–81. Such incompatibilities are known to affect endosperm func- Carroll TS, Liang Z, Salama R, Stark R, de Santiago I. 2014. Impact of artifact removal on ChIP quality metrics in ChIP-seq and ChIP-exo tion (Rebernig et al. 2015) and will also occasionally affect the data. Front Genet. 5:75. regulation of specific TEs and favor new bouts of transposi- Cavrak VV, et al. 2014. How a retrotransposon exploits the plant’s heat tion. Without a global robustness of TE regulation to genome stress response for its activation. PLoS Genet. 10(1):e1004115. merging, hybridization and polyploidization would probably Chae E, et al. 2014. Species-wide genetic incompatibility analysis identifies be much less important for plant diversification (Alix et al. immune genes as hot spots of deleterious epistasis. Cell 159(6):1341–1351. 2017). We note, however, that the impact of hybridization Cheng F, et al. 2016. Epigenetic regulation of subgenome dominance on TE transposition rate has not been examined in this study. following whole genome triplication in Brassica rapa. New Phytol. It thus remains possible that TEs with undetectable expression 211(1):288–299. level transpose at a rate that is higher in F1 interspecific Chester M, et al. 2012. Extensive chromosomal variation in a recently hybrids than within species. Within A. thaliana,elements of formed natural allopolyploid species, Tragopogon miscellus (Asteraceae). Proc Natl Acad Sci U S A. 109(4):1176–1181. about half of the TE families were shown to have recently Cox MP, et al. 2014. An interspecific fungal hybrid reveals cross-kingdom transposed (Quadrana et al. 2016). rules for allopolyploid gene expression patterns. PLoS Genet. 10(3):e1004180. de la Chaux N, Tsuchimatsu T, Shimizu KK, Wagner A. 2012. The pre- Supplementary Material dominantly selfing plant Arabidopsis thaliana experienced a recent Supplementary data areavailableat Genome Biology and reduction in transposable element abundance compared to its out- crossing relative Arabidopsis lyrata. Mobile DNA 3(1):2. Evolution online. De Meaux J, Pop A, Mitchell-Olds T. 2006. Cis-regulatory evolution of Chalcone-Synthase expression in the genus Arabidopsis.Genetics 174(4):2181–2202. Acknowledgments Dobin A, et al. 2013. STAR: ultrafast universal RNA-seq aligner. We thank Franziska Turck (Max Planck Institute for Plant Bioinformatics. 29(1):15–21. Durand S, Bouch e N, Perez Strand E, Loudet O, Camilleri C. 2012. Rapid Breeding Research) for sharing her ChIP-seq protocol and ex- establishment of genetic incompatibility through natural epigenetic perience. This research was supported by the Deutsche variation. Curr Biol. 22(4):326–331. Forschung Gesellschaft (DFG) with grant ME2742/2-1 and Gan ES, Xu Y, Ito T. 2015. Dynamics of H3K27me3 methylation and de- ME2742/6-1 in the realm of the special study program methylation in plant development. Plant Signal Behav. SPP1529, as well as by the European Research Council with 10(9):e1027851. Grant 648617 “AdaptoSCOPE.” The raw data has been de- Ghani MA, et al. 2014. The role of small RNAs in wide hybridisation and allopolyploidisation between Brassica rapa and Brassica nigra.BMC posited in NCBI Sequence Read Archive (SRA) and are acces- Plant Biol. 14:272. sible through SRA study number SRP148881 and SRP148853 Ha M, et al. 2009. Small RNAs serve as a genetic buffer against genomic for RNA-seq and ChIP-seq data, respectively. shock in Arabidopsis interspecific hybrids and allopolyploids. Proc Natl Acad Sci U S A. 106(42):17835–17840. He F, et al. 2012. Widespread interspecific divergence in cis-regulation of Literature Cited transposable elements in the Arabidopsis genus. Mol Biol Evol. Akama S, Shimizu-Inatsugi R, Shimizu KK, Sese J. 2014. Genome-wide 29(3):1081–1091. quantification of homeolog expression ratio revealed nonstochastic He F, et al. 2016. The footprint of polygenic adaptation on stress- gene regulation in synthetic allopolyploid Arabidopsis. Nucleic Acids responsive cis-regulatory divergence in the Arabidopsis genus. Mol Res. 42(6):e46. Biol Evol. 33(8):2088–2101. Alc azar R, et al. 2010. Natural variation at Strubbelig Receptor Kinase 3 Hegarty MJ, et al. 2006. Transcriptome shock after interspecific hybridiza- drives immune-triggered incompatibilities between Arabidopsis thali- tion in senecio is ameliorated by genome duplication. Curr Biol. ana accessions. Nat Genet. 42(12):1135–1139. 16(16):1652–1659. Alix K, et al. 2017. Polyploidy and interspecific hybridization: partners for Henry IM, et al. 2014. The BOY NAMED SUE quantitative trait locus con- adaptation, speciation and evolution in plants. Ann Bot. fers increased meiotic stability to an adapted natural allopolyploid of 120(2):183–194. Arabidopsis. Plant Cell 26(1):181–194. Benjamini Y, Hochberg Y. 1995. Controlling the false discovery rate: a Hohmann N, Wolf EM, Lysak MA, Koch MA. 2015. A time-calibrated road practical and powerful approach to multiple testing. J Roy Stat Soc map of Brassicaceae species radiation and evolutionary history. Plant B. 57(1):289–300. Cell 27(10):2770–2784. 1414 Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Robustness of Transposable Element Regulation GBE Hollister JD, Gaut BS. 2009. Epigenetic silencing of transposable elements: Robinson MD, Oshlack A. 2010. A scaling normalization method for dif- a trade-off between reduced transposition and deleterious effects on ferential expression analysis of RNA-seq data. Genome Biol. 11:R25. neighboring gene expression. Genome Res. 19(8):1419–1428. Senerchia N, Felber F, Parisod C. 2015. Genome reorganization in F1 Hollister JD, et al. 2011. Transposable elements and small RNAs contribute hybrids uncovers the role of retrotransposons in reproductive isolation. to gene expression divergence between Arabidopsis thaliana and Proc Biol Sci. 282(1804):20142874. Arabidopsis lyrata. Proc Natl Acad Sci U S A. 108(6):2322–2327. Seki M, et al. 2002. Monitoring the expression profiles of 7000 Arabidopsis Hu TT, et al. 2011. The Arabidopsis lyrata genome sequence and the basis genes under drought, cold and high-salinity stresses using a full-length of rapid genome size change. Nat Genet. 43(5):476–481. cDNA microarray. Plant J. 31(3):279–292. Jin Y, Tam OH, Paniagua E, Hammell M. 2015. TEtranscripts: a package for Seymour DK, Koenig D, Hagmann J, Becker C, Weigel D. 2014. Evolution including transposable elements in differential expression analysis of of DNA methylation patterns in the Brassicaceae is driven by differ- RNA-seq datasets. Bioinformatics 31(22):3593–3599. ences in genome organization. PloS Genet. 10(11):e1004785. Kapitonov VV, Jurka J. 2004. Harbinger transposons and an ancient Sequeira-Mendes J, et al. 2014. The functional topography of the HARBI1 gene derived from a transposase. DNA Cell Biol. Arabidopsis genome is organized in a reduced number of linear motifs 23(5):311–324. of chromatin states. Plant Cell 26(6):2351–2366. Kharchenko PK, Tolstorukov MY, Park PJ. 2008. Design and analysis of Shi X, Zhang C, Ko DK, Chen ZJ. 2015. Genome-wide dosage-dependent ChIP-seq experiments for DNA-binding proteins. Nat Biotechnol. and -independent regulation contributes to gene expression and evo- 26(12):1351–1359. lutionary novelty in plant polyploids. Mol Biol Evol. 32(9):2351–2366. Kim MY, Zilberman D. 2014. DNA methylation as a system of plant ge- Shivaprasad PV, Dunn RM, Santos BA, Bassett A, Baulcombe DC. 2012. nomic immunity. Trends Plants Sci. 19(5):320–326. Extraordinary transgressive phenotypes of hybrid tomato are influ- Lafon-Placette C, Ko ¨ hler C. 2015. Epigenetic mechanisms of postzygotic enced by epigenetics and small silencing RNAs. EMBO J. reproductive isolation in plants. Curr Opin Plant Biol. 23:39–44. 31(2):257–266. Landt SG, et al. 2012. ChIP-seq guidelines and practices of the ENCODE Soltis PS, Soltis DE. 2009. The role of hybridization in plant speciation. and modENCODE consortia. Genome Res. 22(9):1813–1831. Annu Rev Plant Biol. 60:561–588. Langmead B, Salzberg SL. 2012. Fast gapped-read alignment with Bowtie Spoelhof JP, et al. 2017. Karyotypic variation and pollen stainability in 2. Nat Methods. 9(4):357–359. resynthesized allopolyploids Tragopogon miscellus and T. mirus.Am Li H, et al. 2009. The Sequence Alignment/Map format and SAMtools. J Bot. 104(10):1484–1492. Bioinformatics. 25(16):2078–2079. Todesco M, et al. 2016. Hybridization and extinction. Evol Appl. Lockton S, Gaut BS. 2010. The evolution of transposable elements in 9(7):892–908. natural populations of self-fertilizing Arabidopsis thaliana and its out- Wang H, et al. 2014. Rapid genetic and epigenetic alterations under inter- crossing relative Arabidopsis lyrata. BMC Evol Biol. 10:10. generic genomic shock in newly synthesized Chrysanthemum morifo- Matsui A, et al. 2008. Arabidopsis transcriptome analysis under drought, lium  Leucanthemum paludosum hybrids (Asteraceae). Genome Biol cold, high-salinity and ABA treatment conditions using a tiling array. Evol. 6(1):247–259. Plant Cell Physiol. 49(8):1135–1149. Wittkopp PJ, Haerum BK, Clark AG. 2004. Evolutionary changes in cis and McCarthy DJ, Chen Y, Smyth GK. 2012. Differential expression analysis of trans gene regulation. Nature 430(6995):85–88. multifactor RNA-Seq experiments with respect to biological variation. Wright KM, Lloyd D, Lowry DB, Macnair MR, Willis JH. 2013. Indirect Nucleic Acids Res. 40(10):4288–4297. evolution of hybrid lethality due to linkage with selected locus in McClintock B. 1984. The significance of responses of the genome to chal- Mimulus guttatus. PLoS Biol. 11(2):e1001497. lenge. Science 226(4676):792–801. Wu Y, et al. 2016. Transcriptome shock in an interspecific F1 triploid hybrid Nordberg H, et al. 2014. The genome portal of the department of energy of Oryza revealed by RNA sequencing. J Integr Plant Biol. joint genome institute: 2014 updates. Nucleic Acids Res. 42:D26–D31. 58(2):150–164. Parisod C, et al. 2009. Rapid structural and epigenetic reorganization near Xiong Z, Gaeta RT, Pires JC. 2011. Homoeologous shuffling and chromo- transposable elements in hybrid and allopolyploid genomes in some compensation maintain genome balance in resynthesized allo- Spartina. New Phytol. 184(4):1003–1015. polyploid Brassica napus. Proc Natl Acad Sci U S A. Parisod C, et al. 2010. Impact of transposable elements on the organiza- 108(19):7908–7913. tion and function of allopolyploid genomes. New Phytol. Xu C, et al. 2014. Genome-wide disruption of gene expression in allopo- 186(1):37–45. lyploids but not hybrids of rice subspecies. Mol Biol Evol. Pietzenuk B, et al. 2016. Recurrent evolution of heat-responsiveness in 31:1066–1076. Brassicaceae COPIA elements. Genome Biol. 17(1):209. Yoo MJ, Szadkowski E, Wendel JF. 2013. Homoeolog expression bias and Quadrana L, et al. 2016. The Arabidopsis thaliana mobilome and its impact expression level dominance in allopolyploid cotton. Heredity (Edinb) at the species level. eLife 5:e15716. 110(2):171–180. Rangwala SH, Richards EJ. 2010. The structure, organization and radiation Zhang D, et al. 2016. Genome-wide gene expressions respond differently of Sadhu non-long terminal repeat retroelements in Arabidopsis spe- to A-subgenome origins in Brassica napus synthetic hybrids and nat- cies. Mobile DNA 1(1):10. ural allotetraploid. Front Plant Sci. 7:1508. Rebernig CA, et al. 2015. Non-reciprocal interspecies hybridization barriers Zhu W, et al. 2017. Altered chromatin compaction and histone methyla- in the Capsella genus are established in the endosperm. PLoS Genet. tion drive non-additive gene expression in an interspecific Arabidopsis 11(6):e1005295. hybrid. Genome Biol. 18(1):157. Robinson MD, McCarthy DJ, Smyth GK. 2010. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26(1):139–140. Associate editor: Ellen Pritham Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 1415 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Genome Biology and Evolution Oxford University Press

Robustness of Transposable Element Regulation but No Genomic Shock Observed in Interspecific Arabidopsis Hybrids

Free
13 pages

Loading next page...
 
/lp/ou_press/robustness-of-transposable-element-regulation-but-no-genomic-shock-0ek7pebVa8
Publisher
Oxford University Press
Copyright
© The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution.
ISSN
1759-6653
eISSN
1759-6653
D.O.I.
10.1093/gbe/evy095
Publisher site
See Article on Publisher Site

Abstract

The merging of two divergent genomes in a hybrid is believed to trigger a “genomic shock”, disrupting gene regulation and transposable element (TE) silencing. Here, we tested this expectation by comparing the pattern of expression of transposable elements in their native and hybrid genomic context. For this, we sequenced the transcriptome of the Arabidopsis thaliana genotype Col-0, the A. lyrata genotype MN47 and their F1 hybrid. Contrary to expectations, we observe that the level of TE expression in the hybrid is strongly correlated to levels in the parental species. We detect that at most 1.1% of expressed transposable elements belonging to two specific subfamilies change their expression level upon hybridization. Most of these changes, however, are of small magnitude. We observe that the few hybrid-specific modifications in TE expression are more likely to occur when TE insertions are close to genes. In addition, changes in epigenetic histone marks H3K9me2 and H3K27me3 following hybridization do not coincide with TEs with changed expression. Finally, we further examined TE expression in parents and hybrids exposed to severe dehydration stress. Despite the major reorganization of gene and TE expression by stress, we observe that hybridization does not lead to increased disorganization of TE expression in the hybrid. Although our study did not examine TE transposition activity in hybrids, the examination of the transcriptome shows that TE expression is globally robust to hybridization. The term “genomic shock” is perhaps not appropriate to describe transcrip- tional modification in a viable hybrid merging divergent genomes. Key words: hybridization, genomic shock, transposable element, stress, Arabidopsis. Introduction Yet, a recent study of the homoploid hybrid of Arabidopsis Interspecific hybridization is an important factor in plant thaliana and Arabidopsis lyrata did not report a major disrup- evolution: While fertile allopolyploid hybrids may become tion of the transcriptome and epigenome. Moderate differ- the founders of new species (Soltis and Soltis 2009), the ences in expression of protein-coding genes were reported merging of two divergent genomes in a hybrid can also and the parental pattern of the histone mark H3K27me3 cause mild outbreeding depression or even result in com- appeared to be maintained (Zhu et al. 2017). These species, plete incompatibility (Bomblies and Weigel 2007; Todesco however, differ in both genomic structure and TE content. et al. 2016). At the molecular level, hybridization can lead A recent increase in TE number has been reported for the to a genome-wide misregulation of the transcriptome and A. lyrata genome (Hu et al. 2011). Arabidopsis thaliana TEs epigenome (Lafon-Placette and Ko¨hler 2015). The mobili- are concentrated in pericentromeric regions, rarely venturing zation of up-regulated transposable elements (TE) can fur- in gene rich chromosome arms. By contrast, in A. lyrata TEs ther result in chromosomal rearrangements. Such account for most of the 40% larger genome and are broadly phenomena, first described by Barbara McClintock in distributed, occurring in closer vicinity to expressed genes (Hu maize (McClintock 1984) and termed “genomic shock”, et al. 2011). The relatively low frequency of insertion poly- are thought to provide a postzygotic barrier against gene morphisms within species revealed evolutionary tensions on flow between species. insertions in gene rich regions, where permanent TE silencing The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non- commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com Genome Biol. Evol. 10(6):1403–1415. doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 1403 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Go ¨ bel et al. GBE can have deleterious consequences on the expression of 1-week intervals. Two micrograms of total RNA were used for neighboring genes (Hollister and Gaut 2009; Lockton and library preparation following the TruSeq Illumina RNA Sample Gaut 2010; Hollister et al. 2011). Finally, several studies Preparation v2 Guide. This includes poly-Aþ RNA selection have observed that A. lyrata TEs tended to be more highly and the use of random primers for reverse-transcription. expressed than A. thaliana alleles in hybrids (He et al. 2012). In Sequencing was performed on Illumina HiSeq2000 following view of these differences in number, distribution, and regula- the manufacturer’s protocols and paired-end 100-bp-long tion of transposable elements, we hypothesized that for reads were obtained. In total, 15–20 million paired reads for hybrids between these species, a genomic shock is more likely the parents and 30–40 million reads for the hybrid were pro- to occur at the level of TE expression. To test this hypothesis, duced (supplementary table 1, Supplementary Material on- we quantified the impact of parental versus hybrid genomic line). RNA-seq reads were filtered and trimmed for quality background on TE regulation. controlasin He et al. (2016) and mapped on the hybrid ge- We studied transposable elements in the A. thaliana (Col- nome, a concatenation of the A. thaliana Col-0 reference 0)  A. lyrata (MN47) hybrid and its respective parental lines, genome (TAIR10, www.arabidopsis.org; last accessed may using RNA-Seq expression data and ChIP-seq data of the his- 25, 2018) and the A. lyrata MN47 reference genome tone marks H3K27me3 and H3K9me2. It is difficult to estab- (Araly1, Nordberg et al. [2014]), using STAR with default lish orthology relationships for TEs, because many of them parameters (Dobin et al. 2013). duplicated after the species separation (Hu et al. 2011; de la Chaux et al. 2012). Thus, our analysis focuses on trans-acting effects experienced by TEs of each genome after hybridiza- tion. Contrary to our primary expectations, we find TE silenc- ing to be largely unaffected in the hybrid. No systematic TE and Gene Read Counting relationship could be found between a change of TE expres- sion in the hybrid and the gain or loss of histone marks. We The gene annotations TAIR10 (Col-0, Berardini et al. 2015), further exposed interspecific hybrids and their parents to se- Alyrata_107_v1.0 (MN47, Hu et al. 2011) and TE annotations vere dehydration stress and confirm that TE regulation (e.g., position and TE family memberships) for A. lyrata and A. remains robust to hybridization in stress conditions. thaliana (Pietzenuk et al. 2016) were merged to form the genome annotation file used in this analysis. Aligned reads were filtered as follows: for the MN47 genome, only the main Materials and Methods scaffolds 1–8 were considered. Read or fragment alignments Plant Materials and Growth Conditions shorter than 20 or longer than 700 base pairs (bp) were dis- carded. A minimum alignment score of 33 was required in Seeds of A. thaliana accession Col-0 were obtained from the order to accept an RNA-seq match. We allowed multiple Arabidopsis Biological Resource Center (ABRC, USA). RNA-seq read matches on TEs if they were entirely within a Arabidopsis lyrata ssp. lyrata genotype MN47 was obtained single (super)family, following Pietzenuk et al. (2016), and Jin from the lab of D. Weigel (Max Planck Institute for et al. (2015). A read was assigned to a (super)family F on Development, Tu ¨ bingen, Germany). Parental lines were genome G if its primary match was in F on G, and any sec- crossed by pollinating emasculated A. thaliana flowers with ondary matches on G were in F, too. In order to assign a A. lyrata pollen, as described in De Meaux et al. (2006).We match to a (super)family F, a minimum overlap of 50 bp to obtained approximately 40 seeds per silique, 90% of which a member of F was required. Secondary matches on the al- aborted. We did not use embryo rescue to generate more ternative genome were allowed but not counted. Each read hybrids as in Zhu et al. (2017). Reciprocal crosses using thus either contributed a single count to a (super)family or A. thaliana as pollen donor were not successful. was discarded as multiply matching. If a read matched to TE .. TE of family F, the read was defined to contribute 1 n RNA Sampling and Sequencing in Standard and Stress 1/n counts to each of these TEs. Reads from a region of over- Conditions lap a TE and a gene were not discarded, but counted with the Seeds were stratified for 3 days at 4 C, germinated on soil and TE. These summed fractional per-element counts were the grown for 4 weeks in a growth chamber at 20 C under 14 h final output from the counting procedure. They were either 1 2 light/16 C 10 h dark under dim light (100 mmol s m ). A used directly as counts per individual TE or the counts of the dehydration treatment was applied following Seki et al. members of (super)families were summed up to yield aggre- (2002). Plants were removed from the soil and dehydrated gated per-(super)family counts. We implemented the count- on paper for 0 and 3h, in the same growth chamber. The ing procedure in the R programming language. After the aerial part of the plant was sampled, flash-frozen in liquid counting step, annotated TEs without any read count were nitrogen, and RNA was extracted as previously described excluded from the analysis. Genes were counted with the (He et al. 2016). Four biological replicates were collected in same procedure as the TEs. 1404 Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Robustness of Transposable Element Regulation GBE Read Count Normalization each bin of the reference sample, the same number of TEs as contained in the bin was randomly selected from the corre- Read counts from genes and TEs were simultaneously nor- sponding distance bin of all TEs. The union of the resampled malized. Each count was first multiplied by a factor of 1,000/l, bins then formed one resampled set. The width of the bins where l is the length in base pairs of the respective annotated increased with increasing distance: up to a distance of 1 kb, element (gene, individual TE, or entire TE (super)family). This the bin width was 200 base pairs (bp), between 1 and 5 kb it removes the dependence of the read count on element length was 500 bp, and beyond 5 kb it was 1,000 bp, to avoid having or (super)family size. The total length of a TE (super)family was bins with sparse numbers of TEs. A median, 5% and 95% defined to equal the summed length of all (super)family mem- quantile of each superfamily frequency was extracted from bers considered, discounting any regions of overlap with the 10,000 random samples and compared with the observed other annotated elements. frequency of TEs with hybrid-dependent expression. Differential Expression Analysis The Generalized Linear Models (GLM) support of the R pack- Chromatin Immune-Precipitation Experiments age edgeR (Robinson et al. 2010; McCarthy et al. 2012)was used to test the significance of differences in length- Plants were germinated and grown on germination medium normalized counts between the four biological replicates of containing = Murashige and Skoog salts, 3% sucrose, and hybrid and parental transcriptomes. P values were adjusted 0.8% agar. Seeds were stratified for 5 days at 4 C, and then with the FDR correction (Benjamini and Hochberg 1995). For transferred to a growth chamber under long-day conditions visualization in scatterplots and for the stringent definition of (14 h cool white light supplemented with red light at 20 C, differential TE expression, length-normalized counts were and 10 h dark at 18 C). One week after germination, plants convertedtoCounts per Million(CPM) andscaled were transferred to soil and grown in the same chamber. (Robinson and Oshlack 2010) using the calcNormFactors Finally, leaves larger than 0.5 cm were sampled from function of the R package edgeR. This transformation mirrors 3-week-old Col-0, and 6-week-old MN47 and F1 hybrids what edgeR does internally when computing differential ex- for chromatin immunoprecipitation (ChIP). Leaves from 10 pression. It removes global count biases due to library size to 15 plants were pooled. Sampling was shifted in the hybrid differences, easing the interpretation of the scatterplots and and the A. lyrata individual, in order to sample material of allowing the set-up of a sample-independent cutoff on the similar development (3 cm rosette diameter). absolute expression value at 0 and 10 CPM. ChIP experiments were performed with the MAGnify Chromatin Immunoprecipitation System (49-2024, Thermo Fisher Scientific) according to the manufacturer’s instructions, Analysis of Distance Distribution of TEs to Neighboring with the following modifications. Plant material was fixed in Genes 20 ml Crosslinking Buffer (0.4 M Sucrose, 10 mM Tris–HCl, pH Distances of differential and nondifferential TEs to the nearest 8.0, 10 mM MgCl2, 1% formaldehyde) by applying vacuum gene or gene 5 -end were binned into classes 0–1, 1–2, 2–3, twice for 15 min. Nuclei and chromatin were purified as in He and 3 kb. (TEs overlapping a gene were not considered.) A et al. (2012). Chromatin was sonicated with a BioRuptor de- pseudocount of 1 was added to each bin count. A saturated vice (Diagenode, Lie ` ge, Belgium) for 30 times 30 s at high log-linear model of the counts was computed using the R setting with 30 s intermittent cooling in ice-water. A DNA frag- glm() function (family ¼ “poisson”), using sum-to-zero con- ment size of 300–600 bp was controlled by running an aliquot trasts. The P value for distance d in supplementary table 1, of decrosslinked and purified DNA on 1.5% agarose gels. The Supplementary Material online, is the Pr(>jzj) value of the following antibodies were used in immunoprecipitation: anti- interaction coefficient between the factors distance category rat IgG (R9255, Sigma; St. Louis, MO), as well as anti- and significance of hybrid-dependent TE expression. H3K27me3, anti-H3K9me2 and anti-H3 (07-449, 07-441, and 07-690, respectively, Millipore; Temecula, CA). To reach Random Expectations of Superfamily Frequencies enough precipitated DNA, the ChIP-DNA precipitation proto- To evaluate the enrichment of TEs changing their expression col was performed independently for eight aliquots of the in hybrid background within each superfamily, we computed collected leaf material, pooled and purified by Qiagen a random expectation. For this, we used the set of TEs show- MinElute Reaction Cleanup Kit. Finally, we obtained ChIP- ing differential expression in the hybrid background as a ref- DNA samples from Col-0, MN47, and Col-0xMN47 with anti- erence set and performed 10,000 random resamplings of the bodies targeting H3K27me3 and H3K9me2 marks. In addi- same number of TEs from all TEs with evidence of expression. tion, we generated a ChIP-DNA sample targeting Histone 3 To account for the skewed distribution of these TEs in the in a Col-0xMN47 sample to verify that both genomes could be vicinity of non-TE genes, both the reference set and the set of efficiently immuno-precipitated and assess ChIP-seq quality all TEs were binned by distance to the closest gene 5 -end. For (see supplementary Methods, Supplementary Material online). Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 1405 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Go ¨ bel et al. GBE Table 1 Numbers and Fractions of Differentially Expressed TEs at Two Different False-Discovery Rates (FDR 0.01 or0.05) and Two Different Thresholds on the Normalized Counts Per Million (CPM) (all: 0CPM; 10 CPM) FDR £ 0.01 FDR £ 0.05 All 10 CPM All 10 CPM Arabidopsis thaliana 3 þ 2(0.3%) 1 þ 1(0.1%) 5 þ 5(0.6%) 1 þ 2(0.2%) Arabidopsis lyrata 57 þ 2(0.8%) 14 þ 0(0.2%) 73 þ 6(1.1%) 18 þ 1(0.3%) NOTE.—Numbers are given as “number of up-regulations” þ “number of down-regulations.” Fractions were computed relative to all TEs with evidence of expression in standard conditions, discounting members of genome-specific families. reads were filtered with samtools (Li et al. 2009)for a map- ChIP-Seq Library Preparation and Sequencing ping quality of 3 (-q 3), which discards multimapping reads DNA was sheared with Ion Shear reagents and the Ion Xpress with two equally good mappings but retains those with a Fragment Library Kit (Thermo Fisher Scientific, Waltham, MA) single best mapping location. Duplicate reads from the am- was used for subsequent library preparation. [AQ10] After plification step during library preparation were removed with ligation of the respective Ion Xpress barcode adapters, sam- the MarkDuplicates program from Picard Tools (http://broad- ples were amplified by 17 or 23 PCR cycles. Amplified libraries institute.github.io/picard; last accessed May 25, 2018). ChIP- were quantified on Bioanalyzer DNA High-Sensitivity Chips seq read counts for TEs and non-TE genes were counted fol- (Agilent Technologies, Santa Clara, CA) and diluted to lowing the same pipeline as for RNAseq reads. A first analysis 9 pM. Emulsion PCR was performed on the Ion OneTouch of ChIP-seq count distribution confirmed that the distribution System, followed by enrichment for template positive Ion of marks recovered known chromatin domains but revealed Sphere Particles. Sequencing sample was loaded on an Ion 64 genomic regions giving aberrant ChIP signals (see supple- Proton chip and sequenced on the Ion proton sequencer mentary Methods, Supplementary Material online). These 64 (Thermo Fisher Scientific) according to the manufacturer’s regions were removed from all analyses presented in the instructions. manuscript. ChIP-Seq Read Processing, Mapping and Peak Detection Results Raw data were preprocessed using the TORRENT SUITE ver- Only a Minor Fraction of All Expressed TEs Is Differentially sion 4.0.1 (Thermo Fisher Scientific) to trim adapter sequen- Expressedinthe Hybrid ces. Using the FASTX-Toolkit (http://hannonlab.cshl.edu/ fastx_toolkit/; last accessed May 25, 2018), reads were then Our annotation contains 10863 TEs on the A. thaliana Col- trimmed from the 3 end to remove low quality bases (phred 0 genome and 30,868 elements on the A. lyrata MN47 <15). Reads shorter than 30 bases and reads of poor overall genome. After excluding genome-specific TE families, quality (more than half of the bases with phred< 20) were 10,805 and 29,150 elements remained on the Col-0 and discarded. The mapping of the reads was performed with MN47 genomes, respectively. Of these, 1,592 (15%) of the bowtie2 (Langmead and Salzberg 2012), with the end-to- A. thaliana TEs and 7,045 (24%) of the A. lyrata TEs end parameter, without allowing mismatches in the seed (- showed evidence of expression in at least one replicate. N 0), and allowing for up to two alignments per read (-k 2); We quantified TE expression in normalized Counts Per although then only the best alignment was kept. All samples Million (CPM, see Materials and Methods). TEs contributed were mapped to a hybrid reference genome obtained by but a minor fraction of all reads (on average, only 0.07% in concatenating the genomes of the parent as described above A. thaliana samples and 0.24% in A. lyrata samples). This for RNAseq read mapping. Cross-hybridization between fraction remained comparable in the hybrid transcriptome genomes was negligible (0.5–1% of mapped reads). A sum- (on average 0.06% and 0.30% of A. thaliana and A. lyrata mary of ChIP seq read data are given in supplementary table reads, respectively in the hybrid). Only a small fraction 2, Supplementary Material online. (maximally 2%) of the expressed TEs of both genomes Quality assessment of the assays was performed using were significantly differentially expressed between paren- phantompeaktools (Kharchenko et al. 2008; Landt et al. tal and hybrid backgrounds, using different cutoffs to de- 2012) to calculate NSC and RSC metrics (normalized and rel- fine significance (table 1). Moreover, this percentage ative strand coefficients, respectively), and the bioconductor dropped below 0.5% when the median expression was ChIPQC package (Carroll et al. 2014) to calculate the square required to be above ten CPM in at least one of the two sum of deviation. NSC was low but RSC metrics were well backgrounds (table 1). These observations show that TE within the recommended range (Landt et al. 2012). Mapped transcription is not severely disturbed in the hybrid. 1406 Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Robustness of Transposable Element Regulation GBE FIG.1.—Correlation of TE expression levels, estimated as the median of normalized counts per million over the replicate plants grown in standard conditions, in parent and hybrid for individual TEs, TE families, and superfamilies. TEs or (super)families with significant expression change are marked in red. In A. lyrata, the 59 (0.8%) TEs that respond significantly to Specific Rather than Global Effects Act on TE Expression in hybridization at FDR 0.01 were generally induced. In A. the Hybrid thaliana, the 5 (0.3%) affected elements were either induced In order to investigate the effect of hybridization at the levels or repressed (supplementary fig. S1, Supplementary Material of individual elements, families and superfamilies of TEs, we online). SINEs dominated by Sadhu elements (Rangwala and correlated expression of individual TEs and combined expres- Richards 2010) and Harbingers dominated by the ATIS112A sion levels of TE families between parents and hybrids (fig. 1). family (Kapitonov and Jurka 2004) stood out as enriched Although TE transcription in the two genetic backgrounds among the responding A. lyrata TEs (fig. 2). These superfami- correlated well at all levels, the correlation increased from lies were significantly overrepresented relative to the set of the individual TEs to the family and superfamily levels. expressed TEs (1 RNA-seq count, permutation test, Spearman’s rank correlation of the median expression of rep- FDR< 0.05). The fact that only a small number of TE families licates in the two backgrounds under control conditions was showed enrichment for differential expressed elements indi- 0.73, 0.89, 0.97 on the three levels of individual TE, family, cates that hybridization does not have a global effect on TE and superfamily, respectively in A. thaliana, and 0.65, 0.86, expression. 0.91, respectively, in A. lyrata. Including four replicates in the RNA-seq analysis was necessary to establish these relation- ships with increased confidence. Using subsets of only two Rare TE Expression Changes Coincide with Alterations in replicates, for example, reduced the parent–hybrid correlation Nearby Gene Expression coefficient of TE expression to 0.53, 0.77, 0.87 on the three levels for A. thaliana, and to 0.51, 0.84, 0.89, respectively, for TEs can change expression of nearby genes by influencing A. lyrata. local chromatin accessibility. We therefore asked whether The small subset of TEs with significantly different expres- hybrid-specific expression changes of TEs predict a corre- sion level between parent and hybrid was restricted to lower sponding change in neighboring genes. level units (individual TEs and families), indicating that the Compared with TEs with expression not affected by hybrid- differences are case-specific rather than caused by a global ization, A. lyrata TEs expressed differentially in the hybrid were disruption of TE silencing. depleted in genome regions distant more than 3kb from Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 1407 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Go ¨ bel et al. GBE FIG.2.—Observed and expected superfamily distribution of differentially expressed TEs (parent vs. hybrid). Colored bars show the observed numbers of TEs showing hybrid-dependent expression in each TE superfamily. Each color corresponds to one superfamily as indicated in the legend. Superfamily names are as given in the annotation of (Pietzenuk et al. 2016). Names with trailing question marks were annotated as deviant/uncertain members of the respective superfamily. Transparent bars indicate expected median count for random sampling within each superfamily. Error bars indicate the 5% and 95% quantile of the 10,000 random draws. Single black stars mark superfamilies whose observed count is significantly higher than expected from random expectations. Gray stars indicate depletion (maximally 5% of the universe have counts the observed count of the respective superfamily). Double stars indicate significant enrichment or depletion that is robust to FDR correction. Table 2 Distance Distributions to Nearby Genes of Arabidopsis lyrata TEs with Differential Expression in Parent and Hybrid Under Control Conditions Distance to Closest Gene Closest 5 End of a Gene <1 kb 1–2 kb 2–3 kb 3kb <1 kb 1–2 kb 2–3 kb 3kb Number of TEs Diff. exp. PvsH 15 18 13 13 7 24 12 16 Not Diff. Expr. 1,238 1,675 913 3,160 654 1,455 1,078 3,799 P (frequency in distance bin) 0.239 0.539 0.061 4.0e-4 0.554 0.01 0.517 1.5e-4 NOTE.—The table combines distance information of TEs with a significant change in expression in the hybrid to either the nearest gene (any orientation) or the nearest 5 of any gene. Numbers of significant and nonsignificant TEs are listed for each distance bin (FDR 0.01). In the lower part, P values are given for the contribution of each distance bin to the difference between the distributions of significant and nonsignificant TEs (see Materials and Methods). Bins with significant depletion or enrichment in significant TEs are highlighted in blue and orange, respectively. 0 4 the closest gene 5 ends (log-linear model, P 5  10 , preferential location of hybridization-sensitive TEs in euchro- table 2). In addition, counts of hybridization-sensitive TEs matic, gene-rich regions in A. lyrata.In A. thaliana, there were were enriched in the bins 1–2 kb upstream of gene 5 -ends too few differentially expressed TEs in the hybrid background (log-linear model, P 0.01, table 2). This indicates a to reliably compare the distance distributions. 1408 Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Robustness of Transposable Element Regulation GBE We further analyzed if changed expression of TEs and their neighboring genes in a hybrid background were correlated. Hybrid-specific changes in transcription levels of non-TE genes have already been reported (Zhu et al. 2017). We also quan- tified non-TE gene expression levels in parents and hybrids. Note that in our experiment, the timing of sampling accounted for the developmental differences between the fast developer A. thaliana and both A. lyrata and the hybrid, which are late flowering. Markedly fewer genes were af- fected by hybridization in our experiment (234 in A. thaliana and 395 in A. lyrata), yet, we observed a significant overlap in genes reported to change expression in the hybrid back- ground in Zhu et al. (2017) (Hypergeometric test of overlap in genes changing their expression level between hybrid and parents, P< 10 for all comparisons). In A. lyrata 153 TEs with detectable expression were within 2 kb of a gene with hybridization-sensitive expression (FDR 0.01). Among these, 8 (5%) TEs displayed a change in expression in the hybrid background. This percentage dropped to 1% (71/6,892) for TEs with a gene neighbor showing no significant change in expression in hybrid versus parent (P< 5  10 ,hypergeo- metric test). In the A. thaliana genome, the number of expressed TEs located in the vicinity of genes was too small to calculate correlations. We conclude that hybrid-specific TE expression can coincide with a modification of the expression of nearby protein-coding genes, but this only affects a very small number of expressed loci. Histone Marks Are Not Modified in the Hybrid Background Two epigenetic marks are associated with transcriptional silencing in plants: the plastic mark H3K27me3 and the permanent mark H3K9me2 (Kim and Zilberman 2014; Sequeira-Mendes et al. 2014). We profiled these two histone modifications by Chromatin Immuno-Precipitation followed by sequencing (ChIP-Seq) in hybrids and A. thaliana and A. lyrata parents. Although the ChIP-seq H3K27me3 yielded a comparatively lower number of reads, especially for Col-0 (supplementary tables 2, 5,and 6, Supplementary Material online), the distribution of marks on TEs was in agreement with previous reports (fig. 3A, Seymour et al. 2014) and cap- tured well the chromatin domains that were defined in A. thaliana (Sequeira-Mendes et al. 2014, supplementary Methods, Supplementary Material online). In our data, H3K9me2 marks showed a high correlation of normalized read counts between parental and hybrid genomes, both genome-wide (calculated across 10-kb genomic windows) and for individual TEs, TE families, and superfamilies. The FIG.3—Correlation of histone mark levels on TEs in parent and hybrid. genome-wide Spearman correlation coefficients of (A) Genome-wide distribution of H3K9me2 and H3K27me2 on the pa- H3K9me2 read count along the genomes were 0.66 and rental and hybrid genomes. The circular plots represent the hybrid genome used as reference, with the five A. thaliana chromosomes in green and the eight A. lyrata chromosomes in red. Inner tracks represent gene, TE, FIG.3—Continued H3K9me2 and H3K27me3 densities in 500 kb bins. Top left: Col-0, Top counts. The color code indicates superfamily membership, using the right: MN47, Bottom: hybrid. (B and C) Scatterplots of ChIP-seq read palette of figure 2. (B) H3K9me2; (C) H3K27me3. Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 1409 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Go ¨ bel et al. GBE FIG.4—No change in histone mark occupancy is detected for differentially expressed TEs. Shown are scatterplots of each histone marks on each of the two genomes for TEs with significantly different expression level in parent versus hybrid under standard conditions. Superfamily membership of the TEs is indicated using the palette of figure 2. Symbol shapes indicate the direction of gene expression change in hybrid versus parent (upward triangle: up, downward: down). 0.92 for A. thaliana or A. lyrata, respectively. Spearman’s rank Nevertheless, the correlation of H3K27me3 occupancy in correlation coefficients of H3K9me2 for individual TEs, TE- parent and hybrid remained high even for TEs, in A. lyrata families, and superfamilies were 0.76, 0.94, and 0.96 for A. (Spearman’s rank correlation 0.58, 0.78, 0.92 for individual thaliana and 0.87, 0.97, and 0.99 for A. lyrata (fig. 3B and C). TEs, families, and superfamilies). For the A. thaliana/hybrid Although the high correlation values pointed to a largely ad- comparison, the correlation dropped to 0.2, 0.27, 0.21 for ditive inheritance of the mark in the hybrid on all levels of the individual TEs, families, and superfamilies, respectively annotation, the fact that correlation was lowest on the level (fig. 3C). Note however that the Col-0 H3K27me3 ChIP- of individual elements again underscores the role of element- seq experiments yielded a comparatively lower number of specific responses in the hybrid. reads (supplementary table 2, Supplementary Material While the function of H3K9me2 is to permanently silence online). TEs, H3K27me3 mainly serves to reversibly silence genes We tested then whether changes in histone mark occu- whose expression is restricted to specific conditions or de- pancy could explain TE expression differences between velopmental time windows (Gan et al. 2015). The distribu- parents and hybrids. TEs with a modified expression level in tion of H3K27me3 marks was strongly correlated genome- hybrid background generally showed similar H3K27me3 level wide between parental and hybrid backgrounds (Spearman in parent and hybrid (fig. 4). When a departure in epigenetic correlation coefficients: 0.63 and 0.86 for A. thaliana and A. levels was detected between parents and hybrids, the direc- lyrata, respectively, fig. 3A). H3K27me3 is not a typical si- tion of the change in expression of the TE was not necessarily lencing mark of TEs (Sequeira-Mendes et al. 2014). This was in the expected direction (fig. 4). A. lyrata TEs with a modified most apparent in our data for the LTR superfamily, which expression level in hybrid background also generally showed tended to be associated with the highest level of H3K9me2 similar H3K9me2 mark levels in parents and hybrids. We and the lowest level of H3K27me3 (fig. 3B and C). detected a concordant loss of H3K9me2 marks and increased 1410 Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Robustness of Transposable Element Regulation GBE Table 3 Numbers and Fractions of Differentially Expressed TEs in Stress Conditions (3h) and in Response to Stress (3h specific), at Two Different False-Discovery Rates (FDR 0.01 or 0.05) and Two Different Thresholds on the Normalized Counts Per Million (CPM) (all: 0CPM; 10 CPM) FDR £ 0.01 FDR £ 0.05 All 10 CPM All 10 CPM Arabidopsis thaliana 3h 4 þ 14 (1.1%) 2 þ 4(0.4%) 7 þ 26 (2.1%) 2 þ 5(0.4%) 3h specific 3 þ 12 (0.9%) 2 þ 3(0.3%) 4 þ 23 (1.7%) 2 þ 4(0.4%) Arabidopsis lyrata 3h 33 þ 18 (0.7%) 4 þ 4(0.1%) 54 þ 37 (1.3%) 10 þ 15 (0.4%) 3h specific 9 þ 17 (0.4%) 3 þ 4(0.1%) 12 þ 34 (0.7%) 7 þ 14 (0.3%) NOTE.—Fractions were computed relative to all TEs with evidence of expression under stress conditions, discounting members of genome-specific families. expression in the hybrid for 6 A. lyrata TEs (fig. 4). But there conditions, the correlation of TE expression levels between again, we observed for 2 TEs that the change in H3K9me2 parents and hybrids increased at the family and superfamily occupancy did not reflect the change in expression. levels. The Spearman’s rank correlation of the replicates’ me- Altogether, in our data, a differential deposition of dian expression in the two backgrounds was (0.7, 0.81, 0.83) H3K27me3 and H3K9me2 histone methylation marks in the for A. thaliana, and (0.67, 0.91, 0.94) for A. lyrata, for indi- hybrid is rare and generally does not coincide with hybrid- vidual, family and superfamily levels, respectively. dependent changes in TE expression. Interestingly, we observed that more TEs were down- regulated in the hybrids compared with the parents under stress conditions than under standard conditions (fig. 5A and table 3). Only the Harbinger and LINE_L1 superfamiles Robustness of TE Expression to Hybridization Is Maintained were enriched among expressed A. lyrata TEsinstresscon- under Stress Conditions ditions (fig. 5B; permutation test using TEs with 1 RNA-seq It is assumed that transposition activity is induced by stress count as reference, FDR< 0.05). This further corroborates the conditions (McClintock 1984; Cavrak et al. 2014). conclusion that in our two species, TE regulation is robust to Dehydration stress is associated with a vast transcriptome hybridization even under stress conditions. reprogramming (Matsui et al. 2008). We analyzed the depen- dence of transcript abundance upon interspecific hybridiza- tion for TEs and non-TE genes under strong drought stress Discussion conditions in the same way as described for control condi- No evidence for a genome shock upon hybridization of tions. Having found that TEs with a hybrid-specific expression A. lyrata and A. thaliana change tend to be close to genes, and given the reprogram- The term “genomic shock,” was initially used to describe the ming of the gene transcriptome by the stress, we asked breakage and large scale rearrangement of chromosomes in whether TE regulation was robust to hybridization even under maize (McClintock 1946, 1984). This term is thus appropriate stress conditions. In our data, 30–40% of the genes of both to describe genome instability following interspecific hybridi- species changed expression after 3h of dehydration, with sim- zation, that is, unbalanced segregation of homoelogous chro- ilar proportions in the parental and hybrid backgrounds. The mosomes, chromosome loss or translocations (Xiong et al. fraction of stress responsive TEs was much smaller: 3–4% on 2011; Chester et al. 2012; Henry et al. 2014; Spoelhof the A. thaliana genome and 2–3% for A. lyrata (not shown). et al. 2017). Over the years, the term “genome shock” has In addition, we observed no increase in the variance of TE also been used to describe the anticipated alteration in ge- expression in hybrid context (supplementary fig. 3, nome activity immediately following interspecific hybridiza- Supplementary Material online). Although the fraction of tion or allopolyploidization, although transposon A. thaliana TEs with a significant expression difference be- reactivation has seldom been rigorously documented imme- tween parent and hybrid is higher under dehydration stress diately after genome merging (Parisod et al. 2010). In plants, than under control conditions, both the absolute numbers (4 DNA cytosine methylation (mC) and the histone mark up and 14 down) and the fraction (1.1% of all TEs with read H3K9me2 embed TEs in a chromatin context unfavorable count1) remained small (table 3). The fraction of respond- for transcription (Kim and Zilberman 2014). The establishment ing A. lyrata TEs decreased slightly compared with control and maintenance of these marks is not only inherited through conditions, and again the numbers remained small (33 up, DNA replication, it is also guided by small interfering RNAs 18 down, fraction 0.7%). The percentage dropped below (siRNAs) (see Kim and Zilberman 2014 for a review of plant TE 0.5% if an expression level of10 CPM in at least one of silencing mechanisms). Because of their potential to exert the backgrounds was required (table 3). Like under control Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 1411 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Go ¨ bel et al. GBE genome-wide trans effects via base pair homology of the short 24 nt siRNAs, the population of siRNAs expressed by a genome must be adapted to its gene content to ensure proper gene expression. Interspecific hybridization, therefore, brings together trans-regulators of TE expression that are not mutually adapted and could deeply disorganize gene regula- tion. Gene expression modifications following interspecific hy- bridization or allopolyploidization have been monitored across plant genera as diverse as Gossypium, Senecio, Brassica, Tragopogon, or Oryza (e.g., Hegarty et al. 2006; Buggs et al. 2009;Yoo et al. 2013). Many studies have concluded that epigenome incompatibilities are important contributors to hybrid sterility (Parisod et al. 2009; Ghani et al. 2014; Wang et al. 2014; Xu et al. 2014; Senerchia et al. 2015; Wu et al. 2016). In Aegilops hybrids, for example, parental divergence in TE content was associated with the emergence of new TE copies in hybrids (Senerchia et al. 2015). Allopolyploidy is also believed to be associated with drastic epigenetic reorganiza- tion although evidence for a significant reactivation of TE transposition remains scant (Parisod et al. 2010). Our analysis contradicts the conclusions of these previous reports. Comparing the TE expression and epigenetic profile of the A. lyrata and A. thaliana genomes in parental and hy- brid context, we find no evidence for a major transcriptomic or epigenomic “shock” driven by inappropriate TE silencing and epigenetic mechanisms. Both TE transcription and the decoration of TEs by the silencing histone mark H3K9me2 are globally unaffected in the hybrid. We observe that only a handful of specific TEs show a change in expression in hybrids. In addition, the changes in TE expression are generally of small magnitude. We further show that the distribution of silencing epigenetic marks is strongly correlated in parents and hybrids, and its variation does not match the changes in TE expression. Finally, under severe dehydration stress, a FIG.5.—(A) Correlation of TE expression levels in stress conditions, large number of genes change their transcriptional activity estimated as Median normalized counts per million over the replicates, in (Matsui et al. 2008). This broad genome-wide remodeling parent and hybrid for individual TEs, TE families, and superfamilies. TEs or of transcriptional activity, however, far from inducing a sev- (super)families with significant expression change are marked in red (sig- ered genome shock, tends to increase the robustness of TE nificant only under control conditions), orange (significant both under expression patterns to hybridization. These findings thus sug- standard and stress conditions), or yellow (significant only in stress con- gest that stress does not increase the impact of hybridization, ditions). (B) Observed and expected superfamily distribution of TEs whose further supporting the idea that the control of TE elements is expression is hybrid-dependent specifically in stress conditions. Colored robust to genome merging. Although these observations par- bars show the observed numbers of TEs showing hybrid-dependent ex- tially contradict previous reports, they are not isolated. In re- pression in each TE superfamily. Each color corresponds to one superfamily cent years, an increasing number of studies, mostly as indicated in the legend. Superfamily names are as given in the annota- conducted in the Brassicaceae family, has begun to report tion of Pietzenuk et al. (2016). Names with trailing question marks were annotated as deviant/uncertain members of the respective superfamily. cases where hybridization does not massively alter gene reg- Transparent bars indicate expected median count for random sampling ulation (Akama et al. 2014; Shi et al. 2015; He et al. 2016; within each superfamily. Error bars indicate the 5% and 95% quantile of Zhang et al. 2016; Zhu et al. 2017). In addition, in tomato, the 10,000 random draws. Single black stars mark superfamilies whose transgressive expression of small RNAs was associated with observed count is significantly higher than expected from random expect- the dysregulation of stress responsive genes, but not with that ations. Gray stars indicate depletion (maximally 5% of the universe have of TEs (Shivaprasad et al. 2012). counts the observed count of the respective superfamily). Double stars The contradiction between these studies and previous indicate significant enrichment or depletion that is robust to FDR reports may have various explanations. First, it is striking to correction. note that many of the studies that conclude on the absence of 1412 Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Robustness of Transposable Element Regulation GBE a global gene expression dysregulation were conducted in the A. thaliana allele of FLC, which is up-regulated by species with well-known genomes, where gene expression FRIGIDA, is one of the alleles most strongly impacted by hy- levels can be quantified with the best accuracy. They were bridization (supplementary table 3, Supplementary Material not conducted in polyploids with large and imperfectly assem- online; Zhu et al. 2017). Interspecific differences in trans-act- bled genomes. Second, biological and experimental variance ing factor activity will therefore determine the extent to which between replicates affects the strength of the correlation of gene expression changes upon hybridization. If F1 hybrids are expression in parental versus hybrid backgrounds, and thus fertile, these changes can segregate in subsequent F2 gener- the apparent global impact of hybridization. Finally, it is nec- ations as in, for example, Shivaprasad et al. (2012). essary to contrast diverse sources of information to test Here, focusing on the trans-acting effect of the passage whether hybridization has a broad impact on gene regulation. from a native parental background to a hybrid background on In the allotetraploid A. suecica, modification in siRNA levels TE expression, we show that at most 1–2% of expressed TEs was not associated with nonadditive gene expression in newly have expression that is differentially controlled in trans. formed allopolyploids, suggesting that changes in siRNA ex- Several elements indicate that these trans-acting effects may pression were not the result of epigenetic disorganization in be associated with interspecific differences in mechanisms of the hybrid (Ha et al. 2009). Similarly, in Arabidopsis homo- gene regulation in euchromatic regions. When combined in a ploid hybrids, we observed that hybrid-specific TE expression hybrid, these trans-acting differences modify the regulatory changes do not coincide with obvious epigenetic changes, a environment of a few TE subfamilies. Indeed, the small pro- pattern also reported for protein-coding genes (Zhu et al. portion of TEs with differential expression in hybrids tend to 2017). Taken together, the different molecular data sets be located in the vicinity of genes. Interspecific differences in used in this study show that TE expression and the distribution the regulation of euchromatic regions are however likely to be of epigenetic marks, which control gene expression, are not minor, given the small number of elements affected. massively disorganized in hybrids. Interspecific differences in gene regulation controlled in cis can also create a pattern of genome dominance in hybrids, with a majority of orthologous transcripts showing a transcrip- trans-and cis-Acting Differences between Species Can tion bias towards one parental genome (Yoo et al. 2013; Manifest in F1 Hybrids Cheng et al. 2016; Zhang et al. 2016). For example, in The comparative analysis of gene expression in hybrids and A. thaliana  A. lyrata hybrids, a bias towards preferential their parents reflects the variability of trans-and cis-regulatory expression of the A. lyrata allele has been reported by inde- mechanisms controlling gene expression in the two species pendent studies (He et al. 2012, 2016; Zhu et al. 2017). (Wittkopp et al. 2004). After hybridization, changes in total Cis-effects have not been examined in this study, because expression will simply reflect the number and magnitude of orthologous relationships cannot be established for most trans-regulatory differences between species, whereas the TEs, but for the subset of TEs that are located in orthologous specific cis-regulatory variant of each transcript will result in positions in both species, the A. lyrata TE copy tended to be allelic expression differences. These differences will manifest more highly expressed than the A. thaliana copy in the hybrid in hybrids on a restricted number of genes or elements, which (He et al. 2012). does not imply that there is a major “genomic shock”. Only trans-acting differences will have nonadditive conse- Incompatibilities Can Also Have an Oligogenic Basis quences on expression. In fungi and cotton allopolyploid, >50% of homeologous genes inherit the expression pattern It may seem at first surprising that TE silencing shows so little of the parents (Yoo et al. 2013; Cox et al. 2014). In rice, allelic disruption in interspecific A. thaliana  A. lyrata hybrids, given expression bias in hybrids between Oryza sativa and O. japon- the estimated 6 My of divergence and their vastly different TE ica also reflects mostly interspecific differences. In hybrids be- content (Hu et al. 2011; Hohmann et al. 2015). The proba- tween the diploid O. sativa and the tetraploid O. punctata, bility to observe specific disruptions in gene expression in a only 16% of the genes showed an expression pattern differ- hybrid background, however, does not necessarily scale with entfrom thatofthe parents (Wu et al. 2016). The analysis of evolutionary distance. The proportion of genes with hybrid- protein-coding gene expression in A. thaliana  A.lyrata specific changes in expression is not greater for hybrids be- hybrids reported a larger proportion of nonadditive gene ex- tween the closely related subspecies O. sativa and O. japonica pression changes for the A. thaliana genome compared with than between O. sativa and the more distant species O. punc- the A. lyrata genome (Zhu et al. 2017). This pattern, like in our tata (Wu et al. 2016). A simple oligogenic combination of study, could not be associated with chromatin modifications alleles, in fact, can cause what is known as a Dobzhansky– induced by hybridization or any other genome-wide change Muller incompatibility, even in the absence of large in gene regulatory mechanisms (Zhu et al. 2017). However, evolutionary distances between parents. Dobzhansky–Muller the A. thaliana parental genotype Col-0 is missing an active incompatibilities, which can cause major disruption at the allele of the major developmental regulator FRIGIDA.In fact, phenotypic and, sometimes at the transcriptome level, are Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 1413 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Go ¨ bel et al. GBE Berardini TZ, et al. 2015. The Arabidopsis Information Resource: making not seldom within species, and were even detected within and mining the “gold standard” annotated reference plant genome. populations (Alcazar et al. 2010; Bomblies and Weigel Genesis 53(8):474–485. 2010; Durand et al. 2012; Wright et al. 2013; Chae et al. Bomblies K, Weigel D. 2007. Hybrid necrosis: autoimmunity as a potential 2014). We conclude that the term “shock” is not appropriate gene-flow barrier in plant species. Nat. Rev. Genet. 8(5):382. to describe the transcriptional consequences of the merging Bomblies K, Weigel D. 2010. Arabidopsis and relatives as models for the study of genetic and genomic incompatibilities. Philos Trans R Soc B of divergent genomes in a viable F1 hybrid. Our data rather Biol Sci. 365(1547):1815–1823. supports a less dramatic scenario, where the merging of two Buggs RJ, et al. 2009. Gene loss and silencing in Tragopogon miscellus genomes creates a new trans-acting background, in which (Asteraceae): comparison of natural and synthetic allotetraploids. oligogenic Dobzhansky–Muller incompatibilities can manifest. Heredity 103(1):73–81. Such incompatibilities are known to affect endosperm func- Carroll TS, Liang Z, Salama R, Stark R, de Santiago I. 2014. Impact of artifact removal on ChIP quality metrics in ChIP-seq and ChIP-exo tion (Rebernig et al. 2015) and will also occasionally affect the data. Front Genet. 5:75. regulation of specific TEs and favor new bouts of transposi- Cavrak VV, et al. 2014. How a retrotransposon exploits the plant’s heat tion. Without a global robustness of TE regulation to genome stress response for its activation. PLoS Genet. 10(1):e1004115. merging, hybridization and polyploidization would probably Chae E, et al. 2014. Species-wide genetic incompatibility analysis identifies be much less important for plant diversification (Alix et al. immune genes as hot spots of deleterious epistasis. Cell 159(6):1341–1351. 2017). We note, however, that the impact of hybridization Cheng F, et al. 2016. Epigenetic regulation of subgenome dominance on TE transposition rate has not been examined in this study. following whole genome triplication in Brassica rapa. New Phytol. It thus remains possible that TEs with undetectable expression 211(1):288–299. level transpose at a rate that is higher in F1 interspecific Chester M, et al. 2012. Extensive chromosomal variation in a recently hybrids than within species. Within A. thaliana,elements of formed natural allopolyploid species, Tragopogon miscellus (Asteraceae). Proc Natl Acad Sci U S A. 109(4):1176–1181. about half of the TE families were shown to have recently Cox MP, et al. 2014. An interspecific fungal hybrid reveals cross-kingdom transposed (Quadrana et al. 2016). rules for allopolyploid gene expression patterns. PLoS Genet. 10(3):e1004180. de la Chaux N, Tsuchimatsu T, Shimizu KK, Wagner A. 2012. The pre- Supplementary Material dominantly selfing plant Arabidopsis thaliana experienced a recent Supplementary data areavailableat Genome Biology and reduction in transposable element abundance compared to its out- crossing relative Arabidopsis lyrata. Mobile DNA 3(1):2. Evolution online. De Meaux J, Pop A, Mitchell-Olds T. 2006. Cis-regulatory evolution of Chalcone-Synthase expression in the genus Arabidopsis.Genetics 174(4):2181–2202. Acknowledgments Dobin A, et al. 2013. STAR: ultrafast universal RNA-seq aligner. We thank Franziska Turck (Max Planck Institute for Plant Bioinformatics. 29(1):15–21. Durand S, Bouch e N, Perez Strand E, Loudet O, Camilleri C. 2012. Rapid Breeding Research) for sharing her ChIP-seq protocol and ex- establishment of genetic incompatibility through natural epigenetic perience. This research was supported by the Deutsche variation. Curr Biol. 22(4):326–331. Forschung Gesellschaft (DFG) with grant ME2742/2-1 and Gan ES, Xu Y, Ito T. 2015. Dynamics of H3K27me3 methylation and de- ME2742/6-1 in the realm of the special study program methylation in plant development. Plant Signal Behav. SPP1529, as well as by the European Research Council with 10(9):e1027851. Grant 648617 “AdaptoSCOPE.” The raw data has been de- Ghani MA, et al. 2014. The role of small RNAs in wide hybridisation and allopolyploidisation between Brassica rapa and Brassica nigra.BMC posited in NCBI Sequence Read Archive (SRA) and are acces- Plant Biol. 14:272. sible through SRA study number SRP148881 and SRP148853 Ha M, et al. 2009. Small RNAs serve as a genetic buffer against genomic for RNA-seq and ChIP-seq data, respectively. shock in Arabidopsis interspecific hybrids and allopolyploids. Proc Natl Acad Sci U S A. 106(42):17835–17840. He F, et al. 2012. Widespread interspecific divergence in cis-regulation of Literature Cited transposable elements in the Arabidopsis genus. Mol Biol Evol. Akama S, Shimizu-Inatsugi R, Shimizu KK, Sese J. 2014. Genome-wide 29(3):1081–1091. quantification of homeolog expression ratio revealed nonstochastic He F, et al. 2016. The footprint of polygenic adaptation on stress- gene regulation in synthetic allopolyploid Arabidopsis. Nucleic Acids responsive cis-regulatory divergence in the Arabidopsis genus. Mol Res. 42(6):e46. Biol Evol. 33(8):2088–2101. Alc azar R, et al. 2010. Natural variation at Strubbelig Receptor Kinase 3 Hegarty MJ, et al. 2006. Transcriptome shock after interspecific hybridiza- drives immune-triggered incompatibilities between Arabidopsis thali- tion in senecio is ameliorated by genome duplication. Curr Biol. ana accessions. Nat Genet. 42(12):1135–1139. 16(16):1652–1659. Alix K, et al. 2017. Polyploidy and interspecific hybridization: partners for Henry IM, et al. 2014. The BOY NAMED SUE quantitative trait locus con- adaptation, speciation and evolution in plants. Ann Bot. fers increased meiotic stability to an adapted natural allopolyploid of 120(2):183–194. Arabidopsis. Plant Cell 26(1):181–194. Benjamini Y, Hochberg Y. 1995. Controlling the false discovery rate: a Hohmann N, Wolf EM, Lysak MA, Koch MA. 2015. A time-calibrated road practical and powerful approach to multiple testing. J Roy Stat Soc map of Brassicaceae species radiation and evolutionary history. Plant B. 57(1):289–300. Cell 27(10):2770–2784. 1414 Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Robustness of Transposable Element Regulation GBE Hollister JD, Gaut BS. 2009. Epigenetic silencing of transposable elements: Robinson MD, Oshlack A. 2010. A scaling normalization method for dif- a trade-off between reduced transposition and deleterious effects on ferential expression analysis of RNA-seq data. Genome Biol. 11:R25. neighboring gene expression. Genome Res. 19(8):1419–1428. Senerchia N, Felber F, Parisod C. 2015. Genome reorganization in F1 Hollister JD, et al. 2011. Transposable elements and small RNAs contribute hybrids uncovers the role of retrotransposons in reproductive isolation. to gene expression divergence between Arabidopsis thaliana and Proc Biol Sci. 282(1804):20142874. Arabidopsis lyrata. Proc Natl Acad Sci U S A. 108(6):2322–2327. Seki M, et al. 2002. Monitoring the expression profiles of 7000 Arabidopsis Hu TT, et al. 2011. The Arabidopsis lyrata genome sequence and the basis genes under drought, cold and high-salinity stresses using a full-length of rapid genome size change. Nat Genet. 43(5):476–481. cDNA microarray. Plant J. 31(3):279–292. Jin Y, Tam OH, Paniagua E, Hammell M. 2015. TEtranscripts: a package for Seymour DK, Koenig D, Hagmann J, Becker C, Weigel D. 2014. Evolution including transposable elements in differential expression analysis of of DNA methylation patterns in the Brassicaceae is driven by differ- RNA-seq datasets. Bioinformatics 31(22):3593–3599. ences in genome organization. PloS Genet. 10(11):e1004785. Kapitonov VV, Jurka J. 2004. Harbinger transposons and an ancient Sequeira-Mendes J, et al. 2014. The functional topography of the HARBI1 gene derived from a transposase. DNA Cell Biol. Arabidopsis genome is organized in a reduced number of linear motifs 23(5):311–324. of chromatin states. Plant Cell 26(6):2351–2366. Kharchenko PK, Tolstorukov MY, Park PJ. 2008. Design and analysis of Shi X, Zhang C, Ko DK, Chen ZJ. 2015. Genome-wide dosage-dependent ChIP-seq experiments for DNA-binding proteins. Nat Biotechnol. and -independent regulation contributes to gene expression and evo- 26(12):1351–1359. lutionary novelty in plant polyploids. Mol Biol Evol. 32(9):2351–2366. Kim MY, Zilberman D. 2014. DNA methylation as a system of plant ge- Shivaprasad PV, Dunn RM, Santos BA, Bassett A, Baulcombe DC. 2012. nomic immunity. Trends Plants Sci. 19(5):320–326. Extraordinary transgressive phenotypes of hybrid tomato are influ- Lafon-Placette C, Ko ¨ hler C. 2015. Epigenetic mechanisms of postzygotic enced by epigenetics and small silencing RNAs. EMBO J. reproductive isolation in plants. Curr Opin Plant Biol. 23:39–44. 31(2):257–266. Landt SG, et al. 2012. ChIP-seq guidelines and practices of the ENCODE Soltis PS, Soltis DE. 2009. The role of hybridization in plant speciation. and modENCODE consortia. Genome Res. 22(9):1813–1831. Annu Rev Plant Biol. 60:561–588. Langmead B, Salzberg SL. 2012. Fast gapped-read alignment with Bowtie Spoelhof JP, et al. 2017. Karyotypic variation and pollen stainability in 2. Nat Methods. 9(4):357–359. resynthesized allopolyploids Tragopogon miscellus and T. mirus.Am Li H, et al. 2009. The Sequence Alignment/Map format and SAMtools. J Bot. 104(10):1484–1492. Bioinformatics. 25(16):2078–2079. Todesco M, et al. 2016. Hybridization and extinction. Evol Appl. Lockton S, Gaut BS. 2010. The evolution of transposable elements in 9(7):892–908. natural populations of self-fertilizing Arabidopsis thaliana and its out- Wang H, et al. 2014. Rapid genetic and epigenetic alterations under inter- crossing relative Arabidopsis lyrata. BMC Evol Biol. 10:10. generic genomic shock in newly synthesized Chrysanthemum morifo- Matsui A, et al. 2008. Arabidopsis transcriptome analysis under drought, lium  Leucanthemum paludosum hybrids (Asteraceae). Genome Biol cold, high-salinity and ABA treatment conditions using a tiling array. Evol. 6(1):247–259. Plant Cell Physiol. 49(8):1135–1149. Wittkopp PJ, Haerum BK, Clark AG. 2004. Evolutionary changes in cis and McCarthy DJ, Chen Y, Smyth GK. 2012. Differential expression analysis of trans gene regulation. Nature 430(6995):85–88. multifactor RNA-Seq experiments with respect to biological variation. Wright KM, Lloyd D, Lowry DB, Macnair MR, Willis JH. 2013. Indirect Nucleic Acids Res. 40(10):4288–4297. evolution of hybrid lethality due to linkage with selected locus in McClintock B. 1984. The significance of responses of the genome to chal- Mimulus guttatus. PLoS Biol. 11(2):e1001497. lenge. Science 226(4676):792–801. Wu Y, et al. 2016. Transcriptome shock in an interspecific F1 triploid hybrid Nordberg H, et al. 2014. The genome portal of the department of energy of Oryza revealed by RNA sequencing. J Integr Plant Biol. joint genome institute: 2014 updates. Nucleic Acids Res. 42:D26–D31. 58(2):150–164. Parisod C, et al. 2009. Rapid structural and epigenetic reorganization near Xiong Z, Gaeta RT, Pires JC. 2011. Homoeologous shuffling and chromo- transposable elements in hybrid and allopolyploid genomes in some compensation maintain genome balance in resynthesized allo- Spartina. New Phytol. 184(4):1003–1015. polyploid Brassica napus. Proc Natl Acad Sci U S A. Parisod C, et al. 2010. Impact of transposable elements on the organiza- 108(19):7908–7913. tion and function of allopolyploid genomes. New Phytol. Xu C, et al. 2014. Genome-wide disruption of gene expression in allopo- 186(1):37–45. lyploids but not hybrids of rice subspecies. Mol Biol Evol. Pietzenuk B, et al. 2016. Recurrent evolution of heat-responsiveness in 31:1066–1076. Brassicaceae COPIA elements. Genome Biol. 17(1):209. Yoo MJ, Szadkowski E, Wendel JF. 2013. Homoeolog expression bias and Quadrana L, et al. 2016. The Arabidopsis thaliana mobilome and its impact expression level dominance in allopolyploid cotton. Heredity (Edinb) at the species level. eLife 5:e15716. 110(2):171–180. Rangwala SH, Richards EJ. 2010. The structure, organization and radiation Zhang D, et al. 2016. Genome-wide gene expressions respond differently of Sadhu non-long terminal repeat retroelements in Arabidopsis spe- to A-subgenome origins in Brassica napus synthetic hybrids and nat- cies. Mobile DNA 1(1):10. ural allotetraploid. Front Plant Sci. 7:1508. Rebernig CA, et al. 2015. Non-reciprocal interspecies hybridization barriers Zhu W, et al. 2017. Altered chromatin compaction and histone methyla- in the Capsella genus are established in the endosperm. PLoS Genet. tion drive non-additive gene expression in an interspecific Arabidopsis 11(6):e1005295. hybrid. Genome Biol. 18(1):157. Robinson MD, McCarthy DJ, Smyth GK. 2010. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26(1):139–140. Associate editor: Ellen Pritham Genome Biol. Evol. 10(6):1403–1415 doi:10.1093/gbe/evy095 Advance Access publication May 18, 2018 1415 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1403/4999384 by Ed 'DeepDyve' Gillespie user on 20 June 2018

Journal

Genome Biology and EvolutionOxford University Press

Published: May 18, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off