Identification of Novel Putative Bacterial Feruloyl Esterases From Anaerobic Ecosystems by Use of Whole-Genome Shotgun Metagenomics and Genome Binning

Feruloyl esterases (FAEs) can reduce the recalcitrance of lignocellulosic biomass to enzymatic hydrolysis, thereby enhancing biorefinery potentials or animal feeding values of the biomass. In addition, ferulic acid, a product of FAE activity, has applications in pharmaceutical and food/beverage industries. It is therefore of great interest to identify new FAEs to enhance understanding about this enzyme family. For this purpose, we used whole-genome shotgun metagenomics and genome binning to explore rumens of dairy cows, large intestines of horses, sediments of freshwater and forest topsoils to identify novel prokaryotic FAEs and trace the responsible microorganisms. A number of prokaryotic genomes were recovered of which, genomes of Clostridiales order and Candidatus Rhabdochlamydia genus showed FAE coding capacities. In total, five sequences were deemed as putative FAE. The BLASTP search against non-redundant protein database of NCBI indicated that these putative FAEs represented novel sequences within this enzyme family. The phylogenetic analysis showed that at least three putative sequences shared evolutionary lineage with FAEs of type A and thus could possess specific activities similar to this type of FAEs, something that is not previously found outside fungal kingdom. We nominate Candidatus Rhabdochlamydia genus as a novel FAE producing taxonomic unit.


INTRODUCTION
Production of biofuels and biochemicals from lignocellulosic biomass, a non-food renewable carbon resource, has increasingly become of great importance due to increasing global demands for energy and chemicals, increasing prices of fossil fuels and environmental concerns associated with fossil fuels. Lignocellulosic biomass mainly comprises three structural polymers namely cellulose, hemicellulose and lignin. In the cell walls of monocots (e.g., grass, cereals), lignin and hemicellulose interconnect, forming a matrix that encrusts the cellulose (Wong, 2006;Rubin, 2008). This configuration creates a complex structure, believed to be the main cause of recalcitrance of lignocellulosic biomass to enzymatic hydrolysis (Rubin, 2008;Pu et al., 2013).
The linkage between lignin and hemicellulose is mainly mediated by ferulic acid (FA), forming ester bonds with hemicellulose from the carboxylic side and ether bonds with lignin from the phenolic side of the molecule. These ester bonds in the cell walls of plants can be cleaved with feruloyl esterases (FAEs) (EC 3.1.1.73), member of carboxylic ester hydrolases (EC 3.1.1.-) (Jeske et al., 2019), to reduce complexity of cell wall configuration, thereby enhancing utilization of lignocellulosic biomass (Wong, 2006). Further importance of FAEs is in pharmaceutical and food/beverage industries as FA, a product of FAE activity, has evidently antioxidant properties (Pohl and Lin, 2018) and can also be used to produce vanillin (Chen et al., 2016). In addition, several attempts have already been made to improve digestibility of forages in dairy cattle rations by use of FAE producing lactic acid bacteria (Muck et al., 2018).
Feruloyl esterases are classified into four types (A, B, C, and D) based on substrate specificity against model methyl esters and release of diferulic acid (5-5') from plant cell walls (Crepin et al., 2004). The efficiency of FAEs in breaking lignin-hemicellulose interconnections seems to differ among different FAEs. FAEs-A break these interconnections in the cell walls of cereals at higher rates than FAEs-B (Crepin et al., 2004). Based on phylogenetic analysis, fungal FAEs were classified into seven subfamilies (Benoit et al., 2008) but the phylogeny was further improved in a later attempt, with recognition of 13 subfamilies of fungal FAEs (Dilokpimol et al., 2016). These attempts showed that FAEs did not evolve from a common ancestor (Benoit et al., 2008;Dilokpimol et al., 2016). In a novel approach, protein descriptors, derived from amino acid sequences, were used in conjunction with a machine learning method to classify fungal, bacterial and plant FAEs, which resulted in formation of 12 families of FAEs (Udatha et al., 2011). There is to some extent agreement between the A-D classification and 1-13 subfamily classification as for instance subfamilies 6 and 7 solely include FAEs-B and FAEs-A, respectively. However, the subfamily 1 includes both FAEs-B and FAEs-C and subfamily 5 contains FAEs-A and FAEs-D. It appears that the classification of FAEs can further be improved in the near future when more data is available.
Several fungal and bacterial species are known to produce FAEs, including Aspergillus spp., a number of anaerobic fungal species, Bacillus spp., Lactobacillus spp., etc. (Donaghy et al., 1998;Dilokpimol et al., 2016). Due to industrial significance of FAEs, there is an ever-growing interest to identify new FAEs and new microorganisms with this ability. Potential habitats of FAE producing microorganisms are ecosystems in which, plants are degraded, such as digestive tract of herbivores, soil or aquatic ecosystems. The rapid development of sequencing platforms and metagenomic methodologies has enabled to effectively explore these ecosystems for such purpose. In this work, we explored rumens of dairy cows, large intestines of horses, sediments of freshwater and topsoils of forests by means of wholegenome shotgun metagenomics and genome binning to study prokaryotic capacities for FAE production and potential novelty of the predicted FAEs.

Sampling
Approximately 50 mL rumen content was sampled from four adult Swedish Red and White breed dairy cows through permanent rumen fistula. Cows had been fed standard diets, containing forage and concentrate, based on their production levels. These cows had been fitted with fistula previously, approved by the Uppsala Ethics Committee (C 93/12 and C 142/14) and were maintained at the Livestock Research Centre of the Swedish University of Agricultural Sciences (SLU) for research/education purposes. Horse fecal samples (ca. 75 g) were directly taken from rectum of four adult horses, fed conventional forage-based diet. These horses were maintained at SLU for research/education purposes approved by the Uppsala Ethics Committee (C 148/13). All the animals used were maintained under SLU policy for use of animals in research and education (SLU.ua 2015.1.1.1-4840). Sediment samples were collected from one stream, one river, one lake and one pond from shallow locations in where, water was still and sediment contained dead plant biomass and thus, sampling locations were considered ecologically similar. Four topsoil samples were obtained from four pine-deciduous forests from locations with decaying plant biomass. All samples were collected in the region of Uppsala, Sweden during spring 2017.

Library Preparation and Sequencing
DNA extraction was done with NucleoSpin R soil (MACHEREY-NAGEL, Düren, Germany). DNA quality and quantity were checked with Agilent 2200 TapeStation System (Agilent, Santa Clara, CA, United States) and Qubit R 3.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, United States), respectively by the Science for Life Laboratory (SciLifeLab), Uppsala, Sweden. Library preparation was done with TruSeq DNA PCR-Free kit (Illumina, Inc., San Diego, CA, United States) and pairedend sequencing (2 × 125) was performed using Illumina HiSeq2500 system and v4 sequencing chemistry (Illumina, Inc., San Diego, CA, United States) in one lane by the SciLifeLab.

