1. Introduction

The Sino-Himalayan region is a temperate biodiversity hotspot with a high level of species endemism (Myers et al., 2000). Much of this diversity and endemism is attributed to the heterogeneous topography of the region and the diverse climates it supports (Shrestha et al., 2018). The formation of plant diversity in this region is synchronous with the uplift of the Qinghai-Tibet Plateau and has been intensified by climatic oscillations, especially the glaciale-interglacial cycle during the Quaternary (Muellner-Riehl et al., 2019; Sun et al., 2017; Wen et al., 2014; Wu, 1988). The mechanisms of plant speciation and diversification in this region are complicated and various (Muellner-Riehl et al., 2019; Wen et al., 2014). On the one hand, heterogeneous topography can create geographical barriers, which promote allopatric species divergence (e.g., Luo et al., 2017), and recent climatic oscillations can give rise to habitat fragmentation, which also interrupts the gene flow (e.g., Ma et al., 2017). On the other hand, divergent selection pressures across elevational gradients or heterogeneous habitats can also cause adaptive divergence in sympatry or parapatry (Funk et al., 2016), and even lead to adaptive radiation. This type of ecological divergence has been shown to be an essential speciation process for plants in the Sino-Himalayan region (e.g., Zhao et al., 2016).

Megacodon (Hemsl.) H. Smith is a small genus nested within a polytomy that contains Swertia L. and other genera in the Swertiinae (Chassot et al., 2001; Favre et al., 2016). However, Megacodon is morphologically distinct from other Swertiinae genera by having large campanulate flowers (Ho and Pringle, 1995). Megacodon contains two species: Megacodon stylophorus (C. B. Clark) H. Smith and Megacodon venosus (Hemsl.) H. Smith (Ho and Pringle, 1995). M. stylophorus is found at high elevations between 3000 and 4400 m a.s.l in western Sichuan, north-western Yunnan, southern Tibet of China, also in Bhutan, Nepal, and north-east India. M. venosus is endemic to China and is mainly found in Sichuan, western Hubei, and northern Chongqing at low to middle elevations between 600 and 2400 m a.s.l (Ho and Pringle, 1995). Previous studies of Megacodon have primarily focused on M. stylophorus, its phytochemistry (Liu et al., 2014), genetic diversity (Ge et al., 2005), embryology (Xue and Li, 2005), and breeding system (Meng et al., 2012). In contrast, M. venosus has been scarcely collected or reported. M. venosus was considered extinct until its rediscovery in 2007 by Prof. Zhen-Yu Li (Institute of Botany, the Chinese Academy of Sciences) in Northern Chongqing (Liu et al., 2014). Thus, it is vital to carry out investigations and conservation studies for this genus.

With rapid advances in next-generation sequencing (NGS), RNA-sequencing (RNA-seq) has become a more efficient approach for collecting gene sequences and studying ecological divergence in nonmodel organisms (Hudson, 2008; Strickler et al., 2012). For example, comparative transcriptome analysis can be used to identify candidate genes that underlie species diversification and adaptation, and then link these genes to specific biological processes (Voelckel et al., 2017). Previous studies have used genomeand transcriptome-based approaches to investigate adaptations to the high elevations of the Qinghai-Tibet Plateau region; however, few studies have examined these adaptations in plants, e.g., Crucihimalaya himalaica (Edgew.) Al-Shehbaz, O'Kane & R.A.Price (Qiao et al., 2016; Zhang et al., 2019), Primula poissonii Franchet and Primula wilsonii Dunn (Zhang et al., 2013), Cupressus gigantea W. C. Cheng & L. K. Fu and C. duclouxiana Silba (Ma et al., 2019), Salix brachista C. K. Schneider (Chen et al., 2019), Eutrema heterophyllum (W. W. Smith) H. Hara and E. yunnanense Franchet (Guo et al., 2018), Notopterygium incisum C. C. Ting ex H. T. Chang (Jia et al., 2019). For these studies, Gene Ontologies (GOs) of biotic and abiotic stress often provide novel insights into the local adaptation and species divergence in the extreme environments of the Qinghai-Tibet Plateau region.

