Interspecific and Intraspecific Transcriptomic Variations Unveil the Potential High-Altitude Adaptation Mechanisms of the Parnassius Butterfly Species

Abstract

Parnassius butterflies have significantly advanced our understanding of biogeography, insect-plant interactions, and other fields of ecology and evolutionary biology. However, to date, little is known about the gene expression patterns related to the high-altitude adaptation of Parnassius species. In this study, we obtained high-throughput RNA-seq data of 48 adult Parnassius individuals covering 10 species from 12 localities in China, and deciphered their interspecific and intraspecific expression patterns based on comparative transcriptomic analyses. Though divergent transcriptional patterns among species and populations at different altitudes were found, a series of pathways related to genetic information processing (i.e., recombination, repair, transcription, RNA processing, and ribosome biogenesis), energy metabolism (i.e., oxidative phosphorylation, thermogenesis, and the citrate cycle), and cellular homeostasis were commonly enriched, reflecting similar strategies to cope with the high-altitude environments by activating energy metabolism, enhancing immune defense, and concurrently inhibiting cell growth and development. These findings deepen our understanding about the molecular mechanisms of adaptative evolution to extreme environments, and provide us with some theoretical criteria for the biodiversity conservation of alpine insects.

Keywords: Parnassius butterflies; adaptative evolution; alpine environment; high-throughput sequencing; transcriptome.

Figure 1

The sampling localities of the Parnassius samples in this study.

Figure 2

Multiple regression analysis based on altitude, geographic location, and host plant. (a) Modeling with Jaccard distance matrix; (b) modeling with Bray–Curtis distance matrix.

Figure 3

The number of differentially expressed genes for each pairwise comparison. Bar chart diagrams showing overlaps of DEGs with increased (red) or decreased (blue) transcript abundance in three pairs of comparisons. (a) P. glacialis was used as the reference genome; (b) P. cephalus was used as the reference genome.

Figure 4

(a–c) KEGG enrichment results of the DEGs (P. glacialis was used as the reference genome: left, up-regulated; right, down-regulated) in pairwise comparisons of HA vs. LA, MA vs. LA, and HA vs. MA (the latter group used as the control, similarly hereinafter), respectively. For each comparison, only the top 10 pathways with the most significant enrichment are shown. The shared pathways that appear in the same pairwise comparison are highlighted in red, regardless of the reference genome (similarly hereinafter).

Figure 5

(a–c) KEGG enrichment results of the DEGs (P. cephalus was used as the reference genome: left, up-regulated; right, down-regulated) in pairwise comparisons of HA vs. LA, MA vs. LA, and HA vs. MA, respectively. For each comparison, only the top 10 pathways with the most significant enrichment are shown.

Figure 6

KEGG enrichment results of the DEGs (P. glacialis is used as the reference genome: left, up-regulated; right, down-regulated) in pairwise comparisons of And_el vs. And_my (a), Sim_el vs. Sim_bc (b), Nom_ck vs. Nom_bm (c), and Epa_hs vs. Epa_gg (d), respectively. For each comparison, only the top 10 pathways with the most significant enrichment are shown.

Figure 7

KEGG enrichment results of the DEGs (P. cephalus is used as the reference genome: left, up-regulated; right, down-regulated) in pairwise comparisons of And_el vs. And_my (a), Sim_el vs. Sim_bc (b), Nom_ck vs. Nom_bm (c), and Epa_hs vs. Epa_gg (d), respectively. For each comparison, only the top 10 pathways with the most significant enrichment are shown.

Figure 8

(a) Hierarchical clustering tree (gene dendrogram) showing 11 modules of genes co-expressed by WGCNA using P. glacialis as the reference genome. The major tree branches constitute 11 modules, labeled with different colors. (b) Module–locality relationship (each row represents a module, each column represents a specific sampling locality, and the correlation coefficient between module and locality is represented by the value in each cell at the row–column intersection, with the p-value shown in parentheses. (c,d) KEGG enrichment analyses of the genes in the turquoise and red modules, respectively. For each module, only the top 10 pathways with the most significant enrichment are shown (Table S7).

Figure 9

(a) Hierarchical clustering tree (gene dendrogram) showing nine modules of genes co-expressed by WGCNA analyses using P. cephalus as the reference genome. The major tree branches constitute nine modules, labeled with different colors. (b) Module–locality relationship (each row represents a module, each column represents a specific sampling locality, and the correlation coefficient between module and locality is represented by the value in each cell at the row–column intersection, with the p-value shown in parentheses. (c,d) KEGG enrichment analyses of the genes in the turquoise and brown modules, respectively. For each module, only the top 10 pathways with the most significant enrichment are shown.

Figure 10

Validation expression patterns in P. glacialis and P. cephalus representative samples determined by qPCR. The left ordinate represents the qPCR-based expression levels, and the right ordinate represents the RNA-seq-based expression levels.