Annotation
Annotation of recovered genomes was done with Prokka (Seemann, 2014) (vs. 1.12, default settings) and predicted proteins were further annotated with InterProScan (Jones et al., 2014) (vs. 5.30-69.0, default settings). The InterPro database classifies proteins with similar domains/sites into single entries. The IPR011118 family includes FAE of faeB gene (accession ID: AJ309807), tannase and some other proteins. The IPR010126 family contains some lipases, FAE of faeC gene (accession ID: AJ505939) and acetyl xylan esterase. The IPR034429 family comprises FAEs-C and IPR002921 domain corresponds to a domain in FAEs-A, similar to a domain in fungal lipases.
Proteins classified as members of IPR011118, IPR034429, IPR010126 or IPR002921 entries were scanned with ScanProsite (de Castro et al., 2006) to identify sequence motifs. They were further queried with BLASTP (vs. 2.7.1, default settings) against a set of reviewed FAEs of the UniProt database (2019-10-03) (The UniProt Consortium, 2019). The protocol used to select these reference FAEs was: searching for "ec:3.1.1.73" at the UniProt database and filtering by "Reviewed." This resulted in 44 sequences among which, one sequence was incomplete (UniProt ID: P0CT85) and was thus excluded, resulting in total of 41 fungal sequences and 2 bacterial sequences. In addition, the predicted proteins belonging to the FAE-containing entries of InterPro database were subjected to BLASTP search (vs. 2.7.1, default settings) against non-redundant protein database of NCBI (2019-04-03) to assess their novelty.
Two conditions were opted for annotation as a putative FAE: more than 90% primary sequence similarity to reference FAEs or possession of the serine active site motif 2 . The consensus pattern of this motif is , with square and curly brackets indicating acceptable and unacceptable amino acids in the respective positions, respectively. An overview of the annotation protocol is in Figure 1.

