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

    Comprehensive analysis of coding and non-coding RNA transcriptomes related to hypoxic adaptation in Tibetan chickens

    2021-12-17 11:52:56YingZhangWoyuSuBoZhangYaoLingWooKyunKimandHaoZhang

    Ying Zhang,Woyu Su,Bo Zhang,Yao Ling,Woo Kyun Kimand Hao Zhang*

    Abstract

    Background:Tibetan chickens,a unique native breed in the Qinghai-Tibet Plateau of China,possess a suite of adaptive features that enable them to tolerate the high-altitude hypoxic environment.Increasing evidence suggests that long non-coding RNAs(lncRNAs)and microRNAs(miRNAs)play roles in the hypoxic adaptation of high-altitude animals,although their exact involvement remains unclear.

    Results:This study aimed to elucidate the global landscape of mRNAs,lncRNAs,and miRNAs using transcriptome sequencing to construct a regulatory network of competing endogenous RNAs(ceRNAs)and thus provide insights into the hypoxic adaptation of Tibetan chicken embryos.In total,354 differentially expressed genes(DE genes),389 differentially expressed lncRNAs(DE lncRNAs),and 73 differentially expressed miRNAs(DE miRNAs)were identified between Tibetan chickens(TC)and control Chahua chickens(CH).GO and KEGG enrichment analysis revealed that several important DE miRNAs and their target DE lncRNAs and DE genes are involved in angiogenesis(including blood vessel development and blood circulation)and energy metabolism(including glucose,carbohydrate,and lipid metabolism).The ceRNA network was then constructed with the predicted DE gene-DE miRNA-DE lncRNA interactions,which further revealed the regulatory roles of these differentially expressed RNAs during hypoxic adaptation of Tibetan chickens.

    Conclusions:Analysis of transcriptomic data revealed several key candidate ceRNAs that may play high-priority roles in the hypoxic adaptation of Tibetan chickens by regulating angiogenesis and energy metabolism.These results provide insights into the molecular mechanisms of hypoxic adaptation regulatory networks from the perspective of coding and non-coding RNAs.

    Keywords:ceRNA,Hypoxic adaptation,Non-coding RNA,Tibetan chicken,Transcriptome

    Background

    MicroRNAs(miRNAs)are small non-coding RNA molecules that inhibit gene expression by binding to specific mRNAs.LncRNAs are another type of noncoding RNA,with a length of more than 200 nucleotides(nt),that regulate gene expression through a number of mechanisms,including epigenetic regulation(genetic imprinting and chromatin remodeling),transcriptional regulation(transcription interference),post-transcriptional regulation(splicing),and so on[1].Both of these non-coding RNAs play different roles in various aspects of cellular function,including the regulation of hypoxia-related genes.Through their activity on target transcripts,some critical miRNAs regulate angiogenesis,the cell cycle,as well as lactate transport under hypoxic conditions[2–4].HIF-1 and VEGF are regulated by miR-17-92,miR-20b,and miR-199a,which play roles in tumor development,cancer cell proliferation,and cardiac myocytes[5–7].Various studies have shown that lncRNAs can directly bind to HIF-1 and activate its expression,in turn modulating hypoxic responses similar to the suppression of cancer cell adaptation to hypoxia,which provides important insights into how tumor cells sense and adapt to hypoxic stress[8–10].Although there have been various studies on genes involved in high-altitude adaptation in humans and animals,the regulatory mechanism of non-coding RNAs involved in hypoxic adaptation remains largely unknown.

    Tibetan chickens,which live in Tibetan highlands and have undergone long-term natural selection,exhibit stable genetic adaptation to high-altitude environments,characterized by high hatchability under hypoxic incubation when compared to lowland chickens[11,12].This indicates that Tibetan chickens have special physiological and genetic mechanisms enabling them to adapt to extreme high-elevation environments,and ensuring their normal hatching,reproduction,and growth in plateau environments[13–15].Previous studies using genomic,transcriptomic,and proteomic profiles have identified several functional genes related to hypoxic adaptation in Tibetan chickens[16–18].However,the mechanisms involved in the regulation of hypoxic adaptation by non-coding RNA are unclear.

    Chorioallantoic membrane(CAM),which functions as both a respiratory and circulatory organ for the chicken embryo,contains numerous blood vessels.Tibetan chicken embryos have blunted responses to hypoxia on CAM angiogenesis,which might benefit blood flow and oxygen transportation during hypoxic incubation[18].The characteristics of hypoxic adaptation may involve regulation of expressed mRNAs and non-coding RNAs;however,this remains unclear.In this study,full transcriptome sequencing(RNA-seq)of chicken embryos was performed to explore the regulation of non-coding RNAs and candidate genes responsible for hypoxic adaptation traits in Tibetan chickens as compared to the Chahua chickens.Integration of three types of RNAs(mRNAs,lncRNAs and miRNAs)revealed competing endogenous RNA(ceRNA)interactive network and potential regulatory mechanism to explain hypoxic adaptation in Tibetan chickens.These findings will help us further understand the molecular mechanism of adaptability to hypoxic conditions in Tibetan chickens from a non-coding RNA perspective.

    Methods

    Sample collection and preparation

    The experiments and animal care protocol were approved by the animal welfare committee of the State Key Laboratory for Agro-Biotechnology of the China Agricultural University(approval number,XK257),and all methods were performed in accordance with the relevant guidelines and regulations.

    We collected eggs from Tibetan chickens(TCs)and Chahua chickens(CHs)from the Experimental Station for Poultry Genetic Resources and Breeding,China Agricultural University(CAU).The two groups of eggs were incubated in a hypoxic incubator(13%±0.2% O2)and the temperature and humidity of the incubator were 37.8°C and 60%,respectively.Eighteen CAM samples from Tibetan chickens(n=9)and Chahua chickens(n=9)were collected from female embryos at day 11 of incubation under hypoxic condition,immediately frozen in liquid nitrogen,and stored at ?80°C for RNA extraction and sequencing.

    RNA extraction,library construction,and sequencing

    Total RNA was extracted from the CAM samples using an RNA pure Tissue Kit(Tiangen Biotech Co.Ltd.,Beijing,China)according to the manufacturer’s protocol.Six RNA-seq libraries(TC1,TC2,TC3,CH1,CH2,and CH3)were constructed,three biological duplicates were prepared,and each sample consisted of a mix from three individuals of the same breed.

    For mRNA and lncRNA sequencing,1.5μg of total RNA was used as input material for rRNA removal using the Ribo-Zero rRNA Removal Kit(Epicentre,Madison,WI,USA).Sequencing libraries were generated using the NEBNextRUltra? Directional RNA Library Prep Kit for IlluminaR(NEB,USA)following the manufacturer’s instructions,and 150 paired-end reads were generated.For small RNA sequencing,1.5μg of total RNA per sample was used to ligate the 3′SR and 5′SR Adaptor.Reverse transcription synthetic first chain and PCR amplification were carried out,PAGE gel was used to screen the electrophoresis fragment,then rubber cutting and recycling.The PCR products were purified(AMPure XP system)and library quality was assessed.Clustering of the indexcoded samples was performed using TruSeq PE Cluster Kit v4-cBot-HS(Illumina)according to the manufacturer’s instructions,and 50-nt single-end reads were generated.

    Read mapping and transcriptome assembly

    For mRNA and lncRNA data,using in-house Perl scripts,raw 150bp paired-end reads were filtered by removing adapters and reads containing poly-N and low-quality sequences to obtain clean reads.Clean reads were aligned to the chicken reference genome Gallus-gallus-6.0(GRCg6a-96.Gallus-gallus.GRCg6a-96.genome.fa)using the HISAT2 software[19].The transcriptome was assembled using StringTie and annotated using the gffcompare program.For miRNA data,the 3′linker sequences were removed from the raw 50-nt single-end reads,which were then filtered by removing reads with≥10% base N and shorter than 18 or longer than 30 nucleotides to obtain clean reads.The clean reads were aligned with the Silva,GtRNAdb,Rfam,and Repbase databases to filter other non-coding RNAs,such as ribosomal RNA(rRNA),transfer RNA(tRNA),small nuclear RNA(snRNA),small nucleolar RNA(snoRNA),and repeating sequences to obtain unannotated reads containing miRNAs.Unannotated reads were mapped to the reference genome(Gallus-gallus.GRCg6a-96)using the Bowtie software[20].The Q20,Q30,and GC content of the clean data were calculated for all RNAs.

    Identification of lncRNAs and prediction of their target genes

    The unknown transcripts were screened as putative lncRNAs or protein-coding RNAs.The screening methods were used according to previous reports[21,22],which briefly included length≥200bp,exon number≥2,and FPKM≥0.1,using tools of CPC2(coding potential calculator),CNCI(coding-non-coding index),CPAT(coding potential assessment tool),and pfam protein domain analysis tools[23–26].The transcripts were aligned with lncRNA databases(lncRNAdb and NONCODE)to identify known lncRNAs.Integration of following two methods to predict target genes of the lncRNAs:(1)lncRNA regulation of the expression of its neighboring genes,which is based on the distance(within 100kb)between lncRNA and potential target genes in the genome,and(2)correlations between lncRNA and potential target gene expression in all samples that were significantly positive.

    Identification of miRNAs and prediction of their target genes

    Known miRNAs were identified by mapping the matched mature miRNAs in miRBase(v22)to the reference genome.Novel miRNAs were predicted using miRDeep2 from the unannotated reads[27].Target genes of miRNAs were predicted by miRanda and TargetScan based on the miRNA and gene sequence information[28,29].

    RNA quantification and differential expression analysis

    Expression level of mRNAs and lncRNAs were quantified based on fragments per kilobase of transcript per million fragments mapped(FPKM)using the StringTie software[30].Expression levels of miRNAs were quantified and normalized using the TPM algorithm(TPM=Readcount×1,000,000/mapped reads)[31].The DESeq2 software package was used to detect differentially expressed lncRNAs(DE lncRNAs),differentially expressed genes(DE genes),and differentially expressed miRNAs(DE miRNAs)with criteria of log2|fold change|>1 and P<0.05[32].

    GO and KEGG enrichment analysis of differentially expressed RNA

    Gene Ontology(GO)enrichment analysis of the DE genes was performed using the GOseq R package based on Wallenius’non-central hypergeometric distribution.Kyoto Encyclopedia of Genes and Genomes(KEGG)is a database resource for understanding high-level gene functions and utilities of a biological system,such as the cell,organism,and ecosystem,using molecular-level information,especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies(http://www.genome.jp/kegg/).We used the KOBAS software to test the statistical significance of DE gene enrichment in KEGG pathways.

    Construction of the ceRNA network

    Competitive endogenous RNAs(ceRNAs)are novel transcriptional regulators that can modulate the expression of transcripts through competitive binding.The interactions of miRNA-mRNA and miRNA-lncRNA was predicted using TargetScan and miRanda,and the competitive pairs of mRNA-lncRNA were identified using hypergeometric test and Pearson correlation coefficient.Based on the sharing of miRNA binding sites and competitive pairs between lncRNAs and mRNAs,ceRNA network was constructed as follows:P<0.01 and FDR<0.01 in hypergeometric test[33–35].Visualization of ceRNA using Cytoscape software(http://cytoscape.org/).

    Expression analysis of DE lncRNAs,DE miRNAs,and DE genes using qRT-PCR

    To confirm the differentially expressed RNAs identified through RNA-seq,quantitative real-time PCR(qRT-PCR)was carried out for six DE lncRNAs,six DE miRNAs,and seven DE genes.We used the PrimeScript RT reagent kit with gDNA(Tiangen Biotech Co.Ltd.,Beijing,China)to convert total RNA to cDNA,using random hexamers for mRNA and lncRNAs,and stem-loop RT primers for miRNAs.The qPCR was performed using the SYBR Green PCR Kit(Tiangen Biotech Co.Ltd.,Beijing,China)according to the manufacturer’s instructions.The internal controls included GAPDH(glyceraldehyde-3-phosphate dehydrogenase)for lncRNAs,chicken 5S rRNA for miRNAs,and HPRT(hypoxanthine guanine phosphoribosyl transferase)for mRNAs.All primers used in qRT-PCR are shown in Additional file 1:Table S1.The BioRad CFX96(Bio-Rad,CA,USA)was used to perform qRT-PCR with the SYBR Green PCR Master Mix(Tiangen Biotech Co.Ltd.).Each qRT-PCR experiment was performed in triplicate,relative RNA expression levels were calculated using the 2?ΔΔCTmethod,and ANOVA testing via SPSS 25.0 was applied.The results are expressed as mean±SE,and P≤0.05 was considered statistically significant.Software and associated parameters used are shown in Additional file 2:Table S2.

    Results

    Overview of RNA transcriptomic profiles of chicken CAM tissues

    After low-quality reads were removed,110.3–114.1 M paired-end reads(2×150 bp in length)for the mRNAs and lncRNAs,and 23.8–29.2M single-end reads(50 bp in length)for the miRNAs,were obtained from the CAM of Tibetan and Chahua chicken samples(Additional file 3:Table S3).The majority of the lncRNAs(60.7%)were intergenic long-stranded non-coding RNA (lincRNA),followed by intronic-lncRNA(17.7%),antisense-lncRNA(12.4%),and 9.2% sense-lncRNA(Fig.1a).The length of mRNAs ranged from 51 to 30,0721bp,while that of lncRNAs ranged from 202 to 223,626bp;the majority of mRNAs and lncRNAs were 201–400bp and>3000bp long,respectively(Additional file 4:Fig.S1).In total,20,912 expressed genes,15,523 lncRNAs(3194 known),and 1439 miRNAs were identified in the CAM tissues of Tibetan and Chahua chickens.The 1439 miRNAs included 471 known and 968 novel miRNAs,and most of the miRNA reads were between 20 and 24nt in length(Fig.1b,Additional file 5:Table S4).The values of FPKM of mRNA/lncRNA and TPM of miRNA distribution were similar among the six samples(Fig.1c,d,e).

    Fig.1 Comparison of the features of RNAs in the chorioallantoic membrane of chicken embryos.(a)Statistics of long non-coding RNA,including intergenic non-coding RNAs(lincRNAs),anti-sense lncRNAs,intronic lncRNAs,and sense lncRNAs;(b)Length distribution of 1439 miRNAs;(c)Distribution of mRNA expression values(FPKM)of six samples;(d)Distribution of lncRNA expression values(FPKM)of six samples;(e)Distribution of miRNA expression values(TPM)of six samples.The ordinate shows the log10(FPKM+1)(c,d)and the log10(TPM+1)(e)of each RNA from the six sets of RNA-seq data;the middle line in the box represents the median of FPKM/TPM.TC,Tibetan chicken;CH,Chahua chicken

    Differential expression analysis of RNAs between TC and CH

    We found 389 DE lncRNAs,73 DE miRNAs,and 354 DE genes between the TC and CH groups,among which 194 DE lncRNAs,46 DE miRNAs,and 180 DE genes were upregulated,whereas 195 DE lncRNAs,27 DE miRNAs,and 174 DE genes were downregulated in the Tibetan chickens compared to Chahua chickens(Fig.2,Additional file 6:Table S5).The DE lncRNAs,DE miRNAs,and DE genes were classified into two categories for the six samples,in which the expressed quantity of the DE lncRNAs,DE miRNAs,and DE genes had good repeatability within groups respectively(Additional file 4:Fig.S2).

    Fig.2 Identification of differentially expressed(DE)lncRNAs,DE miRNAs,and DE genes between Tibetan(TC)and Chahua chicken(CH).Number of DE lncRNAs,DE miRNAs,and DE genes(a);Volcano plots displaying DE lncRNAs(b),DE miRNAs(c),and DE genes(d)between Tibetan and Chahua chicken.Upregulated and downregulated genes are shown in red and green,respectively.Black dots represent genes with similar expression levels

    Functional analysis of differentially expressed non-coding RNAs(DE lncRNAs and DE miRNAs)

    For 389 DE lncRNAs,we predicted 722 cis-and transtarget genes that were mainly involved in the GO terms of angiogenesis,blood vessel development,regulation of vasoconstriction,regulation of blood pressure,response to oxygen levels,ATP metabolic process,fatty acid metabolic process,glucuronate metabolic process,enzyme binding,and ATP binding(Fig.3a,Additional file 7:Table S6a).The KEGG analysis revealed the involvement of target genes in carbon metabolism,fatty acid metabolism,vascular smooth muscle contraction,and the Wnt,mTOR,MAPK,VEGF,and calcium signaling pathways(Fig.3b,Additional file 7:Table S6b).

    Fig.3 GO and KEGG enrichment analysis of differentially expressed(DE)lncRNAs.(a)GO enrichment analysis of DE lncRNAs between Tibetan(TCs)and Chahua chickens(CHs).KS represents the statistical significance of enriching the GO term,the smaller the KS value mean the more significant the enrichment;(b)KEGG enrichment analysis of DE lncRNAs between Tibetan and Chahua chickens

    A total of 1766 target genes were predicted for 73 DE miRNAs.The target genes were mainly enriched in the GO terms of vasculogenesis,positive regulation of fatty acid oxidation,glycogen metabolic process,angiogenesis,vasculature development,response to hypoxia,oxygen metabolic process,vasodilation,ATP binding,and enzyme binding(Fig.4a,Additional file 8:Table S7a).Notably,the enriched KEGG pathways were cysteine and methionine metabolism,amino sugar and nucleotide sugar metabolism,apoptosis,vascular smooth muscle contraction,fatty acid metabolism,carbon metabolism,vitamin B6metabolism,and the MAPK,Notch,VEGF,mTOR,and calcium signaling pathways(Fig.4b,Additional file 8:Table S7b).

    Fig.4 GO and KEGG enrichment analysis of differentially expressed(DE)miRNAs.(a)GO enrichment analysis of DE miRNAs between Tibetan(TCs)and Chahua chickens(CHs).The enrichment factor is the ratio of DE miRNA numbers annotated in the GO term to the total gene numbers annotated in the same GO term.(b)KEGG enrichment analysis of DE miRNAs between Tibetan and Chahua chickens

    Comparative parsing of lncRNAs and mRNAs and functional analysis of DE genes

    To compare the exon number,open reading frame(ORF),and expression levels of lncRNA and mRNA,we needed to consider the differences in structure and sequence.The number of exons corresponding to lncRNAs was mainly less than 8,whereas the number of exons in mRNAs was relatively large,ranging from 1 to>30(Additional file 4:Fig.S3a,b).The majority of lncRNA ORFs were mainly between 50 and 250nt long,while the length of the ORFs of mRNAs was mainly between 100 and 1,000nt(Additional file 4:Fig.S3c,d).Expressed mRNAs and lncRNAs were located on all chromosomes of the chicken genome(Additional file 4:Fig.S4).The average and maximum values of mRNA expression were higher compared to those of lncRNA expression(Fig.5a).

    The 354 DE genes were mainly enriched in GO terms of defense response,homeostatic process,regulation of lipid metabolic process,glucose homeostasis,blood microparticle,oxygen transport,blood vessel endothelial cell migration,blood circulation,and respiratory system development.Enriched KEGG pathways included fatty acid metabolism,vascular smooth muscle contraction,as well as the calcium,PPAR and MAPK signaling pathways(Fig.5b,Additional file 9:Table S8).Twelve DE genes(ACTC1,KCNMB4,NCS1,NGFR,ADAM8,CASQ2,IRF4,PTPRZ1,CALML3,ERBB4,ARR3,and NTSR1)were mainly associated with angiogenesis,blood circulation,hematopoiesis,as well as calcium and MAPK signaling pathways.Fifteen DE genes(SSTR5,NR1H4,HTR2C,APOA1,KCNB1,P3H2,CHST8,LYZ,HAO2,ACER1,ACSBG1,ELOVL2,ELOVL3,GBE,and NOX3)were related to glucose and carbohydrate metabolism,fatty acid metabolism,and lactate oxidation(Table 1).

    Fig.5 Comparative parsing of lncRNAs and mRNAs and functional analysis of differentially expressed(DE)genes.(a)A comparison of expression between lncRNAs and mRNAs.The box chart statistics from top to bottom are the maximum value,upper quartile,middle value,lower quartile,and the minimum value;(b)GO and KEGG enrichment analysis of DE genes between Tibetan(TC)and Chahua chicken(CH)

    Construction of the ceRNA network

    We constructed a ceRNA network that yielded 529 pairs of candidate ceRNAs(lncRNA-miRNA-mRNA)from 162 DE lncRNAs and 108 DE genes through 25 miRNA target-mediated relationships(Additional file 4:Fig.S5).Considering that the network contains enormous information,and each relationship cannot be displayed in the figure,we constructed a mini-ceRNA(DE lncRNA-DE miRNA-DE mRNA)network of important RNAs.We identified 10 known mRNAs that were related to hypoxic adaptation and involved in angiogenesis and blood circulation(NGFR,ACTC1,CASQ2,ERBB4,KCNMB4,NCS1,and NTSR1)(Fig.6a),as well as energy metabolism(SSTR5,NR1H4,and GBE)(Fig.6b).A total of 39 lncRNA-miRNA-mRNA interactions were identified in the mini-ceRNA network constructed with these 10 DE genes and their targeted 37 DE lncRNAs and 9 DE miRNAs(Table 1,Additional file 10:Table S9).

    Fig.6 Mini-ceRNA network constructed with 10 differentially expressed(DE)genes and their 37 target DE lncRNAs and 9 target DE miRNAs,involved in angiogenesis(a)and energy metabolism(b).Red color represents the up-regulated and blue color represents down-regulated levels

    Table 1 Potential key differentially expressed mRNAs(DE genes),their targeted differentially expressed miRNAs(DE miRNAs)/differentially expressed lncRNAs(DE lncRNAs),and their functions related to hypoxic adaption in the Tibetan chicken

    Table 1 Potential key differentially expressed mRNAs(DE genes),their targeted differentially expressed miRNAs(DE miRNAs)/differentially expressed lncRNAs(DE lncRNAs),and their functions related to hypoxic adaption in the Tibetan chicken(Continued)

    Expression levels of DE lncRNAs,DE miRNAs,and DE genes using qRT-PCR

    The expression levels of six DE lncRNAs(MSTR G.118362.24,MSTRG.25780.4,MSTRG.19949.13,MSTR G.29590.15,MSTRG.32082.2,and MSTRG.73129.28),six DE miRNAs(gga-miR-726-3p,gga-miR-6608-3p,gga-miR-6606-5p,novel-miR-676,gga-miR-460b-3p,and gga-miR-7442-5p),and seven DE genes(ARR3,NGFR,NTSR1,ADAM8,ACTC1,CASQ2,and CALML3)were measured using qRT-PCR to validate expression differences identified through RNA-seq.The results indicated that the expression levels of all RNAs were significantly different between TC and CH(P<0.05),and fold-changes in expression followed the same trend in qRT-PCR and RNA-seq(Fig.7).

    Fig.7 The expression of six differentially expressed(DE)lncRNAs(a),six DE miRNAs(c)and seven DE genes(e)validated with qRT-PCR in chicken embryo chorioallantoic membrane(CAM)tissue.Note:The Y-axis represents the expression value of RNA in chicken embryo CAM,and the X-axis represents the types of RNAs.Error bars represent SE of expression.*on the bars indicate P<0.05 and**indicate P<0.01 between TC and CH breeds.Tibetan chickens(TCs,n=4–8),and Chahua chickens(CHs,n=4–8).The fold changes of the six DE lncRNAs(b),six DE miRNAs(d)and seven DE genes(f)showed that qRT-PCR results were consistent with the RNA-seq data

    Discussion

    Through evolution,Tibetan chickens have developed genetic adaptability to survive in low-oxygen environments.Previous reports have shown that Tibetan chickens exhibit unique phenotypic and physiological characteristics,including enhanced blood oxygen binding capacity and increased blood circulation[12,17].The ceRNA hypothesis was reported a few years ago and has now been widely accepted,which has improved our understanding of human diseases,as well as led to important advances in animal research[36,37].The recently discovered noncoding RNAs represent a new type of master regulator that affects gene expression and modulates a variety of cellular processes.Based on these theories,we constructed a ceRNA(lncRNA-miRNA-mRNA)regulatory network to clarify the molecular mechanism of Tibetan chicken embryo CAM adaptation to hypoxic environments.

    The CAM,an important respiratory and circulatory organ for chicken embryo development,contains a large number of blood vessels,and hypoxia can induce increased blood vessel density[18,38].From the constructed ceRNA(lncRNA-miRNA-mRNA)network,we identified several miRNAs,including down-regulated novel-miR-589,novelmiR-815,novel-miR-85,novel-miR-669,as well as upregulated miR-6606-5p,novel-miR-676,novel-miR-589,novel-miR-815,novel-miR-85,and novel-miR-567.Among these,miR-6606-5p was significantly up-regulated in the Tibetan chicken,and its target differentially expressed genes(ACTC1,NCS1,and KCNMB4)were found to be involved in vascular development under hypoxic conditions.ACTC1 and NCS1,up-regulated in Tibetan chickens,were enriched in the GO term of blood circulation and regulated calcium ion transport;Ca2+is necessary for the activity of HIF-1,which is a major transcriptional regulator of cells and development in response to hypoxia[39,40].KCNMB4 has been reported to regulate blood pressure and was enriched in the GO terms of regulation of vasoconstriction,ion transport,and vascular smooth muscle contraction pathway[41].Consequently,these findings indicate that Tibetan chickens can improve blood circulation and stimulate the activity of HIF-1 through upregulation of ACTC1 and NCS1,and Tibetan chickens promote vasodilation through downregulation of KCNMB4 in CAM to adapt to hypoxic conditions.It is worth mentioning that with the exception of miR-6606-5p,other differentially expressed miRNAs,including miRNA-155,miRNA-302a,miR-302b-3p,miR-302b-5p,miR-302c-3p,miR-460b-3p,and miR-460b-5p,also play important roles in angiogenesis under hypoxia.The miRNA-155 promotes HIF-1α activity during prolonged hypoxia and participates in the PI3K/AKT pathway [42,43].Previous studies have shown that the miR-302 family(miR-302a/b/c/d)suppresses the proliferation,migration,and angiogenesis of vascular endothelial cells by targeting VEGFA[44–46].The miR-460b-3p and miR-460b-5p were identified in animal models of hypoxic pulmonary hypertension and were involved in the regulation of HIF-1α[47,48].Moreover,other differentially expressed genes(NGFR,ADAM8,CASQ2,IRF4,PTPRZ1,CALML3,ERBB4,ARR3,and NTSR1)were enriched in blood vessel development,angiogenesis,blood circulation,hematopoiesis,response to hypoxia,oxygen transport,vascular smooth muscle contraction,calcium,and MAPK signaling pathways.The regulatory mechanism of these genes in hypoxic adaptation of Tibetan chickens is the same as that of ACTC1,NCS1 and KCNMB4 genes.In general,they regulate angiogenesis through the up or downregulation of genes in their pathways.For example,NGFR is the receptor of NGF,which can induce chick CAM neovascularization[49].Studies have shown that ADAM8 is significantly induced by hypoxia and plays a role in the proliferation and migration of endothelial cells during angiogenesis[50,51].The explicit regulatory mechanism of these genes in hypoxic adaption of Tibetan chickens needs to be further studied.Therefore,we speculate that the enhanced tolerance of Tibetan chicken embryo CAM under hypoxic conditions is ascribed to miRNA-mediated modulation of the related target mRNAs,further regulating HIF,which enables Tibetan chicken to maintain hypoxic adaptation via the promotion of angiogenesis and blood circulation.

    Glucose uptake and carbohydrate metabolism are the basis for the maintenance of normal physiological functions in humans and animals.Under hypoxic conditions,oxygen and carbon dioxide metabolism mainly depend on mitochondrial respiration and make use of adenosine triphosphate(ATP).Recently,a number of studies have shown the key regulatory role of energy metabolism,including glucose,carbohydrate,and lipid metabolic processes,among others,during hypoxicadaptation[52–55].In the current research,from the constructed ceRNA network,we identified seven differentially expressed miRNAs (novel-miR-819,novel-miR-676,novel-miR-85,novel-miR-693,novelmiR-775,novel-miR-669,and novel-miR-867)and sixteen differentially expressed genes related to energy metabolism including glycosylation,glucose metabolism process,carbohydrate metabolism process,fatty acid metabolism,and ATP binding.For example,SSTR5 and NR1H4 are essential for glucose homeostasis and play a pivotal role in glucose metabolism in animals[56,57].Cortisol regulates the metabolism of mouse adipose cells through the serotonin receptor gene HTR2C,and genetic variation of the APOA1 gene is linked to lipid metabolism and cardiovascular disease risk[58,59].It should be considered that more differentially expressed genes(KCNB1,P3H2,CHST8,LYZ,HAO2,ACER1,ACSBG1,ELOVL2,ELOVL3,GBE,and NOX3)are targeted by key miRNAs,and several studies have shown that these mRNAs are also involved in energy metabolism.Therefore,we speculate that Tibetan chickens have enhanced energy metabolism due to the function of these RNAs,allowing for adaptation to hypoxic conditions.

    Accumulating evidence indicates that lncRNAs play the role of ceRNAs(or miRNA sponges)in a variety of biological processes,including high-altitude adaptation.Such as polymorphisms of LINC-PINT and LINC00599 are associated with high-altitude pulmonary edema in Chinese populations[60].Another study reported the expression profiles of lncRNAs in mice with highaltitude hypoxia-induced brain injury and provided new insights into the molecular mechanism of its treatment[61].In addition,a previous study reported the expression profiles of lncRNAs responsible for fatness and fatty acid composition traits in Tibetan pigs[62].In the current work,we found that the aforementioned differentially expressed genes involved in angiogenesis and energy metabolism were targeted by 37 differentially expressed lncRNAs in the ceRNA network,suggesting that these lncRNAs may also function as miRNA sponges and may play a role in the hypoxic adaptation of Tibetan chicken embryos with regard to angiogenesis and energy metabolism.

    Based on the ceRNA theory and the ceRNA network constructed in this study,we propose a mode of action of differentially expressed lncRNAs,miRNAs,and mRNAs during hypoxic adaptation of Tibetan chicken embryos(Fig.8).Under hypoxic conditions,miRNAs act as key regulators to modulate the up or down-regulation of important differentially expressed genes.As a consequence,the angiogenesis/blood circulation of chorioallantoic capillaries and energy metabolism,such as glucose/carbohydrate metabolism,are stimulated,leading to enhanced hypoxia adaptability of Tibetan chicken embryos.In this process,some lncRNAs act as ceRNAs to competitively bind the miRNA response element of miRNAs,which may indirectly affect the expression of mRNA.

    Fig.8 Proposed model of hypoxic adaptation in Tibetan chicken embryos.Under hypoxic conditions,miRNAs act as key regulators that directly target important differentially expressed(DE)genes.In this process,some lncRNAs can act as ceRNAs to competitively bind the miRNA response element(MRE)of miRNAs,indirectly affecting the expression of mRNAs

    Conclusions

    In conclusion,389 DE lncRNAs,73 DE miRNAs,and 354 DE genes were identified between Tibetan and Chahua chickens.The transcriptomic data revealed several key candidate ceRNAs (DE lncRNAs-DE miRNAs-DE genes)that may play high-priority roles in the hypoxic adaptation of Tibetan chickens by regulating angiogenesis and energy metabolism.These results provide insights into the molecular mechanisms of hypoxic adaptation regulatory networks from the perspective of coding and non-coding RNAs.

    Abbreviations

    lncRNAs:Long non-coding RNA;miRNA:MicroRNA;TC:Tibetan chicken;CH:Chahua chicken;CAM:Chorioallantoic membrane;DE:Differentially expressed gene;DE lncRNA:Differentially expressed lncRNA;DE miRNA:Differentially expressed miRNA;ceRNA:Competing endogenous RNA;GO:Gene Ontology;KEGG:Kyoto Encyclopedia of Genes and Genomes;qRTPCR:Quantitative real-time PCR;ORF:Open reading frame

    Supplementary Information

    The online version contains supplementary material available at https://doi.org/10.1186/s40104-021-00582-2.

    Additional file 1.Table S1.The qRT-PCR primers of six DE lncRNAs,six DE miRNAs and seven DE genes for identification by RNA-seq.

    Additional file 2:Table S2.List of software and the parameters used in this study.

    Additional file 3:Table S3.Summary of sequencing reads aligned with the Gallus gallus genome.

    Additional file 4:Figure S1.Lengths of mRNAs and lncRNAs distribution.Figure S2.Figure of hierarchical clustering of 389 DE lncRNAs(a),73 DE miRNAs(b),and 354 DE genes(c)in Tibetan and Chahua chickens.Figure S3.Comparative analysis of exon and ORF lengths of lncRNAs and mRNAs.Figure S4.Distribution of the DE lncRNAs and DE mRNAs on the different chromosomes.Figure S5.Network of ceRNA(lncRNA-miRNA-mRNA)including 162 DE lncRNAs,25 DE miRNAs,and 108 DE genes.

    Additional file 5:Table S4.The 3194 known lncRNAs,471 known miRNAs,and 968 novel miRNAs.

    Additional file 6:Table S5.The 389 differentially expressed lncRNAs(DE lncRNAs),73 DE miRNAs,and 354 DE genes.

    Additional file 7:Table S6.GO and KEGG pathway enrichment of differentially expressed lncRNA(DE lncRNA)target genes.

    Additional file 8:Table S7.GO and KEGG pathway enrichment of differentially expressed miRNA(DE miRNA)target genes.

    Additional file 9:Table S8.GO and KEGG enrichment of differentially expressed mRNAs(DE genes).

    Additional file 1.:Table S9.Important difference in circRNA-miRNA-mRNA interactions based on the competing endogenous RNA(ceRNA)network.

    Acknowledgements

    The authors would like to thank BioMarker Technology(Beijing,China)for the sequencing services and bioinformatic analysis.

    Authors’contributions

    Y.Z.,and H.Z.designed the research;Y.Z.,W.S.,Y.L.,and B.Z.collected samples,extracted RNAs,and performed qRT-PCR verification experiments;Y.Z.performed most of the bioinformatics studies;Y.Z.,H.Z.,and W.K.wrote and revised the manuscript;and all authors reviewed the manuscript.The authors read and approved the final manuscript.

    Funding

    This work was supported by the National Natural Science Foundation of China(31972532),the China Agricultural Research System(CARS-40-K05),and the Innovation Base Cultivation and Development Project(Z171100002217072).

    Availability of data and materials

    All RNA-seq data were deposited in the Gene Expression Omnibus under accession number GSE160324.

    Declarations

    Ethics approval and consent to participate

    All animal work was conducted following the guidelines for the care and use of experimental animals.The experiments and animal care protocol were approved by the animal welfare committee of the State Key Laboratory for Agro-Biotechnology of the China Agricultural University(approval number,XK257).

    Consent for publication

    Not applicable.

    Competing interests

    The authors declare that they have no competing interests.

    Received:14 November 2020 Accepted:8 March 2021

    50天的宝宝边吃奶边哭怎么回事| 免费在线观看影片大全网站| 99国产精品一区二区三区| 精品久久久久久久久久免费视频| 国产人伦9x9x在线观看| 色婷婷久久久亚洲欧美| 精品电影一区二区在线| 深夜精品福利| a在线观看视频网站| e午夜精品久久久久久久| 人妻丰满熟妇av一区二区三区| 757午夜福利合集在线观看| 久久久久亚洲av毛片大全| 99热只有精品国产| 亚洲成av人片免费观看| 人妻丰满熟妇av一区二区三区| 亚洲五月色婷婷综合| 久久久国产欧美日韩av| 国产亚洲精品久久久久5区| 可以在线观看的亚洲视频| 久久久久久久久免费视频了| 欧美日韩黄片免| 精品国产美女av久久久久小说| 国产成人精品久久二区二区91| e午夜精品久久久久久久| 国产1区2区3区精品| 亚洲中文字幕日韩| 99精品欧美一区二区三区四区| 亚洲精品久久国产高清桃花| 成人18禁高潮啪啪吃奶动态图| 美女国产高潮福利片在线看| 亚洲国产看品久久| 亚洲五月天丁香| 成人亚洲精品av一区二区| 美国免费a级毛片| 韩国精品一区二区三区| 天天添夜夜摸| 国产精品亚洲av一区麻豆| 国产成人精品久久二区二区免费| 亚洲精品在线观看二区| 制服丝袜大香蕉在线| 国产麻豆成人av免费视频| 国内精品久久久久精免费| 欧美成人午夜精品| av中文乱码字幕在线| 亚洲成人国产一区在线观看| 久久久国产精品麻豆| 黄频高清免费视频| 久久亚洲真实| 国产精品久久视频播放| 亚洲美女黄片视频| 国产在线观看jvid| or卡值多少钱| 久久香蕉激情| 91大片在线观看| 亚洲专区字幕在线| 老熟妇仑乱视频hdxx| 日韩 欧美 亚洲 中文字幕| 视频区欧美日本亚洲| 夜夜爽天天搞| 制服诱惑二区| av片东京热男人的天堂| 黄色丝袜av网址大全| 好看av亚洲va欧美ⅴa在| 少妇粗大呻吟视频| 日本五十路高清| 99热只有精品国产| 免费搜索国产男女视频| 久久久久久亚洲精品国产蜜桃av| 久久久国产成人免费| 一边摸一边抽搐一进一小说| 国产成人欧美在线观看| 一级毛片女人18水好多| 麻豆成人av在线观看| 婷婷精品国产亚洲av在线| 男女做爰动态图高潮gif福利片 | 久久人人精品亚洲av| 啦啦啦免费观看视频1| 免费av毛片视频| 嫩草影院精品99| 国产精品久久久久久人妻精品电影| 波多野结衣一区麻豆| 日本在线视频免费播放| 搡老熟女国产l中国老女人| 久久伊人香网站| 国产国语露脸激情在线看| 性欧美人与动物交配| 精品一区二区三区av网在线观看| 精品久久久精品久久久| 亚洲电影在线观看av| 国产三级在线视频| 中文字幕人妻丝袜一区二区| 免费女性裸体啪啪无遮挡网站| 午夜福利视频1000在线观看 | 日韩中文字幕欧美一区二区| 国产成人系列免费观看| 色综合站精品国产| 一进一出抽搐动态| 亚洲 欧美一区二区三区| 亚洲熟女毛片儿| 男女做爰动态图高潮gif福利片 | 欧美成人免费av一区二区三区| 成人精品一区二区免费| 免费在线观看影片大全网站| 亚洲视频免费观看视频| 淫秽高清视频在线观看| 在线观看日韩欧美| 啦啦啦观看免费观看视频高清 | 国产成人精品久久二区二区91| 麻豆av在线久日| 国产国语露脸激情在线看| 色综合欧美亚洲国产小说| 黑人欧美特级aaaaaa片| 女人爽到高潮嗷嗷叫在线视频| 麻豆成人av在线观看| 无遮挡黄片免费观看| 国产精品,欧美在线| 久久人妻福利社区极品人妻图片| 如日韩欧美国产精品一区二区三区| 久久精品国产99精品国产亚洲性色 | 国产99久久九九免费精品| 女人被躁到高潮嗷嗷叫费观| 国内毛片毛片毛片毛片毛片| 日本三级黄在线观看| 女同久久另类99精品国产91| 欧美日本中文国产一区发布| 最新在线观看一区二区三区| 老熟妇乱子伦视频在线观看| АⅤ资源中文在线天堂| 精品久久蜜臀av无| 亚洲 欧美 日韩 在线 免费| 日韩视频一区二区在线观看| 成人18禁在线播放| 欧美日韩黄片免| 久久精品人人爽人人爽视色| 成人亚洲精品av一区二区| 久久国产精品影院| 久久久久久大精品| 亚洲国产精品久久男人天堂| 老熟妇仑乱视频hdxx| 黑丝袜美女国产一区| 性欧美人与动物交配| 亚洲午夜精品一区,二区,三区| 丰满的人妻完整版| 一进一出抽搐动态| 国产精品一区二区精品视频观看| 国产三级黄色录像| 美女扒开内裤让男人捅视频| 午夜视频精品福利| 国产伦人伦偷精品视频| 国产在线精品亚洲第一网站| 涩涩av久久男人的天堂| 黄片大片在线免费观看| 丁香六月欧美| avwww免费| 久久久久九九精品影院| 91在线观看av| 在线免费观看的www视频| 国产欧美日韩综合在线一区二区| 色老头精品视频在线观看| 久久天堂一区二区三区四区| e午夜精品久久久久久久| 一边摸一边抽搐一进一出视频| 宅男免费午夜| 久久久国产成人精品二区| 黄色成人免费大全| 天天添夜夜摸| 成人18禁高潮啪啪吃奶动态图| 亚洲全国av大片| 黄频高清免费视频| 怎么达到女性高潮| 久久国产精品男人的天堂亚洲| 视频在线观看一区二区三区| 色综合婷婷激情| 午夜福利欧美成人| bbb黄色大片| 熟妇人妻久久中文字幕3abv| 国产私拍福利视频在线观看| 桃色一区二区三区在线观看| 日本黄色视频三级网站网址| 丝袜美腿诱惑在线| 精品国产乱子伦一区二区三区| 99国产极品粉嫩在线观看| 大型av网站在线播放| 黄片播放在线免费| 国产成人系列免费观看| 在线观看日韩欧美| 欧美 亚洲 国产 日韩一| 国产av一区二区精品久久| 中文字幕人妻熟女乱码| 国产精品98久久久久久宅男小说| 亚洲狠狠婷婷综合久久图片| 最新美女视频免费是黄的| 大型黄色视频在线免费观看| 午夜福利在线观看吧| 窝窝影院91人妻| 午夜福利18| 亚洲五月天丁香| x7x7x7水蜜桃| 亚洲成人免费电影在线观看| 天天添夜夜摸| 欧美一级a爱片免费观看看 | 好男人电影高清在线观看| 欧美日韩乱码在线| 免费在线观看日本一区| 欧美+亚洲+日韩+国产| 国产真人三级小视频在线观看| 亚洲成人久久性| 男女床上黄色一级片免费看| bbb黄色大片| 亚洲国产日韩欧美精品在线观看 | 国产精品久久电影中文字幕| 国产精品秋霞免费鲁丝片| 日韩成人在线观看一区二区三区| 美女大奶头视频| 老汉色av国产亚洲站长工具| 午夜免费激情av| 国产精品 欧美亚洲| 日韩精品中文字幕看吧| 国产成人系列免费观看| 淫秽高清视频在线观看| 日韩大码丰满熟妇| 久久热在线av| 精品不卡国产一区二区三区| 国产精品久久久久久精品电影 | 最好的美女福利视频网| 制服丝袜大香蕉在线| 精品久久久久久成人av| 在线观看66精品国产| 亚洲激情在线av| 桃红色精品国产亚洲av| 国产精品一区二区免费欧美| 悠悠久久av| 岛国视频午夜一区免费看| 制服丝袜大香蕉在线| 深夜精品福利| 老鸭窝网址在线观看| 国产av一区在线观看免费| 老汉色∧v一级毛片| 久久香蕉国产精品| 国产成人欧美在线观看| 国产精品免费视频内射| 99在线人妻在线中文字幕| 很黄的视频免费| 91在线观看av| 久久婷婷人人爽人人干人人爱 | 丰满人妻熟妇乱又伦精品不卡| 悠悠久久av| 午夜亚洲福利在线播放| 丰满人妻熟妇乱又伦精品不卡| 不卡一级毛片| 美女高潮到喷水免费观看| 最新在线观看一区二区三区| 一边摸一边抽搐一进一出视频| 丰满的人妻完整版| 欧美日韩一级在线毛片| 色播在线永久视频| 男女下面进入的视频免费午夜 | 国产一卡二卡三卡精品| 一边摸一边抽搐一进一小说| 两个人看的免费小视频| 亚洲欧美激情在线| videosex国产| 国产成人av教育| 亚洲中文日韩欧美视频| 亚洲自拍偷在线| 高清黄色对白视频在线免费看| 最新美女视频免费是黄的| 一二三四在线观看免费中文在| 美女高潮喷水抽搐中文字幕| 日韩大码丰满熟妇| 成人欧美大片| 日本五十路高清| 亚洲精品中文字幕一二三四区| www.熟女人妻精品国产| 国产欧美日韩一区二区三区在线| 亚洲欧美一区二区三区黑人| 老司机午夜福利在线观看视频| 一级a爱视频在线免费观看| 欧美午夜高清在线| 一级作爱视频免费观看| 久久精品国产综合久久久| 日韩欧美一区视频在线观看| 成年人黄色毛片网站| 搞女人的毛片| 国产xxxxx性猛交| 国产av又大| 亚洲av片天天在线观看| 午夜福利18| 91国产中文字幕| 成人18禁在线播放| 曰老女人黄片| 天堂动漫精品| 人妻久久中文字幕网| 亚洲中文字幕日韩| 欧美日韩中文字幕国产精品一区二区三区 | 亚洲aⅴ乱码一区二区在线播放 | 婷婷精品国产亚洲av在线| 一边摸一边抽搐一进一出视频| 国产精品久久久久久精品电影 | 欧美日韩一级在线毛片| 免费一级毛片在线播放高清视频 | 国产精品亚洲av一区麻豆| 麻豆成人av在线观看| 91成年电影在线观看| 亚洲avbb在线观看| 999精品在线视频| 亚洲九九香蕉| 日本在线视频免费播放| 亚洲 欧美 日韩 在线 免费| 亚洲精品美女久久av网站| 黄色片一级片一级黄色片| 91九色精品人成在线观看| 一级a爱视频在线免费观看| 国产精品一区二区三区四区久久 | 最近最新中文字幕大全免费视频| 亚洲欧洲精品一区二区精品久久久| 日韩欧美国产一区二区入口| 午夜福利在线观看吧| 两性午夜刺激爽爽歪歪视频在线观看 | 黄色成人免费大全| 18禁美女被吸乳视频| 在线观看舔阴道视频| 黄色视频不卡| 我的亚洲天堂| 欧美黄色片欧美黄色片| 国产亚洲欧美精品永久| 97人妻精品一区二区三区麻豆 | bbb黄色大片| 欧美日韩黄片免| 丝袜人妻中文字幕| 黑人巨大精品欧美一区二区mp4| 高清黄色对白视频在线免费看| 女人爽到高潮嗷嗷叫在线视频| 美女高潮到喷水免费观看| 高潮久久久久久久久久久不卡| 91麻豆精品激情在线观看国产| 婷婷丁香在线五月| 国产不卡一卡二| www国产在线视频色| av中文乱码字幕在线| 18禁黄网站禁片午夜丰满| 老司机靠b影院| 大香蕉久久成人网| 一级黄色大片毛片| 一本久久中文字幕| 激情视频va一区二区三区| 又大又爽又粗| 午夜两性在线视频| 禁无遮挡网站| 身体一侧抽搐| 99在线视频只有这里精品首页| 亚洲电影在线观看av| 日韩高清综合在线| 亚洲熟女毛片儿| 麻豆av在线久日| 777久久人妻少妇嫩草av网站| 国产精品av久久久久免费| 国产成人免费无遮挡视频| 男女床上黄色一级片免费看| 精品国产乱子伦一区二区三区| 美女高潮喷水抽搐中文字幕| 怎么达到女性高潮| 亚洲中文av在线| 脱女人内裤的视频| 欧美激情极品国产一区二区三区| 成人欧美大片| 国产私拍福利视频在线观看| 极品人妻少妇av视频| 久久国产精品影院| 国产精品,欧美在线| www.999成人在线观看| 自拍欧美九色日韩亚洲蝌蚪91| 身体一侧抽搐| а√天堂www在线а√下载| 午夜免费鲁丝| 午夜a级毛片| 很黄的视频免费| 亚洲专区国产一区二区| 亚洲无线在线观看| 精品久久久久久成人av| 少妇粗大呻吟视频| 两个人视频免费观看高清| 国产av一区在线观看免费| 久久久久久免费高清国产稀缺| 欧美激情极品国产一区二区三区| 丁香欧美五月| 夜夜看夜夜爽夜夜摸| 99香蕉大伊视频| 国产单亲对白刺激| 日本精品一区二区三区蜜桃| 可以在线观看的亚洲视频| 亚洲自偷自拍图片 自拍| 婷婷精品国产亚洲av在线| 这个男人来自地球电影免费观看| 看黄色毛片网站| 国产精品香港三级国产av潘金莲| 国产亚洲av嫩草精品影院| 日韩精品中文字幕看吧| 久久久国产精品麻豆| 国内久久婷婷六月综合欲色啪| 亚洲全国av大片| 1024视频免费在线观看| 亚洲欧美日韩高清在线视频| 亚洲伊人色综图| 操出白浆在线播放| 亚洲欧美日韩另类电影网站| 免费观看人在逋| 午夜免费鲁丝| 在线观看免费视频网站a站| 十八禁人妻一区二区| 国产精品乱码一区二三区的特点 | 香蕉国产在线看| 一边摸一边抽搐一进一出视频| 天堂√8在线中文| 亚洲第一电影网av| 久久久国产欧美日韩av| 中文字幕人妻熟女乱码| 久热爱精品视频在线9| 国产av精品麻豆| 中文字幕色久视频| 级片在线观看| 午夜两性在线视频| 中亚洲国语对白在线视频| 久久国产亚洲av麻豆专区| 久久久国产欧美日韩av| 制服诱惑二区| 少妇熟女aⅴ在线视频| 国产一区二区在线av高清观看| 级片在线观看| 久久人人爽av亚洲精品天堂| a级毛片在线看网站| 久久久精品国产亚洲av高清涩受| 午夜两性在线视频| 老鸭窝网址在线观看| 色综合婷婷激情| 久热这里只有精品99| 国产97色在线日韩免费| 久久精品91蜜桃| 亚洲无线在线观看| 国内久久婷婷六月综合欲色啪| 天天躁夜夜躁狠狠躁躁| 国产亚洲av嫩草精品影院| 国产亚洲精品久久久久5区| 满18在线观看网站| 中出人妻视频一区二区| 精品国产超薄肉色丝袜足j| 久久天堂一区二区三区四区| 中亚洲国语对白在线视频| 亚洲熟妇中文字幕五十中出| 欧美日韩瑟瑟在线播放| 男人舔女人下体高潮全视频| 午夜精品久久久久久毛片777| 9热在线视频观看99| 一夜夜www| 一区二区三区国产精品乱码| АⅤ资源中文在线天堂| 成人18禁高潮啪啪吃奶动态图| 国产91精品成人一区二区三区| 制服人妻中文乱码| 久久婷婷人人爽人人干人人爱 | 午夜激情av网站| 18禁美女被吸乳视频| 日韩精品中文字幕看吧| 午夜福利一区二区在线看| 成人亚洲精品av一区二区| av在线播放免费不卡| 91麻豆av在线| 波多野结衣av一区二区av| 国产高清有码在线观看视频 | 亚洲熟妇熟女久久| 国产成人啪精品午夜网站| 亚洲 国产 在线| 91国产中文字幕| 丁香六月欧美| 好看av亚洲va欧美ⅴa在| 免费高清在线观看日韩| 男男h啪啪无遮挡| 1024香蕉在线观看| 国产精品久久电影中文字幕| 亚洲国产毛片av蜜桃av| 亚洲一区高清亚洲精品| 国产99白浆流出| 色播亚洲综合网| 看黄色毛片网站| 欧美激情极品国产一区二区三区| 高潮久久久久久久久久久不卡| 国产欧美日韩一区二区三| 国产精品 欧美亚洲| avwww免费| 亚洲av电影不卡..在线观看| 国产精品久久久av美女十八| 人人妻人人澡人人看| 亚洲精品在线美女| 青草久久国产| 国产午夜福利久久久久久| 如日韩欧美国产精品一区二区三区| 亚洲五月色婷婷综合| 麻豆av在线久日| 国产精品1区2区在线观看.| 丝袜美腿诱惑在线| 9191精品国产免费久久| 亚洲狠狠婷婷综合久久图片| 99久久99久久久精品蜜桃| 两性午夜刺激爽爽歪歪视频在线观看 | 久久人妻熟女aⅴ| 国产成+人综合+亚洲专区| 国内毛片毛片毛片毛片毛片| 精品久久久久久久人妻蜜臀av | 午夜福利视频1000在线观看 | 国产精品野战在线观看| 美女免费视频网站| 国内精品久久久久久久电影| 一区二区三区高清视频在线| 成年人黄色毛片网站| 精品一区二区三区视频在线观看免费| 免费看美女性在线毛片视频| 黄色视频,在线免费观看| 男男h啪啪无遮挡| 亚洲精品美女久久av网站| 亚洲欧洲精品一区二区精品久久久| 亚洲美女黄片视频| 麻豆国产av国片精品| 69精品国产乱码久久久| 午夜激情av网站| 男女下面插进去视频免费观看| 久久青草综合色| 国产精品二区激情视频| 久久人妻av系列| av在线播放免费不卡| 欧美成人免费av一区二区三区| 久久欧美精品欧美久久欧美| 成人国语在线视频| 午夜福利,免费看| 中文字幕久久专区| 国产成人精品久久二区二区91| 999精品在线视频| 久久久国产成人精品二区| 黄色女人牲交| 精品卡一卡二卡四卡免费| 久久精品91无色码中文字幕| 亚洲精品国产区一区二| 精品日产1卡2卡| www.999成人在线观看| 黑人操中国人逼视频| 欧美日韩乱码在线| 色老头精品视频在线观看| 性色av乱码一区二区三区2| 黄色a级毛片大全视频| 黄色a级毛片大全视频| 黄片小视频在线播放| 久久这里只有精品19| 丝袜在线中文字幕| 性欧美人与动物交配| 每晚都被弄得嗷嗷叫到高潮| 又大又爽又粗| 久久精品91蜜桃| 国产精品日韩av在线免费观看 | 国产极品粉嫩免费观看在线| 搞女人的毛片| 亚洲九九香蕉| 嫁个100分男人电影在线观看| 亚洲一区中文字幕在线| 香蕉丝袜av| 国产成人精品久久二区二区免费| 在线天堂中文资源库| 国产熟女午夜一区二区三区| 激情视频va一区二区三区| 久久午夜亚洲精品久久| 亚洲五月色婷婷综合| 国产成人av教育| 亚洲激情在线av| 成人精品一区二区免费| 在线观看午夜福利视频| 亚洲成人精品中文字幕电影| 9色porny在线观看| 欧美激情高清一区二区三区| 亚洲成人久久性| 亚洲男人的天堂狠狠| 美女高潮喷水抽搐中文字幕| 少妇熟女aⅴ在线视频| 在线观看免费午夜福利视频| 老熟妇仑乱视频hdxx| e午夜精品久久久久久久| 欧美黑人欧美精品刺激| 亚洲午夜精品一区,二区,三区| 精品欧美国产一区二区三| 亚洲国产精品sss在线观看| 人人澡人人妻人| 久久国产精品男人的天堂亚洲| 欧美中文综合在线视频| 亚洲精品国产精品久久久不卡| 9热在线视频观看99| 热99re8久久精品国产| 男人舔女人下体高潮全视频| 一级毛片女人18水好多| 极品教师在线免费播放| 国产伦一二天堂av在线观看| 日韩大码丰满熟妇| 搡老岳熟女国产| 精品一区二区三区av网在线观看| cao死你这个sao货| 99国产极品粉嫩在线观看| 国产精品亚洲av一区麻豆| 久久精品国产亚洲av香蕉五月| 丝袜美足系列| 给我免费播放毛片高清在线观看| 日本免费一区二区三区高清不卡 | 国产私拍福利视频在线观看| 日日摸夜夜添夜夜添小说| 久久久久国内视频| 亚洲国产精品999在线| 精品国内亚洲2022精品成人| 国产高清有码在线观看视频 | 天堂影院成人在线观看|