During field investigation in 2017, we found a putative new species of Megacodon in Lushui county of north-west Yunnan at 1700 m a.s.l. This Lushui Megacodon has similar morphology and habitat as that of M. venosus; however, the Lushui Megacodon and M. venosus have extraordinarily disjunctive distributions. The new record of Megacodon from Lushui implies that the evolutionary history of this genus may be more complicated than previously thought. Two species co-occur in the Nushan Mountains, yet at different elevations. Furthermore, these observations suggest that the divergence of these Megacodon species may have been the result of both ecological divergence and allopatric speciation.

In this study, we described a new species of Megacodon from Lushui (Fig. 1) and determined its phylogenetic relationships within Megacodon. To identify genetic signatures of adaptive evolution to different elevations, we compared the transcriptome of M. stylophorus, which is distributed at high elevations, with that of the new species Megacodon lushuiensis, which is distributed at low elevations.

2. Materials and methods 2.1. Material collection

We sampled a total of 12 samples of M. stylophorus from different sites in Tibet and Yunnan, two samples of M. venosus from Chongqing, and two samples of the putative new species from Lushui, Yunnan (Fig. 1 and Table S1). Leaves were collected and dried in silica gel for DNA extraction. Specimens were deposited at the herbarium of Kunming Institute of Botany, Chinese Academy of Sciences (KUN). Flower organs of the Megacodon species from Lushui and M. stylophorus from Daxueshan Mountain were sampled by liquid nitrogen and stored at -80 ℃ for RNA-seq. Each sample was obtained from one individual.

2.2. DNA extraction, PCR amplification, and sequencing

The Plant Genomic DNA Kit (Tiangen Biotech, Beijing, China) was used to extract total genomic DNA following the manufacturer's protocol. The nuclear ribosomal internal transcribed spacer (nrITS) and two chloroplast markers, atpB-rbcL, and trnL-trnF, were amplified and sequenced following the protocols of previous studies (Favre et al., 2016; Matuszak et al., 2016). Primers used in this study for atpB-rbcL, trnL-trnF, and the nrITS region were the same as those used in previous studies (see Hoot et al., 1995; Taberlet et al., 1991; Grudinski et al., 2014 for the nrITS region). The 25 μL volume polymerase chain reactions (PCRs) contained 1.0 μL of genomic DNA (20 ng/μL), 0.3 μL of each primer (10 μM), 12.5 μL Ex Taq Version 2.0 (Takara, Dalian, China), and 10.9 μL of double-distilled H₂O. PCRs were performed in a C1000 Touch^TM Thermal Cycler (Bio-Rad, Hercules, CA, USA), with initiation at 95 ℃ for 2 min followed by 35 cycles of denaturation at 95 ℃ for 1 min, primer annealing at 53-56 ℃ for 1 min, primer extension at 72 ℃ for 1 min and a final extension step at 72 ℃ for 10 min. Alignments and subsequent manual adjustment were performed in BioEdit 7.2.5 (Hall, 1999) by using the multiple alignment software ClustalW (Thompson et al., 1994).

2.3. Phylogenetic analyses

Maximum likelihood (ML) and Bayesian inference (BI) approaches were used for phylogenetic tree reconstruction with combined sequences from nrITS, atpB-rbcL, and trnL-trnF. Three Potalia Aubl. and three Lisianthius P. Browne species were set as outgroups (Table S2). In addition, 11 Gentiana L. species and 29 Swertia species were included as closely related genera (Table S2). ML analyses were performed on the CIPRES Science Gateway (www.phylo.org; Miller et al., 2010) through RAxML-HPC2 v.8.2.8 (Stamatakis, 2014) with 1000 bootstrap replicates. The best-fitting nucleotide substitution model for BI analyses was GTR + G, selected using the Akaike information criterion (AIC) in jModeltest 2.1.7 (Darriba et al., 2012). BI analyses were also performed on CIPRES using MrBayes v.3.2.6 (Ronquist et al., 2012) on XSEDE with the following settings: ngen = 10, 000, 000, samplefreq = 1000, printfreq = 1000 and sumt burnin = 2500.