Community Analysis and Taxonomic Classification
Prokaryotic community composition in pooled datasets was estimated by means of taxonomic classification of reads using Kaiju (Menzel et al., 2016) (vs. 1.7.2, default settings) and non-redundant protein database of NCBI (2019-06-25). Recovered genomes were assigned taxon with phylophlan (Segata et al., 2013) (vs. 0.99, default settings), using predicted proteins from Prokka annotation as input. The taxonomic assignments were further evaluated with CheckM and METAXA2 (Bengtsson-Palme et al., 2015) (vs. 2.2 beta 9, default settings). In case of agreement among predictions, the lowest taxonomic rank given by any of the software was reported and in case of disagreement, the lowest common taxonomic rank was assigned.

Assembly and Binning
The Cow, Horse, Sediment, and Soil datasets had 2 × 85,285,247, 2 × 46,961,631, 2 × 79,657,128, and 2 × 72,876,571 reads, respectively. Assembly statistics is shown in Table 1. The longest contig was assembled in the Soil dataset with a length of 514,904 bp. More contigs were assembled in the Cow and Horse datasets than in the Sediment and Soil datasets and the Sediment dataset had the poorest assembly statistics. Binning resulted in formation of 87, 83, 15, and 10 binned genomes in the Cow, Horse, Sediment, and Soil datasets, respectively. There was no genome redundancy based on ANI. FIGURE 2 | Prokaryotic composition (family level) of Cow, Horse, Sediment, and Soil ecosystems based on classification of reads by Kaiju. Families shown make up at least 75% of the community (together with the unclassified and unassigned sequences). "unclassified" sequences have no classification by Kaiju and "unassigned" sequences are those not assigned to a non-viral species.
FIGURE 3 | Number of binned genomes per taxon (order level) recovered from Cow, Horse, Sediment, and Soil datasets.
Frontiers in Microbiology | www.frontiersin.org Only the highest scores are shown. Proteins annotated as putative feruloyl esterase (FAE), based on the possession of serine active site (PS00120), are in bold (see Supplementary Sequence File S1 for putative protein sequences). a Bitscore, E.value, and identical match were respectively used to report subject sequences matched.

Community Composition and Recovered Taxa
The prokaryotic community composition was notably similar between the Cow and Horse samples, with the dominance of Prevotellaceae, Lachnospiraceae, Ruminococcaceae, and Clostridiaceae (Figure 2). There were also similarities between the Sediment and Soil samples, with the community comprising a wide range of taxa from Acidobacteria and Actinobacteria to Alphaproteobacteria and Betaproteobacteria (Figure 2). Viruses and archaea accounted for a combined total of at most 1.2% in each sample (data not shown).

Annotation and Phylogenetic Relationship
In total, 35 hypothetical proteins were classified as members of IPR011118, IPR010126, IPR002921 entries ( Table 2). None of the predicted proteins belonged to the IPR034429 entry. The BLASTP bitscores against reference FAEs were generally low, as were the sequence identities which ranged from 21 to 46% ( Table 2). In the BLASTP search against the full non-redundant protein database, the scores were higher, with the sequence identities ranging from 27 to 100%. One of the predicted proteins (Horse.16: FOA763) identically matched to a bacterial FAE and two others (Horse.14: NAH160; Soil.2: DAH257) had slight similarities (32 and 50%, respectively) to bacterial FAEs. Five proteins of IPR002921 entry contained the serine active site motif and were thus considered as putative FAEs (Table 3). Taxonomic classifications of binned genomes with FAE coding capacities are in Table 4. Genomes of the Clostridiales order in the Cow and Horse datasets coded for FAEs, as did a genome of Candidatus Rhabdochlamydia genus in the Sediment dataset.
Overall, the topologies of ML ( Figure 4A) and NJ ( Figure 4B) trees were similar with formation of three main clades. In the ML tree, the putative FAEs and FAEs-A formed a clade that had a moderate support (0.443). The other two clades had high supports and collectively included FAEs-B, FAEs-C and a FAE-D (Q7RWX8) that was included in our reference dataset. The main difference in the NJ tree was that the putative FAEs of Horse.7 and Horse.24 were not placed with the other putative sequences and FAEs-A in one clade but were placed basal to the other two main clades.

