mRNA-Seq characteristics of meiotic recombination pathway mutants
Previous reports from our lab and others revealed that the increased meiotic crossover frequencies observed in Arabidopsis grown at elevated temperatures are mediated by the Type I interference sensitive pathway [2, 3]. To identify genes potentially involved in the hyper-recombinant phenotype, we collected stage 4–7 anthers that contain male meiocytes [14] from WT (Col-0), mus81, and msh4 plants. We collected two biological replicates for each genotype grown at 20 °C and 28 °C, then sequenced the 12 mRNA libraries. A total of 690,371,408 raw reads were retrieved from the 12 mRNA-Seq datasets with an average of 53 million mapped reads per dataset and a 93% average mapping rate (Supplementary Table 1). The biological replicates for each genotype-temperature combination had a minimum 0.92 correlation coefficient (Supplementary Table 2) indicating high reproducibility among our datasets. We examined the number of expressed genes in each condition with a threshold equal to or greater than 1 TPM (Transcripts Per Kilobase Million). At 20 °C, WT, mus81 and msh4 have 20,489, 20,755 and 20,116 expressed genes, respectively (Fig. 1a). After 28 °C heat stress treatment, the number of expressed genes decreased in all three genotypes (Fig. 1a). Both WT and mus81 have 1.6% fewer expressed genes, while msh4 has 3.3% fewer. We also compared the average level of gene expression in two biological replicates for each genotype at 20 °C and 28 °C. The TPM values do not differ significantly for WT (18.4 and 19.4 at 20 °C verus 18.8 and 20.1 at 28 °C, Mann–Whitney test, P value = 0.05) and mus81 (18.4 and 18.1 at 20 °C verus 18.5 and 18.4 at 28 °C, Mann–Whitney test, P value = 0.18), but did for msh4 (17.7 and 17.9 at 20 °C versus 20.0 and 20.2 at 28 °C; Mann–Whitney test, P value < 2.2e-16) (Fig. 1b), suggesting that heat stress decreases the number of expressed genes in msh4 but increases the average strength of gene expression. We also compared the number of expressed genes in seven TPM value groups from low to high and found that all treatments and genotypes are similarly distributed, and that the group with TPM values between 1–10 is the largest (Fig. 1c). We then examined the overlap of expressed gene in WT, mus81, and msh4 at 20 °C and 28 °C. At least 90% (18,398) of expressed genes are shared between all datasets, and 63 genes are shared by only WT and mus81 at 28 °C (Supplementary Fig. 1a). Only one enriched GO term, phospholipase C activity from Molecular Function (GO: 0004629, FDR = 0.02, by AT3G03530 and AT3G03540) was detected among these 63 genes, no meiotic recombination related genes were found, and most of the 63 genes have relatively low gene expression (Supplementary Fig. 1b).
WT and mus81 have similar differential gene expression profiles
By comparing expression levels of genes at 28 °C to those at 20 °C in each genotype, we found 2,922, 2,363 and 4,009 differentially expressed genes (DEGs) from WT, mus81 and msh4, respectively. Those include 1,196 and 1,727 up- and down-regulated genes in WT, 1,015 and 1,348 in mus81, and 835 and 3,174 in msh4. We analyzed the intersections of the up- and down-regulated DEGs (Fig. 2a, b) to determine whether any of the genotypes had similar profiles. On a proportional basis, WT and mus81 share 22% of their up-regulated DEGs and 25% of their down-regulated DEGs, WT and msh4 share 19% and 19% respectively, and mus81 and msh4 share 13% and 12% respectively. WT and mus81 share significantly more up-regulated DEGs compared to mus81 and msh4 (χ2 = 36.78, P value = 1.32E-09), and not significantly more to WT and msh4 (χ2 = 3.56, P value = 0.06). Similarly, WT and mus81 share significantly more down-regulated DEGs compared to WT and msh4 (χ2 = 24.51, P value = 7.38E-07), as well as mus81 and msh4 (χ2 = 158.02, P value < 2.2E-16). It should be noted that disruption of the Type I CO pathway in msh4 compromises but does not abolish pollen development and fertility, while WT and mus81 are fertile, and that these differences in phenotype likely influence the differential gene expression patterns. Taken together, these data suggest that WT and mus81 have the most similar DEG profiles of the genotypes examined (Fig. 2a).
We performed GO enrichment analysis to explore the characteristics of the DEGs from WT, mus81 and msh4 grown at 20 °C and 28 °C. Among the up-regulated DEGs, 45 Biological Process (BP), 7 Molecular Function (MF) and 3 Cellular Component (CC) GO terms were enriched in WT; 71, 19 and 7 in mus81; and 36, 0 and 0 in msh4 (Supplementary Data 1). Several GO terms enriched among the 1,196 DEGs from WT and 1,015 from mus81 are not enriched in the 835 DEGs from msh4 including: 22 BP GO terms including cellular response to decreased oxygen levels (GO:0036294), response to heat (GO:0009408), protein folding (GO:0006457), homologous chromosome segregation (GO:0045143), and chromosome organization involved in meiotic cell cycle (GO:0070192), 4 MF GO terms including unfolded protein binding (GO:0051082) and 2 CC GO terms including anchored component of membrane (GO:0031225) (Fig. 3 and Supplementary Data 1). Among the down-regulated DEGs, there were 26 BP, 4 MF and 0 CC GO terms enriched in WT; 34, 18 and 8 in mus81; and 31, 13 and 9 in msh4 (Supplementary Data 1). Only 6 BP GO terms, including response to absence of light (GO:0009646) and glucose import (GO:0046323), and 1 MF (carboxylic ester hydrolase activity; GO:0052689) are enriched in the down-regulated DEGs from WT and mus81 (Supplementary Data 1 and Supplementary Fig. 2). We noticed that msh4 has the most down-regulated DEGs, and enriched GO terms not shared with WT or mus81 that relate to male gametogenesis, such as pollen wall assembly (GO:0010208), pollen sperm cell differentiation (GO:0048235) and pollen exine formation (GO:0010584) (Supplementary Fig. 2A). As mentioned above, this may be a consequence of the compromised pollen development phenotype observed in msh4.
To determine whether the up- and down-regulated DEGs from WT, mus81 and msh4 display any pattern along the 5 Arabidopsis chromosomes we plotted their density as a function of physical position (Fig. 4a, c). We did not detect any global patterns in distribution: up-regulated DEGs from WT and mus81 have a correlation coefficient of 0.46, while WT and msh4 have a correlation coefficient of 0.42. Similarly, the down-regulated DEGs have a correlation coefficient of 0.51 between WT and mus81 and 0.43 between WT and msh4. To examine the local distribution patterns of DEGs, we divided each chromosome into 12 euchromatic segments of equal length and one heterochromatic segment and tested whether the number of WT DEGs in each segment differs significantly from those of mus81 or msh4. The number of up-regulated DEGs from mus81 and msh4 do not differ significantly from WT (Fig. 4b). Interestingly, two segments on the north arms of chromosome 1 and 4 show significantly more down-regulated DEGs in msh4 than in WT (Fig. 4d). Overall, the distribution of DEGs does not differ dramatically among the genotypes but is consistent with WT and mus81 being more similar to one another than msh4.
GO terms enriched in the DEGs shared by WT and mus81
Exposure to elevated temperature results in an increase in Type I (MSH4-dependent) COs in Arabidopsis [2, 3]. To determine if that phenotype is correlated with changes in gene expression, we looked for DEGs shared between WT and mus81 (in which the Type II pathway is inactive) but not msh4. We found 278 and 384 down-regulated DEGs that may be implicated in modulating the activity of the two pathways (Fig. 2a, b). GO analysis of the 278 up-regulated DEGs reveals enrichment of 3 BP GO terms: response to heat (GO:0009408), protein folding (GO:0006457) and response to toxic substance (GO:0009636), and 1 MF GO term: unfolded protein binding (GO:0051082) (Fig. 3d and Supplementary Data 2). Among the 384 down-regulated DEGs, 14 BP GO terms were enriched including: cellular response to decreased oxygen levels (GO:0036294), response to absence of light (GO:0009646), response to water deprivation (GO:0009414), and circadian rhythm (GO:0007623) (Fig. 3d and Supplementary Data 2).
ASY1 is exclusively up-regulated in WT and mus81 but not in msh4
We were interested in whether meiosis related genes are differentially expressed under heat stress. To this end, we examined the expression patterns of 148 genes that have been previously reported to have meiotic function (Supplementary Data 3). A heatmap comparison of these genes does not reveal any significant global difference in expression levels at 28 °C compared to 20 °C in WT, mus81 or msh4 (Fig. 5). To our surprise, ASY1 and ASY2 are the only consistently up-regulated DEGs exclusively shared by WT and mus81 but not msh4 (Fig. 2c, d; Supplementary Fig. 3A and 8C) (ASY1 and ASY2 show an increase below the 2-fold cut-off at 28 °C compared to 20 °C in msh4, in contrast to 2.9 and 2.3-fold increases in WT, and 4.6 and 3.5-fold increases in mus81). ASY1 encodes a component of the chromosome axis that forms along the length of replicated sister chromatids during meiosis [15], and is required for proper interhomolog interaction including chromosome pairing, synapsis and recombination [16], as well as ensuring crossover interference [17]. ASY2 is a putative functional homolog of ASY1 [15]. No meiosis-related DEGs were found among the 384 down-regulated DEGs shared by WT and mus81 (Supplementary Data 3). Interestingly, 4 genes are up-regulated in all three genotypes: the meiosis-specific cyclin-like gene SDS [18], meiosis-specific recombinase DMC1 [19], meiosis-specific cohesion complex component SYN1/REC8 [20], and SWI1 which functions in meiotic sister chromatid cohesion [21] (Supplementary Data 3). However, none of the 14 known meiotic genes involved in either the Type I or Type II CO pathways were differentially expressed (Table 1). The apparent up-regulation of MSH4 in the msh4 mutant background at 28 °C is most likely the result of the truncated msh4 being driven by the 35S promoter carried by the T-DNA transgene insertion in this allele (Supplementary Fig. 3B). Interestingly, we found one of the 21 Arabidopsis Heat Stress transcription Factor (HSF) genes, AtHSFA2 (AT2G26150) [22], is among the 278 DEGs up-regulated in WT and mus81 but not msh4 (Supplementary Data 4). In addition, AtHSFB2a (AT5G62020) and AtHSFA3 (AT5G03720) are up-regulated in WT, and AtHSFB1 (AT4G36990) is up-regulated in mus81. Surprisingly, no HSF transcription factors are up-regulated in msh4 (Supplementary Data 4). Further experimental work is needed to understand the potential interaction of the heat shock transcriptional regulatory networks and the meiotic recombination machinery.
Table 1
Expression value (TPM) of genes in meiotic recombination type I and type II pathways.
Gene name | Alias | WT 20C rep1 | WT 20C rep2 | WT 28C rep1 | WT 28C rep2 | mus81 20C rep1 | mus81 20C rep2 | mus81 28C rep1 | mus81 28C rep2 | msh4 20C rep1 | msh4 20C rep2 | msh4 28C rep1 | msh4 28C rep2 |
Type I CO pathway regulators |
AT4G17380 | MSH4 | 7.3 | 9.8 | 9.0 | 10.4 | 5.3 | 4.6 | 5.8 | 7.4 | 8.6 | 9.6 | 38.0 | 50.1 |
AT3G20475 | MSH5 | 20.4 | 25.2 | 20.4 | 24.5 | 14.6 | 12.3 | 19.0 | 20.3 | 23.3 | 20.5 | 16.5 | 16.5 |
AT5G48390 | ZIP4 | 4.8 | 5.7 | 7.2 | 7.9 | 3.1 | 2.9 | 4.8 | 5.9 | 5.3 | 4.9 | 6.6 | 6.4 |
AT3G27730 | RCK/MER3 | 4.8 | 8.1 | 6.3 | 8.2 | 3.0 | 2.7 | 4.6 | 5.1 | 7.6 | 5.5 | 8.0 | 7.8 |
AT1G12790 | PTD | 15.1 | 15.4 | 14.5 | 17.2 | 8.9 | 12.3 | 9.9 | 11.8 | 14.3 | 9.7 | 12.9 | 13.4 |
AT5G52290 | SHOC1/ZIP2 | 1.3 | 2.0 | 2.5 | 3.4 | 1.0 | 1.0 | 1.2 | 1.6 | 1.8 | 1.6 | 3.1 | 3.1 |
AT1G53490 | HEI10 | 6.9 | 9.2 | 9.7 | 10.7 | 26.5 | 28.1 | 27.6 | 22.9 | 7.4 | 6.3 | 12.6 | 11.5 |
AT4G09140 | MLH1 | 26.5 | 26.8 | 26.3 | 29.6 | 21.2 | 21.0 | 23.0 | 21.5 | 24.8 | 26.4 | 27.2 | 26.2 |
AT4G35520 | MLH3 | 6.8 | 7.5 | 6.8 | 9.1 | 5.8 | 5.9 | 9.0 | 6.8 | 7.1 | 8.1 | 4.8 | 5.9 |
AT1G63770 | MPA1 | 122.1 | 115.0 | 116.4 | 128.4 | 105.0 | 96.5 | 111.5 | 114.4 | 117.6 | 120.3 | 124.8 | 119.7 |
Type II CO pathway regulators |
AT4G30870 | MUS81 | 22.6 | 24.0 | 25.0 | 30.2 | 4.5 | 5.4 | 5.0 | 4.9 | 21.0 | 20.2 | 17.2 | 19.1 |
AT2G21800 | EME1A | 4.4 | 5.6 | 7.5 | 7.2 | 4.7 | 5.7 | 7.0 | 5.6 | 4.7 | 5.5 | 8.4 | 9.2 |
AT2G22140 | EME1B | 3.1 | 4.2 | 6.0 | 5.6 | 3.8 | 3.0 | 2.9 | 3.2 | 3.7 | 3.4 | 5.8 | 5.4 |
AT4G14970 | FANCD2 | 5.9 | 8.3 | 10.7 | 12.1 | 5.2 | 4.6 | 7.1 | 7.9 | 7.8 | 7.6 | 12.9 | 14.1 |
Novel anther-specific transcripts associated with A-rich promoter motifs
Our RNA-seq analysis of anthers identified 51 novel transcripts that had not been previously annotated in the Arabidopsis genome (Supplementary Data 5). 39% of the novel transcripts have relatively low TPM values between 1–5, 14% have TPMs of 5–10 (Fig. 6a) similar to previously annotated genes, 35% of which have TPMs of 1–10 (Fig. 1c). The novel transcripts are short with a median length of 809 bp compared to 1,787 bp from the 38,194 annotated genes (Fig. 6b). However, 26 (46%) of the novel transcripts have 2 exons, 14 (25%) have 1, and 8 (14%) have 4 exons, compared to annotated genes in which the most abundant group (21%) has only one exon (Fig. 6c). Comparison of the transcripts expressed at 20 °C and 28 °C in the different genotypes reveals Differentially Expressed Novel Transcripts (DENTs) which include 5, 3 and 7 up-regulated and 6, 2 and 6 down-regulated DENTs from WT, mus81 and msh4, respectively (Fig. 6d). Noticeably, all DENTs in mus81 are shared by WT (Fig. 6d), which is consistent with the previous observations that the transcriptome profiles of WT and mus81 are more similar compared to msh4.
The loci encoding the novel transcripts are not preferentially partitioned into euchromatin or heterochromatin (Fisher's Exact Test, P value = 0.18) (Fig. 6e). However, we wondered whether they are expressed in a tissue-specific manner. To address this question, we generated two RNA-seq datasets from WT leaf tissue and compared them to the anther datasets (Supplementary Table 1). We found that 20 (40%) of the 51 novel transcripts are anther-specific and 5 more are preferentially expressed in the anther when comparing 20 °C anther datasets to WT leaves (Fig. 6e; Supplementary Data 5). Like the larger novel transcript set, the 20 anther-specific transcripts are not preferentially partitioned into euchromatin or heterochromatin (Fisher's Exact Test, P value = 0.29) (Fig. 6e). To test whether these anther-specific novel transcripts are regulated by similar cis-regulatory elements we used MEME to search for conserved DNA motifs [23]. 43 (84%) of the novel transcripts had an upstream A-rich motif in their putative promoter regions, including 90% (18 from 20) of the anther-specific novel transcripts (E value = 8.8e-047) (Fig. 6f). These A-rich motifs have been reported to be associated with open chromatin structure and meiotic recombination hotspots [24–27], suggesting that the anther-specific novel transcripts may be expressed in response to a change in chromatin state associated with meiotic recombination. Moreover, of the 51 novel transcripts we identified 9 can potentially code for proteins (Supplementary Data 6). Of those MSTRG.16271, resembles a Danio rerio gene that encodes the guanine nucleotide exchange factor subunit RIC1 (Fig. 6g).