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

    Long-Term Cultivation and Meta-Omics Reveal Methylotrophic Methanogenesis in Hydrocarbon-Impacted Habitats

    2023-11-14 08:03:06YiFnLiuJingChenZhongLinLiuZhoWeiHouBoLingLiYingWngLeiZhouLiBinShouDnDnLinShiZhongYngJinFengLiuXioLinWuJiDongGuBoZhongMu
    Engineering 2023年5期

    Yi-Fn Liu, Jing Chen, Zhong-Lin Liu, Zho-Wei Hou, Bo Ling,e, Li-Ying Wng,e, Lei Zhou,e,Li-Bin Shou,e, Dn-Dn Lin,e, Shi-Zhong Yng,e, Jin-Feng Liu,e, Xio-Lin Wu,*, Ji-Dong Gu,Bo-Zhong Mu,e,*

    a State Key Laboratory of Bioreactor Engineering and School of Chemistry and Molecular Engineering, East China University of Science and Technology, Shanghai 200237, China

    b Exploration and Development Research Institute of Daqing Oilfield Company Limited, PetroChina, Daqing 163712, China

    c Environmental Engineering and Engineering Group, Guangdong Technion Israel Institute of Technology, Shantou 515063, China

    d Shanghai Institute of Pollution Control and Ecological Security, Shanghai 200092, China

    e Engineering Research Center of Microbial Enhanced Oil Recovery (MEOR), East China University of Science and Technology, Shanghai 200237, China

    Keywords:Methanogenic hydrocarbon degradation Oily sludge Bioremediation Alkanes

    ARTICLEINFO The microbial conversion of alkanes to methane in hydrocarbon-contaminated environments is an intrinsic bioremediation strategy under anoxic conditions.However, the mechanism of microbial methanogenic alkane degradation is currently unclear.Under ten years of continuous efforts, we obtained a methanogenic n-alkane-degrading (C15–C20) enrichment culture that exhibited sustained improvements in the kinetic properties of methane production.The integrated metagenomic and metatranscriptomic analyses revealed that n-alkanes were mainly attacked by members of Desulfosarcinaceae, Firmicutes,and Synergistetes using the fumarate addition strategy, and were then further degraded in a common effort by Tepidiphilus members.Meanwhile, the abundant members of Anaerolineaceae were mainly responsible for cell debris recycling.However, according to the metatranscriptomic analyses, methane was predicted to be produced mainly via H2-dependent methylotrophic methanogenesis,primarily from necromass-derived trimethylamine mediated by Methanomethyliaceae within the candidate phylum Verstraetearchaeota.These findings reveal that H2-dependent methylotrophic methanogens, as well as methylotrophic methanogens, may play important ecological roles in the carbon cycle of hydrocarbon enriched subsurface ecosystems.

    1.Introduction

    The worldwide use and distribution of crude oil and its derivatives has led to hydrocarbon contamination in various environments [1].Alkanes, which comprise the major fraction of crude oils, are commonly found in hydrocarbon-contaminated habitats[2].Environmental contamination with high organic loads usually results in the development of anoxic conditions,under which alkanes can be naturally converted to methane via microorganisms[3]after nitrate,iron(III),and sulfate are depleted as electron acceptors[4].Despite having been studied for two decades, beginning with the first definitive enrichment culture that produced methane from n-hexadecane [5], the underlying mechanism of methanogenic alkane degradation has not yet been fully resolved.

    Several studies have shown that methanogenic alkane degradation requires the concerted efforts of metabolically different types of microorganisms [6], including an initial attack on the inert C–H bond of alkane compounds by syntrophic/fermentative bacteria,subsequent conversion to methanogenic substrates (such as formate, H2, and acetate), and methanogenic archaea that consume these substrates to a very low level to keep the overall biodegradation process thermodynamically favorable [7].

    Among the bacterial members, Deltaproteobacteria, especially the genus Smithella,are continuously detected and often dominant in low- to high-temperature methanogenic alkane degradation enrichment cultures from different sources [5,8–12].Increasing studies based on single-cell and metagenomic sequencing have revealed that Smithella-like genomes contain genes encoding canonical alkylsuccinate synthase (ASS) [9,10], which catalyzes alkane activation via the addition of fumarate.Meanwhile,studies based on metatranscriptomic and reverse transcriptionpolymerase chain reaction(PCR)analyses have definitively demonstrated that ASS gene(assA)affiliated with Smithella are expressed during methanogenic alkane degradation [9,10].Therefore, it has been proposed that members of Smithella are major‘‘alkane breakers,”via fumarate addition in methanogenic alkane-degrading cultures [11,12].

    Archaeal communities in methanogenic alkane-degrading cultures, however, are often made up of hydrogenotrophic (e.g.,Methanocalculus and Methanoculleus) and acetoclastic methanogens(e.g.,Methanosaeta)[6,10,13].In a study considering a methanogenic hexadecane-degrading enrichment culture [5], more than 2/3 of the methane was predicted to be produced from acetate,and thermodynamic calculations indicated that hexadecane conversion to methane is driven primarily via acetoclastic methanogenesis [14].Meanwhile, by modeling isotopic signatures for gaseous CO2and methane collected from heavily degraded oils in the mesothermic Peace River Oil Sands, it was suggested that hydrogenotrophic methanogenesis likely predominates in subsurface degraded oil reservoirs where the most easily degradable compounds are straight-chain n-alkanes [15].A similar result was obtained from a DNA-stable isotope probing(SIP)-based study of a methanogenic13C-hexadecane degrading consortium, where hydrogenotrophic methanogens were found to be dominant in the heavy DNA fractions [16].Another radiotracer experiment using14C-labeled bicarbonate and acetate as substrates revealed that the rates of hydrogenotrophic methanogenesis were 42–68 times greater than those of acetoclastic methanogenesis.Furthermore, a metatranscriptomic analysis of alkane-degrading enrichment consortia [9,17] and production water samples [18]collected from oil reservoirs suggested that hydrogenotrophic and acetoclastic methanogenesis are the main pathways for converting alkane-derived substrates into methane.These results collectively support the exclusive role of hydrogenotrophic CO2-reducing and acetoclastic methanogenesis in the decomposition of oil alkanes in anoxic hydrocarbon-associated environments[8].However, alternative methanogenesis pathways, such as methylotrophic methanogenesis and H2-dependent methylotrophic methanogenesis, are often overlooked, despite the fact that methanogens(such as the Methanohalophilus,Methanolobulus,and Methanosarcina species) that can use methyl compounds to produce methane were commonly detected in hydrocarbonassociated environments [19,20].Although methane generation from methylated C1compounds is possible, it is difficult to envision those methylated compounds can be produced during the degradation of alkanes or other hydrocarbon molecules [21].

    In a previous study, Wang et al.[22] obtained methanogenic enrichment cultures from oily sludge(OS),and the following transfer was amended with n-alkanes (C15–C20) as the sole carbon and energy source.Subsequent transfer by Liang et al.[23] into fresh medium resulted in a methanogenic alkane-degrading enrichment culture after more than 1300 days.A phylogenetic analysis based on the 16S ribosomal RNA (rRNA) gene and methyl coenzyme M reductase gene (mcrA) amplicon libraries showed that members of Anaerolineaceae and Methanosaeta dominated the bacterial and archaeal communities, respectively [19].Hence, Anaerolineaceae was proposed to play an important role in fermentation and alkane oxidation, and Methanosaeta was mainly responsible for methane production via the acetoclastic pathway [19].However, because of a lack of cultured representatives and genomic information, the exact functions of Anaerolineaceae and other organisms in the culture remain largely unknown.

    To study the mechanism of methanogenic alkane degradation and resolve the roles of microorganisms in the community, we conducted a third transfer from a long-term n-alkane-degrading(C15–C20) enrichment culture.The cells were harvested at the late exponential phase of methane production and investigated for microbial metabolic potential and transcriptional activity using an integrated metagenomic and metatranscriptomic method.

    2.Material and methods

    2.1.Establishing enrichment cultures

    The initial enrichment culture was established in a previous study using OS collected from the Shanghai Oil Refinery (China)[22].Enrichment cultures in this study (the third transfer) were established by transferring approximately 10 mL (20% v/v) of the second transfer culture with a mixture of n-alkanes(C15–C20)incubated long-term as previously described[23,24]in a 120 mL autoclaved serum bottles that included 40 mL sterile basal medium(0.20 g·L-1NaCl, 1.20 g·L-1MgCl2·6H2O, 0.10 g·L-1CaCl2·2H2O,0.25 g·L-1NH4Cl, 0.20 g·L-1KH2PO4, 1.30 g·L-1KCl, 2.50 g·L-1NaHCO3, and 0.0001 g·L-1rezasurin), and 1.0 mL·L-1of vitamin stock solution and trace element stock solution [20].The basal medium was reduced by adding 0.50 g·L-1Na2S·9H2O to the medium and adjusted to a final pH of 7.2,after which it was sealed with butyl rubber stoppers and aluminum crimps (Bellco Glass, Inc.,USA).Each sample (three replicates) was supplemented with 15 μmol of each n-alkane (C15, C16, C17, C18, C19, and C20) (≥99%;Sigma-Aldrich, USA) as described previously [23], and incubated at 37 °C in the dark for more than three years.

    2.2.Methane detection

    Methane accumulation was periodically monitored via gas chromatography (GC; Agilent 6890; Agilent Technologies, Inc.,USA) using the standard calibration curve for methane(R2= 0.994, n = 6).Specifically, 200 mL of headspace gas was injected into the GC using a sterile micro-syringe.The procedures was as follows: the initial column temperature was set to 60 °C for 12 min, increased to 200 °C at a rate of 15 °C·min-1, and then maintained at 200°C for 24 min.The temperatures of the injector,flame ionization detector, and thermal conductivity detector were maintained at 200 °C.

    2.3.Kinetic model of methane production

    The modified Gompertz model(Eq.(1))is one of the most suitable models for fitting the growth of microorganisms [25].Therefore, the modified Gompertz model was used to process the experimental data for the lag phase, maximum methane production rate, and methane yield potential under the assumption that the rate of biogas production is proportional to the microbial activity in the bioreactor.The modified Gompertz model is as follows:

    where Y represents the cumulative biogas production(μmol)in the bioreactor at an incubation time, t (d), Amis the maximum biogas yield potential (μmol), Vmis the maximum biogas production rate(μmol·d-1), tlagis the lag phase (d), and e is the natural constant(2.7183).Note that the Am, Vm, and tlagconstants were estimated based on the triplicate datasets using Origin8.0(Originlab Corporation, USA).

    2.4.PCR amplification

    The primers for the canonical assA genes have been described previously [19].The PCR conditions were as follows: 95 °C for 5 min, followed by 38 cycles of 95 °C for 30 s, 52–58 °C for 45 s,72 °C for 60–90 s, and a final elongation step at 72 °C for 10 min.The annealing temperature and extension time for each primer were modified according to the manufacturer’s instructions(Appendix A Table S1).

    2.5.Chemical analysis

    Residual n-alkanes and intermediate fatty acids in the enrichment culture were detected at the end of incubation using gas chromatography–mass spectrograph (GC–MS; Agilent 6890 GC)fitted with an HP-5MS capillary column (30 m × 0.25 mm × 0.25 μm),in which the 5975 MS was run in the full-scan mode(Agilent Technologies, Inc.) [19].To extract the residual alkanes, 50.0 μL of cetyl chloride was added to the serum bottles as a surrogate standard, after which the entire serum bottle was extracted with nhexane and subsequently dried over anhydrous Na2SO4.The pooled organic layer was then transferred to a new vial and concentrated under a stream of N2.The extracts were injected into the GC–MS with the following program settings: initial oven temperature of 120°C for 3 min,and then increased to 260°C at a rate of 8 °C·min-1, at which the temperature was maintained for 10 min.The GC peak area of the surrogate standard(cetyl chloride)and alkanes ranging from C15to C20were integrated to obtain the n-alkane-to-standard ratio.Finally, the quantity of each n-alkane was determined based on the peak area ratio.

    Long-chain fatty acids (LCFAs) and volatile fatty acids (VFAs)were detected by taking 5 mL of culture aliquot from the triplicate samples, adding ammonia water to raise the pH to more than 12,and drying in an oven at 110°C.This was followed by esterification at 90°C for 60 min by adding 0.5 mL 10%butanol/sulfate solution.Then, the esterified derivatives were extracted with 0.5 mL nhexane and 0.5 mL n-dodecane.After, the extracts were injected into the GC–MS for detection.For the LCFA analysis,the oven temperature was maintained at 120°C for 3 min and increased at a rate of 8 °C·min-1to 260°C,where it was maintained for 10 min.Conversely,for the VFA analysis,the oven temperature was maintained at 60 °C for 1 min and then increased to 130 °C at a rate of 15 °C·min-1.

    The theoretical methane yield was calculated using the Buswell equation[25],which represents the complete oxidation of alkanes to methane and CO2(Eq.(2)).The equation was modified as previously described [23] (Eq.(3)) by considering the intermediate or dead-end metabolites detected from the incomplete oxidation of alkanes.This calculation was further modified by subtracting the amount of carbon in the intermediate metabolites(including LCFAs and VFAs) from the total consumed alkanes.

    To measure the cations and anions in the culture,2 mL of medium from the triplicate cultures was extracted and centrifuged before the supernatant was further filtered.The anion concentrations of Cl-, NO3-, PO43-, CO32-, SO42-, and S2O32-were measured using an ion chromatography model ICS-1100 (Dionex, USA)equipped with an IonPac AS11-HC(Thermo Scientific,USA)column(temperature: 30 °C; eluent: KOH elution; detector: conductivity detector and ampere detector; suppressor: anion suppressor SHY-A-6;Shine Technology,China).Meanwhile,the cation concentrations of K+, Mg2+, Ca2+, methylamine, dimethylamine, and trimethylamine were measured using ion chromatography model IC25A (Dionex) equipped with an IonPac CS17 (Thermo Scientific)column (temperature: 35 °C; eluent: 20 mmol·L-1methylsulfonic acid isocratic elution; detector: conductivity detector;suppressor:cation suppressor, CERS-500, Thermo Scientific).For Na+and ammonia detection, an ion chromatography model CIC-D120(Shine Technology) equipped with an SH-CC-3(Shine Technology)column was used.

    2.6.DNA/RNA extraction and 16S rRNA high-throughput sequencing

    After more than 1400 days of anaerobic incubation, the total DNA and RNA were extracted from the pellets of one of the biological triplicates using a PowerSoil DNA Isolation Kit (MO BIO, USA)and a PowerMicrobiomeTMRNA Isolation Kit(MO BIO),respectively.To enrich the messenger RNA(mRNA),rRNA was depleted from the total RNA using a Ribo-Zero rRNA Removal Kit (Bacteria) (Illumina?, USA) according to the manufacturer’s instructions.The complementary DNA (cDNA) was synthesized as previously described [9].The bacterial and archaeal 16S rRNA genes were PCR-amplified from metagenomic DNA using primer sets 515F/907R [26] and 344F/915R [27], respectively.All 16S rRNA gene amplicon libraries were sequenced on an Illumina Hiseq 2500 platform, generating 250 base pairs (bp) paired-end reads.Quality filtering, paired-end read merging, adapter removal, and primer trimming were performed as described previously [11].The high-quality sequences obtained using the Illumina Miseq platform (Majorbio, China) were clustered into operational taxonomic units at a 97% similarity level based on the SILVA database(release 123; Germany) [28].

    2.7.Metagenomic and metatranscriptomic sequencing

    Metagenomic sequencing and metatranscriptomic libraries were generated using the NEB Next?UltraTMDNA Library Prep Kit for Illumina?(New England Biolabs, USA), in accordance with the manufacturer’s recommendations, to which the index codes were added.The library quality was assessed using a Qubit 3.0 Fluorometer (Life Technologies, USA) and Agilent Technologies, Inc.4200(Agilent Technologies, Inc.) system.Finally, the library was sequenced on an Illumina Hiseq X-ten platform,generating a total of 150 bp paired-end reads.

    2.8.Metagenomic assembly and binning

    Raw reads were quality-filtered, and adapters were removed using Trim_galore?? http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/.under the default settings.Filtered reads of the metagenome dataset were assembled using metaSPAdes v3.7.0 with‘‘k-mer=127”[29],and all contigs shorter than 500 bp were removed.Read mapping was performed using Bowtie 2 [30], after which the mapped reads were filtered by removing the reads with a mapping quality of less than 3 [31] in order to remove ambiguously mapped reads.The resulting sequence alignment/map (SAM) file was sorted and converted into a binary alignment/map (BAM) file using SAMtools v1.6[32].Scaffolds were binned into population genomes using three different tools based on different clustering algorithms and features: ①MaxBin 2 v2.2.3 [33], which is based on an expectation–maximization algorithm,uses tetranucleotides,differential coverage,and marker genes,and was used with marker gene sets of 107 and 40 for bacterial and archaeal lineages,respectively;②MetaBAT2 v0.32.4[34], which applies a k-medoid clustering on tetranucleotide frequencies and differential coverage, was used with its default settings; and ③CONCOCT v0.4.0 [35] uses Gaussian mixture models and tetranucleotide frequencies with differential coverages.The resulting bins were supplied to the DAS Tool (v1.1.1) for consensus binning and de-replication, generating an optimized set of metagenome-assembled genomes (MAGs) [36].The resulting MAGs were further refined by removing outlier scaffolds by setting the percentile of divergent guanine and cytosine,sequencing coverage,and tetranucleotide frequencies to 95 in RefineM v0.0.23[37].Subsequent manual checks were performed using Anvi’o v5 [38].Briefly, the MAGs were split into major clusters based on their split coverage and guanine and cytosine content.These clusters were merged if the separated clusters were parallel in the same genome tree position.

    2.9.Taxonomic and functional annotation

    Genes from the entire assembly were labeled by Prodigal v2.6[39] in ‘‘meta” mode.Then, the amino acid files were submitted to the GhostKOALA server[40]in the prokaryotic species database to assign the KEGG orthology (KO) numbers.Additionally, the carbohydrate-active genes were identified in the dbCAN web server[41]with cutoffs used in a previous study[42](coverage >0.40;Evalue <1 × 10-18; identity >30%).Genes involved in anaerobic hydrocarbon degradation were identified using Blastp via a comparison with a previously reported database [42] (coverage >0.40;E-value <1×10-20;identity >30%).Hydrogenases were first identified using the HMMER(v3)search function[43](E-value <1×10-20)with custom HMM models(Appendix A File S1),and hits were further classified into different hydrogenase groups on the HydDB web server [44].We also uploaded the nucleotide MAG sequences to the RAST server to confirm the gene annotations using subsystem technology [45].Taxonomic assignments of the unbinned genes were performed using Kaiju [46] against the national center for biotechnology information (NCBI) nr_euk database (2018-7-20) in‘‘Greedy” mode.Finally, the taxonomic classification of the MAGs was mainly based on GTDBtk (release 86) prediction?? https://github.com/Ecogenomics/GTDBTk.? https://inkscape.org/., which uses recently described relative evolutionary distance metrics to predict the divergence of newly binned clades[47].The classifications were further verified using a polygenetic analysis of the conserved protein sequences and 16S rRNA gene sequences.

    2.10.Genome coverage and gene FPKM calculation

    To estimate the genome sequencing coverage, we pooled the scaffold files of the final MAGs, the corresponding read mapping files generated from the metagenomic dataset, and the metatranscriptomic dataset into Anvi’o v5 [38] after filtering to remove reads with a mapping quality of less than 3, as previously described.The mean coverage for the inner quartiles (Q2 and Q3)for each MAG was estimated by summarizing the metagenome and metatranscriptome profiles.We used eXpress v1.5.1 [48] to calculate the fragments/kilobase of transcript per million fragments mapped (FPKM) of individual genes, in which the entire assembly was used as an input file.The FPKM values of the metagenome and metatranscriptome datasets (referred to as FPKM in this study) were individually calculated using sorted BAM files derived from the read mapping of the metagenome and metatranscriptome datasets, respectively.

    2.11.Phylogenetic analysis

    The small subunit rRNA genes in the MAGs were identified using checkM v1.0.16 [37] with the ‘‘ssu_finder” function, and sequences longer than 400 bp were used for the phylogenetic analysis.The selected sequences were aligned using MAFFT [49] with iterative refinement methods ‘‘G-INS-i,” before they were manually refined (retained columns with <10% gaps).To construct the genome tree, all reference genomes (listed in Appendix A Table S2) and MAGs were pooled into PhyloPhlAn v0.99, which extracted and aligned 400 conserved protein sequences [50].The concatenated alignment file was then extracted for phylogenomic tree building.For both alignments,maximum likelihood trees were reconstructed using IQ-tree v1.6.7,under standard model selection with 1000 ultrafast bootstrap conditions.The trees were further decorated using Inkscape v0.92.3?? https://github.com/Ecogenomics/GTDBTk.? https://inkscape.org/.and visualized using iTOL [51].

    2.12.Data availability

    The raw sequence data used in this study were submitted to the NCBI database under BioProject PRJNA792315.

    3.Results

    3.1.Ten-year successive enrichment

    In this study, we set up an enrichment culture with OS as the inoculum and an n-alkane mixture (C15–C20) as the sole carbon and energy source [19,22].After ten years of cultivation and three consecutive transfers,the culture was found to improve the kinetic properties of methane production.This was concluded by considering the period of the lag phase, the maximum methane concentration, and the maximum methane production rate (Fig.1 and Appendix A Fig.S1).In the third transfer, the culture produced a total of 352 μmol of methane 1413 days after a lag phase of 267.1 days (predicted by fitting the data into the modified Gompertz model [24,52]) (Fig.1 and Fig.S1).After incubation, almost all n-alkanes were depleted (Appendix A Table S3), and putative degrading intermediate metabolites,including octadecanoate,hexadecanoate, tetradecanoate, dodecanoic, butyrate, propionate,acetate, and formate, were detected in the culture (Appendix A Table S4).The detected methane in the third transfer accounted for 63.4%–71.5% of the theoretical methane yield.

    3.2.Microbial community succession

    The biomass in the third transfer was harvested during the lateexponential methane production phase for nucleotide sequencing(Fig.1).Based on the 16S rRNA gene survey results, the microbial communities varied greatly over the incubation period (Fig.2).Compared with the second transfer, the four major bacterial lineage (Enterobacteriaceae, Unclassified_Cloacimonadales, Ignavibacteriaceae, and Leptospiraceae) concentrations decreased by more than 50% in the third transfer (Fig.2).Subsequently,Hydrogenophilaceae, Anaerolineaceae, and Desulfosarcinaceae,which accounted for more than 80% of the bacterial community composition,became the dominant bacteria(Fig.2).Lineages such as Syntrophaceae,Synergistaceae,Thermodesulfovibrionaceae,and Kosmotogaceae, which were found to be occasionally detected in OS samples, the initial culture, or the first transfer,existed at concentrations below the detection limit in the second transfer sample, but were detected in low abundance in the third transfer(Fig.2).This variation might be because of the extremely low abundance of these minor members in the samples, which could be overlooked when applying low-throughput amplicon sequencing from clone libraries [19].

    Fig.1.Methane production values of the n-alkane-degrading enrichment cultures during three sequential transfers.Methane data from this study are marked with dotted borders.The curve with black diamonds represents the initial culture and subsequent enrichment cultures amended with the n-alkanes mixture(C15–C20)as the sole carbon and energy source(three replications),and the curve with red squares represents the n-alkanes-free control group(three replications).Error bars indicate the standard error of the mean of the biological triplicates.Each time point for transfer to the fresh medium is marked with a green arrow, and the time point for metagenomic and metatranscriptomic sequencing is marked with a yellow arrow.

    The archaeal communities in all the samples from each stage were mainly composed of methanogenic lineages, including Methanomicrobiaceae, Methanosaetaceae, Methanosarcinaceae,Methanoregulaceae, and Methanocorpusculaceae.Compared with the second transfer, however, the dominant archaea members in the third transfer shifted from Methanosaetaceae (95.9% in the second transfer) to Methanosarcinaceae (51.2% in the third transfer).Methanomethyliaceae(22.4%),which is in the candidate phylum Verstraetearchaeota was detected for the first time in every stage of the third transfer.Note that the different PCR primers and sequencing techniques used for 16S rRNA gene amplicon sequencing could potentially bias the microbial composition[19,22].

    The assembly and binning of metagenome data resulted in a total of 23 MAGs with >50% completeness and <10% contamination (estimated by CheckM [37]).Reconstructed MAGs comprised taxonomically diverse members from a total of 13 bacterial and 2 archaeal phyla (Fig.3 [50] and Appendix A Table S5).Within the bacteria domain, members of the phylum Chloroflexi, especially the class Anaerolineae, were highly represented (6 of 23 MAGs with 7.8× to 198.5× mean sequencing coverage, as compared with the median coverage of 58.2× all MAGs, Appendix A Fig.S2).Within the archaea domain, two MAGs representing Methanosarcinaceae exhibited medium to high sequencing coverage(54.9×and 102.9×for Methanosarcinaceae archaeon bin5 and Methanosarcina sp.bin9, respectively; Fig.S2).We also found one Methanomethyliaceae-like MAG with a slightly lower sequencing coverage (46.7×) that represented the candidate phylum, Verstraetearchaeota.Overall,the pattern of MAG sequencing coverage was mostly in accordance with the microbial composition based on the 16S rRNA gene amplicon data.In total,18 MAGs were classified into three classes, four orders, three families, and six genus level uncultured lineages that have not been described before (Fig.3 and Table S5).

    3.3.Anaerobic alkane degradation

    To explore the metabolic potential of the third transfer,we analyzed the MAGs and the unbinned assembly fragments that could not be classified into any MAGs for specific functions (see Section 2).We used the FPKM values to demonstrate the relative abundances of the genes in the metagenome data, in which the FPKM values that ranked in the top 10% of all genes (i.e.,FPKM >13.0) were considered to be relatively high.

    In this defined medium, the organic carbon source is derived from either inoculated n-alkanes (C15–C20) or detrital biomass[42].To identify the potential for the microbial degradation of hydrocarbons, we searched for functional marker genes encoding enzymes that initiate anaerobic hydrocarbon biodegradation via the activation of C–H bonds [53].The metagenomic results indicated a possible oxygen-independent C–H activation pathway,wherein there is n-alkane addition to fumarate (known as the fumarate addition pathway) [54] by glycyl-radical enzymes.Putative assA, which may mediate n-alkane addition to fumarate,were found in four MAGs (Desulfosarcinaceae bacterium bin16, Chloroflexi bacterium bin26,Phycisphaerae bacterium bin27,and Spirochaetes bacterium bin3), as well as several unbinned fragments representing Chloroflexi, Deltaproteobacteria, Alphaproteobacteria, Synergistetes, Firmicutes, Spirochaetes, Actinobacteria, and Candidatus (Ca.) Aminicenantes (Appendix A Fig.S3), with a high relative abundance(most FPKM values >13.0;Fig.4 and Appendix A Table S6).The assA-like sequences identified were phylogenetically distant from the canonical Smithella-like ASS[55](amino acid identities <45%), but formed a common clade with the putative ASS, which has been proposed to mediate alkane activation in anaerobic alkane degraders Vallitalea guaymasensis L81 [56], Thermococcus sibiricus MM 739 [57], and Archaeoglobus fulgidus VC-16[18,58](Fig.S3).Meanwhile,no Smithella-like MAGs or operational taxonomic units (OTUs) were generated from the metagenome or 16S rRNA gene clone library.The lack of canonical assA was further verified via PCR amplification using previously designed canonical assA/adenosine triphosphate-binding cassette (bssA)-targeting primers(Table S6),as no amplicon sequences were retrieved from the samples.

    Fig.2.Relative abundance of major microbial lineages (>5% in at least one sample) at family level represented by 16S rRNA gene amplicon sequences detected in OS and enrichment cultures.16S rRNA gene sequences from clone libraries of OS, initial culture (Initial), and two transfers were retrieved from previous studies [19,22], and the relative abundance based on the high-throughput 16S-tag sequencing of the third transfer is shaded in gray.The total sequences in the bacterial and archaeal libraries are also shown.Ca.: Candidatus; No.: numbers.

    Following the initial fumarate addition of alkanes, alkylsuccinates are converted to fatty acids through several steps, including acetyl-coenzyme A (CoA) activation, carbon skeleton rearrangement,and decarboxylation[59].The genes involved in these steps,namely CoA-synthetase/ligase gene (assK), putative methylmalonyl-CoA mutase gene (mcmLS), and putative propionyl-CoA carboxylase gene (pccAB), were detected in several MAGs (Appendix A Table S7).However, no MAGs containing the putative assA gene contained the complete gene set required for alkylsuccinate metabolism, but some possessed a high degree of completeness (for example, the completeness of the Desulfosarcinaceae bacterium bin16 was 96.3%).Nevertheless, genes related to the pathway downstream of the n-alkane fumarate addition were found in several MAGs that did not have the assA(Table S7), suggesting an interplay between these MAGs and the assA-containing MAGs when degrading n-alkanes.

    3.4.Necromass recycling and fermentation

    In addition, necromass, which is mainly composed of complex carbohydrates, proteins, and lipids, could serve as a substrate for many community members in the third transfer[42].On this basis,we searched for genes associated with peptide/amino acid degradation,complex carbohydrate degradation,and fatty acid degradation.Most of the MAGs(19 out of 23;Fig.5 and Table S7)contained secreted carbohydrate-active enzymes and a near-complete Embden–Meyerhof pathway (glycolysis), suggesting the ability to break down complex carbohydrates from cell debris.Meanwhile,several genes encoding adenosine triphosphate-binding cassette(ABC)-type peptide transporters and amino acid transporters were observed in 11 MAGs (Fig.5 and Table S7).Genes encoding peptidases, aminotransferases, 2-oxoacid-ferredoxin oxidoreductases,and aldehyde-ferredoxin oxidoreductases were also detected in these MAGs.Therefore, we inferred that these MAGs are capable of importing and degrading peptides step-wise into amino acids,then to 2-keto acids, and finally to acetyl-CoA [60].Further, the genes required for fatty acid degradation(β-oxidation)were abundant in the metagenome data (Fig.4), and a complete β-oxidation pathway was observed in seven bacterial MAGs (Fig.5).

    As no exogenous electron acceptors (e.g., oxygen, oxidized metal cations, oxidized sulfur compounds, or oxidized nitrogen compounds) were supplied to the enrichment culture, the heterotrophic microorganisms were usually disposed of catabolismderived reducing equivalents by either producing H2or fermentation.Genes encoding H2evolving-only hydrogenases([FeFe]group A/B and [NiFe] group 4a/b/d/g) were found in 17 MAGs, whereas hydrogenases for H2utilization were found in 13 MAGs (Fig.5).Various MAGs and assembled fragments contained genes that indicated the capability for the fermentative production of lactate (6 MAGs) and ethanol (12 MAGs) (Figs.4 and 5).Moreover, acetogenic CO2reduction via the Wood–Ljungdahl pathway could be an alternative pathway for the disposal of electrons derived from hydrocarbon/necromass oxidation [61,62].Nearly all the bacterial MAGs contained orthologous genes that are responsible for acetate formation.However, complete gene sets for the Wood–Ljungdahl pathway were only found in four bacterial MAGs(Fig.5),of which the Wood–Ljungdahl pathway found in the Methanosarcinaceae MAGs was most likely associated with methanogenesis or biosynthesis [63].

    Fig.3.Phylogenetic placement of reconstructed metagenome-assembled genomes.(a)Phylogenetic affiliation of the metagenome-assembled genomes(MAGs)and reference microorganisms based on their 16S rRNA genes,shaded in blue.(b)Phylogenomic affiliation of 23 MAGs based on 400 conserved protein sequences[50].Both alignments are based on MAFFT and then filtered with columns containing >90%gaps.The trees were built using IQ-Tree under the standard model selection with 1000 bootstrap replicates.Bootstrap values >80% are indicated.The MAGs assembled in this study are marked in red for bacteria and green for archaea.Scale bars indicate the average number of substitutions/site.

    3.5.Methanogenesis

    To obtain the maximum free energy available from organic matter in the absence of external electron acceptors, microbial communities will continue their carbon turnover reactions until they reach the most thermodynamically stable states, that is, methane and CO2[64].Genes specific for methanogenesis were only found in members of Methanosarcinales (Methanosarcinaceae archaeon bin5 and Methanosarcina sp.bin9) and Methanomethyliales(Methanomethyliaceae archaeon bin15) (Fig.4).Methylcoenzyme M reductase genes in these MAGs were clustered with MCR enzymes that catalyze methane metabolism, rather than MCR-like complexes associated with anaerobic methane/alkane oxidation (Appendix A Fig.S4) [65,66].Therefore, a putative role for methane production was proposed for these organisms.As all these methanogenic signature genes had a high abundance(FPKM >13.0)(Fig.4),they had a high potential for methanogenesis in the sample.Furthermore, the complete gene sets for hydrogenotrophic, acetoclastic, and methylotrophic methanogenesis from methanol and tri/di/methylamine (T/D/MMA) were found in the MAGs of Methanosarcina sp.bin9 and Methanosarcinaceae archaeon bin5, indicating a versatile metabolic potential (Fig.5).Conversely, the genome analysis revealed a strict H2-dependent methylotrophic methanogenic lifestyle for Methanomethyliaceae archaeon bin15 (Fig.5), which is consistent with the previous description of this lineage [67].

    With respect to the ubiquitous distribution of methylotrophic methanogenic capability in all three archaeal MAGs, we searched for genes responsible for methyl compound generation.Marker genes for converting choline, trimethylamine N-oxide, betaine,and L-carnitine (cut, tor, grd, and cnt) to trimethylamine (TMA),and genes for glutamate conversion(mgnm)and methanol generation (mdh1) were observed in the metagenome data at varying abundances (Fig.4), implying that these substrates are readily available.

    Fig.4.Reconstruction of the potential carbon cycle in the third transfer.The FPKM values of the genes were summed according to their taxonomic assignment.Gene names:fwd,formylmethanofuran dehydrogenase gene;ftr,formylmethanofuran–tetrahydromethanopterin N-formyltransferase gene;mer,5,10-methylenetetrahydromethanopterin reductase gene; mtd, methylenetetrahydromethanopterin dehydrogenase gene; mch, methenyltetrahydromethanopterin cyclohydrolase gene; mcr, methyl-coenzyme M reductase alpha subunit gene; mtr, N5-methyltetrahydromethanopterin:coenzyme M methyltransferase; acs, acetyl-coenzyme A (CoA) synthetase (EC: 6.2.1.1) gene; acd,acetate-CoA ligase (adenosine diphosphate (ADP)-forming) gene; cdh, acetyl-CoA decarbonylase/synthase; mta, methanol-5-hydroxybenzimidazolylcobamide Comethyltransferase gene; mtm, methylamine-corrinoid protein co-methyltransferase; mtb, dimethylamine-corrinoid protein Co-methyltransferase; mtt, trimethylaminecorrinoid protein co-methyltransferas gene; acox/fadE/acd/bcd, acyl–acyl carrier protein (ACP) dehydrogenase gene; paaF/crt/fadJ/fadB/ech, enoyl-CoA hydratase gene; fadJ/fadB/fadN/ech/paaH/EHHADH,3-hydroxyacyl-CoA dehydrogenase gene;fadA,acetyl-CoA acyltransferase gene.TMA,trimethylamine;cut,choline trimethylamine-lyase gene;tor, trimethylamine-N-oxide reductase (cytochrome c) gene; grd, betaine reductase gene; cnt, carnitine monooxygenase subunit gene; MGNM, methylamine–glutamate Nmethyltransferase gene; mdh, methanol dehydrogenase gene; ALDH, aldehyde dehydrogenase gene; CoA, coenzyme A; por, pyruvate ferredoxin oxidoreductase gene; metF,methylenetetrahydrofolate reductase gene;mut,methylmalonyl-CoA mutase gene;acyP,acylphosphatase gene;ADA,acetaldehyde dehydrogenase gene;kor,2-oxoglutarate/2-oxoacid ferredoxin oxidoreductase gene; ldh,L-lactate dehydrogenase gene;fol,methenyltetrahydrofolate cyclohydrolase gene; fhs,formate-tetrahydrofolate ligase gene;cooF, anaerobic carbon-monoxide dehydrogenase iron sulfur subunit gene.For clarity, the Embden–Meyerhof pathway genes were omitted, and more details are listed in Table S6.

    3.6.Metatranscriptomic sequencing coverages

    Based on the metatranscriptomic sequencing coverage, several MAGs (10 out of 23) exhibited zero transcriptional activity at the time of sampling (Fig.S2).Meanwhile, all the MAGs with metagenomic sequencing coverage above the median had transcripts,with the exceptions of Methanosarcina sp.bin9 and Syntrophorhabdaceae bacterium bin2.Specifically, Tepidiphilus sp.bin11 had the highest mean coverage (43.2×), followed by Mesotoga sp.bin14(8.8×), Methanomethyliaceae archaeon bin15 (8.4×), Anaerolineaceae bacterium bin23 (7.2×), Anaerolineaceae bacterium bin21(5.5×),Desulfosarcinaceae bacterium bin16(4.8×),Ca.Hydrogenedentota bacterium bin7 (4.1×), Bellilinea sp.bin22 (3.4×), Enterobacter bacterium bin6 (2.8×), Thermodesulfovibrio sp.bin10(2.6×), Coriobacteriia bacterium bin12 (1.2×), Deferribacterales bacterium bin8 (0.6×), and Ca.Sumerlaeota bacterium bin4(0.6×) (Fig.S2).

    3.7.Gene transcriptomic pattern

    To analyze the specific functional activity in the metatranscriptome of the third transfer, the gene-specific transcriptional levels were calculated using the FPKM values obtained from the metatranscriptome data.As less than 5% of genes had transcripts, we defined the FPKM values higher than the median value of all transcribed genes (i.e., FPKM >10.2) as relatively high.

    Despite having been detected in several MAGs and unbinned fragments, only the assA-like genes in the Desulfosarcinaceae bacterium bin16(FPKM=38.0)and several fragments classified within the phyla Synergistetes and Firmicutes (FPKM = 3.3 and 8.1,respectively) were transcribed at the sampling time (Fig.5 and Table S8 in Appendix A), which suggests that the members of Deltaproteobacteria (Desulfosarcinaceae bacterium bin16), Synergistetes, and Firmicutes are responsible for n-alkane activation in the culture.In addition to assA,the assK gene in Desulfosarcinaceae bacterium bin16 also had high-level transcripts(FPKM=28.1).Further,the genes associated with pathways downstream of n-alkane fumarate addition in the Tepidiphilus sp.bin11 were all highly transcribed.

    Fig.5.Heatmap showing the transcription level of functional genes or pathways present in each MAG.MAGs with specific functional genes or pathways present are marked with color,and the transcriptional levels are expressed as the log-transformed FPKM(fragments/kilobase of transcript per million fragments mapped)values of the functional genes or the average value of all the detected genes in the pathways.Pathways are indicated as being present if ≥75% of the genes exist in the genome.TMAO:trimethylamine N-oxide.

    Regarding necromass recycling,the genes for peptide transport and degradation were transcribed in Tepidiphilus sp.bin11,Enterobacter bacterium bin6,Anaerolineaceae bacterium bin21 and bin23, Bellilinea sp.bin22, and Mesotoga sp.bin14 (Fig.5).Meanwhile, the genes for intracellular and extracellular carbohydrateactive enzymes,as well as those involved in the Embden–Meyerhof pathway (glycolysis), were transcribed in Desulfosarcinaceae bacterium bin16, Enterobacter bacterium bin6, Anaerolineaceae bacterium bin23, Ca.Hydrogenedentota bacterium bin7, Ca.Sumerlaeota bacterium bin4, Thermodesulfovibrio sp.bin10, and Mesotoga sp.bin14(Fig.5).The transcriptional activity of fatty acid degradation was found in four MAGs: Desulfosarcinaceae bacterium bin16, Tepidiphilus sp.bin11, Ca.Hydrogenedentota bacterium bin7, and Enterobacter bacterium bin6 (Fig.5).

    Moreover, the genes for methanogenesis in Methanosarcina sp.bin9 and Methanosarcinaceae archaeon bin5 were both below the detection limit(Fig.S2),indicating a dormant state.Conversely,the genes related to methanogenesis, as well as the trimethylamine methyltransferase gene in Methanomethyliaceae archaeon bin15,were found to be highly transcribed(e.g.,the FPKM values of mcrA and mttB were 20.3 and 13.7, respectively) with the transcripts of the H2-utilizing gene.This suggests the H2-dependent methylotrophic methanogenic activity of this organism.Regarding TMA production, the cut and grd genes encoded in unbinned fragments exhibited transcriptional activity (FPKM values of 3.2 and 14.9, respectively), whereas the genes generating other methyl compounds (methyl-/dimethyl-amine (M/DMA) and methanol)were inactive in the transcriptome.

    4.Discussion

    The succession of microbial communities during long-term incubation is considered to be the result of adapting to the cultivation environment.Thus, the change in microbial communities likely contributes to sustained improvements in the kinetic properties of methane production in the enrichment cultures.The dominant lineages in this culture were primary lineages that have been consistently detected in hydrocarbon-polluted environments,such as petroleum-contaminated soils [68], hydrocarbon-contaminated aquifers [1], and oil fields [69], thereby supporting the theory of the in situ microbial conversion of hydrocarbons to methane as a bioremediation strategy in hydrocarbon-polluted environments.However, exogenous nutrient supply (such as phosphate) may be required to establish an active methanogenic alkane-degrading microbial community, as hydrocarbon-polluted environments are usually limited to specific nutrients [70].

    In this study, the metagenomic analysis revealed that the third transfer contained a wide range of uncultured lineages not observed in other habitats, which warrants an investigation into their metabolic potential.As revealed by the integrated metagenomic and metatranscriptomic analyses,the major organic carbon and energy sources in the cultures, n-alkanes, were utilized by anaerobic microorganisms via the fumarate addition pathway,wherein they were mediated by members of Deltaproteobacteria,Firmicutes, and Synergistetes at the time of sampling.Furthermore,the pathway reconstruction from 23 MAGs recovered from the metagenome indicated incomplete fumarate addition in Desulfosarcinaceae bacterium bin16.The metatranscriptomic analysis suggested that n-alkane degradation via fumarate addition could be a joint achievement of Desulfosarcinaceae bacterium bin16 and Tepidiphilus sp.bin11.This expands the current understanding of Smithella, which is considered to be the major alkane breaker in many other low- to high-temperature alkanedegrading enrichment cultures [9,11,12,71,72].In addition to nalkanes,the microorganisms in the third transfer actively recycled cell debris, peptides/amino acids, carbohydrates, and fatty acids.We also propose that fermentation (ethanol) and acetogenic and syntrophic interactions (H2) were major electron sinks for heterotrophic microorganisms,thereby driving microbial interactions.To date, microorganisms affiliated with Anaerolineaceae have been consistently detected in methanogenic cultures [10,71–73], and in a previous study, they were found to be the dominant bacteria in the second transfer [19].However, the functional roles of these organisms in the n-alkane methanogenic are largely unknown due to the lack of genomic representatives [19,74].In this study, we demonstrated that members of Anaerolineaceae in the third transfer were mainly responsible for the recycling of necromass.Consistent with this finding, a previous study examined an independent enrichment culture amended with13C-labeled hexadecane that used the same inoculum employed in this study and revealed that Anaerolinea and Tepidiphilus phylotypes have significantly higher relative abundances in the heavy fractions of DNA, as compared with the light fractions [75].This result indicates the active involvement of Tepidiphilus in the transformation and biodegradation of alkanes and suggests the cross-labeling of Anaerolinea by incorporating13C from necromass [76] when considering the long-term incubation (920 days) [75].

    In previous studies, hydrogenotrophic and acetoclastic methanogenesis have been considered to be the major processes downstream methanogenic hydrocarbon degradation [5,6,20,77].However,in this study,metatranscriptomic analysis demonstrated the exclusive transcriptional activity of methanogenesis from methyl compounds.It has been proven that methanogenesis from methyl compounds is common where these compounds are abundant, such as marine sediments [78,79].In this study, the methyl compounds derived from the necromass were readily available substrates for methanogenesis.Note that the Methanosarcinales MAGs constructed in this study have the metabolic potential to utilize hydrogen and methyl compounds for methane production(Fig.4),but they were underrepresented in the metatranscriptome.This might be partially explained by the fact that methanogens without cytochromes(such as Methanomethyliaceae)usually outcompete methanogens with cytochromes (species within Methanosarcinales)for H2because they have a lower H2threshold concentration [63].However, the exact relationship between the members of Methanomethyliaceae and Methanosarcinales will only be resolved with additional research, in which a pure culture is studied to isolate a representative of Methanomethyliaceae.

    Overall,we do not intend to rule out the significant contribution of the hydrogenotrophic and acetoclastic methanogenic pathways,as they have been repeatedly observed to be involved in methanogenic hydrocarbon degradation in both enrichment cultures and environmental samples.However, this study expands the current knowledge that, at least in the late-exponential phase, wherein necromass-derived methyl-compounds are available, H2-dependent methylotrophic methanogenesis and methylotrophic methanogenesis may also play major roles in the conversion of alkanes to methane.H2-dependent methylotrophic and methylotrophic methanogenesis may have been overlooked in previous enrichment cultures because they were harvested earlier than the late-exponential phase [19,20,23].Another possible reason for this difference is that the initial inoculums for previous cultures(mostly oil reservoir production waters and ditch sediments) are different from that used in this study (OS from oil refinery)[5,23].Based on the results of this study, as well as the fact that obligate methylotrophic methanogens (such as Methanohalophilus and Methanolobulus) and facultative methylotrophic methanogens(such as Methanosarcina)are ubiquitous in hydrocarbon-associated environments [21,69], these organisms may be assigned the ecological role of methylotrophic methanogenesis from necromassderived methyl-compounds.

    Herein, we found that the transcripts of several MAGs were below the detection limit, which was unexpected, as some possessed a high relative abundance in the metagenome.Nevertheless,the extremely low transcriptional levels were consistent with the previous descriptions of metatranscriptomes from biodegraded oil fields [18].A possible reason for this may be the extremely small growth rate of methanogenic alkane-degrading cultures,which implies that very little gene expression occurred.The metagenomic analysis demonstrated that the transcriptionally underrepresented members have a metabolic potential associated with methanogenic alkane degradation and may actively participate in this process at a time before sampling.Note that the results of this metatranscriptomic analysis are only a snapshot of the microbial activity occurring in the late-exponential phase of methane production.In a future study, increased sampling times are required to capture the overall microbial activity profile at different phases, including the lag, early exponential methane production, mid-exponential methane production, late-exponential,and stationary phases.In addition, H2-dependent methylotrophic methanogenesis during anaerobic alkane degradation requires further verification, probably by another transfer cultivation combined with DNA-SIP labeling technology.

    5.Conclusions

    This study evaluated the mechanism of microbial methanogenic alkane degradation, revealing that H2-dependent methylotrophic methanogens, as well as methylotrophic methanogens, may play an important ecological role in the carbon cycle in hydrocarbonassociated environments.Additionally, we achieved sustained improvements in the kinetic properties of methane production in the alkane-degrading consortium through several enrichment cultivations, thereby demonstrating potential remediation of hydrocarbon-contaminated anoxic environments.Overall, this study will help to inform technologies that the microbially enhanced energy recovery of petroleum from marginally productive oil fields may convert approximately 30% of the total crude oil remaining underground as residues into methane gas, which can then be easily recovered.

    Acknowledgments

    This work was supported by grants of National Natural Science Foundation of China (42061134011, 52074129, and 42173076);the Natural Science Foundation of Shanghai (21ZR1417400); the Shanghai International Collaboration Program (18230743300);the Fundamental Research Funds for the Central Universities(JKJ012016028) to Yi-Fan Liu, Li-Bin Shou, Shi-Zhong Yang,Jin-Feng Liu, and Bo-Zhong Mu; and the NSFC/RGC Joint Research Fund (41161160560) to Ji-Dong Gu.

    Authors’ contribution

    Zhao-Wei Hou and Xiao-Lin Wu provided the samples.Li-Ying Wang, Bo Liang and Jing Chen set up the enrichment cultures.Jing Chen monitors the methane production during incubation.Zhong-Lin Liu performed the PCR experiments.Dan-Dan Lin and Li-Bin Shou were committed to all the experiments.Yi-Fan Liu performed the statistic and metagenomic and metatranscriptomic analysis.Yi-Fan Liu,Ji-Dong Gu,Lei Zhou,and Bo-Zhong Mu wrote the paper.Bo-Zhong Mu oversaw this project.All authors read and approved the final manuscript.

    Compliance with ethics guidelines

    Yi-Fan Liu, Jing Chen, Zhong-Lin Liu, Zhao-Wei Hou, Bo Liang,Li-Ying Wang,Lei Zhou,Li-Bin Shou,Dan-Dan Lin,Shi-Zhong Yang,Jin-Feng Liu, Xiao-Lin Wu, Ji-Dong Gu, and Bo-Zhong Mu declare that they have no conflict of interest or financial conflicts to disclose.

    Appendix A.Supplementary data

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

    午夜免费鲁丝| 精品国产乱子伦一区二区三区| 可以免费在线观看a视频的电影网站| 高清毛片免费观看视频网站 | 露出奶头的视频| 色精品久久人妻99蜜桃| 中出人妻视频一区二区| 久久久久久久久中文| 久久国产精品影院| 久久国产精品男人的天堂亚洲| 黄色视频,在线免费观看| 桃红色精品国产亚洲av| 亚洲欧洲精品一区二区精品久久久| 美国免费a级毛片| 亚洲 国产 在线| 久久精品国产清高在天天线| 精品福利永久在线观看| 精品无人区乱码1区二区| 很黄的视频免费| 亚洲男人的天堂狠狠| 中文字幕最新亚洲高清| 久久国产精品影院| 国产av一区二区精品久久| 天堂中文最新版在线下载| 黄片大片在线免费观看| 熟女少妇亚洲综合色aaa.| 韩国av一区二区三区四区| 女人被躁到高潮嗷嗷叫费观| 国产蜜桃级精品一区二区三区| 香蕉丝袜av| 国产成人av教育| 久久青草综合色| 水蜜桃什么品种好| 国产高清激情床上av| 欧美日韩乱码在线| 性欧美人与动物交配| 亚洲av成人一区二区三| 大型黄色视频在线免费观看| 午夜福利免费观看在线| 国产色视频综合| 国产片内射在线| 日韩精品中文字幕看吧| 亚洲国产精品999在线| 91成年电影在线观看| 99在线视频只有这里精品首页| 久久热在线av| 欧美日韩亚洲高清精品| 啦啦啦 在线观看视频| 午夜免费激情av| 欧美日韩亚洲国产一区二区在线观看| 亚洲欧美精品综合一区二区三区| 在线免费观看的www视频| 最近最新中文字幕大全免费视频| 国产成人精品久久二区二区91| 无人区码免费观看不卡| 在线观看舔阴道视频| 成人国产一区最新在线观看| 91av网站免费观看| 成熟少妇高潮喷水视频| 亚洲欧美精品综合一区二区三区| 欧美激情 高清一区二区三区| 最近最新免费中文字幕在线| 国产精品一区二区免费欧美| 欧美一区二区精品小视频在线| 免费高清视频大片| 后天国语完整版免费观看| 精品国产乱子伦一区二区三区| 国产成人欧美| 成年人免费黄色播放视频| 国产av又大| 最新美女视频免费是黄的| 亚洲专区国产一区二区| 亚洲久久久国产精品| 成年人黄色毛片网站| 黑人欧美特级aaaaaa片| 悠悠久久av| 欧美日本中文国产一区发布| 国产欧美日韩综合在线一区二区| av电影中文网址| 亚洲成人免费av在线播放| 国内毛片毛片毛片毛片毛片| 狂野欧美激情性xxxx| 超碰成人久久| 他把我摸到了高潮在线观看| 男人的好看免费观看在线视频 | 亚洲精品粉嫩美女一区| x7x7x7水蜜桃| 黑人操中国人逼视频| 久久久久亚洲av毛片大全| 露出奶头的视频| 国产区一区二久久| 两性午夜刺激爽爽歪歪视频在线观看 | av在线播放免费不卡| 国产成人精品无人区| 久久久国产一区二区| 成年版毛片免费区| 岛国视频午夜一区免费看| 一区在线观看完整版| 欧美最黄视频在线播放免费 | 一区二区三区精品91| 久久热在线av| 新久久久久国产一级毛片| 久久中文字幕一级| 校园春色视频在线观看| 国产精品98久久久久久宅男小说| 岛国视频午夜一区免费看| 中文字幕人妻熟女乱码| 国产成+人综合+亚洲专区| 亚洲国产精品sss在线观看 | 免费在线观看视频国产中文字幕亚洲| 另类亚洲欧美激情| 亚洲人成电影免费在线| 国产免费男女视频| 亚洲视频免费观看视频| 国产色视频综合| 免费一级毛片在线播放高清视频 | 亚洲性夜色夜夜综合| 如日韩欧美国产精品一区二区三区| svipshipincom国产片| 一级黄色大片毛片| 国产三级黄色录像| 国产成人欧美在线观看| 黄色成人免费大全| 午夜福利,免费看| 99香蕉大伊视频| 精品高清国产在线一区| 国产不卡一卡二| 午夜福利在线观看吧| 欧美久久黑人一区二区| 99国产极品粉嫩在线观看| 极品教师在线免费播放| 亚洲欧美精品综合久久99| 老熟妇乱子伦视频在线观看| 午夜福利一区二区在线看| 大香蕉久久成人网| 黑人操中国人逼视频| 精品一区二区三卡| 午夜免费激情av| 免费一级毛片在线播放高清视频 | 亚洲av第一区精品v没综合| 中出人妻视频一区二区| 亚洲熟女毛片儿| 亚洲熟妇熟女久久| 大码成人一级视频| 超碰97精品在线观看| 亚洲av电影在线进入| www国产在线视频色| 侵犯人妻中文字幕一二三四区| 88av欧美| 精品一区二区三区四区五区乱码| 757午夜福利合集在线观看| 另类亚洲欧美激情| 操美女的视频在线观看| 亚洲一区二区三区不卡视频| 欧美不卡视频在线免费观看 | 成人永久免费在线观看视频| 在线观看66精品国产| 亚洲成人精品中文字幕电影 | 别揉我奶头~嗯~啊~动态视频| 成人18禁在线播放| 午夜免费鲁丝| 亚洲专区国产一区二区| 国产精品九九99| 最新在线观看一区二区三区| 日本黄色视频三级网站网址| 免费高清在线观看日韩| 成熟少妇高潮喷水视频| 国产欧美日韩一区二区三区在线| 亚洲欧美一区二区三区黑人| 日韩大尺度精品在线看网址 | 国产精品av久久久久免费| 一级片'在线观看视频| 亚洲熟女毛片儿| 久久久久久大精品| 国产激情久久老熟女| 国产精品免费视频内射| 两个人免费观看高清视频| 琪琪午夜伦伦电影理论片6080| 日韩人妻精品一区2区三区| 久久香蕉激情| av片东京热男人的天堂| 久久国产亚洲av麻豆专区| 久热爱精品视频在线9| 久久人人精品亚洲av| 中文字幕人妻丝袜一区二区| 亚洲av电影在线进入| 久久99一区二区三区| 成年人黄色毛片网站| 伦理电影免费视频| 涩涩av久久男人的天堂| 中亚洲国语对白在线视频| xxxhd国产人妻xxx| 精品少妇一区二区三区视频日本电影| 一级a爱片免费观看的视频| 99热只有精品国产| 一夜夜www| 亚洲人成电影免费在线| 九色亚洲精品在线播放| 黑丝袜美女国产一区| 国产高清videossex| 老熟妇仑乱视频hdxx| 高清黄色对白视频在线免费看| 国产精品一区二区免费欧美| 在线观看一区二区三区| 免费不卡黄色视频| 亚洲欧美日韩另类电影网站| 久久精品国产亚洲av高清一级| 精品第一国产精品| 涩涩av久久男人的天堂| 成年人黄色毛片网站| 欧美 亚洲 国产 日韩一| 国产精品免费一区二区三区在线| 婷婷丁香在线五月| 亚洲 国产 在线| 一级,二级,三级黄色视频| 久久 成人 亚洲| 在线观看一区二区三区| 国产区一区二久久| 久久久久久免费高清国产稀缺| 丝袜美腿诱惑在线| 国产成人av激情在线播放| 国产av在哪里看| 一级片免费观看大全| 成人精品一区二区免费| 他把我摸到了高潮在线观看| 亚洲欧美激情综合另类| 亚洲精品中文字幕在线视频| 亚洲国产欧美一区二区综合| 搡老熟女国产l中国老女人| 日韩欧美三级三区| 国产亚洲欧美精品永久| 国产av在哪里看| 这个男人来自地球电影免费观看| 日韩国内少妇激情av| 欧美老熟妇乱子伦牲交| 成人三级做爰电影| 日本五十路高清| 91国产中文字幕| 国产免费男女视频| 免费高清在线观看日韩| 91老司机精品| 午夜福利一区二区在线看| aaaaa片日本免费| 亚洲熟女毛片儿| а√天堂www在线а√下载| 欧美一级毛片孕妇| 亚洲少妇的诱惑av| 无人区码免费观看不卡| 悠悠久久av| 最近最新免费中文字幕在线| 免费日韩欧美在线观看| 久久人妻福利社区极品人妻图片| 啦啦啦在线免费观看视频4| 国产av在哪里看| 欧美日韩视频精品一区| 亚洲五月色婷婷综合| 极品教师在线免费播放| 男男h啪啪无遮挡| 少妇被粗大的猛进出69影院| 婷婷丁香在线五月| 色在线成人网| 日本wwww免费看| 欧美大码av| 波多野结衣av一区二区av| 欧美日韩视频精品一区| 777久久人妻少妇嫩草av网站| 国产免费现黄频在线看| 免费观看精品视频网站| 啦啦啦免费观看视频1| 女警被强在线播放| 欧美日韩黄片免| 大型av网站在线播放| 久久国产亚洲av麻豆专区| 宅男免费午夜| 久久草成人影院| 国产精品综合久久久久久久免费 | 极品教师在线免费播放| 看黄色毛片网站| 亚洲第一青青草原| 亚洲avbb在线观看| 成人三级黄色视频| 色婷婷av一区二区三区视频| 99久久人妻综合| 中文欧美无线码| 精品福利观看| 欧美成人性av电影在线观看| 精品卡一卡二卡四卡免费| 亚洲av成人av| 久久人人爽av亚洲精品天堂| 他把我摸到了高潮在线观看| 法律面前人人平等表现在哪些方面| 久久人人精品亚洲av| 18禁黄网站禁片午夜丰满| 动漫黄色视频在线观看| 午夜免费激情av| 精品福利观看| 精品免费久久久久久久清纯| 午夜日韩欧美国产| 五月开心婷婷网| 高潮久久久久久久久久久不卡| 国产精品久久久人人做人人爽| 国产97色在线日韩免费| 婷婷精品国产亚洲av在线| 男人的好看免费观看在线视频 | 99精品久久久久人妻精品| 午夜两性在线视频| 久久久久九九精品影院| 久久国产精品影院| 一边摸一边抽搐一进一小说| 又紧又爽又黄一区二区| 夜夜爽天天搞| 亚洲欧美精品综合一区二区三区| 久久久久久久精品吃奶| www日本在线高清视频| 亚洲人成77777在线视频| 亚洲精品久久午夜乱码| 久久久水蜜桃国产精品网| 久久国产精品人妻蜜桃| 亚洲男人天堂网一区| 久久久久精品国产欧美久久久| av在线播放免费不卡| 亚洲一区高清亚洲精品| 亚洲欧美日韩无卡精品| 国产麻豆69| 午夜精品国产一区二区电影| 99香蕉大伊视频| 亚洲精品中文字幕在线视频| 欧美人与性动交α欧美精品济南到| 侵犯人妻中文字幕一二三四区| 性色av乱码一区二区三区2| 人人妻人人添人人爽欧美一区卜| 18禁国产床啪视频网站| 最近最新中文字幕大全电影3 | 老熟妇乱子伦视频在线观看| 满18在线观看网站| 香蕉久久夜色| 亚洲中文av在线| 50天的宝宝边吃奶边哭怎么回事| 深夜精品福利| 国产99久久九九免费精品| 亚洲成av片中文字幕在线观看| 国产aⅴ精品一区二区三区波| 午夜a级毛片| 欧美日韩亚洲国产一区二区在线观看| 亚洲片人在线观看| 高清在线国产一区| 国产精品免费一区二区三区在线| 国产高清videossex| 80岁老熟妇乱子伦牲交| 久久午夜综合久久蜜桃| 日韩欧美三级三区| 一级毛片女人18水好多| 深夜精品福利| 国产99久久九九免费精品| 免费高清视频大片| 国产成人免费无遮挡视频| 亚洲精品粉嫩美女一区| 精品高清国产在线一区| 99国产综合亚洲精品| 欧美日韩国产mv在线观看视频| 免费高清在线观看日韩| 高清欧美精品videossex| 国产精品爽爽va在线观看网站 | 久久久久久久久久久久大奶| 精品卡一卡二卡四卡免费| 一二三四社区在线视频社区8| 久久影院123| 乱人伦中国视频| 日韩高清综合在线| 身体一侧抽搐| 在线观看日韩欧美| 亚洲成人久久性| 真人一进一出gif抽搐免费| 色哟哟哟哟哟哟| 一区二区三区精品91| 天堂俺去俺来也www色官网| 男女之事视频高清在线观看| 国产在线精品亚洲第一网站| 人人妻,人人澡人人爽秒播| 精品免费久久久久久久清纯| 欧美日韩亚洲国产一区二区在线观看| 成人av一区二区三区在线看| 国产伦人伦偷精品视频| 精品国产乱子伦一区二区三区| 久久狼人影院| 久久这里只有精品19| 亚洲三区欧美一区| 丝袜美腿诱惑在线| 亚洲 欧美一区二区三区| 一区二区日韩欧美中文字幕| 99久久精品国产亚洲精品| 亚洲精品在线美女| 亚洲av第一区精品v没综合| av国产精品久久久久影院| 80岁老熟妇乱子伦牲交| 免费不卡黄色视频| 亚洲一区二区三区色噜噜 | 午夜久久久在线观看| 自线自在国产av| 欧美精品一区二区免费开放| 99久久久亚洲精品蜜臀av| 99热国产这里只有精品6| 99久久人妻综合| 久久久久久亚洲精品国产蜜桃av| 国产精品香港三级国产av潘金莲| 国产高清国产精品国产三级| 天堂动漫精品| 国产亚洲精品综合一区在线观看 | 国产精品免费视频内射| 国产精品成人在线| 亚洲国产欧美日韩在线播放| 99久久国产精品久久久| 国产深夜福利视频在线观看| 99riav亚洲国产免费| 99久久99久久久精品蜜桃| 天天影视国产精品| 别揉我奶头~嗯~啊~动态视频| 麻豆成人av在线观看| 亚洲av电影在线进入| 天堂影院成人在线观看| 亚洲国产毛片av蜜桃av| 久久香蕉激情| 侵犯人妻中文字幕一二三四区| 免费女性裸体啪啪无遮挡网站| 国产精品电影一区二区三区| 久久天堂一区二区三区四区| 中文字幕色久视频| av天堂在线播放| 性色av乱码一区二区三区2| 日韩三级视频一区二区三区| 大码成人一级视频| 新久久久久国产一级毛片| 亚洲五月天丁香| 欧美黄色淫秽网站| www.自偷自拍.com| 国产极品粉嫩免费观看在线| 亚洲专区中文字幕在线| 制服诱惑二区| 日韩中文字幕欧美一区二区| av超薄肉色丝袜交足视频| 又黄又爽又免费观看的视频| 性色av乱码一区二区三区2| 国产一区二区三区视频了| 亚洲av美国av| 多毛熟女@视频| 免费少妇av软件| 99久久国产精品久久久| 丝袜在线中文字幕| 青草久久国产| 精品少妇一区二区三区视频日本电影| 男女下面插进去视频免费观看| 日日夜夜操网爽| 岛国视频午夜一区免费看| 精品久久久久久成人av| 国产亚洲精品第一综合不卡| 午夜精品久久久久久毛片777| 亚洲人成网站在线播放欧美日韩| 18禁美女被吸乳视频| 久久精品国产综合久久久| 国产欧美日韩一区二区三| 久久精品国产99精品国产亚洲性色 | av视频免费观看在线观看| 久久久精品国产亚洲av高清涩受| 亚洲一区二区三区色噜噜 | 国产片内射在线| 日韩精品免费视频一区二区三区| 久久久国产成人精品二区 | 免费在线观看亚洲国产| 国产99白浆流出| 亚洲av熟女| 男男h啪啪无遮挡| 久久精品aⅴ一区二区三区四区| 国产麻豆69| 妹子高潮喷水视频| 国产成人精品无人区| 女人爽到高潮嗷嗷叫在线视频| 亚洲精品国产精品久久久不卡| 91老司机精品| 亚洲成人国产一区在线观看| 亚洲avbb在线观看| 亚洲国产毛片av蜜桃av| 女人精品久久久久毛片| 男人舔女人的私密视频| 在线观看日韩欧美| 女人被躁到高潮嗷嗷叫费观| 成人手机av| 国产亚洲精品久久久久久毛片| 久久精品国产亚洲av高清一级| 在线国产一区二区在线| 亚洲av美国av| 在线播放国产精品三级| 一a级毛片在线观看| 欧美黄色淫秽网站| 日本五十路高清| 电影成人av| 亚洲一区二区三区不卡视频| 精品久久久久久久毛片微露脸| 免费久久久久久久精品成人欧美视频| 亚洲精品国产区一区二| 日韩人妻精品一区2区三区| 日韩免费高清中文字幕av| 久久中文看片网| 欧美丝袜亚洲另类 | 欧美av亚洲av综合av国产av| 国产国语露脸激情在线看| 亚洲精品中文字幕在线视频| 男人的好看免费观看在线视频 | 在线永久观看黄色视频| av网站在线播放免费| 乱人伦中国视频| 国产精品秋霞免费鲁丝片| 日韩欧美国产一区二区入口| 国产麻豆69| 新久久久久国产一级毛片| e午夜精品久久久久久久| 午夜影院日韩av| 久久久精品国产亚洲av高清涩受| 亚洲一区二区三区色噜噜 | 大香蕉久久成人网| 欧美性长视频在线观看| 亚洲中文av在线| 亚洲情色 制服丝袜| 国产av一区二区精品久久| 精品午夜福利视频在线观看一区| 新久久久久国产一级毛片| 两性夫妻黄色片| 亚洲精品中文字幕在线视频| 真人一进一出gif抽搐免费| 久久国产乱子伦精品免费另类| 黄色丝袜av网址大全| 国产精品爽爽va在线观看网站 | 人妻久久中文字幕网| 亚洲精品久久成人aⅴ小说| 精品日产1卡2卡| 男女下面进入的视频免费午夜 | 天堂俺去俺来也www色官网| 两性夫妻黄色片| 女人被躁到高潮嗷嗷叫费观| 久久香蕉激情| 久久久久久人人人人人| 人人妻人人澡人人看| 成人国产一区最新在线观看| 在线观看免费午夜福利视频| 乱人伦中国视频| 十分钟在线观看高清视频www| 波多野结衣av一区二区av| 中出人妻视频一区二区| 久久狼人影院| 欧美大码av| 深夜精品福利| 久久国产精品人妻蜜桃| 午夜a级毛片| 多毛熟女@视频| 欧美午夜高清在线| √禁漫天堂资源中文www| 亚洲精品久久成人aⅴ小说| 美女扒开内裤让男人捅视频| 啦啦啦免费观看视频1| 成人亚洲精品av一区二区 | 一区在线观看完整版| 91字幕亚洲| 性欧美人与动物交配| 一a级毛片在线观看| 级片在线观看| 亚洲精品美女久久av网站| 亚洲色图综合在线观看| 国产一区二区三区在线臀色熟女 | 18禁美女被吸乳视频| 黄色成人免费大全| 午夜精品在线福利| 波多野结衣高清无吗| 大码成人一级视频| 丰满的人妻完整版| 欧美久久黑人一区二区| 国产三级黄色录像| 亚洲专区字幕在线| 国产精品 国内视频| 久久久久久免费高清国产稀缺| bbb黄色大片| 美女大奶头视频| 成人影院久久| 亚洲av成人av| 精品卡一卡二卡四卡免费| 久久人妻熟女aⅴ| 99国产精品一区二区三区| 国产成人精品无人区| 精品少妇一区二区三区视频日本电影| 久9热在线精品视频| 女人高潮潮喷娇喘18禁视频| 欧美成人免费av一区二区三区| 精品一品国产午夜福利视频| 黄色丝袜av网址大全| www.999成人在线观看| 咕卡用的链子| 男人操女人黄网站| 两个人免费观看高清视频| av电影中文网址| 99精品久久久久人妻精品| 日韩欧美一区二区三区在线观看| 亚洲欧美一区二区三区久久| 欧美成狂野欧美在线观看| 日本 av在线| 精品国产国语对白av| 国产av又大| 国产av一区二区精品久久| 久久精品亚洲精品国产色婷小说| 久久香蕉激情| 欧美午夜高清在线| 国产精品乱码一区二三区的特点 | 两个人看的免费小视频| 久久久久久久精品吃奶| 日本五十路高清| 脱女人内裤的视频| 一边摸一边抽搐一进一小说| 日本欧美视频一区|