DISCUSSION
In this study, we used whole-genome shotgun metagenomics combined with de novo assembly and genome binning to study prokaryotic FAEs of anaerobic (cow rumen, large intestine of horse and sediment of fresh water) and microaerobic (topsoil) ecosystems. The lower assembly quality in the Sediment and Soil datasets ( Table 1) suggests that there was insufficient coverage of microbial genomes in these samples, likely due to a high microbial diversity in these two ecosystems, something evident from Figure 2.
Although several proteins from the binned genomes matched the FAE-containing entries of the InterPro database, these proteins showed very low primary sequence similarities to our reference FAEs. We therefore, explored the reference FAEs to identify sequence motifs of this enzyme family to enable a functional annotation of our predicted proteins. Surprisingly, only FAEs-A consistently contained a motif, i.e., the serine active site (PS00120), a signature of some lipases 4 . It was previously reported that FAEs-A have sequence similarities to lipases (Crepin et al., 2004) and therefore, this finding may not be entirely unexpected. The PS00120 motif is also detected in FAEs of Lactobacillus spp. (Xu et al., 2017). Our finding here indicates that attempts should be made to identify sequence features unique to FAEs to facilitate functional annotation of novel FAEs.
The PS00120 motif was only found in the protein sequences classified as members of IPR002921 domain but not in all of them (only 5 out of 16 sequences). As the IPR002921 entry describes a domain in FAEs-A, we considered these five sequences as putative FAE. The results from BLASTP search against non-redundant proteins (Table 2) suggest that these putative FAEs represent novel sequences in this enzyme family.
Species belonging to the Clostridiales order were previously reported to produce FAEs, including Butyrivibrio fibrisolvens (Dalrymple et al., 1996), B. proteoclasticus (Goldstone et al., 2010), Ruminococcus albus, and R. flavefaciens (McSweeney et al., 1998). To the best of our knowledge, this is the first report about FAE coding capacity within the Candidatus Rhabdochlamydia genus. It should be pointed out that although the Sediment dataset was constructed from a mix of different sources, the genome quality analysis showed that the considered genome had a very low contamination (∼1%, Table 4), indicating that our sampling strategy was adequate. The association between specific activities of FAEs, summarized as the A-D classification scheme (Crepin et al., 2004), and phylogenetic relationships of these enzymes is not straightforward (Dilokpimol et al., 2016), something also evident from our phylogenetic analysis. In both ML and NJ trees, two clusters were formed with each comprising FAEs of mixed specific activities. It is possible that ecological niches and specific needs of individual species largely determine the specific activity of FAEs, also pointed out by Benoit et al. (2008). This was however not the case for the FAEs-A, as in both trees the FAEs-A formed a distinct cluster, not showing close evolutionary relationships with other types of FAEs. The different evolutionary lineage of FAEs-A is further evidenced from that the PS00120 motif was only found in this type of FAEs. Interestingly, three putative sequences were consistently clustered with FAEs-A in both phylogenetic trees, suggesting that these putative FAEs may have specific activities similar to this type of FAEs, something that should be verified experimentally. Production of FAEs-A is until now only found in fungi and in particular in Aspergillus spp.

CONCLUSION
In total, 31, 44, 7, and 6 prokaryotic genomes were reconstructed from the Cow, Horse, Sediment, and Soil datasets, respectively, and were explored for FAE coding capacities. Four genomes of Clostridiales order in the Cow and Horse datasets and one genome of Candidatus Rhabdochlamydia genus in the Sediment dataset were found to have such capacity. In total, five FAEs were predicted. The results from BLASTP against non-redundant protein database of NCBI suggested that these putative FAEs are novel. Phylogenetic analysis suggested that at least three putative sequences might have specific activities similar to FAEs-A.

DATA AVAILABILITY STATEMENT
Raw data is deposited at the Sequence Read Archive database under PRJNA543979 accession number.

AUTHOR CONTRIBUTIONS
KM designed the study and did the sampling, DNA extraction, bioinformatic analyses and preparation of the first draft of the manuscript. JS provided bioinformatic expertise and contributed to the manuscript preparation.

ACKNOWLEDGMENTS
We acknowledge Prof. Anna Schnürer from the SLU, Sweden, and Dr. Sabine Kleinsteuber from the Helmholtz Centre for Environmental Research-UFZ, Germany, for their valuable inputs on experimental setup. Sequencing was performed by the SNP and SEQ Technology Platform in Uppsala. The facility is part of the National Genomics Infrastructure (NGI) Sweden and SciLifeLab. The SNP and SEQ Platform is also supported by the Swedish Research Council and the Knut and Alice Wallenberg Foundation. The computations were performed on resources provided by SNIC through Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) under Project SNIC 2017/7-306.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb. 2019.02673/full#supplementary-material TABLE S1 | Genome quality and taxonomic classifications of recovered genomes searched for feruloyl esterase coding capacities.