2.4. Climatic data analysis

Our fieldwork identified 15 precise locations where Megacodon is distributed. We used the 'raster' package to extract 19 bioclimatic variables for these sites from WorldClim (http://www.worldclim.org; Hijmans et al., 2005). Principal component analysis (PCA) of these environmental factors was performed by the dudi.pca function of the package 'ade4' in R 3.5.1 (R Core Team, 2018). We conducted paired two-tailed t-tests in Microsoft® Excel® 2016 (Microsoft Company, USA) to characterize factors that diverged between M. stylophorus and the other two Megacodon species (including the suspected new species from Lushui).

2.5. Transcriptome sequencing and functional annotation

TRIzol reagent (Sigmae-Aldrich) was used to extract total RNA of samples following the manufacturer's instructions. A total amount of 3 μg of RNA per sample was used to construct a cDNA library with a fragment length range of 250-300 bp. The cDNA libraries were created with NEBNext® Ultra^TM RNA Library Prep Kit for Illumina® (NEB, USA), following the manufacturer's recommendations. Illumina Hiseq 2500 was used for sequencing the library preparations, and paired-end reads were generated. Library preparation and Illumina sequencing were performed at Novogene Bioinformatics Technology Co., Ltd., Beijing, China (www.novogene.cn). Clean data was obtained by trimming read adapter sequences, and filtering low-quality reads from raw data. De novo assembly of each individual was accomplished using the Trinity program (Raychowdhury et al., 2011) with min_kmer_cov set to 2, and all other parameters set to default.

Functions of all assembled unigenes were annotated based on the following databases: Nr (NCBI non-redundant protein sequences); Nt (NCBI non-redundant nucleotide sequences); Pfam (Protein family); KOG/COG (Clusters of Orthologous Groups of proteins); Swiss-Prot (A manually annotated and reviewed protein sequence database); KO (KEGG Ortholog database); and GO (Gene Ontology).

2.6. Homologous gene identified and positive selection analysis

We used OrthoFinder software v. 2.2.7 (Emms and Kelly, 2018, 2015) to identify homologous gene clusters (orthogroups). OrthoFinder was run with multiple sequences alignment (MSA) for gene tree inference and RAxML for tree inference method (-T ramxl -M msa). Orthogroups of single-copy genes (1:1 orthologs) were obtained for substitution rates estimation. Each orthogroup was aligned by Prank (Löytynoja and Goldman, 2005) at the condon level (with options "-F -codon"). To estimate a general evolutionary pattern of selective pressure for the two species, codeml program from PAML 4.9 (Yang, 2007) was used to calculate the d_N, d_S and d_N/ d_S with the free ratio model (model = 2, NSsites = 0) and ML pairwise alignment (runmode = -2). To reduce false-positive findings in subsequent evolutionary analyses, genes with length less than 150 bp, d_S < 0.001 (to avoid a very large d_N/d_S ratio) and d_S > 0.1 (to avoid suspicious paralogs) were discarded (Bustamante et al., 2005; Sun et al., 2018; Zhang et al., 2013).

Positive selection genes (d_N/d_S > 1) were used as the test set for the GO enrichment analysis, which was performed on web-based KOBAS software (Xie et al., 2011) by Fisher's exact test, while whole single-copy genes were selected as the background set. Significant enriched GO terms (P < 0.05) were selected and further analyzed using REVIGO (http://revigo.irb.hr/; Supek et al., 2011) to investigate in-depth GO similarities.

3. Results 3.1. Morphological comparisons and molecular variation

Comparative analysis of morphological traits indicates that the Megacodon species from Lushui is closely related to M. venosus (Fig. 2 and Table 1). The Megacodon species from Lushui is similar to M. venosus by its narrow leaf blades shape, calyx shape and the existence of glands at the base of calyx, but differs from M. venosus by its taller and thinner plant, lower stem leaf blades ellipticspathulate to elliptic-obovate (vs. elliptic-lanceolate), middle to upper stem leaf blades veins 3-5 (vs. 5-7), 11-15-flowered thyrses (vs. 7-11 flowered), and dense glands at the 1/3 base of the calyx (vs. sparse glands at the 1/2 base of the calyx). Two individuals from Lushui and two individuals from M. venosus differed by a total of eleven nucleotide substitutions and six insertions/deletions in the alignment dataset (including ITS, atpB-rbcL, trnL-trnF, Table S3).

3.2. Phylogenetic relationships of Megacodon

The data set of phylogenetic analyses comprised 62 accessions with 2383 characters. Because the topologies of phylogenetic trees from ML and BI analyses were mostly identical, only the ML tree is shown with the posterior probabilities (PP) from BI analysis (Fig. 3). Phylogenetic inference showed that Megacodon is most likely a monophyletic group (ML-BS = 100/BI-PP = 1.00 in Fig. 3). Twelve accessions of M. stylophorus from Tibet and Yunnan formed a wellsupported clade (96/1.00). In comparison, two accessions of M. venosus from Chongqing and two accessions from Lushui formed the other clade of Megacodon. The population of M. venosus and the species from Lushui each formed a monophyletic clade in this branch.

3.3. Taxonomic treatment

According to both morphological and molecular evidence, as well as the extraordinary disjunctive distribution with M. venosus, we propose the Megacodon plants from Lushui as a new species here.

M. lushuiensis J. C. Peng & H. Sun sp. nov.(Fig. 2A-E, and Fig. 4).

Holotype. China, Yunnan Province: Laowo Xiang, Lushui County, on the hillside, 1700 m a.s.l, 25.85°N, 99.03°E, flowering, 21 September 2017, J. C. Peng & L. Sun PM001 (KUN).

Paratype. The same locality, flowering and fruiting, 7 November 2017, J. C. Peng, L. Sun & Z. Chen PM002 (KUN).

Description. Plant perennial, herbaceous ca. 130-240 cm tall. Basal to middle stems erect, terete, 1.1-1.8 cm in diam at base; upper stems slim, 3-9 mm in diam. Lower stem leaves ellipticspathulate to elliptic-obovate, base cuneate and decurrent, veins 5-7 and arcuate; middle stem leaves elliptic-lanceolate to lanceolate, apex acuminate, 15-45 × 5-14 cm, veins 3-5 and arcuate; upper stem leaves linear-lanceolate to linear, 5-10 × 1.1-1.5 cm, veins 3-5 and arcuate. Thyrses 11-15-flowered. Pedicel 3-9 cm. Calyx campanulate, 1.7-3.6 cm, tube 4-6 mm; lobes ellipticlanceolate, apex acuminate, dense pustulous glands at the base ca. 1/3. Corolla white to pale yellow, with green veins, ca. 5-7 × 4-5 cm, tube 8-10 mm; lobes oblong-spatulate, 4-5.8 cm, apex rounded. Filaments flat and linear, 2.2-2.8 cm, alternate with lobes; anthers oblong-quadrate, 8-11 mm; Capsule ovoid-ellipsoid, 1.7-2.0 cm.

Phenology. Flowering and fruiting from Sep. to Nov.

Distribution and habitat. M. lushuiensis was found in Laowo Xiang of Lushui County, north-west Yunnan, south-west China (Fig. 1), growing on hillsides (north-west slope) and under shrub at ca. 1700 m a.s.l.

Conservation status. M. lushuiensis is only known from one population in Lushui county (criterion B2a), with an area of occupancy estimated to be less than 500 km² (criterion B2b (i, ii)). The habitat of this population is near to cultivated field (less than 2 km) and suffering from invasive plants such as Ageratina adenophora (criterion B2b (iii)). The new species should be regarded as "Endangered" based on the IUCN (2012) criteria.

Etymology. The specific epithet refers to the type locality, Lushui County of Yunnan Province, in China.

3.4. Climatic niche

PCA of 19 bioclimatic variables indicated clear climate niche differentiation between 13 M. stylophorus sites, one M. venosus site, and one M. lushuiensis site (Fig. 5). The first two PCs identified two climatic groups and explained 59% and 26.2% of the total variation, respectively (Fig. 5). The climatic niche of M. lushuiensis is closer to the population of M. venosus from Chongqing, whereas 13 sites of M. stylophorus formed another group. Among these 19 bioclimatic variables, 12 bioclimatic variables (BIO1, BIO5, BIO6, BIO8-BIO12, BIO14, BIO15, BIO17, and BIO19) showed significant differences between M. stylophorus and the other two species (P < 0.01; Table S4).

3.5. Overview of RNA-seq analysis

We generated 55.2 and 55.3 million clean reads of RNA-seq data for M. lushuiensis and M. stylophorus, with a total size of 8.28 and 8.30 Gb, respectively (NCBI GenBank accession number: SRR11743766, SRR11743765 and Table 2). The median unigene length for M. lushuiensis after de novo assembly was 981 bp, with an N50 length of 2043 bp. The median unigene length for M. stylophorus was 1018 bp, with an N50 length of 2058 bp. A total of 24, 086 unigenes from M. lushuiensis and 23, 262 unigenes from M. stylophorus were annotated in at least one database with a significant match (e-value < 10^-5).

3.6. Orthogroups identified and positive selection analysis

In this study, a total of 20, 679 and 20, 471 unigenes were detected in M. lushuiensis and M. stylophorus, respectively, and 23, 611 orthogroups were identified. Among these, 8926 orthogroups contained putative one-to-one single-copy genes between two species. To detect genes that are involved in ecological adaptation to high elevations, d_N, d_S, and d_N/d_S ratio of single-copy genes were calculated by PAML. A peak of d_S distribution between M. stylophorus and M. lushuiensis was observed at 0.0411 ± 0.0210 (Fig. 6A). The d_N/d_S of most genes was less than one (Fig. 6B), which indicates purifying selection. Only 370 single-copy genes were identified as positively selected genes (PSGs, ω > 1) between two species.

In enrichment analyses, the 370 positively selected genes were designated as the test sets, and 7749 single-copy genes (0.001 < d_S < 0.1) were designated as the background set. Twentyfive GO categories were annotated to 49 pairs of positively selected genes, which were overrepresented (Fisher's exact test, p < 0.05) in the test sets (Table 3) by GO enrichment analysis in the KOBAS website. Then based on GO similarities, five clusters of biological processes were shown in the REVIGO treemap (Fig. 6C). A majority of biological processes in treemap were found to be related to specific adaptation traits, such as response to water (11 PSGs), regulation of response to osmatic (3 PSGs), cellular response to nutrient levels (7 PSGs), and cellular response to external stimulus (7 PSGs) (Fig. 6C, Table 3 and Table S5).

4. Discussion 4.1. New insights into the evolutionary history of Megacodon

In this study, we described a new species of Megacodon, M. lushuiensis, and determined the phylogenetic relationships among Megacodon species. The three existing Megacodon species formed a well-supported monophyletic clade (Fig. 3). All samples of M. stylophorus formed a monophyletic clade. Furthermore, in our phylogenetic tree, M. lushuiensis is closely related to M. venosus, which is consistent with the morphological similarity between these two species. The phylogenetic relationship between Megacodon and four other lineages (Potalia, Lisianthius Gentiana, and Swertia) is congruent with previous studies (Favre et al., 2016). Although sharing many morphological similarities, M. lushuiensis can be easily distinguished from M. venosus by its leaf blade, inflorescence, extrafloral nectaries, and more. Importantly, the extrafloral nectaries in the basal portion of the calyx were neglected in previous morphological descriptions of M. venosus. The extrafloral nectaries of M. venosus and M. lushuiensis not only secrete nectar but also cover the petal and make it sticky (Fig. 2D, E, I, and J). We speculate that in addition to attracting ants to participate in ant-plant mutualisms, the secretions from M. venosus and M. lushuiensis extrafloral nectaries may directly participate in the protection of the flowers (Marazzi et al., 2013).

The exact distribution of M. lushuiensis in north-west Yunnan indicates that the most recent common ancestor of M. lushuiensis and M. venosus may occupy a more extensive area, from west Yunnan to central China. In the past, Megacodon populations from the Nujiang River valley may have been connected to M. venosus populations from central China through the eastern edge of the Qinghai-Tibet Plateau region or the Yunnan-Guizhou Plateau. The current fragmented distribution of this species pair may be the result of a vicariance event closely related to the climatic changes of the Sino-Himalayan region and its low reproductive capacity. Because of its humid climate and large elevational gradient, the Nujiang River valley harbors many distinct species that are identical to or sister to species in central or south China, such as Ichtyoselmis macrantha (Oliver) Lidén, Ypsilandra Franch., Asteropyrum peltatum (Franch.) Drumm. ex Hutch. The large elevational gradient and complex climate conditions in Nushan Mountain should provide more ecological opportunities for the maintenance of relict populations or species. Despite this, the fragmented distribution pattern of M. venosus and M. lushuiensis is still rare and may have important significance in biogeographical studies.

The occurrence of M. lushuiensis in north-west Yunnan also alters our understanding of the divergence of M. stylophorus and related species. Speciation of M. venosus and M. stylophorus was previously thought to have occurred through allopatric speciation, similar to that of Pleurospermum hookeri C. B. Clarke and P. giraldii Diels, which have a fragmented distribution between central China and the Qinghai-Tibet Plateau region (Bai et al., 2015). However, the fragmented distribution of M. lushuiensis and M. stylophorus occur at different elevations along the Nushan Mountain, suggesting that ecological differentiation may have played an important role in their divergence.

4.2. Functions under positive selection between M. lushuiensis and M. stylophorus

To further investigate why M. lushuiensis and M. stylophorus occupy distinct ecological niches along an elevational gradient, we used comparative transcriptome analysis. The majority of genes identified by GO enrichment of M. lushuiensis and M. stylophorus function in response to water deprivation (GO:0009414) and other external stimuli (e.g., GO:0047484 and GO:0071496) (Fig. 6C). These findings suggest that genes that respond to the environment may have been triggered by natural selection during ecological adaptation (Table 3). However, unlike many other alpine plants, especially alpine scree plants (Chen et al., 2019; Guo et al., 2018; Jia et al., 2019; Zhang et al., 2019), M. stylophorus does not appear to have adapted to high radiation. This may be because M. stylophorus grows in the forest or scrub habitat. Analysis of climatic variables revealed that M. stylophorus favors the colder and higher climatic drought conditions of high elevations (Fig. 5 and Table S1). In contrast, field observations indicate that M. lushuiensis prefers to grow in drier soil conditions such as semi-humid hillsides at low elevations (Fig. 2A). Habitat differentiation between these two species may have led some populations to experience water deficiency at different periods, which can influence plant growth and development. Our transcriptome analysis identified several genes (e.g., PGGT-1, ALDH10A9, and EDL3) under positive selection that is directly involved in water deprivation and hormone signaling pathways (Johnson et al., 2005; Koops et al., 2011; Missihoun et al., 2014). DRIP1 is a negative regulator in drought toleranceresponsive gene expression that was also found to be under positive selection (Qin et al., 2008). Our study also found that the following genes were under positive selection: LEW1, which response to endoplasmic reticulum stress, drought, and darkinduced senescence by catalyzing dolichol biosynthesis (Zhang et al., 2008); RCD1, which is involved in the regulation of stomatal conductance and the flavonoid accumulation pathway (Ahlfors et al., 2004; Suvorov et al., 2017); and CSDP1, which contains a cold shock domain and is involved in cold and salt stress responses (Yang and Karlson, 2013). The functions of these homologous genes emphasize that the genetic pathways that regulate efficient water use may have been important drivers of the divergence of M. lushuiensis and M. stylophorus.

Under a short growing season in the alpine region, M. stylophorus also require proper phenological responses. Our transcriptome analysis identified several genes that may play roles in these phenological responses. For example, in addition to its response to water deprivation, CSDP1 can also regulate seed germination (Yang and Karlson, 2013). Furthermore, bZIP53 is an enhancer for seed maturation gene transcription (Alonso et al., 2009), while BZOZH3 (also known as bZIP63) is involved in the regulation of a circadian oscillator gene PRR7 (Frank et al., 2018). The correct circadian clock guarantees that plants have photosynthetic advantages and can accurately anticipate environmental changes (Dodd et al., 2014, 2005). These homologous genes may strengthen the ability of M. stylophorus to respond to habitatspecific variability. Similarly, EX1 and FRD3, which were identified in our analysis, may help M. stylophorus tolerate higher oxidative stress and iron deficiency that result from radical-initiating factors at high elevations (Dogra et al., 2019; Lei et al., 2014), especially at the end of the growing period (Sakata et al., 2006). Overall, adapting to extreme environments (such as climatic drought, oxidative stress, and low temperature) in Sino-Himalayan regions may play an important role in the speciation of Megacodon.

5. Conclusion

In this study, we used field observations in conjunction with morphological and molecular data to propose a new Megacodon species from the Nujiang River valley, M. lushuiensis. The discovery of this species calls attention to the importance of ecological divergence between M. stylophorus and low-elevation Megacodon species. We also presented the phylogenetic position M. venosus for the first time and described a neglected trait in this species (sparse glands at 1/2 base of the calyx). Although population genetic and genomic studies are essential to understanding the demographic and divergence history of Megacodon, our analysis of the transcriptomes of M. stylophorus and M. lushuiensis provide valuable resources for the future genomic research and general insight into how ecological divergence occurred between these two Megacodon species. M. venosus and M. lushuiensis are both extremely rare in the field; we, therefore, recommend that conservation efforts be carried out timely and properly for these two species.

Author contributions

HS, XGM and YHW conceived and designed experiments. JCP performed field investigations, analyzed the data, and wrote the manuscript. XGM and HS revised the manuscript.

Declaration of Competing Interest

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

Acknowledgments

This study was supported by the Key Projects of the Joint Fund of the National Natural Science Foundation of China (U1802232), the Second Tibetan Plateau Scientific Expedition and Research (STEP) program (2019QZKK0502). Thanks to Si-Rong Yi (Chongqing Three Gorges Medical College) and Wei Zhang for providing valuable distribution information, Sha Wen for her line drawing of the new species, Ti-Cao Zhang (Yunnan University of Chinese Medicine) and Yan-Ting Hu for useful advice on RNA-seq analysis. We also appreciate the help of Jian Wen Zhang, Zhuo Zhou, Yang Niu, Yi Yang, Lu Sun, Zhe Chen, Bo Xu (Southwest Forestry University), Hong-Liang Chen, Li-Shen Qian, and Yong-Zeng Zhang during fieldwork.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2020.05.003.