Transcriptome-wide analysis of introgression-resistant regions reveals genetic divergence genes under positive selection in Populus trichocarpa

Abstract

Comparing gene expression patterns and genetic polymorphisms between populations is of central importance for understanding the origin and maintenance of biodiversity. Based on population-specific gene expression levels and allele frequency differences, we sought to identify population divergence (PD) genes across the introgression-resistant genomic regions of Populus trichocarpa. Genes containing highly diverged loci [i.e., genetic divergence (GD)] or showing expression divergence (ED) between populations were widely distributed in the genome and substantially enriched in functional categories related to stress responses, disease resistance, timing of flowering, cell cycle regulation, plant growth, and development. Nine genomic regions showing evidence of strong positive selection were overlapped with GD genes, which had significant differences between Oregon (a southernmost peripheral deme) and the other demes. However, we did not find evidence that genes under positive selection show an enrichment for ED. PD genes and genes under selection pertained to the same gene classes, such as SERINE/CYSTEINE PROTEASE, ABC TRANSPORTER, GLYCOSYLTRANSFERASE and other transferases. Our analysis also revealed that GD genes were polymorphic within the species (41.9 ± 3.66 biallelic variants per gene), as previously reported in herbaceous plants. By contrast, ED genes contained less genetic variants (10.73 ± 1.14) and were likely highly expressed. In addition, we found that trans- rather than cis-acting variants considerably contribute to the evolution of >90% PD genes. Overall, this study elucidates that cohorts of PD genes agree with the general attributes of known speciation genes and GD genes will provide substrates for positive selection to operate on.

Fig. 1. Introgression-resistant transcriptome-wide portrait of population divergence at genetic and expression levels.

Color coding for each deme is the same throughout this study. a Multidimensional scaling (MDS) plot based on pairwise identity-by-state distance for the introgression-resistant transcriptome. An inset graph shows drainages colored by deme based on Supplementary Fig. S1. b Linkage disequilibrium (LD) difference between demes. LD decays quickly within 10–15 kb for all demes but South-Oregon. All the four demes but South-Oregon have similar LD decay patterns. A horizontal dashed line marks r² dropping to 0.2. c, d nucleotide diversity (π) and Tajima’s D for each deme (shown by chromosome in Supplementary Figs. S3 and S4, respectively). Average values displayed in boxes of the two plots. *** for significance at P < 0.0001 by Wilcoxon test. e BEAST phylogenetic tree for the transcriptome along the introgression-resistant regions. f Percentage of differentially expressed (DE) genes between demes were based on fold-changes > 2 and FDR < 5% using edgeR. Top-left and bottom-right triangle matrix for leaves and xylem data, respectively.

Fig. 2. Integrated results of transcriptomes, selective scans, and selected SNPs.

In graph (a), tracks from outside to inside are: a Chromosome ideogram for the P. trichocarpa genes (construction detailed in Supplementary Note S2). b and c Gene expression in stem xylem and leaves, respectively. d–f Results of selective sweeps averaged over the demes using three approaches, namely, SweepFinder2, OmegaPlus, and OutFLANK, respectively. Significant loci shown in red. g Introgressed genomic regions from a previous study (see the Methods section). h Link ribbons showing the genes (Supplementary Table S2) with r² > 0.45 and pairwise distances between 250 and 500 kb. Graph (b) delineates the distribution of SNPs within genes ± 800 bp (N = 172,408 SNPs) across the transcribed introgression-resistant regions.

Fig. 3. Partial representation of mean individuals for each deme and groups representation based on the first six principal components (Dim 1–6) of the MFA.

In a, the balanced representation of each deme is the barycenter of the points summarizing partial points of view linked through colored lines. The projection of partial representations for leaf gene expression and SNPs of the five demes onto Dim 1 are very close, which means that leaf gene expression and SNPs define similar structures for demes on Dim 1; but this projection for xylem gene expression of the five demes is not the case: North (N) and North-South (N_S) are negative (coordinate close to −2) while the other demes (S, S_O, and O for South, South-Oregon, and Oregon, respectively) are positive (close to 1). Dim 1 is, therefore, specific to the xylem gene expression point of view that differentiates North and North-South from the other demes. On the Dim 2, the mean South-Oregon individuals from the partial xylem gene expression representation is positive (around 2 in the coordinate) while such a representation for the other demes is negative (0 to −0.5). Dim 2 is also specific to the xylem gene expression allowing to distinguish South-Oregon from the other demes. This is confirmed by analyzing the groups representation in b: only xylem gene expression has a coordinate value close to 1 (max. value) on Dim 1 and 2. Likewise, Dim 3–6 are specific to SNPs, which provide a partition of Oregon (and South-Oregon) with the other demes.

Fig. 4. eQTL, clustering, and GO enrichment of ED genes.

a cis- and trans-eQTLs identified in all genes (N = 26,328 and 27,548 for xylem and leaves, respectively) and ED genes. b Heatmap of the expression of ED genes based on Dim 1, 2, 4, and 6 of the MFA. Ward’s method for maximum distance used for clustering. c Adjusted P-value based on the MIC test against fold-change in gene expression for two-group stratified demes on each Dim shown in a. Red dots are ED genes of each Dim and horizontal dashed lines represent significance at α = 0.05. d Top five GO terms identified by the weight algorithm for down-weighting genes in the GO enrichment. Boxes indicate the five most important terms, where filled color represents the relative significance with dark red and light yellow for most and least significant, respectively. No color ellipses for the rest of relatively less important GO terms. Black arrows indicate is-a relationships. Note that all P-values > 0.05 suggest no significant difference among the GO terms for candidate ED genes.

Fig. 5. Merge of loci showing strong positive selection and PD with those distinct in Oregon.

a Loci identified by Dim 3 and 5 of the MFA distinguish Oregon from the other demes. Red asterisks indicate loci under strong positive selection only for Oregon and/or South-Oregon (more detail in Table 1). Regions with very low π and high F_ST in Oregon relative to the other demes (Supplementary Figs. S3 and S5) shaded in light gray in chr01, 02, and 15. Loci showing PD across demes or in Oregon based on the HMM are marked at the top (in black) or bottom (in gray) of each panel, respectively. PD genes were also marked in black dots. b P-values were calculated using the MIC test for allele frequencies between Oregon and the other demes. The horizontal red line marks the significance at α = 0.05.

Fig. 6. Local ancestry analysis and comparisons of gene expression (mean ± SE) and allele frequency between pure and hybrid individuals.

a The ancestry assignment is based on a whole-genome local ancestry analysis for these individuals. b According to the ancestry analysis, 22 pure P. trichocarpa individuals (N = 151) were identified, in which 10 and 12 individuals were originated from North and South demes, respectively. c, d Using the MIC test, we probed 10 genes with significant differential expression (adjusted P < 0.05 and absolute fold-change > 2) and six genes each with ≥10 SNPs having significant difference in allele frequency (P < 0.05). See Supplementary Fig. S22 for details. Gene functions are briefly described beneath each panel. Allele frequency ⊂ {0, 1, 2}.