Evolutionary Genomics Of Parnassius Glacialis (Papilionidae:Parnassiinae) Dispersal Out Of The Qinghai-Tibet Plateau And Local Adaptation

Global environmental change and biological evolution are currently hot issues in the study areas of evolutionary biology.As far as we know,many organisms have dispersed from the Qinghai-Tibetan Plateau（QTP）in to Asia,Europe and North America during the Quaternary Ice Age,which has made an important impact on the formation of biodiversity for the Northern Hemisphere.However,little is known about the underlying mechanisms of genome evolution and local adaptation for these alpine organisms out of the QTP.Parnassius is a genus of typical alpine butterflies belonging to the subfamily Parnassiinae（Lepidoptera:Papilionidae）,and in recent years,we have known that they have the largest genomes compared with other butterfly groups.Previous studies have shown that Parnassius was originated at the early Miocene period（about 21～17 million years ago（Ma））under the complex geological and environmental change of the QTP’s orogenic movement and later rapidly radiated and dispersed into other areas of the Northern Hemisphere.During the long evolutionary process,most of the Parnassius butterfly species stayed endemically in the high-altitude QTP,while a few Parnassius species that have spread into the middle-and low-altitude areas of Europe,North America,and East Asia（mainly central-eastern China）.Among these"out of the QTP"species,Parnassius glacialis is the only one that has spread from northwest to southeast China,reaching up to the southern areas of the Yangtze River in the Quaternary Ice Age,which provides an ideal species of the"out of the QTP"organisms for studying the underlying mechanisms of genome evolution and population differentiation.In this study,we assembled a high-quality chromosome-level genome of the P.glacialis through Denovo genome sequencing and resequenced the genome sequences of their 41 individuals from 9 populations at different elevations ranging from 300 to 1,800m.Based on these sequence data,we analyzed the pattern of large genome formation,gene family expansion,population genetic structure,historical population dynamics,and selective sweeps of the P.glacialis through comparative genomics and population genetics methods.Furthermore,we explored the impact of transposable element（TE）,the key factor in shaping and reorganizing the genome,on genome composition,gene duplication,population genetic differentiation and selective evolution in P.glacialis.This work will provide new insights into the evolutionary adaptation of"Out of the QTP"alpine organism during the Quaternary period.Main results of this study are listed as follows:1.Genome size assessment of 10 Parnassius species at different elevations in ChinaWe assessed the genome size of 10 Parnassius species（P.glacialis,P.stubbendorfi,P.apollo,P.apollonius,P.acdestis,P.andreji,P.simo,P.imperator P.orleans and P.nomion）at different elevations（300～5000 m）based on genome survey sequencing.The results showed that the genome sizes of the Parnassius butterflies mainly range from 1.00to 1.90 Gb,and the genome size is significantly correlated with the elevation of altitudinal distribution.Among them,the endemic species distributed in the high-altitude region（3000～5000 m）of the QTP,such as P.acdestis,P.andreji,P.simo,P.imperator and P.orleans,generally have the smaller genome sizes ranged from 1.0 to 1.22 Gb,whereas the species at middle-and low-altitudes（300～3000 m）,such as P.nomion,P.stubbendorfi,P.apollo,P.apollonius and P.glacialis,harbored the relatively larger genome sizes ranged from 1.27 to 1.90 Gb.In conclusion,the Parnassius species outside the QTP generally have larger genome size than the endemic species on the QTP.2.Genome assembly and genome size formation in P.glacialisWe completed the chromosome-level genome assembly and annotation of the P.glacialis butterfly based on Pac Bio and Hi-C sequencing.The genome assembly results showed that the genome of the P.glacialis butterfly was about 1.35 Gb in size,composing of 29 hypothetical chromosomes,with a Contig N50 length of 6.45 Mb,a Scaffold N50length of 49 Mb,and a BUSCO integrity of 95.8%.Gene annotation results showed that P.glacialis probably harbor total 17,080 protein-coding genes,with an average transcript length of about 21.08 kb,an average intron length of 4.52 kb,and an average number of5.37 exons.Compared with the Papilio butterflies,the genes of P.glacialis have longer transcripts,longer introns and fewer exons.Repeat annotation showed that the P.glacialis genome contains about 971 Mb of repetitive sequences（71.8%）,which is one of the highest in the known butterfly species.Among these repetitive sequences,LINE,DNA,LTR,SINE and Helitron transposons accounted for about 38.9%,10.64%,6.28%,3.40%and 4.38%in the P.glacialis genome,respectively.Combined with the published genomes of P.orleans and P.apollo,we further analyzed the Kimura distribution and insertion time of transposons in 3 Parnassius butterflies（P.glacialis,P.apollo and P.orleans）,and assessed the unequal recombination rate of LTR retrotransposons（Solo-LTRs/complete LTRs）.The results revealed that transposon presented a shared outbreak during the Quaternary period（～3 Ma to the present）,and this transposon outbreak event mainly involved a few types of transposons,including LINE/RTE,LINE/CR1,LTR/Gypsy and LTR/Pao.We speculate that that the glacial-interglacial cycles in the Quaternary period may evoke the transposon activity to format the larger genome of the Parnassius butterflies.Meanwhile,the unequal recombination rates of LTR retrotransposons were about 1.15 and 0.08 in P.glacialis and P.apollo outside the QTP,which were much smaller than 2.67 of P.orleans on the QTP.It suggests that the Parnassius species at high altitude inhibits the increase of genome size through the stronger unequal recombination,whereas the Parnassius species at the low and middle altitudes cannot inhibit effectively,which could explain the pattern of increased transposons and larger genome size for the Parnassius species out of the QTP.3.Mechanisms of gene family expansion and evolutionary adaptation in P.glacialisThe phylogenetic relationships and divergence times of 9 lepidopteran species were constructed based on 8 butterflies（Maniola hyperantus,Danaus plexippus,Pieris rapae,Papilio xuthus,Papilio bianor,P.glacialis,P.apollo and P.orleans）as well as Bombyx mori.The results showed that the divergence time between genera Parnassius and Papilio was about 63.6 Ma,and that between P.glacialis and P.orleans was about 5.91 Ma respectively.In addition,703 expanded gene families mainly enriched in the ribosome metabolism pathway,especially the RPLP2 genes,were detected in P.glacialis through comparative genomic analysis.Among these RPLP2 genes,2 and 434 were detected in P.apollo and P.glacialis,respectively,while only one was detected in other butterflies.Moreover,only one of the 434 RPLP2 genes in P.glacialis was a complete functional gene containing the introns（Pglac-RPLP2）,while the other 433 showed the characterization of processing pseudogenes such as no intron,promoter deletion or premature codon termination.Further analysis indicated that the functional Pglac-RPLP2mediated the duplications of other RPLP2 pseudogenes by forming a complex chimera as 5’LTR-Pseudogene-IN-3’LTR.Meanwhile,most of the RPLP2 pseudogenes in P.glacialis were likely to occur from about 3 Ma to the present,which coincided with the timing of the transposon outbreak（about～3 Ma to the present）mainly during the Quaternary period.The average interspecific Ka/Ks value of functional RPLP2 genes between the P.glacialis and other butterflies was about 0.08,while the average intraspecific value of the P.glacialis’s RPLP2 pseudogenes reached up to 0.36.This case indicated that the RPLP2pseudogenes of P.glacialis have faster evolutionary rate than the functional genes.Additionally,transcriptomic analysis showed that most of the RPLP2 pseudogenes in P.glacialis had no expression,but 111 and 67 pseudogenes presented a variety of expression levels in the high-altitude population（BQXL）and the low-altitude population（BQTM）,respectively.Furthermore,4 pseudogenes specifically expressed in the low-altitude population（BQTM）have evolved putative start codons.Considering that the RPLP2 gene is related to protein synthesis,antioxidant and inhibition of cellular senescence,we proposed that these transposon-mediated RPLP2 pseudogenes have undergone rapid evolution to obtain the ability of functional genes,possibly as a response of P.glacialis to the low-altitude environments outside the QTP,such as physical development,oxidative damage and pathogenic microorganism invasions.4.Mechanisms of population genetic and local adaptation in P.glacialisThe genomic resequencing data of 41 samples from nine different altitude populations（300～1800 m）were used to analyze the population diversity,genetic structure and historical population dynamics for P.glacialis.The results showed that the genetic diversity of the low-altitude populations in eastern China was generally lower than that of the middle-and high-altitude populations in central or western China.The average nucleotide polymorphismθπvalues of the three eastern low-altitude populations（BQKY,BQLS,BQTM）were only 0.00156,whereas the averageθπvalues of the three central mid-altitude populations（BQLJ,BQTT,BQTA）and the three western high-elevation populations（BQXL,BQHD,BQSN）were 0.00172 and 0.00178,respectively.The results of genetic structure analysis indicated that the southeastern low-altitude population BQTM（～300 m）was one cluster,while the other eight groups formed another cluster.At the same time,BQTM was found to harbor high genetic differentiation(F_ST)with other populations.The historical population dynamics analysis showed that the nine populations had similar effective population sizes and dynamics before the Last Glacial Maximum（LGM）,with a differential size increase in western populations,and a sharp size decrease in eastern populations（especially BQTM）during the global warming after LGM.In conclusion,the low-altitude populations of the P.glacialis was likely to be exerted stronger environmental stresses,and thus displaying relatively smaller population size,lower nucleotide polymorphism（θπ）and higher genetic differentiation(F_ST).Moreover,we identified selective sweep signals in low-altitude population BQTM（～300 m）using four methods(θπratio,F_ST,XPEHH and Tajima’s D)and identified 454candidate genes with selective signatures.Among them,9 genes（LAP3 and 8 Anpep）were enriched in the glutathione metabolic pathway,6 genes（Ds,Ed,Lgl,Mer,Pat J and Sdt）were enriched to the Hippo signaling pathway,6 genes（PSH,Toll,P38,Dl/Dif and2 PGRP）were enriched to the Toll and Imd signaling pathways,and 4 genes（2 STE24and 2 FDPS）were enriched to terpene skeleton biosynthesis.4 genes（THSD4,FST,INHBA and INHBB）enriched to TGF-βsignaling pathway,and 3 genes（NAGA and 2HEXA＿B）enriched to glycosphingolipid biosynthesis pathway.These genes with selective signatures are mainly associated with the functions of antioxidant,growth and development,and immune responses in insects,which could contribute to local adaptation of the new low-altitude habitat for P.glacialis.In addition,we identified 140,666 transposon-mediated structural variants（TE-SVs）between the high-altitude population（BQXL）and low-altitude population（BQTM）to further investigate the effects of transposons on population genetic differentiation and selective sweep.The results showed that the average F_ST of regions（10 kb window）with TE-SVs were significantly larger than those of regions without TE-SVs in P.glacialis.Furthermore,we compared the number of TE-SVs located in the genes as well as their upstream and downstream regions in P.glacialis,and the results showed that the average number of TE-SVs in 454 genes with selective signatures is about 6.5,which is significantly larger than that（～5.0）in the other genes without selective signatures.In conclusion,our results support that transposons play an important role in population genetic differentiation and local adaptation of the P.glacialis butterfly.