Protein Domain Analysis of Genomic Sequence Data Reveals Regulation of LRR Related Domains in Plant Transpiration in Ficus

Predicting protein domains is essential for understanding a protein’s function at the molecular level. However, up till now, there has been no direct and straightforward method for predicting protein domains in species without a reference genome sequence. In this study, we developed a functionality with a set of programs that can predict protein domains directly from genomic sequence data without a reference genome. Using whole genome sequence data, the programming functionality mainly comprised DNA assembly in combination with next-generation sequencing (NGS) assembly methods and traditional methods, peptide prediction and protein domain prediction. The proposed new functionality avoids problems associated with de novo assembly due to micro reads and small single repeats. Furthermore, we applied our functionality for the prediction of leucine rich repeat (LRR) domains in four species of Ficus with no reference genome, based on NGS genomic data. We found that the LRRNT_2 and LRR_8 domains are related to plant transpiration efficiency, as indicated by the stomata index, in the four species of Ficus. The programming functionality established in this study provides new insights for protein domain prediction, which is particularly timely in the current age of NGS data expansion.


Introduction
With the advent of next-generation sequencing (NGS) technology, a massive amount of DNA data is currently being produced in both model and non-model species. However, there are many problems associated with de novo assembly, i.e., when there is no reference genome on which to map reads, especially when the genome structure is complex with large parts of repetitive elements, as it is often the case in plant species [1]. In such cases, the DNA reads can only be assembled to scaffold or contig level [2]. Thus, methods based on an analysis of the fragments are needed.
A protein domain is a conserved part of a protein sequence which has a specific structure and function. The typical length of a protein domain is from about 25 to 500 amino acids. For some protein domain analysis, the whole protein sequence is not required [3]. Hence, some of the problems associated with fulllength assembly without a reference genome can be avoided by protein domain analysis.
In the present study, fig trees belonging to the Ficus genus of the Moraceae family were examined to verify the above hypothesis. The Ficus genus has been found to have great diversity in tropical and subtropical areas, which is linked to geographical evolution within the genus [4,5]. Ficus altissima Blume, Ficus tinctoria G.
Forst, Ficus langkokensis Drake and Ficus fistulosa Reinw. ex Blume usually have overlapping distributions. However, their ecological niches are different due to their physiology. F. altissima and F. tinctoria are semi-epiphytic and their leaves are coriaceous. As a result, they can tolerate environments with drought episodes [6]. In contrast, F. langkokensis and F. fistulosa grow in relatively humid habitats, such as waterside rocks, and their leaves are thin coriaceous [7]. The ecological differences in the growing areas of these different Ficus species might thus exert different types of drought stress pressures, leading to different responses in stomatal development and morphology [8]. Hence, it would be valuable to develop a model that predicts the peptide domains of proteins for genes potentially involved in responses to drought stress, using genomic data.
One of the strategies used by plants to respond to drought stress events is plant transpiration efficiency. In the model plant Arabidopsis, plant transpiration efficiency is a quantitative trait, which has been shown to be controlled by several genes based on quantitative trait loci (QTLs) mapping studies [9]. To date, only a few contributing genes have been identified, one of which is the ERECTA gene, which explains 21-46% of the total phenotypic variation in D(leaf carbon isotopic discrimination) [9]. In Arabidopsis, ERECTA is one of the best studied receptor like kinases (RLKs) with leucine rich repeat (LRR) domains, which not only participates in plant transpiration efficiency but also regulates aerial architecture, stomatal patterning and confers resistance to the pathogenic bacteria Ralstonia solanacearum, the necrotrophic fungi Plectosphaerella cucumerina and Pythium irregulare [10,11]. Structurally, the protein encoded by the ERECTA gene in Arabidopsis has one LRRNT_2 protein domain at the N-terminal, two LRR_8 protein domains in the middle part, and one Pkinase domain at the C-terminal (Fig. 1A). The LRR_8 domains form the hydrophobic core of the proteins, and they are frequently involved in the formation of protein-protein interactions [11,12]. The LRRNT_2 domain of the protein encoded by ERECTA in Arabidopsis has LRRs flanked by cysteine rich sequences (Fig. 1B).
In contrast to model species, the molecular mechanism of plant transpiration efficiency still remains unclear in many plant and tree species, especially those without reference genomes. Improving functional annotation of assembled data obtained from NGS technology may provide new insights into genes potentially involved in this important trait. In this study, our first objective was to develop a method for obtaining high quality contigs from low coverage NGS data. Secondly, we attempted to predict protein domains from contigs obtained via the above method. Finally, we utilized our programming functionality to predict LRR domains homologous to those from the Arabidopsis ERECTA gene in four Ficus species that respond differently to drought environments and examined the relationship between LRR domain numbers and plant transpiration efficiency.

