• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Phylogenetic relationships of the zokor genus Eospalax(Mammalia, Rodentia, Spalacidae) inferred from wholegenome analyses, with description of a new species endemic to Hengduan Mountains

    2022-06-07 10:50:06TaoZhangMengLongLeiHaoZhouZhongZhengChenPengShi
    Zoological Research 2022年3期
    關(guān)鍵詞:鼢鼠木里

    Tao Zhang, Meng-Long Lei, Hao Zhou, Zhong-Zheng Chen,4,*, Peng Shi,5,6,*

    1 State Key Laboratory of Genetic Resources and Evolution, Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming,Yunnan 650204, China

    2 Jiangsu Key Laboratory of Neuropsychiatric Diseases and College of Pharmaceutical Sciences, Soochow University, Suzhou, Jiangsu 215123, China

    3 Joint Laboratory of Animal Models for Human Diseases and Drug Development, Soochow University and Kunming Institute of Zoology,Chinese Academy of Sciences, Kunming, Yunnan 650223, China

    4 Collaborative Innovation Center of Recovery and Reconstruction of Degraded Ecosystem in Wanjiang Basin Co-founded by Anhui Province and Ministry of Education, School of Ecology and Environment, Anhui Normal University, Wuhu, Anhui 241000, China

    5 School of Future Technology, University of Chinese Academy of Sciences, Beijing 100049, China

    6 Center for Excellence in Animal Evolution and Genetics, Chinese Academy of Sciences, Kunming, Yunnan 650223, China

    ABSTRACT Zokors in the genus Eospalax, which are endemic to northern and western China, are subterranean rodents that inhabit various niches, including grasslands, high-altitude meadows, forests, and farmlands.Six species in Eospalax were described a century ago but their taxonomy and phylogeny remain controversial.In this study, we performed high-depth whole-genome sequencing of 47 zokor samples, comprising all six previously described species.Genomic analyses revealed a reliable and robust phylogeny of Eospalax and supported the validity of the six named species.According to the inferred phylogenetic relationships, Eospalax first divergent into two clades in the early Pliocene (ca.4.68 million years ago (Ma)), one inhabiting the highaltitude Qinghai-Xizang (Tibet) Plateau (QTP) and adjacent regions, and the another inhabiting the lowaltitude Loess Plateau and Qinling-Daba Mountains.The most recent divergences occurred between E.baileyi and E.smithii and between E.rufescens and E.rothschildi in the late Pliocene (ca.2.09 and 2.19 Ma, respectively).We also collected specimens of zokors in the southern Hengduan Mountains (Muli County, Sichuan Province), far from the known distributions of all other zokors.Morphological andmolecular analyses strongly suggested that the specimens represent a new species, formally described here as Eospalax muliensis sp.nov.The new species belongs to the high-altitude clade and diverged from closely related species (ca.4.22 Ma)shortly after the first divergence in Eospalax.Interestingly, Eospalax muliensis sp.nov.possesses more supposedly plesiomorphic characters, suggesting a possible origin of the genus in the Hengduan Mountains.

    Keywords: Zokor; Eospalax; Phylogenomic analyses; New species; Hengduan Mountains

    INTRODUCTION

    Zokors (subfamily Myospalacinae) are a group of strictly subterranean rodents endemic to East Asia, which belong to the family Spalacidae and are closely related to the Rhizomyinae and Spalacinae subfamilies (Bouckaert et al.,2019; Guo et al., 2021; Lin et al., 2014; Norris et al., 2004).Zokors are cylindrical in shape, with stocky and strong forelimbs and hardened nose pads to excavate their complex burrow systems.Based on morphological and molecular evidence, extant zokors are classified into two genera: i.e.,Myospalax, which includes species with flat occiputs, andEospalax, which includes species with convex occiputs (Su et al., 2014; Zheng, 1994; Zhou & Zhou, 2008).Myospalaxspecies are mainly distributed in the grasslands of northeast China, Mongolia, and Russia, whileEospalaxspecies inhabit various niches in northern and western China (Figure 1A) (Fan& Shi, 1982), including the high-altitude Qinghai-Xizang (Tibet)Plateau (QTP).The phenotypic and genetic adaptations of the plateau zokor (E.baileyi) to the hypoxic conditions on the QTP have been well researched (Wei et al., 2006; Xu et al., 2021;Zhang et al., 2021).

    However, despite studies on the physiology, ecology, and genetics of certain species, the taxonomy and phylogeny ofEospalaxremain controversial.Six species ofEospalaxwere described in the late 19thand early 20thcentury, but the taxonomic status ofEospalaxspecies has changed with time,especially forE.baileyi,E.cansus, andE.rufescens(Table 1).Eospalax fontanierii,E.rothschildi, andE.smithiiwere recognized as distinct species by Allen (1940), with this classification adopted by many monographs until recently(Smith & Xie, 2008; Wilson & Reeder, 2005).Based on comprehensive morphological comparisons, Fan & Shi (1982)found thatE.cansusis stably smaller thanE.fontanieriiin body size and skull andE.baileyiis distinct fromE.fontanieriiin nose pad shape and presence of a postero-external lobe on M3.Thus, they elevatedE.baileyiandE.cansusto species level.Subsequent studies using mitochondrial DNA markers have generally validated the status of all sixEospalaxspecies(He et al., 2012; Su et al., 2014; Zhou & Zhou, 2008).However, the classification ofE.baileyihas remained somewhat controversial due to insufficient molecular markers and incomplete sampling (Jiang et al., 2017; Wei et al., 2021).In addition, the potential existence of new undescribed species within theEospalaxgenus may complicate established taxonomy.Here, based on a literature review of species distribution, we found reports of zokors in the southern Hengduan Mountains (Muli County, SichuanProvince) (Liu et al., 2007), far from the known distribution of all other zokors (Figure 1A).Thus, in the current study, we collected 14Eospalaxspecimens from Muli County between October 2020 and October 2021.Externally, these zokors differed considerably from all other known species due to their small size and long tail (Figure 1B), indicating that they may represent a distinct taxon, which we name herein asEospalax muliensissp.nov.

    Table 1 Previous and proposed classifications of Eospalax genus

    Figure 1 Sampling of species in Eospalax

    In addition to taxonomic controversies, the phylogeny ofEospalaxremains incomplete and inconsistent.Conflicting phylogenies amongEospalaxspecies have been reported in previous studies based on both mitochondrial DNA markers(Cai et al., 2020; Su et al., 2014; Zhou & Zhou, 2008; Zou et al., 2020) and morphological traits (Li & Chen, 1986; Li &Wang, 1996).Notably, phylogenetic relationships inferred from mitochondrial markers are unstable when different markers or tree-construction methods are applied (Su et al., 2014; Zhou &Zhou, 2008).Therefore, the ambiguous and conflicting phylogenetic relationships reported inEospalaxhamper studies of this remarkable lineage and underscore the need for a reliable and robust phylogeny.

    In the current study, we aimed to clarify the taxonomy and phylogeny ofEospalax.All six known species ofEospalaxwere sampled and whole genomes were sequenced to high depth.Whole-genome analysis revealed reliable and robust phylogenetic relationships and validated the status of all six previously described species.Furthermore, by integrating morphological and genomic evidence, the zokor specimens collected from Muli County were found to represent a new species.Thus, the phylogeny and diversification of the new species were also studied.

    MATERIALS AND METHODS

    Sample collection

    Fieldwork was carried out from April 2017 to October 2021.Skeletal muscle was collected for DNA extraction.All collections followed the animal use protocols approved by the Animal Care and Ethics Committee of the Kunming Institute of Zoology (KIZ), Chinese Academy of Sciences (CAS)(Approval No.SMKX-SQ-2021-067).

    Sampling sites of each species were located in their generally acknowledged distribution areas (Figure 1A;Supplementary Table S1).The seven sampling sites ofE.baileyiwere reported in our previous work (Zhang et al.,2021).Because the distribution ofE.cansusis wide and adjacent to several other species inEospalax, three sampling sites ofE.cansuswere chosen across Shaanxi, Ningxia, and Gansu provinces to ensure representativeness.For other species, one or two sampling sites were selected.Two zokors collected near Hulun Lake (Inner Mongolia Autonomous Region, China) were recognized asMyospalax psilurusaccording to their skull and hairless tails.Thus, these two zokors were used as an outgroup.

    Whole-genome sequencing and de novo assembly-based variant calling

    In this study, 47 zokor samples were collected and whole genomes were sequenced to a high depth, with an average of 137.2 Gb of data for each sample (Figure 1A; Supplementary Table S1).Combined with previously reported whole-genome sequences of 38E.baileyisamples (Zhang et al., 2021), our dataset covered all six previously describedEospalaxspecies,as well asEospalax muliensissp.nov.and two samples of closely relatedM.psilurus.

    DNA extraction, resequencing library construction, and quality control of raw sequencing reads were performed following Zhang et al (2021).Briefly, the resequencing library was constructed following standard protocols and sequenced using Illumina HiSeq or BGISEQ paired-end mode.Raw sequencing reads were trimmed using Trimmomatic v0.36(Bolger et al., 2014) to remove low-quality reads.Although aE.baileyigenome assembly was recently published (Zhang et al., 2021), all other zokor genomes are currently unavailable.Thus, we sequenced all samples to high depth and appliedde novoassembly-based variant calling following the FermiKit pipeline (Li, 2015b).For comparability, variant calling of allE.baileyisamples was also applied using this pipeline.First,clean reads of each sample were error corrected using BFC v.r181 (Li, 2015a).Second, Fermi2 v.r178 (Li, 2012) was used tode novoassemble the genome of each sample.Third, the assembly of each sample was mapped to the plateau zokor genome using minimap2 v.2.17 (Li, 2018) with parameters “-ax asm5/asm10/asm20 -k 23 -w 11 -f 1 000 -n 2 -s 100 --endbonus 5”.According to the genetic distance between samples and plateau zokors, different parameters of -ax were applied to balance the sensitivity and specificity of mapping.Specifically, asm5, asm10, and asm20 were applied forsamples ofE.baileyi,E.smithii, and all other species,respectively.After mapping, all sample variants were called jointly using the “pileup” functions in HTSbox v.r310 with parameters “-d -V 0.15 -S 300 -q 20 -Q 3 -s 5”.Raw genotypes were first filtered according to allele balance.For homozygous genotypes, at least 90% of reads should support these genotypes.For heterozygous genotypes, the two types of alleles should have 25%-75% of read-support simultaneously.If the above criteria were not satisfied, the genotypes were set as missing.After allelic balance filtering,we only retained those biallelic single nucleotide variants(SNVs) with a missing rate of no more than 20% using VCFtools v.0.1.17 (Danecek et al., 2011).

    Phylogenomic analyses

    Phylogenetic relationships among all zokor samples were inferred using whole-genome SNVs and maximum-likelihood(ML) algorithms.Considering the computational resources required, we pruned high-quality SNVs by randomly choosing one SNV in each 100 bp window across the genome.ML analyses were conducted in RAxML v.8.2.12 (Stamatakis,2014) with the nucleotide substitution model of GTRGAMMA.To test the reliability of inferred phylogenetic relationships, 100 bootstraps were applied.

    In addition to the construction of phylogenetic trees based on whole-genome pruned SNVs, we also applied a summary of sliding-window ML trees.First, we split the whole genomes into over 50 000 non-overlapping windows (each 50 kb long).For each window, ML trees were constructed using RAxML based on all high-quality SNVs.Weighting and summary of all sliding-window ML trees across whole genomes were explored using topology weighting by iterative sampling of subtrees (Twisst) (Martin & van Belleghem, 2017).

    Divergence time among all species in the genusEospalaxwas determined using SNAPP v2.6.2 (Bryant et al., 2012) in BEAST2 (Bouckaert et al., 2019).Secondary calibration was applied following He et al (2020), which estimated that the divergence betweenMyospalaxandEospalaxoccurred 8.8 million years ago (Ma), with a 95% highest posterior density(HPD) interval from 7.43 to 10.23 Ma.The estimate from He et al (2020) is congruent with fossil records of the genusMyotalpavus, which occurred in the Late Miocene and may be the ancestor of bothEospalaxandMyospalax(Zheng, 1994).Thus, we constrained the divergence time betweenMyospalaxandEospalaxwith a lognormal distribution centered at 8.8 Ma and a standard deviation of 0.085 in SNAPP analysis.For each species, the individual with the lowest proportion of missing data was used to represent that species.The SNVs were first pruned as described previously, with those identified as monomorphic or containing missing data then filtered.Finally, 100 000 randomly selected SNVs were used, and 10 million Markov Chain Monte Carlo (MCMC) iterations were set for SNAPP analysis.Convergence was assessed using Tracer v1.6 (Rambaut et al., 2018) with 10% burn-in.DensiTree v.2.6.2 (https://www.cs.auckland.ac.nz/~remco/DensiTree/)was used to visualize all posterior trees sampled by SNAPP.TreeAnnotator v2.6.2 (https://beast.community/treeannotator)was used to quantify the posterior probabilities of clades as node support in a maximum clade credibility tree with 10%burn-in.Divergence time and phylogenetic relationships were visualized in FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).SNAPP analysis was implemented three times independently and obtained nearly identical results.

    Absolute genetic distance (DXY) among species was estimated in 50 kb windows with a 50 kb step size along the genome using popgenWindows.py in genomics_general(https://github.com/simonhmartin/genomics_general).We excluded windows located in scaffolds <1 Mb due to possible poor assembly of short scaffolds.

    Mitochondrial genome (mitogenome) analysis

    For all whole-genome sequenced samples, we assembled their mitogenomes.We randomly selected 5 Gb of resequencing reads for each sample.Assembly was performed using MitoZ v.2.3 (Meng et al., 2019) or GetOrganelle v.1.7.2b (Jin et al., 2020) with default parameters.The mitogenome of a publishedM.psilurussample (JX014234) (Li et al., 2016) served as a reference in GetOrganelle analysis.Annotation of the assembled mitogenomes was performed with MitoZ.The quality of each mitogenome was carefully checked based on length (>16 kb),gene count (37 genes), and circularity.

    In addition to the mitogenomes assembled in this study, we downloaded all other mitogenomes inEospalaxandMyospalaxreported in previous studies.Thirteen proteincoding genes and two ribosomal RNA (rRNA) genes were included in the following analyses.Gene alignments were performed by CLUSTAL v.2.1 (Larkin et al., 2007) and then concatenated.The best partition schemes for the dataset were identified using PartitionFinder v.2.1.1 (Lanfear et al., 2017)with recommended corrected Akaike Information Criterion(AICc).This analysis partitioned 15 mitochondrial genes into eight subsets.Mitochondrial phylogenomic analyses were inferred using RAxML with 1 000 bootstrap replicates.

    Morphological data and analyses

    A total of 62 specimens ofEospalaxwere collected.All specimens and tissue samples were deposited in KIZ, CAS.Specimens were identified following Fan & Shi (1982) and Luo et al (2000).Forty-eight specimens were assigned to the six recognized species, includingE.baileyi(n=34),E.cansus(n=5),E.fontanierii(n=1),E.rothschildi(n=2),E.rufescens(n=2), andE.smithii(n=4).Fourteen specimens collected from Muli County, Sichuan Province, were assigned to the putative new species (Eospalaxmuliensissp.nov.).

    The body weight (BW) of each specimen was weighed up to 0.01 g using an electronic scale.External measurements,including head and body length (HB), tail length (TL), and hindfoot length (HF), of some individuals were taken in the field using a ruler to the nearest 0.1 mm.Twenty-one cranial and dental measurements were taken using a digital caliper graduated to 0.01 mm following Pan et al (2007).These measurements included greatest length of skull (GLS),condylobasal length (CBL), basal length (BL), palatal length(PL), zygomatic width (ZMW), interorbital breadth (IOB),foramen infraorbital breadth (FIB), rostrum width (RSW),mastoid width (MTW), maximum width across upper second molars (M2-M2), length of upper tooth row (LUTR), length of upper molars (LUM), nasal length (NSL), braincase height(BCH), greatest breadth of foramen magnum (GBFM), length of auditory bulla (LAB), distance between auditory bulla(DAB), mandibular length (MDL), length of below toothrow(LBTR), length of lower molars (LLM), and length of lower incisor (LLI).We performed principal component analysis(PCA) to analyze morphometric variation of these specimens in SPSS v19.0 (SPSS Inc., USA).PCA was only based on the 21 craniodental measurements of adult specimens.All craniodental variables were log10-transformed before PCA.We considered any component with an eigenvalue exceeding 1.0 to be meaningfully interpretable.We further examined specimen morphology as per Fan & Shi (1982), Song (1986),and Luo et al (2000), and followed their terminologies for morphological descriptions.

    RESULTS

    Phylogenetic relationship and molecular dating

    A total of 85 whole-genome sequenced samples were included in this study (Figure 1A; Supplementary Table S1),covering all six previously describedEospalaxspecies as well asEospalax muliensissp.nov.and two samples of closely relatedMyospalax psilurus.

    After stringent quality control of genome alignment and variant calling, ~227 million high-quality SNVs were identified.Based on the genome-wide SNVs, ML analysis yielded a highly resolved phylogeny.All nodes that separated species within the genusEospalaxexhibited 100% bootstrap support(Figure 2A).Accordingly, species inEospalaxcould be divided into two major clades (i.e., high- and low-altitude clades).The high-altitude clade, comprised ofE.baileyi,E.smithii, andEospalax muliensissp.nov., generally inhabited the highaltitude QTP and adjacent regions (2 700-4 300 m a.s.l.),while the low-altitude clade, comprised of the other fourEospalaxspecies, generally inhabited the relatively lowaltitude Loess Plateau and Qinling-Daba Mountains.Within the high-altitude clade,Eospalax muliensissp.nov.was found in a basal position and inhabited extremely high altitudes (3 700 m a.s.l.).Within the second clade,E.fontanieriisplit first from the other species.Notably, threeE.cansuspopulations located remotely from each other formed a tight monophyletic cluster.

    The genetic distances between theM.psilurusoutgroup and species inEospalaxranged from 3.5% to 3.7%(Figure 2B).Within the genusEospalax, the genetic distances among species ranged from 1.0% to 2.2% (Figure 2B).For each species pair separated by the same internal node, their genetic distances were comparable (Figure 2B),demonstrating the reliability of the inferred phylogenetic relationships and that the applied methods enabled phylogenomic analysis.

    Phylogenetic relationships inferred from the mitogenomes were also constructed (Supplementary Figure S1).Consistent with nuclear genomic analyses,Eospalax muliensissp.nov.was also identified as a distinct genetic lineage based on mitogenomic analyses (Supplementary Figure S1).However,significant incongruences among phylogenies inferred from nuclear genomes and mitogenomes were found, especially in the positions ofEospalax muliensissp.nov.andE.fontanierii(Figure 2A; Supplementary Figure S1) (see Discussion).

    To confirm the inferred phylogenetic relationships and investigate the alternative topology, we split the genomes into>50 000 windows (50 kb in size) and constructed window trees based on the corresponding SNVs in each window.Most window trees (75.8%) were consistent with the whole-genome-based phylogenetic tree (Figure 3A, B).In addition,four alternative topologies were supported by at least 1% of windows (Figure 3A).The top two alternative topologies differed from the majority tree in the position ofEospalax muliensissp.nov.(Figure 3C,D).The other two lesssupported topologies revealed an uncertain relative position ofE.fontanieriiandE.cansuswithin the second clade(Supplementary Figure S2).

    Figure 2 Phylogenetic relationships among species of Eospalax

    Figure 3 Alternative topologies and divergence times of genus Eospalax

    In divergence time analysis, all posterior trees sampled by SNAPP were identical to the whole-genome-based phylogenetic tree (Supplementary Figure S3).The two major high- and low-altitude clades diverged from each other ~4.68 Ma (95% CI=3.94-5.48 Ma) (Figure 3E).Shortly after,Eospalax muliensissp.nov.diverged fromE.baileyiandE.smithii(~4.22 Ma; 95% CI=3.57-4.97 Ma).Within the lowaltitude clade,E.fontanieriifirst diverged from other species~3.87 Ma (95% CI=3.26-4.54 Ma).Thus, the speciation events within the genusEospalaxgenerally occurred in the Pliocene (Figure 3E).

    Morphometric comparison and morphological diagnosis

    The external and craniodental measurements of each species are given in Table 2.A total of 57 intact skulls of adult specimens were used for PCA.The first three principal components (PC1-3) had eigenvalues greater than 1.0, and together accounted for >79.78% of total variance(Supplementary Table S2).PC1 accounted for 44.08% of the variation (eigenvalue=9.26) and was strongly loaded with most variables, except IOB, GBFM, DAB, and M2-M2, indicating it may be a marker of size.PC2 accounted for 21.57% of the variation (eigenvalue=4.53) and was strongly dominated by RSW, IOB, FIB, and LLI (loading>0.62), while PC3 accounted for 14.13% of the variation (eigenvalue=2.97) and was dominated by GBFM, DAB, and M2-M2.

    Table 2 Body weight and external and cranial measurements (mm) (mean±SD and ranges) of Eospalax specimens (n) examined in this study

    The PC1-2 plot (Figure 4) showed thatEospalax muliensissp.nov.occupied the lower-left corner and was well separated from the other species.Specimens ofEospalax muliensissp.nov.were plotted on the negative field of both PC1 and PC2, reflecting its generally smaller size, narrower RSW, IOB, and FIB, and shorter LLI.Specimens ofE.baileyioccupied a large area, from -1.34 to 1.45 in PC1 and -0.49 to 1.70 in PC2, indicating great variation in size, with a relatively wide RSW, IOB, and FIB, and long LLI.Specimens ofE.cansus,E.fontanierii,E.rufescens, andE.smithiiwere mostly plotted on the positive region of PC1, indicating relatively large skulls.

    Both molecular and morphological evidence supportedEospalax muliensissp.nov.as a new species, which is formally described below.

    Taxonomic account

    Eospalax muliensis Zhang, Chen & Shi, sp.nov.

    Suggested common name: Muli zokor; 木里鼢鼠(Muli Fenshu)

    Holotype:KIZ 040324 (field No.: SC 2110385).Adult female collected by Tao Zhang, Meng-Long Lei, Zhong-Zheng Chen,and Hao Zhou on 21 April 2021.Dried skin, cleaned skull, and alcohol-preserved carcass are deposited in KIZ.Measurements (mm): BW=133.00 g, HB=172.0, TL=56.5,HF=28.0, GLS=39.41, CB=37.46, BL=35.40, PL=27.27,ZMW=25.58, IOB=7.27, FIB=6.12, RSW=8.84, MTW=22.95,M2-M2=9.30, LUTR=23.93, LUM=9.11, NSL=14.65,BCH=15.84, GBFM=6.25, LAB=8.33, DAB=4.12, MDL=24.75,LBTR=20.36, LLM=9.63, LLI=10.51.

    Type locality:Kangwu Ranch (N28.135°, E101.196°), Muli County, Sichuan Province, China, adjacent to Shangri-La Lake at an altitude of ~3 700 m a.s.l.

    Paratypes:Seven specimens were collected from the type locality.Four specimens (KIZ 040325-040328) were deposited in KIZ, CAS; three specimens (SC 2110384,2110386, and 2110451) were deposited in the Biological Museum of Anhui Normal University.All specimens were prepared as dried skins with cleaned skulls and alcoholpreserved carcasses.

    Etymology:The species namemuliensisis derived from MuliCounty, the type locality of the new species in Sichuan Province, China, and the Latin adjectival suffix -ensismeans“belonging to”.

    Diagnosis:Eospalax muliensissp.nov.can be distinguished from other described species ofEospalaxby a combination of the following characters: size small (mean BW=155.64 g;mean GLS=40.52 mm); tail relatively long, densely hairy;nose-pad trifoliate; rostrum weak, nearly rectangular, nasals small, posterior border of nasals transverse, anterior halves of premaxilla nearly aligned with nasals; braincase well-domed,temporal ridges not conspicuous, parallel in front, lambdoid crests only present on side of skull; occipital shield with developed occipital ridges, extending well posteriorly, almost forming plane with occipital condyle; about 1/3 of incisive foramina included within maxillae, remaining 2/3 of incisive foramina included within premaxillae; palate and pterygopalatine fossa large, M2-M2almost equal to LUM; M3with two reentrant angles on outer side.

    Description:Small-sized zokor (Table 2), with head and body length <175 mm (range: 145.0-175.0 mm) and skull length<42.20 mm (range: 38.01-42.17 mm).Only one elderly male(KIZ 0403327) significantly larger, with head and body length of 204 mm and skull length of 45.52 mm.Dorsal pelage dark grayish brown, with slightly cinnamon-colored tips, ventral pelage slightly paler; most individuals without white blaze on forehead, small when present.Lips and muzzle white,surrounded by short white hairs.Nose-pad trifoliate.Eyes very small and external ears absent (Figure 1B).Tail relatively long, densely hairy, proximal 2/3 grayish-brown, and distal 1/3 white.Forepaws strong and powerful; third claw longest and stoutest; second and fourth claws almost equal in length,about 2/3 of third claw; fifth claw stout, about half length of fourth finger; thumb very small.Dorsal surfaces of hands and feet covered with short white hair.

    Skull relatively weaker than other zokors (Figure 5).Rostrum relatively weak, thin, and nearly rectangular; nasals relatively small and trapezoidal, extending far beyond incisors,posterior border of nasals transverse.Premaxillaundeveloped, anterior half not expanded, nearly aligned with nasals.Zygomatic arches slightly expanded and thin;zygomatic plates narrow and simple.Interorbital region slightly constricted, almost same width as rostrum.Braincase welldomed, rounded ventrally.Temporal ridges not conspicuous(but much more developed in old male KIZ 0403327), parallel in front.Lambdoid crests sharp, but only on side of skull;occipital shield and occipital ridges developed, extending well posteriorly, almost forming plane with occipital condyle; small medial occipital ridge present between occipital ridges.Incisive foramina small, about 2/3 included within premaxilla,posterior 1/3 within maxillaries.Palate and pterygopalatine fossa broad, relatively wider than otherEospalaxspecies.Body of mandible curved.Coronoid process weak, pointed,and curved to posterior.Condyloid process stout and squared.Angular process small, curved upward almost joining with ramus.

    Figure 4 Plot of principal components 1 and 2 from analysis of 21 craniodental measurements of Eospalax genus

    Figure 5 Dorsal, ventral, and lateral views of skull and mandible of Eospalax muliensis sp.nov.(KIZ 040324, holotype) (A), E.baileyi (DF 004) (B), E.cansus (LM 001) (C), E.fontanierii (FS 001) (D), E.rothschildi (ZB 001) (E), E.rufescens (FP 001) (F), and E.smithii (XZ 001) (G)

    Teeth relatively light.Incisor teeth stout with tips curving slightly backward, outer side orange.Upper molar tooth rows slightly divergent to posterior; M1with two reentrant angles on outer and inner side; M2and M3with two reentrant angles on outer side, and one inside, no postero-external lobe on M3;lower incisor short, slightly curved, outer side orange; M1with two reentrant angles on outer side and three on inner side; M2and M3with two reentrant angles on outer and inner sides.Comparison: According to the measurements,Eospalax muliensissp.nov.is similar in size toE.rothschildi, but much smaller than other species in all dimensions (Table 2;Supplementary Table S3).The new species has some features unique in the genus.The skull ofEospalax muliensissp.nov.is obviously weaker than that of the otherEospalaxspecies, especially in the rostrum region (Figure 5).Therefore,the rostrum-related measurements (i.e., RSW, FIB, and NSL)ofEospalax muliensissp.nov.are the smallest amongEospalaxspecies (Table 2).The nasals of the new species are very short and narrow, the smallest in the genus.The premaxilla ofEospalax muliensissp.nov.is depressed and very thin, with the anterior half parts nearly aligned with the nasals, in contrast with other species, in which the premaxillae are well-developed and broader than the nasals.The palate and pterygopalatine fossa are relatively wider inEospalax muliensissp.nov than in other species, thus M2-M2is almost equal to LUM in the new species but is much smaller than LUM in other species.The braincase ofEospalax muliensissp.nov.are well-domed and temporal ridges being ambiguous, characteristics only seen in some individuals ofE.cansusand different from the flattened braincase with developed temporal ridges seen in other species.

    Eospalax muliensissp.nov.can be distinguished from small-sizedE.rothschildiby its strong and powerful forepaws(vs.slender); longer tail TL=45-70 mm (vs.short tail TL=23-38 mm, Supplementary Table S3); and proximal 2/3 of tail grayish-brown and distal 1/3 white (vs.gray above white below).Eospalax muliensissp.nov.can be distinguished fromE.fontanieriiandE.cansusbased on trifoliate nose-pad and densely hairy tail (vs.oval nose-pad and nearly naked tail).Eospalax muliensissp.nov.can be clearly distinguished fromE.baileyiandE.smithiibased on incisive foramina about 1/3 expanding premaxillary boundaries (vs.incisive foramina wholly contained in premaxillary).Eospalax muliensissp.nov.can be distinguished fromE.baileyiby the postero-external lobe on M3; the parallel temporal ridges ofE.muliensissp.nov.is distinct from the temporal ridges meeting in the median line ofE.smithii.The zygomatic arches ofE.rufescensare much stouter and flare wider thanEospalax muliensissp.nov.

    Distribution and habitat:Eospalax muliensissp.nov.is currently known only from the Kangwu Ranch, Muli County,Sichuan Province, China (elevational range ~3 700 m a.s.l.).The habitat is alpine meadow, surrounded by shrubs(Supplementary Figure S4).This species may be adapted to high-elevation habitats, and thus may occur in other alpine meadows in the southern Hengduan Mountains.

    DISCUSSION

    Taxonomic implications

    Both morphological and molecular analyses strongly suggested that the zokors inhabiting Muli County represent a new species (Eospalax muliensissp.nov.) ofEospalax.The new species can be diagnosed from all others withinEospalaxby the combination of small size, long tail, weak rostrum and nasal, and large palate and pterygopalatine fossa.Genomic analyses also revealed thatEospalax muliensissp.nov.diverged from other recognized zokor species ~4.22 Ma, with high genetic distances in nuclear and mitochondrial genomes.

    Our results also strongly supported the status of all six previously describedEospalaxspecies.The most recent divergences inEospalaxoccurred betweenE.baileyiandE.smithii, and betweenE.rufescensandE.rothschildi.Evidence revealed substantial differences among the two pairs of sister species and supported their species status.First, based on whole-genome analysis, the two species pairs diverged 2.09 Ma and 2.19 Ma (Figure 3E), with relatively high genetic distances of 1.05% and 1.02% (Figure 2B), respectively.Genetic distances of the mitochondrial cytbgene between sister species of these pairs reached 11.9% and 12.5%,respectively (Supplementary Table S4) (Bradley & Baker,2001).Second, the morphological characteristics that distinguish these sister species have been long recognized.For example,E.smithiican be distinguished fromE.baileyibased on temporal ridges meeting at the median line(Figure 5) (Fan & Shi, 1982); andE.rothschildican be distinguished fromE.rufescensbased on temporal ridges conspicuous and raised (Figure 5) (Fan & Shi, 1982).Third,the distribution of habitats differs among these sister species.Eospalax baileyiis dominant on the QTP (altitudinal range of 2 700-4 300 m a.s.l.), whileE.smithiiinhabits a small area around the Min Mountains (eastern edge of QTP) at an altitude of 2 700 m a.s.l.and possibly around the Liupanshan Mountains (within the Loess Plateau) at much lower altitudes(Fan & Shi, 1982).In the second sister species pair,E.rufescensoccurs in the Qinling Mountains, whileE.rothschildioccurs in the Daba Mountains and probably the Min Mountains.Although the forest habitats of these two species are similar, their local climate and flora differ, and they are likely geographically isolated by the Hanjiang River.In summary, genomic, morphological, and biogeographicalevidence well supports the delimitation of these two pairs of recently diverged sister species.Furthermore, the status of other species withinEospalaxthat were divergent longer before should also be retained.

    Phylogenetic implications

    Utilizing whole-genome sequencing data, our study revealed reliable and robust phylogenetic relationships among species in the genusEospalax.Before discussing the implications of the whole-genome-based phylogenetic relationships, we first explored the significant incongruences among phylogenies inferred from nuclear genomes, mitogenomes, and previous studies.Regarding phylogenies inferred from mitogenomes and mitochondrial genes, our results and previous studies show general similarity in the close relationship betweenE.baileyiandE.smithiiand the basal position ofE.fontanierii(Supplementary Figure S1) (Cai et al., 2020; He et al., 2012,2020; Su et al., 2014).However, extensive inconsistencies appear in mitochondrial analyses, which may be due to the limited informative sites of mitochondrial genes and complex speciation and post-speciation processes in the genusEospalax.For the latter, some clues may be present in our mitogenomic analysis that integrates the mitogenomes in our and previously published data.Our mitogenomic analyses revealed a robust phylogeny but embedded oneE.cansus(KC514112) individual tightly withinE.baileyi(Supplementary Figure S1).This individual was reportedly sampled from Lintan County, Gansu Province (around Min Mountains) (Su et al.,2013), which contains the type localities ofE.cansus,E.rothschildi, andE.smithii.Furthermore,E.baileyihabitats are also adjacent to this region.Besides the possible misidentification of the specimen or habitat expansion ofE.baileyiin Lintan, another plausible explanation may be postspeciation gene flow betweenE.cansusandE.baileyiaround the Min Mountains.If this hypothesis is true, the evolutionary history of zokors around the Min Mountains may be more complex, which requires further research.

    Several incongruences in the phylogenetic topologies inferred from nuclear and mitochondrial genomes were also found.Notably, nuclear genomic analysis placedE.fontanieriiwithin the low-altitude clade (Figure 2A), whereas mitogenomic analysis and most previous studies placedE.fontanieriiin a basal position, whileEospalax muliensissp.nov.were put aside (Supplementary Figure S1).In addition,nuclear genomic analysis placedEospalax muliensissp.nov.within the high-altitude clade (Figure 2A), while mitogenomic analysis placed it in the basal position of the genus(Supplementary Figure S1).Nuclear and mitochondrial discordance appears to be prevalent in phylogenetic studies and the possible roles of rapid speciation and incomplete lineage sorting have been discussed extensively (He et al.,2010; Mckay & Zink, 2010).Rapid speciation and incomplete lineage sorting may have contributed to the inconsistent relationships found amongE.cansus,E.rufescens, andE.rothschildi, given the rapid cladogenesis of these three species (Figures 2A, 3E).Other factors may also have contributed to the inconsistent positions ofE.fontanieriiandEospalax muliensissp.nov.For example, the clear pattern of mitochondrial phylogenomic analyses is the relevance between geographical distributions and phylogenetic positions of the species inEospalax.Eospalax fontanieriiandEospalax muliensissp.nov.are distributed in the remote northeast and southwest regions (Figure 1A), and their phylogenetic positions are basal.Furthermore, other species are geographically and phylogenetically adjacent.The possible explanation for this pattern may be ancient cross-species hybridization between geographically adjacent species or some other geographically relevant factors (Toews &Brelsford, 2012).However, these hypotheses lack solid evidence and require further comprehensive investigations.Regardless of the actual source of incongruence, the phylogeny inferred from nuclear genomes was more robust and thus serves as our current ‘‘best” estimate (Figure 2A).

    From the inferred nuclear whole-genome-based phylogeny,Eospalaxwas first divided into high- and low-altitude clades.Because uplift of the QTP may predate the separation of these two clades (Spicer et al., 2021) and no other species in the subfamily Myospalacinae or family Spalacidae inhabit the QTP, adaptation to high altitude may have occurred in the common ancestor ofEospalax muliensissp.nov.,E.baileyi,andE.smithii.Nowadays,E.smithiiinhabits areas at relatively low altitudes, implying possible re-adaptation to low elevations.Notably, given the possible existence of other unrecognized species in the Hengduan Mountains, adaptation to high altitude could have occurred in the common ancestor of allEospalaxspecies.Within the low-altitude clade,successful habitation of forest niches appeared inE.rufescensandE.rothschildi.Because these two species are closely related and separated more recently (Figures 2A, 3E),forest adaptation likely occurred in their common ancestor.Thus, although the dispersal capacity of zokors is extremely limited (Wei et al., 1997), species inEospalaxsuccessfully and rapidly adapted to diverse ecological niches in the Pliocene, and further investigations on their adaptive mechanisms are warranted.

    Implications of systematics and macroevolution from Muli zokors

    The discovery ofEospalax muliensissp.nov.in the southern Hengduan Mountains is unexpected and of great interest because the habitats ofEospalax muliensissp.nov.are southernmost of extant zokors and it is the second extant zokor species that inhabit extremely high altitudes (~3 700 m a.s.l.).The Hengduan Mountains are considered a global biodiversity hot spot (Myers et al., 2000).The extremely diverse ecological habitats in three dimensions of Hengduan Mountains (Sky islands) protected species from climate fluctuations and are known inhabited by some relict species(Chen et al., 2017; Chen et al., 2021; He & Jiang, 2014).The discovery ofEospalax muliensissp.nov.in such a biodiverse area suggests the possible existence of other unknown zokor species in the Hengduan Mountains.On the other hand, the discovery could also provide clues for deciphering the origin and evolutionary history of the genusEospalax.

    Based on the paleontological rules and evolutionary trajectories of zokor fossils, Li & Chen (1986) proposed nine plesiomorphic characters of ancientEospalaxspecies.Strikingly, of the nine characters, six are found inEospalaxmuliensissp.nov., in contrast to the five or fewer in other species inEospalax.For example, body size trajectory in zokor fossils is from small- to large-sized (Zheng, 1994),whereasEospalax muliensissp.nov.is among the smaller ones inEospalax.Eospalax muliensissp.nov.also possesses other plesiomorphic characters such as: long tail;weak rostrum and nasal; slightly expanded and thin zygomatic arch; ambiguous temporal ridge; and separated temporal ridge.Absent ofEospalax muliensissp.nov., previous studies proposed the Loess Plateau or Qinling-Daba Mountains as the possible origin center ofEospalax(Li & Chen, 1986; Li &Wang, 1996).Considering the plesiomorphic characters and relatively basal phylogenetic position ofEospalax muliensissp.nov., it is possible thatEospalaxmay have originated from the Hengduan Mountains.This hypothesis seems to contradict fossil evidence, as the earliest knownEospalaxfossils are recorded from the Nihewan Basin, Hebei Province (Zheng,1994) and are dated to ~2 Ma (MN18) (Li et al., 1984).However, the estimated diversification ofEospalaxspecies occurred ~4.68 Ma, and thus the contradiction may be caused by an incompleteEospalaxfossil record, especially in the Hengduan Mountains.Therefore, the origin and evolutionary history ofEospalaxremain obscure and require further interdisciplinary evidence.

    SCIENTIFIC FIELD SURVEY PERMISSION INFORMATION

    All collections followed the animal use protocols approved by the Animal Care and Ethics Committee of the Kunming Institute of Zoology, Chinese Academy of Sciences (Approval No.SMKX-SQ-2021-067).Permission for field sampling in each sampling sites was granted by the local Forestry Bureau.

    DATA AVAILABILITY

    All sequences reported in this study were deposited in the Genome Sequence Archive database (http://gsa.big.ac.cn/)under Accession ID (CRA005933) and in Science Data Bank(https://www.scidb.cn/) under DOI: 10.11922/sciencedb.o00023.00002.Mitogenomes were deposited in GenBank under Accession Nos.OM329951-OM330035.

    SUPPLEMENTARY DATA

    Supplementary data to this article can be found online.

    COMPETING INTERESTS

    The authors declare that they have no competing interests.

    AUTHORS’ CONTRIBUTIONS

    P.S.conceived and supervised the study.T.Z., M.L.L., H.Z.,and Z.Z.C.collected samples.Z.Z.C.and T.Z.performed morphometric comparison and morphological diagnosis.T.Z.analyzed whole-genome resequencing data.T.Z., Z.Z.C., and P.S.wrote the paper.All authors read and approved the final version of the manuscript.

    ACKNOWLEDGMENTS

    We are grateful to Bao-Guo Li (Northwest University) for sharing sample sites and helpful discussion.We are grateful to Jia-Tang Li (Chengdu Institute of Biology, Chinese Academy of Sciences), Zheng-Jun He (Animal Husbandry Bureau of Tibetan Autonomous County of Muli, Sichuan Province), and Shi-Yan Ding (Forestry Bureau of Zhenba County, Shaanxi Province) for help in field sampling.

    猜你喜歡
    鼢鼠木里
    青海木里煤田腐敗案背后的“保護(hù)傘”
    公民與法治(2022年7期)2022-07-22 07:11:54
    以沙為墨指尖畫夢
    小隴山林區(qū)中華鼢鼠的生活習(xí)性及防控措施
    正寧縣中華鼢鼠的為害及有效防治措施
    生態(tài)木里 綠色未來
    防火保安全我們在行動木里縣第一國有林場
    綠色天府(2021年5期)2021-06-24 16:40:58
    防火保安全 我們在行動
    綠色天府(2021年3期)2021-06-24 12:41:58
    中華鼢鼠的防治方法
    小隴山林區(qū)鼢鼠防治研究
    中華鼢鼠的生活習(xí)性及防治探析
    最近最新中文字幕大全电影3| 内地一区二区视频在线| 性色av乱码一区二区三区2| 欧美一区二区亚洲| 精品一区二区三区视频在线观看免费| 亚洲精品在线美女| 欧美zozozo另类| 91狼人影院| 亚洲人与动物交配视频| 亚洲精品粉嫩美女一区| 91久久精品电影网| 日韩欧美在线乱码| 国产综合懂色| 久久久久久久亚洲中文字幕 | 国产大屁股一区二区在线视频| 午夜福利成人在线免费观看| 精品日产1卡2卡| 日本撒尿小便嘘嘘汇集6| 最新在线观看一区二区三区| 国产精品嫩草影院av在线观看 | 精品人妻偷拍中文字幕| 色哟哟·www| 久久久久国内视频| 精品99又大又爽又粗少妇毛片 | 999久久久精品免费观看国产| 一a级毛片在线观看| 俄罗斯特黄特色一大片| 在线看三级毛片| 国内精品久久久久精免费| 欧美日韩瑟瑟在线播放| 一a级毛片在线观看| www.999成人在线观看| 一个人免费在线观看电影| 少妇人妻精品综合一区二区 | 亚洲真实伦在线观看| 亚洲av免费在线观看| 免费搜索国产男女视频| 最后的刺客免费高清国语| 99热这里只有是精品在线观看 | 观看美女的网站| 成人性生交大片免费视频hd| 老女人水多毛片| 亚洲国产精品合色在线| 亚洲国产精品999在线| 国产精品,欧美在线| 听说在线观看完整版免费高清| 久99久视频精品免费| 99热这里只有是精品50| 日本黄色视频三级网站网址| www.999成人在线观看| 男插女下体视频免费在线播放| 级片在线观看| 九九在线视频观看精品| 国产高清视频在线观看网站| 此物有八面人人有两片| 搞女人的毛片| 99精品在免费线老司机午夜| 波多野结衣巨乳人妻| 免费观看精品视频网站| 欧美一级a爱片免费观看看| 成人欧美大片| 午夜久久久久精精品| 又黄又爽又刺激的免费视频.| 国产精品久久电影中文字幕| 成年免费大片在线观看| 国产成人aa在线观看| 久久午夜亚洲精品久久| 人妻久久中文字幕网| 久久精品国产亚洲av涩爱 | 色精品久久人妻99蜜桃| 国产高清三级在线| 身体一侧抽搐| 欧美在线黄色| 一本久久中文字幕| 高潮久久久久久久久久久不卡| 永久网站在线| 天堂动漫精品| 最近在线观看免费完整版| 国产色爽女视频免费观看| 夜夜看夜夜爽夜夜摸| avwww免费| 超碰av人人做人人爽久久| 国产淫片久久久久久久久 | 好男人电影高清在线观看| 男女床上黄色一级片免费看| bbb黄色大片| 精品久久久久久久人妻蜜臀av| 大型黄色视频在线免费观看| 久久久久久久久大av| 亚洲精华国产精华精| 如何舔出高潮| 99视频精品全部免费 在线| 他把我摸到了高潮在线观看| 欧美色视频一区免费| 亚洲av成人精品一区久久| 亚洲第一区二区三区不卡| 久久久久久久久久黄片| 真人一进一出gif抽搐免费| 欧美三级亚洲精品| 搡老岳熟女国产| 久久久久久久久久黄片| 午夜福利在线观看免费完整高清在 | 亚洲精品一区av在线观看| 色综合亚洲欧美另类图片| 欧美日韩乱码在线| 国产欧美日韩一区二区精品| 91麻豆av在线| 国产高清三级在线| 真实男女啪啪啪动态图| 黄色丝袜av网址大全| 国产高清三级在线| 赤兔流量卡办理| 国产亚洲精品综合一区在线观看| 可以在线观看毛片的网站| 日日摸夜夜添夜夜添av毛片 | 国产黄片美女视频| 亚洲成人免费电影在线观看| 免费av不卡在线播放| 国内久久婷婷六月综合欲色啪| 精品午夜福利视频在线观看一区| 国产人妻一区二区三区在| 在线免费观看不下载黄p国产 | 极品教师在线视频| 中文字幕精品亚洲无线码一区| 欧美成人一区二区免费高清观看| 偷拍熟女少妇极品色| 国产私拍福利视频在线观看| av国产免费在线观看| av国产免费在线观看| 中文字幕人成人乱码亚洲影| 成人三级黄色视频| 国产成年人精品一区二区| 琪琪午夜伦伦电影理论片6080| 啪啪无遮挡十八禁网站| www.www免费av| 亚洲欧美激情综合另类| 校园春色视频在线观看| 亚洲天堂国产精品一区在线| 老女人水多毛片| 老女人水多毛片| 国产av一区在线观看免费| 可以在线观看毛片的网站| 搞女人的毛片| 一区二区三区激情视频| 欧美bdsm另类| or卡值多少钱| 美女cb高潮喷水在线观看| 夜夜爽天天搞| 一本久久中文字幕| 非洲黑人性xxxx精品又粗又长| 国产单亲对白刺激| 亚洲美女搞黄在线观看 | 中文字幕av成人在线电影| 亚洲国产高清在线一区二区三| 午夜两性在线视频| 亚洲精品在线美女| 成人国产一区最新在线观看| 色综合亚洲欧美另类图片| 成年女人永久免费观看视频| 久久天躁狠狠躁夜夜2o2o| 久久久久精品国产欧美久久久| 日本黄色视频三级网站网址| 亚洲av一区综合| 欧美成人免费av一区二区三区| 国产伦一二天堂av在线观看| 久久久久久久久久黄片| 国产美女午夜福利| 老熟妇乱子伦视频在线观看| 琪琪午夜伦伦电影理论片6080| 能在线免费观看的黄片| 夜夜躁狠狠躁天天躁| 国产乱人伦免费视频| 日韩中文字幕欧美一区二区| 身体一侧抽搐| 亚洲专区中文字幕在线| 精品国产三级普通话版| 丰满乱子伦码专区| 午夜福利视频1000在线观看| 丰满人妻熟妇乱又伦精品不卡| 亚洲成人久久性| 一个人免费在线观看电影| 老女人水多毛片| 久久欧美精品欧美久久欧美| 日本一二三区视频观看| 欧美精品啪啪一区二区三区| 五月玫瑰六月丁香| 国产视频内射| 国产免费一级a男人的天堂| 在线观看美女被高潮喷水网站 | 国产三级在线视频| 精品欧美国产一区二区三| 国产精品人妻久久久久久| 国产不卡一卡二| 久久午夜亚洲精品久久| 欧美激情国产日韩精品一区| 亚洲精品粉嫩美女一区| av黄色大香蕉| 午夜免费成人在线视频| 欧美成狂野欧美在线观看| 国产亚洲精品av在线| 久久久成人免费电影| 国产精品美女特级片免费视频播放器| 午夜精品在线福利| 日韩有码中文字幕| 午夜日韩欧美国产| 动漫黄色视频在线观看| 日日摸夜夜添夜夜添小说| 国产亚洲欧美98| 一本综合久久免费| 丁香欧美五月| 99久久无色码亚洲精品果冻| 中国美女看黄片| 色综合欧美亚洲国产小说| 久久久久久久久中文| 在线观看66精品国产| 国产真实伦视频高清在线观看 | АⅤ资源中文在线天堂| 一本一本综合久久| 午夜精品一区二区三区免费看| 三级国产精品欧美在线观看| 久久国产乱子伦精品免费另类| 国产伦精品一区二区三区视频9| а√天堂www在线а√下载| 国内精品久久久久久久电影| 久久伊人香网站| 亚洲天堂国产精品一区在线| 少妇熟女aⅴ在线视频| 亚洲一区二区三区色噜噜| 色综合婷婷激情| av天堂在线播放| 99国产综合亚洲精品| 蜜桃亚洲精品一区二区三区| 看片在线看免费视频| 国产一级毛片七仙女欲春2| 亚洲,欧美精品.| 欧美乱妇无乱码| 天堂影院成人在线观看| 亚洲经典国产精华液单 | 国产午夜精品论理片| 国产精品一及| 好男人在线观看高清免费视频| 一本精品99久久精品77| 在线天堂最新版资源| 国产主播在线观看一区二区| 亚洲欧美清纯卡通| av专区在线播放| 一进一出抽搐gif免费好疼| 99在线人妻在线中文字幕| 亚洲人与动物交配视频| 特级一级黄色大片| 午夜a级毛片| 一本久久中文字幕| 久久人妻av系列| 亚洲第一区二区三区不卡| 国产av不卡久久| 热99re8久久精品国产| 在线播放无遮挡| 少妇被粗大猛烈的视频| 好男人在线观看高清免费视频| 伦理电影大哥的女人| 色综合站精品国产| 成人美女网站在线观看视频| 在线a可以看的网站| 欧美三级亚洲精品| 两个人视频免费观看高清| 看片在线看免费视频| 高清日韩中文字幕在线| 免费黄网站久久成人精品 | 亚洲欧美精品综合久久99| 久久久久久九九精品二区国产| 久久天躁狠狠躁夜夜2o2o| 久久久久免费精品人妻一区二区| 久久人人爽人人爽人人片va | 日韩欧美 国产精品| 日本 欧美在线| 亚洲成av人片在线播放无| 亚洲国产欧洲综合997久久,| 国产 一区 欧美 日韩| 日韩精品青青久久久久久| 男插女下体视频免费在线播放| 亚洲欧美日韩东京热| 久久久久久国产a免费观看| 级片在线观看| 久久99热这里只有精品18| 国产69精品久久久久777片| 一区二区三区高清视频在线| 在线观看午夜福利视频| 搞女人的毛片| 久久久精品欧美日韩精品| 午夜精品一区二区三区免费看| 国产精品久久久久久亚洲av鲁大| 两个人视频免费观看高清| 午夜视频国产福利| 最好的美女福利视频网| 女人十人毛片免费观看3o分钟| 欧美高清成人免费视频www| 国产精品一区二区三区四区久久| 久久国产精品人妻蜜桃| 18禁黄网站禁片免费观看直播| 一级黄片播放器| 精品一区二区三区人妻视频| 久久性视频一级片| 三级男女做爰猛烈吃奶摸视频| 少妇裸体淫交视频免费看高清| a级一级毛片免费在线观看| 亚洲av美国av| 国产精品一区二区免费欧美| 亚洲美女黄片视频| x7x7x7水蜜桃| 精品一区二区三区视频在线观看免费| 久久久久九九精品影院| 99视频精品全部免费 在线| 69av精品久久久久久| 亚洲一区二区三区不卡视频| 久久人人爽人人爽人人片va | 亚洲精品成人久久久久久| 国产中年淑女户外野战色| 麻豆久久精品国产亚洲av| 精品久久久久久久末码| 亚洲va日本ⅴa欧美va伊人久久| a级毛片a级免费在线| 中文字幕久久专区| 美女高潮喷水抽搐中文字幕| 乱码一卡2卡4卡精品| 最新中文字幕久久久久| 成人亚洲精品av一区二区| 午夜福利18| 一进一出抽搐动态| 此物有八面人人有两片| 很黄的视频免费| 亚洲专区中文字幕在线| 亚洲人成电影免费在线| 亚洲精品在线美女| 亚洲成av人片免费观看| 国产精品一区二区性色av| 国产亚洲精品久久久com| 在线观看免费视频日本深夜| 日韩高清综合在线| 全区人妻精品视频| 国产亚洲精品久久久com| 三级男女做爰猛烈吃奶摸视频| 91久久精品国产一区二区成人| 一个人看视频在线观看www免费| 听说在线观看完整版免费高清| 免费一级毛片在线播放高清视频| 长腿黑丝高跟| 免费在线观看日本一区| 长腿黑丝高跟| 亚洲在线观看片| 免费电影在线观看免费观看| 国产aⅴ精品一区二区三区波| 久久这里只有精品中国| 久久久久精品国产欧美久久久| 热99在线观看视频| 免费黄网站久久成人精品 | 精品久久久久久久久久久久久| 久久亚洲精品不卡| 欧美极品一区二区三区四区| 亚洲美女视频黄频| 十八禁国产超污无遮挡网站| 日本成人三级电影网站| 精品久久国产蜜桃| 亚洲久久久久久中文字幕| 亚洲专区国产一区二区| 91狼人影院| 亚洲av美国av| 国产成人av教育| 亚洲人成网站在线播放欧美日韩| 精品久久久久久久久亚洲 | 俺也久久电影网| 久久久久免费精品人妻一区二区| 国产欧美日韩一区二区三| 一级av片app| 久久香蕉精品热| 搡老妇女老女人老熟妇| 99精品在免费线老司机午夜| 热99在线观看视频| 欧美丝袜亚洲另类 | 欧美性猛交╳xxx乱大交人| 成人一区二区视频在线观看| 香蕉av资源在线| 一卡2卡三卡四卡精品乱码亚洲| 国产精品久久久久久亚洲av鲁大| av女优亚洲男人天堂| 在线观看免费视频日本深夜| 欧美zozozo另类| 国产伦在线观看视频一区| 在线免费观看的www视频| 久久草成人影院| ponron亚洲| 国产亚洲精品久久久com| h日本视频在线播放| 精品免费久久久久久久清纯| 国产av一区在线观看免费| eeuss影院久久| 国产精品伦人一区二区| 国产野战对白在线观看| 欧美色视频一区免费| 在线天堂最新版资源| 久久天躁狠狠躁夜夜2o2o| 日本熟妇午夜| 制服丝袜大香蕉在线| 亚洲午夜理论影院| 久久人人爽人人爽人人片va | 三级男女做爰猛烈吃奶摸视频| 黄色视频,在线免费观看| 中文字幕久久专区| 日本五十路高清| 毛片一级片免费看久久久久 | 日韩 亚洲 欧美在线| 观看免费一级毛片| 无遮挡黄片免费观看| 午夜老司机福利剧场| 国产成+人综合+亚洲专区| 国产精品日韩av在线免费观看| 夜夜躁狠狠躁天天躁| 亚洲av免费在线观看| 一本精品99久久精品77| 亚洲欧美日韩无卡精品| 色视频www国产| 97碰自拍视频| 99在线视频只有这里精品首页| 全区人妻精品视频| 精品久久久久久久末码| 国产白丝娇喘喷水9色精品| 69av精品久久久久久| 久久久久国产精品人妻aⅴ院| 一个人免费在线观看电影| 91在线观看av| 久久久久久久午夜电影| 又爽又黄a免费视频| www日本黄色视频网| 国产精品亚洲av一区麻豆| 国产三级中文精品| 亚洲欧美日韩高清专用| 全区人妻精品视频| 欧美黑人巨大hd| 国产毛片a区久久久久| 最近中文字幕高清免费大全6 | 久久久久精品国产欧美久久久| 美女免费视频网站| 最新在线观看一区二区三区| 亚洲精品亚洲一区二区| 亚洲av免费高清在线观看| 欧美最新免费一区二区三区 | 琪琪午夜伦伦电影理论片6080| 国产成人aa在线观看| 国语自产精品视频在线第100页| 舔av片在线| 成人美女网站在线观看视频| 无遮挡黄片免费观看| 99在线视频只有这里精品首页| 亚洲av二区三区四区| 亚洲黑人精品在线| 中文资源天堂在线| 色吧在线观看| 直男gayav资源| 搡女人真爽免费视频火全软件 | 99国产精品一区二区三区| 亚洲自偷自拍三级| 欧美另类亚洲清纯唯美| 动漫黄色视频在线观看| 精品免费久久久久久久清纯| 看片在线看免费视频| 久久久久久久久久黄片| 日本三级黄在线观看| 最新中文字幕久久久久| 亚洲一区高清亚洲精品| 午夜免费成人在线视频| 在现免费观看毛片| 三级男女做爰猛烈吃奶摸视频| 欧美性猛交╳xxx乱大交人| 国产在线精品亚洲第一网站| 精品久久久久久久久亚洲 | 能在线免费观看的黄片| 国产精品1区2区在线观看.| 91久久精品国产一区二区成人| 欧美激情久久久久久爽电影| 搡老妇女老女人老熟妇| 亚洲无线观看免费| 99国产综合亚洲精品| 免费av观看视频| 亚洲无线在线观看| 亚洲人成网站高清观看| 美女黄网站色视频| 精品久久久久久久久久久久久| 国产精品免费一区二区三区在线| av在线老鸭窝| 男女床上黄色一级片免费看| 99久久成人亚洲精品观看| 亚洲美女视频黄频| 久久草成人影院| 国产真实伦视频高清在线观看 | 久9热在线精品视频| 国产亚洲精品久久久久久毛片| 91在线观看av| 乱人视频在线观看| 韩国av一区二区三区四区| 十八禁网站免费在线| 啦啦啦观看免费观看视频高清| 国产精品,欧美在线| 天堂影院成人在线观看| 精品乱码久久久久久99久播| 十八禁国产超污无遮挡网站| 丰满人妻一区二区三区视频av| 亚洲av不卡在线观看| 在线观看一区二区三区| 男人狂女人下面高潮的视频| 深夜a级毛片| 精品午夜福利视频在线观看一区| 中文字幕人成人乱码亚洲影| 亚洲人成网站在线播| av国产免费在线观看| 最新在线观看一区二区三区| 国产三级中文精品| 一进一出抽搐动态| 日本免费一区二区三区高清不卡| 亚洲欧美日韩高清专用| 成人特级黄色片久久久久久久| 桃红色精品国产亚洲av| 不卡一级毛片| 简卡轻食公司| 日韩精品青青久久久久久| 精品无人区乱码1区二区| 亚洲精品在线观看二区| 美女高潮的动态| 亚洲国产精品sss在线观看| 国产精品久久久久久久久免 | or卡值多少钱| 91九色精品人成在线观看| 欧美日韩乱码在线| 免费av观看视频| 九色成人免费人妻av| 最近最新中文字幕大全电影3| 欧美黄色片欧美黄色片| 亚洲人成网站高清观看| 日韩欧美免费精品| 成年女人看的毛片在线观看| 国产成人欧美在线观看| 久久久久久久久久黄片| 免费看光身美女| 男人和女人高潮做爰伦理| 国产午夜精品久久久久久一区二区三区 | 国产高清激情床上av| a级一级毛片免费在线观看| 久久亚洲真实| 欧美又色又爽又黄视频| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 中文字幕人妻熟人妻熟丝袜美| 噜噜噜噜噜久久久久久91| 国产视频一区二区在线看| 日本黄色片子视频| 精品午夜福利在线看| 少妇熟女aⅴ在线视频| 别揉我奶头~嗯~啊~动态视频| 久久精品国产自在天天线| 免费一级毛片在线播放高清视频| 亚洲精华国产精华精| 日韩欧美国产一区二区入口| 又爽又黄a免费视频| 99国产综合亚洲精品| 一区二区三区激情视频| 亚洲第一欧美日韩一区二区三区| АⅤ资源中文在线天堂| 少妇高潮的动态图| 午夜激情欧美在线| 亚洲av.av天堂| 一本综合久久免费| 精品欧美国产一区二区三| 精华霜和精华液先用哪个| av专区在线播放| 亚洲五月天丁香| 精品人妻1区二区| 十八禁国产超污无遮挡网站| 男女之事视频高清在线观看| 国产高清视频在线播放一区| 亚洲美女搞黄在线观看 | avwww免费| 日韩欧美在线乱码| 桃红色精品国产亚洲av| 久久国产乱子免费精品| 窝窝影院91人妻| 中亚洲国语对白在线视频| 日本与韩国留学比较| 夜夜躁狠狠躁天天躁| 欧美zozozo另类| 嫩草影视91久久| 中文字幕人妻熟人妻熟丝袜美| 极品教师在线免费播放| 亚洲中文字幕一区二区三区有码在线看| 亚洲欧美精品综合久久99| 久久久久久国产a免费观看| 在线观看66精品国产| 麻豆av噜噜一区二区三区| 亚洲熟妇中文字幕五十中出| 国产精品乱码一区二三区的特点| 18禁黄网站禁片午夜丰满| 欧美午夜高清在线| 老司机深夜福利视频在线观看| av女优亚洲男人天堂| 久久99热这里只有精品18| 欧美激情国产日韩精品一区| 一级作爱视频免费观看| 日韩欧美一区二区三区在线观看| 成人永久免费在线观看视频| 精品一区二区三区视频在线| 他把我摸到了高潮在线观看| 一个人看视频在线观看www免费| 精品一区二区三区视频在线| eeuss影院久久| 日韩欧美在线乱码| 禁无遮挡网站| 搡女人真爽免费视频火全软件 | xxxwww97欧美| 国产欧美日韩精品亚洲av|