DNA extraction and genome sequence
Leaf material of four species, F. altissima, F. tinctoria, F. langkokensis and F. fistulosa, was collected from the Xishuangbanna Tropical Botanical Garden, Yunnan Province, P. R. China (101u259E, 21u419N) in April 2013 and stored in a paper bag with silica gel until DNA extraction. The four species had been  . The proposed programming functionality for predicting protein domains directly from genomic sequence data without a reference genome. The Illumina reads were first trimmed with quality control methods. Then, assembly software ABySS, SOAPdenovo and Velvet were used separately to obtain original contigs. Next, length control methods were used to select contigs larger than 250 base pairs. Afterwards, the assembly software Phrap was used to obtain final contigs and Genscan was used to predict peptides from these contigs. Finally, Hmmsearch was used to predict protein domains. doi:10.1371/journal.pone.0108719.g002

Quality control methods
The raw data from the Illumina Hiseq2000 were trimmed with two programs for performing quality control written in the Practical Extraction and Report Language (PERL). The first program was used to remove nucleotides with a Phred score lower than 20 (Script S1). The second program was used to delete fastq reads with length less than 20 base pairs as well as ''orphanage'' reads (single reads not in a pair) created by the first program (Script S2).

Sequence assembly
To generate a better genome assembly, we used a combination of four popular assembly software packages: ABySS, SOAPdenovo, Velvet and Phrap. ABySS, SOAPdenovo and Velvet were used to align the trimmed Illumina fastq reads to obtain contigs. These contigs were then aligned again with Phrap to improve the alignment.
First of all, ABySS (http://www.bcgsc.ca/platform/bioinfo/ software/abyss) which allows de novo, parallel, paired-end sequence assembly for short sequence reads was used to construct alignments [13][14][15] on our Ficus genome data. We employed 25 as the k-mer length and 10 as the minimum number of pairs needed for two joined contigs.
Secondly, SOAPdenovo (http://soap.genomics.org.cn/ soapdenovo.html) [16], which is particularly designed to assemble Illumina GA short reads, was used for building the contigs. The detailed parameter set was as follows: k-mer length 25; average insert size 250; cutoff value of pair number for a reliable connection between two contigs of pre-scaffolds 3; and minimum alignment length between a read and contig required for a reliable read location 32.
Thirdly, Velvet (http://www.ebi.ac.uk/,zerbino/velvet/) [17], which is a sequence assembler for very short sequence reads, was also applied for the sequence alignment. We set the k-mer length as 25 and the average insert size as 250.
Finally, Phrap (http://www.phrap.org/) [18], which is a program for assembling shotgun DNA sequence data was further applied on the sequence to increase the maximum length and remove redundancy. We analyzed the results of ABySS, SOAPdenovo and Velvet by Phrap (for parameters see Table S1 and some connection Script S3).
Gene structure identification GENSCAN (http://genes.mit.edu/GENSCANinfo.html) was used to identify complete gene structures in genomic DNA. It is a GHMM-based program that can be used to predict the location of genes and their exon-intron boundaries in genomic sequences are from a variety of organisms. The ''Arabidopsis.smat'' file was downloaded and used as parameter file for the Ficus genome data [19].
Protein domain prediction HMMER 3.0 (http://hmmer.janelia.org/) was used for searching sequence databases for homologs of protein sequences and making protein sequence alignments [20]. It employs methods using probabilistic models called profile hidden Markov models (profile HMMs). We used hmmscan to predict protein domains in the gene ERECTA, which were then predicted in Ficus and Populus by hmmsearch.

Experimental analysis
Stomata index evaluation. The study was conducted in the Xishuangbanna Tropical Botanical Garden in Yunnan Province, P. R. China (101u259E, 21u419N) in August 2013. Four to six trees of each species showing good growth performance were sampled. We collected three mature and well-exposed leaves from each tree. To obtain a better view of the stomata, we removed the main vein of leaves and then boiled them in hot alkaline buffer to remove the mesophyll. Treated leaves were examined under a light microscope (DM2500, Leica, Germany). The numbers of stomata and epidermal cells were counted using ImageJ (National Institutes of Health, Bethesda, MD, USA, http://rsb.info.nih.gov/ij/index. html). We used the formula (1006stomatal density)/(stomatal density + epidermal cell density) to calculate the stomata index as in Hara et al., (2009) [21]. We repeated the experiments six times.

Results and Discussion
Genomic contigs assembled from Illumina genomic sequence data The Illumina Hiseq2000 genomic sequence data in fastq format for each of F. altissima, F. tinctoria, F. langkokensis and F. fistulosa was about 15 Gigabytes, and thus the coverage was about 15x [22]. The quality control methods improved the quality of the Illumina genomic sequence data for all four species, where the number of reads was about 135 million for each species. For the gene and protein domain prediction, it was most appropriate to use the contigs that had a base pair number .250 as the input for Phrap as the typical number of nucleotides in a DNA fragment encoding a protein domain of the LRR family is about 250 (Fig. 2). The maximum lengths of contigs predicted by Phrap in F. altissima, F. tinctoria, F. langkokensis and F. fistulosa were 6914 base pairs, 10755 base pairs, 6665 base pairs and 9203 base pairs, respectively. These numbers showed that the genomic contigs finally assembled could be used for the gene prediction (Table 1).

Peptides and protein domains predicted from genomic contig data
Genscan was used to predict peptide sequences from the genomic contigs assembled by Phrap for F. altissima, F. tinctoria, F. langkokensis and F. fistulosa (Fig. 2). The maximum lengths of peptides predicted by Genscan for the four species (880, 467, 550 and 1576 amino acids, respectively) were employed since longer peptides enable better protein domain prediction. These numbers indicated that the predicted peptides could be used to predict protein domains (Table 1).
HMMER was used to predict protein domains from the peptides predicted from Genscan for the four species (Fig. 2). The numbers of LRRNT_2 domains predicted from the contigs assembled by Phrap were 19,9,3, and 2 in F. altissima, F. tinctoria, F. langkokensis and F. fistulosa, respectively, whereas the numbers of LRR_8 domains were 132, 62, 45, and 19, respectively (Table 1).

Phrap can improve assembly of contigs and remove identical segments of genomic sequences
Phrap cannot directly work on Illumina fastq reads. However, it can increase the maximum length of the contigs assembled through ABySS, SOAPdenovo and Velvet, which can directly work on Illumina fastq reads (Table 1). Thus, the maximum length of the peptides predicted by Genscan may also be increased (Fig. 3).
Assembly through Phrap can remove redundancy of the contigs, including the identical genomic sequence segments predicted by ABySS, SOAPdenovo and Velvet. The percentage of redundancy removed from F. altissima, F. tinctoria, F. langkokensis and F. fistulosa was 36.51, 45.91, 18.64 and 46.91 respectively ( Table 2). The Phrap software can be used to combine identical DNA fragments into one sequence, thus avoid the effect of gene expression difference produced by NGS methods. Programs like CAP3 and TIGR Assembler may also offer similar functions. However, it is important to choose correct parameters for different species. When applying on the contigs which are assembled by NGS assembly methods, we found that Phrap has more suitable

The numbers of LRRNT_2 and LRR_8 domains in Ficus correlate with stomata index values
To test whether the LRRNT_2 domains and LRR_8 domains are related to transpiration efficiency, we used our programming functionality to predict their numbers in the four species of Ficus. The mean values of the stomata index for F. altissima, F. tinctoria, F. langkokensis and F. fistulosa were 5.3, 10.9, 13.6 and 15.9, respectively ( Table 3). As the stomata index values increased in these species the numbers of LRRNT_2 and LRR_8 domains decreased accordingly (Fig. 4). To eliminate the contingency in protein domain selection, we used the actin domain from actin1 protein in Arabidopsis thaliana (NCBI accession number NP_850284.1) for control analysis. Actin is a house-keeping protein expressed in every plant cell as a component of the cytoskeleton [23], and thus provides a good control. Among all the peptides predicted for the four Ficus species, one actin domain was found to be longer than 100 amino acids and another was shorter than 50 amino acids (Fig. 4). These results suggest that the transpiration efficiency could be related to the LRRNT_2 and LRR_8 domains in Ficus.
The ERECTA gene has not only a positive regulatory role on respiration in drought conditions but also benefits plants in the absence of water shortage [9]. Therefore, the protein domains in the ERECTA gene might show a cumulative positive evolution. The LRR_8 domain has more LRRs than the LRRNT_2 domain, and thus may have a more important role in protein-protein interactions (Fig. 1). Hence, this could explain why the number of LRR_8 domains was more than that of LRRNT_2 domains (Fig. 4).

Conclusion
The programming functionality in this study was proved to be a useful tool in biological studies by showing that the LRRNT_2 and LRR_8 domains were potentially related to plant transpiration efficiency, as we can see from the stomata index in F. altissima, F. tinctoria, F. langkokensis, and F. fistulosa. The main benefit of the functionality is that it overcomes many of the complex problems associated with de novo assembly. However, with the increasing read lengths produced by NGS and improvements in thirdgeneration sequencing, such problems may also be solved with the rapid developments of de novo assembly methods. The main limitation of the functionality is GENSCAN prediction step, which requires a suitable model. In addition, it is hard for some species to choose a perfect model to predict the gene structure. Confronting with this situation, researchers normally prefer to pick a widely used model which turns out to have more or less shortage. Nevertheless, methods of whole genome protein domain analysis will still help researchers to better understand some mechanisms of biological function from the perspective of genetic sequence, if combined with a large amount of NGS data. Script S1 Perl program used to remove the nucleotides which have Phred score lower than a specific value.

(DOCX)
Script S2 Perl program used to delete the fastq reads which have length less than a specific value as well as to erase the ''orphange'' reads (single reads without pair).

(DOCX)
Script S3 Perl programs used for dealing with the assembly files which were created by Phrap as well as for making statistic analysis. (DOCX)