Background
==========

Horseradish peroxidases (HRPs) originating from the perennial herb *Armoracia rusticana* (*Brassicaceae*) are heme-containing monomeric glycoproteins belonging to the class III plant peroxidase subfamily \[[@B1]\]. These versatile enzymes have traditionally been utilized as reporters in various diagnostic assays and histochemical stainings but have also gained increasing interest in other life science and biotechnological applications ranging from cancer therapeutics \[[@B2]\], protein engineering \[[@B3]\] and transgenics \[[@B1]\] to bioremediation \[[@B4]\], biosensors \[[@B5]\] and biocatalysis \[[@B6]\]. The *in vivo* functions of HRPs have not been fully elucidated owing to the estimated large number of isoenzymes (up to 42) \[[@B7]\], but are known to be very diverse \[[@B8]\] thus offering a wide range of substrate specificities and applications. Although HRP has been studied for decades and in spite of the large diversity of this enzyme family, protein engineering and heterologous expression have mainly been focused on one single isoenzyme C1A, thus neglecting the potential of all others. This was largely due to the lack of sequence information and the low efficiency of HRP expression in heterologous hosts. Commercial preparations are still extracted from the roots of *A. rusticana* and therefore consist of a mixture of various isoenzymes. The quality of these preparations varies greatly and depends on several biotic and abiotic factors, such as seasonal change or origin. Chromatographic purification is needed to isolate highly enriched isoenzyme fractions. Only very few purified isoenzymes were accessible so far and their substrate specificities and enzymatic properties are poorly characterized.

The sequences of only eight isoenzymes are currently known: six nucleotide sequences (C1A, C1B, C1C, C2, C3, N) have previously been published \[[@B9]-[@B11]\] and two amino acid sequences A2 (P80679) \[[@B12]\] and E5 (P59121) \[[@B13]\] have been determined. However, for decades, HRP has been known as a large group of enzymes with versatile physical properties. Already in 1966, Shannon *et al.*\[[@B14]\] described the physical properties of seven isoenzymes. Aibara *et al*. \[[@B15],[@B16]\] characterized five neutral isoenzymes (B1, B2, B3, C1 and C2) and six basic isoenzymes (E1-E6) in 1982 and 1981, respectively. A total number of 42 peroxidase isoenzymes have been identified by isoelectric focusing in commercially available HRP preparations \[[@B7]\], without knowing whether these isoenzymes differ in amino acid sequence or due to different post-translational modifications.

This study has two major goals. First, we resolve the transcriptome sequence of *A. rusticana*, a perennial plant of industrial and medical importance. The availability of a large expressed sequence tag (EST) collection is crucial to support annotation in possible future *A. rusticana* genome sequencing projects. Secondly, we demonstrate the use of an efficient enzyme discovery pipeline including new generation cDNA sequencing technologies, *in silico* isoenzyme discovery and experimental sequence verification, gene synthesis and enzyme production and secretion by *Pichia pastoris* (*Komagataella pastoris*) as a straightforward approach to discover and characterize new isoenzymes from plants or other eukaryotes. Transcriptomes deliver all sequences of expressed genes, at the same time avoiding sequencing introns and providing information about all expressed exons and alternative exon junctions. Studies in *Arabidopsis thaliana*, a model species of the *Brassicaceae* family, have demonstrated the power of massively parallel transcriptome sequencing in providing high-quality representation of transcripts needed for gene discovery \[[@B17]\]. Similar transcriptome sequencing approaches have previously been applied, for example, in marker development, population genomics \[[@B18]\] and predictions of biosynthetic pathways \[[@B19]-[@B23]\]. Over the course of the study, the new 'third-generation' and 'fourth-generation' sequencing techniques have revolutionized the speed and cost of the transcriptome sequencing projects \[[@B24]\]. Hundreds of transcriptomes have been sequenced and annotated, especially by the large sequencing projects 1KP \[[@B25]-[@B27]\], PhytoMetaSyn \[[@B28],[@B29]\] and Medicinal Plant Genomics Resource \[[@B30],[@B31]\]. However, the concept of novel isoenzyme discovery from the large bulk of sequences generated by next generation sequencing (NGS) technologies needs a full pipeline of efficient tools. This is the first study using NGS transcriptome sequencing to discover, discriminate and characterize large numbers of sequence verified isoenzymes of non-model plant origin. Although a similar study was recently performed to identify fungal cellulases \[[@B32]\], the method described utilized combined secretome and transcriptome analyses and was only aiming to show cellulose activity of the cloned cDNA without the need of the full verified sequences to make all discovered isoenzymes available by recombinant expression. The method described in this study can be widely applied for the replenishment of the sequence data in any eukaryotic organism including fungi, plants and animal cell lines or tissues when detailed sequence and gene structure information of enzymes and isoenzymes is needed.

Results
=======

Sample preparation and cDNA library generation
----------------------------------------------

The high quality (RNA integrity number RIN 9.4-9.8) of the RNA samples was confirmed with an Agilent 2100 bioanalyzer. In order to include the majority of all encoded isoenzymes, a mixture of RNA from diverse plant parts, including leaves, roots, sprouts and stems, was chosen for mRNA isolation. cDNA synthesis, normalization, size selection and cloning was performed by LGC Genomics (Berlin, Germany). The normalized cDNA library was subjected to quality control experiments before using it for 454 pyrosequencing: a cDNA fragment size of over 800 bp was ensured and the normalization efficiency was verified by sequencing 96 randomly selected clones (LGC Genomics).

Sequencing and *de novo* assembly
---------------------------------

The normalized cDNA library was sequenced on half a picotiterplate run on the GS FLX using Roche 454/Titanium chemistry. A total of 592507 sequence reads with an average read length of 353 ± 122 nucleotides were obtained. A total of 12798 (2.16%) clonal reads (exact, 3' or 5') were detected. Prior to assembly, the sequence reads were screened for the linker sequence used for concatenation, the linker sequences were clipped and the reads were quality checked (LGC Genomics). The resulting 556269 reads with an average length of 343 bp (Figure [1](#F1){ref-type="fig"}A) were further filtered by Newbler sequence filtering to ensure consistent high quality of the reads used in the assembly. 490285 reads were aligned to individual transcripts using Newbler version 2.5.3 with default settings. The *de novo* assembly generated 18511 contigs with an average length of 718 bp (Figure [1](#F1){ref-type="fig"}B) and an average coverage of 11.4-fold (Figure [1](#F1){ref-type="fig"}C). The contigs were further processed to 14871 isotigs with an average length of 1133 bp (Figure [1](#F1){ref-type="fig"}D). 35950 reads were left as singletons. A detailed summary of the alignment and assembly process is described in Table [1](#T1){ref-type="table"}.

![**Overview of the sequencing and assembly of the*A. rusticana*transcriptome. (A)** Size distribution of the quality-filtered reads. Total number of reads: 556269, average/median length 342.9/379.0 **(B)** Length distribution of the 18511 contigs. Average/median length of the contigs: 717.7/667.0 **(C)** Coverage distribution of the 18511 contigs. Average/median coverage of the contigs: 11.4/7.8 **(D)** Length distribution of the14871 isotigs. Average/median length of the isotigs: 1133.3/1113.0.](1471-2164-15-227-1){#F1}

###### 

Summary of the transcriptome sequencing, assembly and enzyme discovery results

                                                          **Number of sequences**   **% of total**
  ----------------------------------------------------- ------------------------- ----------------
  **Transcriptome sequencing and assembly**                                       
  Total number of reads                                                   592,507               \-
  Clipped reads                                                           556,269               \-
  Reads after Newbler quality control                                     543,439            100.0
  Reads aligned                                                           490,285             90.2
  Reads assembled                                                         433,179             88.4
  Reads partially assembled                                                57,064             10.5
  Singletons                                                               35,950              7.3
  Contigs                                                                  18,511               \-
  Isogroups                                                                10,619               \-
  Isotigs                                                                  14,871               \-
  Average isotig length                                                     1,133               \-
  Largest isotig length                                                     3,659               \-
  Average contig coverage                                                    11.4               \-
  **HRP isoenzyme discovery and verification**                                    
  \# of isotigs with a secretory peroxidase domain                             18               \-
  \# of full length peroxidase genes                                           18               \-
  \# of peroxidase genes after manual revision                                 20               \-
  Total \# of isoenzymes (including allelic variants)                          28              100
  Successfully verified from gDNA                                              26             89.7
  **Enzyme production in*Pichia pastoris***                                       
  Synthetic genes for production                                               26             89.7
  Successful production of an active isoenzyme                                 22             75.9

Plant peroxidase search and manual validation of assembled transcripts
----------------------------------------------------------------------

The position-specific scoring matrix (PSSM) corresponding to known horseradish peroxidases (cd00693) was used in a tblastn search in the assembled transcripts. The settings allowed for a very permissive filtering of putative peroxidases, thus including many false positives, but avoiding the loss of valuable data for further analyses. This search yielded hits in 91 transcripts, which were classified in secretory peroxidases, ascorbate peroxidases, glutathione peroxidases and peroxidase-like proteins by definition of the Conserved Domains Database (CDD, NCBI \[[@B33]\]). All previously known HRPs were classified as secretory peroxidases, so only the contigs comprising a secretory peroxidase conserved domain were kept. The horseradish transcriptome contigs of the 18 resulting secretory peroxidases were manually reviewed. In this process, the coding sequences of four contigs were extended with available assembly data, three contigs were split because of strongly conflicting reads, and two more contigs were discarded because of only a partial domain match that could not be resolved into a full-length sequence. In total, 20 *HRP* genes, with allelic variants corresponding to 28 peroxidase isoenzymes, were identified in the transcriptome of *A. rusticana.* This includes isoenzymes C1A and C3 which could be partially retrieved from the raw reads although they did not form a full-length contig. No read - even partially - corresponding to the previously published "neutral isoenzyme" N (Q42517) \[[@B11]\] was found.

Sanger sequencing and genome walking
------------------------------------

Sequences yielded by the transcriptome assemblies are not necessarily error free but can include incorrect information either caused by the transcription and RNA editing machineries of the plant \[[@B34],[@B35]\], introduced in the sequencing process or resulting from misassemblies. The sequences of all peroxidase genes detected in the isotigs were verified on genome level by Sanger sequencing of amplified genomic DNA. In addition, the sequences of the isoenzymes C1A and C3 available in the databases and partially also found in the raw reads were revised. In case of five full-length contigs where no sufficient read data from untranslated regions was available to enable amplification and sequencing of the complete gene, a genome walking approach was successfully performed in order to verify the 5' and/or 3' regions of the respective gene. Divergent coding sequence information was observed for ten genes in the form of possible allelic variants. PCR artifacts were ruled out by repeated experiments to increase the coverage of the positions. Thus, conflicting sequence information was postulated not to be due to sequencing errors, but rather due to the high sequence similarity of the HRP isoenzymes and the permissive settings used in the assembly. Supporting this assumption, putative allelic sequences could be found as separate raw reads. For altogether nine positions in contigs 22684, 6117, 17517 and 23190 (Table [2](#T2){ref-type="table"}) the nucleotide present in the transcriptome reads could not be found in the genomic DNA sequenced. The sequence of the previously published "neutral isoenzyme" N (Q42517) not present in the transcriptome sequences could not be amplified from genomic DNA (gDNA) either. Therefore, it was not included in the following analyses or experiments. The sequences of the transcript \# 22489 (see Table [2](#T2){ref-type="table"}) could not be verified.

###### 

Comparison of the HRP isoenzyme sequences between GenBank, UniProt, transcriptome and verified genome sequences

    **HRP**    **Sanger sequence**   **Transcriptome sequence**   **GenBank sequence**   **UniProt sequence**                    
  ----------- --------------------- ---------------------------- ---------------------- ---------------------- -------- -------- -----
      C1A              TA                       Y37                        \-                     \-              AT      I37     Y37
    109-110          109-110                                                                                                     
     C1159           intron                      \-                      intron                 G1159           intron     \-    
      C1B            T/C253                    intron                      \-                   intron           T253    intron   \-
    T/C859           intron                      \-                      intron                  C859           intron     \-    
      C1C              ss                        ss                        ss                     ss              \*       \*     \*
    nt1-60           aa1-20                    nt1-60                    aa1-20                                                  
     C178              R60                      C178                      R60                    A118            S40      S40    
    A/T1335          intron                      \-                      intron                   \-              \-       \-    
    A/G1888          T/A165                    A/G493                    T/A165                  G433            A145     A145   
     C1889            A165                    C/T1889                    A/V165                                                  
    C/G1921          Q/E176                    C/G526                    Q/E176                  G466            E156     E156   
      C2               CT                      intron                      \-                   intron            \*     intron   \-
   1250-1251                                                                                                                     
     A1334           intron                      \-                      intron                   \*            intron     \-    
      C3             G/T1294                   intron                      \-                   intron          G1294    intron   \-
    A/T1323          intron                      \-                      intron                 A1323           intron     \-    
    T/C1484           L231                       \-                        \-                   T1484            L231     L231   
    C/T1541           F250                       \-                        \-                   C1541            F250     F250   
      A2               ss                        ss                        ss                     ss              \-       \-     \*
    nt1-93           aa1-31                     1-93                     aa1-31                                                  
      AAT              N78                      AAT                       N47                     \-              \-      D47    
    231-234          231-234                                                                                                     
      GGA             G220                      GGA                       G220                    \-              \-      N189   
    996-998          661-663                                                                                                     
      AAT             N221                      AAT                       N221                    \-              \-      G190   
   999-1001          664-666                                                                                                     
      ACG             T284                      ACG                       T284                    \-              \-      L253   
   1185-1187         850-852                                                                                                     
    G/A1203          A/T290                     G868                      A290                    \-              \-      A259   
      AAT             N334                      AAT                       N334                    \-              \-      D303   
   1335-1337        999-1002                                                                                                     
      E5               ss                        ss                        ss                     ss              \-       \-     \*
    nt1-81           aa1-27                    nt1-81                    aa1-27                                                  
     T419              L82                     C/T246                     L82                     \-              \-      L55    
     C422              D83                     T/C249                     D83                     \-              \-      D56    
     C545             C124                     T/C372                     C124                    \-              \-      C97    
     01805            None                      none                      none                   none             \-       \-     \-
     22684            G1611                     R337                     A1010                   K337             \-       \-     \-
      TGA             D343                      CGG                       G343                    \-              \-       \-    
   1627-1629        1026-1028                                                                                                    
     01350            None                      none                      none                   none             \-       \-     \-
     02021            None                      none                      none                   none             \-       \-     \-
     03523            None                      none                      none                   none             \-       \-     \-
     06117             T30                      V10                      C/T30                   V10              \-       \-     \-
     C1088            I269                      T807                      I269                    \-              \-       \-    
     17517            T190                      Y64                       C190                   H64              \-       \-     \-
     C1157            G282                      T846                      G282                    \-              \-       \-    
     A1232            K307                      G921                      K307                    \-              \-       \-    
    08562.1           None                      none                      none                   none             \-       \-     \-
    08562.4           None                      none                      none                   none             \-       \-     \-
     23190            T1345                     S109                     G1345                   S109             \-       \-     \-
     C1423            G135                     T1423                      G135                    \-              \-       \-    
     T1842            S222                    T/C1842                    S/P222                                                  
     C1850            T224                    A/C1850                     T224                    \-              \-       \-    
     A2221            E348                    T/A2221                    V/E348                   \-              \-       \-    
     04663            None                      none                      none                   none             \-       \-     \-
     06351            None                      none                      none                   none             \-       \-     \-
     05508           G/A346                    A/T116                    G/A346                 A/T116            \-       \-     \-
     22489             \-                        \-                      G/A597                  T199             \-       \-     \-
       .                .                      G/T715                    A/S239                   \-              \-       \-    

All nucleotide (nt) and amino acid (aa) positions were calculated from the start ATG. Variations between the nt positions of the transcriptome sequence compared to other sequence sources are either due to deletions or intronic sequences in the other sources. "-" indicates that no sequence information is available from the respective source. "ss" indicates a putative signal sequence. Deletions or missing sequences are marked with "\*". "N/NX" indicates a variation at position X. "None" indicates that no polymorphisms or differences between transcriptome and genome sequence were found. The isoenzymes C1A and C3 were detected in the transcriptome raw reads with only partial/low coverage (0-2x), thus no consensus sequence was formed.

GC content and codon usage
--------------------------

The average GC content of all 14871 isotigs was calculated to be 42.7% (range 28%-62%) (Figure [2](#F2){ref-type="fig"}), almost identical to the average GC content, 42.5%, reported for *A. thaliana*\[[@B36]\]. This is lower than the average GC content of 45.1% reported by the Codon Usage Database (CUD) \[[@B37]\], suggesting that the small set of available *A. rusticana* genes (14 coding sequences) used in CUD for the calculation is not representative for the whole species. The GC content of the HRP isoenzymes was observed to vary between 42.9% (C2, contig \#04627) and 51.0% (contig \#22489), with an average GC content of 47.1%. The results of the analysis are described in more detail in Table [3](#T3){ref-type="table"}.

![**GC content distribution of the*A. rusticana*isotigs.** GC content distribution of the *A. rusticana* isotigs varies from 28% (min) to 62% (max) with a range of 35%. The average GC content of all isotigs is 42.72%. Mode (x-axis value) 43, mode value (y-axis value) 2376.](1471-2164-15-227-2){#F2}

###### 

Summary of the horseradish peroxidase isoenzymes and associated data produced during this study

   **Contig number**   **Name (\* novel)**   **Length nt**   **GC content%**   **Accession \#**   **CAI calculated**   **Intron \#**   **Signal sequence length**                               **Disulfide bridges (EDBCP prediction)**                               **pI mature protein**
  ------------------- --------------------- --------------- ----------------- ------------------ -------------------- --------------- ---------------------------- -------------------------------------------------------------------------------------------------- -----------------------
          \-                   C1A               1062             43.69            HE963800             0.812                3                     30                                        \[11--91\] \[44-49\] \[97--301\] \[177--209\]                                     5.59
         15901                 C1B               1056             43.94            HE963801             0.808                3                     28                                        \[11--91\] \[44-49\] \[97--301\] \[177--209\]                                     5.84
         25148                C1C\*              1059             45.14            HE963802             0.809                3                     29                                        \[11--91\] \[44-49\] \[97--301\] \[177--209\]                                     6.49
        25148_2               C1D\*              1059             45.04            HE963803             0.810                3                     29                                        \[11--91\] \[44-49\] \[97--301\] \[177--209\]                                     7.04
         04627                 C2                1044             42.91            HE963804             0.800                3                     24                                        \[11--91\] \[44-49\] \[97--301\] \[177--209\]                                     8.56
          \-                   C3                1050             46.76            HE963805             0.781                3                     29                                        \[11--91\] \[44-49\] \[97--300\] \[177--209\]                                     7.71
    Manual assembly           A2A\*              1011             46.79            HE963806             0.761                3                     31                                        \[11--91\] \[44-49\] \[97--299\] \[176--208\]                                     4.93
    Manual assembly           A2B\*              1011             46.69            HE963807             0.761                3                     31                                        \[11--91\] \[44-49\] \[97--299\] \[176--208\]                                     4.93
         04382                 E5                1044             46.07            HE963808             0.771                3                     27                                        \[11--91\] \[44-49\] \[97--300\] \[177--209\]                                     8.84
         01805               1805\*              1065             44.41            HE963809             0.797                3                     31                                        \[11--91\] \[44-49\] \[97--301\] \[177--209\]                                     5.97
         22684              22684.1\*            1050             46.76            HE963810             0.770                3                     29                                        \[11--91\] \[44-49\] \[97--300\] \[177--209\]                                     6.98
        22684_2             22684.2\*            1050             46.67            HE963811             0.772                3                     29                                        \[11--91\] \[44-49\] \[97--300\] \[177--209\]                                     6.37
         01350               1350\*               975             50.97            HE963812             0.707                3                     28                                        \[11--91\] \[44-49\] \[97--292\] \[176--201\]                                     8.67
         02021               2021\*               996             46.08            HE963813             0.788                3                     29                                               \[11--89\] \[44-49\] \[95--297\]                                           9.46
         23190              23190.1\*            1080             49.26            HE963817             0.724                2                   31 42              \[11--92\] \[44-49\] \[98--293\] \[178--205\] or \[22--103\] \[55-60\] \[109--304\] \[189--216\]         8.40 6.58
        23190_2             23190.2\*            1080             49.17            HE963817             0.722                2                   31 42              \[11--92\] \[44-49\] \[98--293\] \[178--205\] or \[22--103\] \[55-60\] \[109--304\] \[189--216\]         8.60 7.09
         04663               4663\*              1077             47.82            HE963814             0.748                3                     31                                        \[11--91\] \[44-49\] \[97--299\] \[176--208\]                                     4.48
         06351               6351\*               945             43.39            HE963816             0.786                3                     18                                        \[17--96\] \[50-55\] \[102--292\] \[180--206\]                                    6.37
         03523               3523\*               960             44.58            HE963820             0.761                0                     22                                        \[11--92\] \[44-49\] \[98--293\] \[177--203\]                                     8.99
         05508              5508.1\*              966             49.28            HE963815             0.735                2                   24 30               \[11-87\] \[44-49\] \[93--287\] \[171--198\] or \[17--93\] \[50-55\] \[99--293\] \[177--204\]           8.49 8.47
        05508_2             5508.2\*              966             49.38            HE963815             0.735                2                   24 30               \[11-87\] \[44-49\] \[93--287\] \[171--198\] or \[17--93\] \[50-55\] \[99--293\] \[177--204\]           8.49 8.47
        22489_1             22489.1\*             978             51.02            HE963818             0.726                2                   23 34               \[22--98\] \[55-60\] \[104--298\] \[182--209\] or \[11-87\] \[44-49\] \[93--287\] \[171--198\]          8.93 8.51
        22489_2             22489.2\*             978             50.82            HE963819             0.727                2                   23 34               \[22--98\] \[55-60\] \[104--298\] \[182--209\] or \[11-87\] \[44-49\] \[93--287\] \[171--198\]          8.93 8.51
         06117               6117\*              1008             47.02            HE963821             0.802                3                   22 32              \[11--91\] \[44-49\] \[97--298\] \[176--208\] or \[21--101\] \[54-59\] \[107--308\] \[186--218\]         5.52 6.16
        17517_1             17517.1\*             972             48.46            HE963822             0.737                2                   23 24               \[12--88\] \[45-50\] \[94--296\] \[171--203\] or \[11-87\] \[44-49\] \[93--295\] \[170--202\]           9.49 9.39
        17517_2             17517.2\*             972             48.56            HE963823             0.739                2                   23 24               \[12--88\] \[45-50\] \[94--296\] \[171--203\] or \[11-87\] \[44-49\] \[93--295\] \[170--202\]           9.52 9.41
        08562_1             08562.1\*             996             47.09            HE963824             0.779                3                   22 28              \[11--91\] \[44-49\] \[97--298\] \[176--208\] or \[17--97\] \[50-55\] \[103--304\] \[182--214\]          9.01 9.03
        08562_4             08562.2\*             996             47.79            HE963825             0.788                3                   22 28              \[11--91\] \[44-49\] \[97--298\] \[176--208\] or \[17--97\] \[50-55\] \[103--304\] \[182--214\]          9.00 9.02

The nucleotide sequences of 28 isoenzymes were submitted to EMBL. "\*" indicates novel. "CAI" = codon adaptation index. If signal sequence predictions gave more than one alternative result, both signal sequence lengths are shown, separated by "/". Disulfide bridges were predicted using both alternatives of the mature protein (EDBCP = Ensemble-based Disulfide Bonding Connectivity Pattern prediction server).

The codon usages of the *HRP* genes and the previously known *A. thaliana* peroxidases were compared in the form of heatmaps (Additional file [1](#S1){ref-type="supplementary-material"}), depicting the fold change of the codon usage frequencies compared to the expected (1/64) frequency (ΔRSCU). The clustering of the isoenzymes according to their codon usage frequencies situated newly discovered isoenzymes with most divergent sequence and gene structure (*HRP*\_3523, *HRP*\_5508, *HRP*\_22489, *HRP*\_17517, *HRP*\_23190) also furthest away from the previously known group C isoenzymes.

Gene structure and phylogenetic analyses
----------------------------------------

Phylogenetic relationships of the HRP isoenzymes are shown in Figure [3](#F3){ref-type="fig"} and Additional file [1](#S1){ref-type="supplementary-material"}. Interestingly, the previously known isoenzymes seem to be closely related to each other, while most of the new isoenzymes discovered in the transcriptome seem to share higher evolutionary distance to them. From the 20 peroxidase gene loci, 15 were confirmed to have three introns by comparing either transcript data or protein sequence data to the verified gDNA sequence. Further four genes (5508, 22489, 17517, 23190) were noted to have only two introns and one gene (3523) no introns (Table [3](#T3){ref-type="table"}). The number of introns correlates with the evolutionary distance so that genes having aberrant intron numbers were situated in separate branches close to each other in the phylogenetic topology. With the information obtained from the reads, no alternative splicing could be shown. According to SignalP, all of the isoenzymes have an N-terminal signal sequence varying in length from 18 amino acids to 31 amino acids. The lengths of the signal sequences are described in Table [3](#T3){ref-type="table"}, and an alignment of the amino acid sequences of the HRP isoenzymes is shown in Additional file [2](#S2){ref-type="supplementary-material"}.

![**Cladogram of all isoenzymes known and discovered during this study.** The dendrogram was cut at a branch length of 0.21 and the resulting sub-trees were colored. All previously described isoenzymes are located in the red and light green trees, respectively; whereas all novel isoenzymes (except for 01805, 22684.1, 22684.2, and 04663) cluster in distinct sub-trees (black, blue, dark green and orange) indicating a larger sequence diversity.](1471-2164-15-227-3){#F3}

Heterologous protein production in *Pichia pastoris*
----------------------------------------------------

Twenty-six *HRP* sequences including allelic variants were codon optimized for *P. pastoris* expression and ordered as synthetic fragments. Twenty-two of them showed activity with at least one of the substrates tested, thus verifying a successful expression in *P. pastoris*. The peroxidase activities of the produced isoenzymes were detected with four substrates having variable assay pH optima. The substrate or pH-specific performance (Table [4](#T4){ref-type="table"}) of the isoenzymes suggests a wide range of possible applications for this versatile group of peroxidases.

###### 

Summary of isoenzyme expression and characterization

   **Contig number**   **Name (\* previously unknown)**   **Activity ABTS**   **Activity TMB**   **Activity guaiacol**   **Activity pyrogallol**
  ------------------- ---------------------------------- ------------------- ------------------ ----------------------- -------------------------
          \-                         C1A                         \+                  \+                   \+                       \+
         15901                       C1B                         \-                  \+                   \-                       \-
         25148                      C1C\*                        \+                  \+                   (+)                      \+
        25148_2                     C1D\*                        \+                  \+                   (+)                      \+
         04627                        C2                         \+                  \+                   \+                       \+
          \-                          C3                         \+                  \+                   \+                       \+
    Manual assembly                 A2A\*                        \+                  \+                   \+                       \+
    Manual assembly                 A2B\*                        \+                  \+                   (+)                      (+)
         04382                        E5                         \+                  \+                   \+                       \+
         01805                  01805\* (B1?)                    \+                  \+                   \-                       \-
         22684                 22684.1\* (B2A?)                  \+                  \+                   \-                       \-
        22684_2                22684.2\* (B2B?)                  \+                  \+                   \-                       \-
         01350                     01350\*                       \+                  \+                   \-                       \-
         02021                     02021\*                       \-                  \-                   \-                       \-
         23190                    23190.1\*                      \-                  \-                   \-                       \-
        23190_2                   23190.2\*                     n.d.                n.d.                 n.d.                     n.d.
         04663                    04663.1\*                      \+                 (+)                   \-                       \-
         06351                     06351\*                       \+                  \+                   \+                       \+
         03523                     03523\*                       \-                  \-                   \-                       \-
         05508                    05508.1\*                      \+                  \+                   \+                       \+
        05508_2                   05508.2\*                     n.d.                n.d.                 n.d.                     n.d.
        22489_1                   22489.1\*                      \+                  \+                   (+)                      \+
        22489_2                   22489.2\*                      \+                  \+                   (+)                      \+
         06117                     06117\*                       \-                  \-                   \-                       \-
        17517_1                   17517.1\*                      \+                  \-                   \-                       \-
        17517_2                   17517.2\*                      \+                  \+                   \+                       \+
        08562_1                   08562.1\*                      \+                  \+                   \-                       \-
        08562_4                   08562.2\*                      \+                  \+                   \-                       \+

Isoenzymes showing obvious peroxidase activity with the assay used are marked with "+". Isoenzymes showing very low but detectable peroxidase activity with the assay used are marked with "(+)". Isoenzymes with no activity detected during an observation period of 2 h are marked with "-". Allelic variants not produced heterologously are marked with n.d. (no data available). Isoenzymes discovered during this study are marked with "\*".

Discussion
==========

Although high throughput sequencing technologies and bioinformatics tools to handle the enormous amount of data generated have been rapidly developing in the recent past, the expressed sequence data of many organisms of wide importance are still not available. Our study demonstrates how NGS technologies can provide a rapid, low-cost basis for the discovery of isoenzymes required for specific industrial, medical or biological applications. Below we discuss the reliability of the approach in identifying and characterizing an important group of isoenzymes, challenges provided by the library generation, sequencing and assembly methods, and suitability of the data obtained from the pipeline for heterologous protein production without laborious manual verification of the sequences.

A normalized cDNA collection originating from multiple *A. rusticana* tissues was sequenced with 454 pyrosequencing technology. In comparison to the other NGS methods, 454 pyrosequencing produces long reads. This is of advantage when sequencing genes coding for isoenzymes and identifying exonic variants. Longer reads increase the probability to uniquely align a given read, which might be problematic with the short reads produced by other NGS methods \[[@B38]\]. Studies in *A. thaliana* suggested that with the high coverage attainable by massively parallel sequencing, all transcripts can be well represented in the sequence data regardless of expression levels \[[@B17]\]. However, the benefits of normalization in non-model species have not been well characterized, and a recent study by Cirulli *et al*. reports low identification rates of exonic SNVs (single nucleotide variations) in non-normalized transcriptome, if the genes of interest are not well-expressed in the source tissue \[[@B38]\]. Since also the genome and transcriptome sizes of *A. rusticana* and the sequence data required to reach also isoenzymes with low expression levels are unknown, the cDNA library sequenced in this study was normalized. Normalization of the cDNA has been reported to be especially important in gene discovery when the cDNA used for sequencing is pooled from many tissues or individuals \[[@B39]-[@B41]\], and to considerably reduce the frequency of abundant transcripts thus increasing the possibility to reach also unique transcripts of isoenzymes with low expression levels. Half a plate run on a GS FLX platform resulted in over 500000 high-quality reads, corresponding to a relatively high average coverage (\>10×) of the assembled 14871 isotigs with a high average length of 1133 bp. Comparable to previous transcriptome sequencing studies \[[@B40],[@B41]\], 88% of the reads could be assembled into contigs. Many of the remaining singletons were of high quality and also represented an important source of information \[[@B18]\]. Singletons could either result from 454 sequencing errors or contaminants from plant parasites \[[@B42]\], they could be caused by over-efficient normalization methods, or simply represent rare transcripts with thin coverage despite the normalized cDNA pool used for sequencing.

The 454 pyrosequencing technique generates relatively long reads including very few technical errors (mainly related to homopolymer runs), and is therefore well-suited for applications such as *de novo* transcriptome sequencing. Although the sequence length achieved by 454 Titanium FLX platform is still clearly shorter than by traditional Sanger sequencing, it has been reported to be adequate for reconstructing full length transcripts \[[@B17]\] and validated to be comparable in accuracy to Sanger sequencing \[[@B43],[@B44]\]. The high average isotig length of 1133 bp and the assembly of 90% (18/20) full length peroxidase genes reached in this study supports the statement of adequate read length and coverage to detect complete transcripts. Only isoenzymes C1A and C3 were not assembled into a contig due to low sequence coverage.

Although the error rates associated with NGS methods have been reported to be low, they could still cause problems in reliable sequence polymorphism detection. The requirement of \>90% match used in this study, combined with a minimal match length of 40 bp was expected to provide a very high number of contigs without collapsing and joining similar isoenzymes into a single contig \[[@B18]\]. However, the isoenzyme sequences were known to be partially almost identical. To validate the assembly and ORF (open reading frame) prediction correctness, and the existence of allelic variants, the isoenzyme sequences were amplified from genomic DNA. The combination of transcriptome sequencing and Sanger sequencing of amplified genomic DNA revealed 38 variable positions in the coding sequences of 11 genes. For nine positions thereof (in four transcripts) the corresponding nucleotide could not be found in the genomic DNA. Manual confirmation of the transcriptome reads revealed that the coverage of all positions was high and that the reads agreed. Closer analysis of the positions also ruled out false variations caused by too high coverage. At relatively low or moderate levels of coverage, sequencing mistakes are not pushed over. Too high sequencing depth increases noise and could have introduced, without very strict quality control, false variations \[[@B38]\]. Therefore, sequencing errors could be excluded as the source of the differences. The differences could be caused either by mismatches introduced by the reverse transcriptase enzyme during library generation, as a result of RNA editing events \[[@B34],[@B35]\] or reflect the natural variation between individual plants and plant parts. The general error rate of the sequencing technique was noted to be very low, but isoenzymes coded by less common allelic variants would have been missed without manual sequence verification. Twenty-six of the resulting coding sequences, including allelic variants missed in the transcriptome sequencing and assembly processes, were codon optimized, synthesized and transformed to *P. pastoris* for expression.

This study reveals that for the coverage of all isoenzymes including allelic variants represented by the cDNA library sequenced, manual work to verify the resulting transcript sequences cannot be avoided. However, the allelic variants represent only a minor part of the newly discovered enzymes. The quality of the sequences is very high and differences to genomic DNA minimal, confirming that the enzyme discovery method described in this study for high-throughput applications would not necessarily require manual verification of the sequencing by laborious Sanger sequencing of amplified genomic DNA. However, manual curation of the contigs of interest and splitting of the data in contigs with clear assembly conflicts can be done and could be worthwhile, as especially the additive effects of amino acid changes in collapsed contigs could cause problematic changes in the enzyme structure. The primary enzyme discovery pipeline utilized in this study provides a functional approach to find proteins of interest for heterologous production, giving an example of an affordable standardized sequencing project. Despite the large amount of useful data produced by the NGS approach, our study showed that sequence confirmation and data validation should not be neglected.

For the heterologous secretory production of single isoenzymes in *P. pastoris*, the codon usages of the coding sequences were optimized for efficient translation, and fragments corresponding to the predicted mature isoenzymes were produced synthetically. When the signal peptide prediction with SignalP led to two alternative signal peptide junctions, the mature peptides corresponding to the longer signal peptide variants were ordered as synthetic fragments, and the signal peptide variants with shorter mature peptide were successfully amplified via PCR. Correct cloning of all genes into *P. pastoris* expression vectors was verified by Sanger sequencing. Since the used *P. pastoris* expression vector already contains the signal sequence of the *S. cerevisiae* mating factor α, the isoenzymes were produced without the predicted natural signal sequences. Twenty-two isoenzymes showed peroxidase activity with at least one of the substrates used. All activities were measured with four different assays over a pH range from 4.5 to 7. Interestingly, for each assay (Table [4](#T4){ref-type="table"}) a different isoenzyme showed the highest activity thus suggesting variable substrate specificities or pH optima. This observation further emphasizes the importance of the availability of a large group of individually produced pure isoenzymes to be able to comprehensively respond to the need of variable performance parameters including substrate specificity, activity, stability, and operating pH optimum. Four of the isoenzymes did not show peroxidase activity with the assays used. This could either be due to very low yields of active enzyme, totally inactive enzyme or unsuitable assay conditions.

The isoenzyme C1A was reported to be the most abundant isoenzyme in *A. rusticana*\[[@B1]\] and was thus expected to be found in the transcriptome. However, only two raw reads covering a minor part of the coding sequence could be detected. This might either suggest over-normalization of the cDNA library decreasing the total counts of the putatively most abundant transcripts to almost zero \[[@B45]\], or happen due to naturally occurring genetic variation with phenotypic correlation to adaptations to natural environments ranging from pathogens, light conditions or abiotic stress to a variety of other environmental perturbations \[[@B8],[@B46],[@B47]\]. Although mRNA originating from all available plant parts was used to reach genes activated at diverse stages of *A. rusticana* growth, some developmental stages were not present and the absence of certain isoenzymes due to missing tissues cannot be ruled out. This finding might illustrate the high variance of *HRP* expression in *A. rusticana* plants and consequently the variance in the commercial HRP preparations, thus underlining the clear need for a reliable heterologous expression system that enables a consistent isoenzyme quality.

Peroxidase isoenzymes have been suggested to have multiple roles in the plant and thus also be variedly expressed depending on both biotic and abiotic factors \[[@B8]\]. To roughly estimate the expression levels of the newly discovered peroxidase genes, their GC contents and codon usages were assessed. Genes that are highly expressed have been suggested to possess a higher GC content and a more biased codon usage than genes with low expression levels \[[@B48]\]. A majority of both eukaryotic and prokaryotic species with large population sizes have been reported to have non-random codon usage mainly due to Darwinian selection between synonymous codons \[[@B49]-[@B51]\]. Highly expressed genes have been reported to use a restricted set of codons to ensure optimal translational efficiency \[[@B52],[@B53]\]. In addition to gene expression levels, GC content has also been connected to gene regulation \[[@B54]-[@B57]\] and correlated with genomic features including methylation pattern \[[@B58]\], short intron length \[[@B59]\] and gene density \[[@B60]\] thus suggesting possible functional relevance. The codon usages and GC contents calculated using the verified coding sequences of the isoenzymes are described in Table [3](#T3){ref-type="table"}. As expected, large variation between isoenzymes exist. These findings could suggest a spatial and temporal distribution of the isoenzymes in cellular processes \[[@B61]\].

Phylogenetic relationships of the HRP isoenzymes are shown in Figure [3](#F3){ref-type="fig"}. Whereas the previously known isoenzymes are closely related to each other, most of the new isoenzymes discovered in the transcriptome share higher evolutionary distance to the previously known HRPs. BLASTX analysis to the peroxidases of *A. thaliana* (Additional file [3](#S3){ref-type="supplementary-material"}) revealed that the *A. rusticana* peroxidases share 81% to 95% sequence similarity to the most similar isoenzyme of *A. thaliana*. Evolutionary distance does not necessarily correlate with altered substrate specificity, specific activity or optimal reaction conditions, but the discovery of new evolutionary branches with higher structural diversity does offer optimal conditions for the generation of an enzyme assortment with diverse properties for a wide variety of biomedical and industrial applications.

A combination of cDNA sequencing and gDNA verification in this study also provided valuable information of the intron-exon boundaries of the *HRP* genes. The number of exons in the isoenzymes was noted to vary from one to four, corresponding to zero to three introns. A large majority (75% of the peroxidase loci) of the isoenzymes were found to have four exons and three introns. Intron numbers have been reported to be highly conserved, but total intron length (total sum of the sizes of all introns within a gene) rather correlated to the GC-content of the gene \[[@B62]\]. Thus, intron number could be informative in terms of evolutionary origin and distance of the enzyme. In this study, intron numbers were found to correlate with the phylogenetic relationships of the amino acid sequences. Contigs with an unusual number of introns (none or two, Table [3](#T3){ref-type="table"}) were situated in close proximity to each other furthest away from the previously known isoenzymes and clustered together when comparing the codon usage frequencies. With the information obtained from the reads, no alternative splicing could be observed.

The well-characterized isoenzyme HRP C1A has been reported to have a signal peptide consisting of 30 amino acids, and a carboxy-terminal extension suggested to target the protein to the vacuoles \[[@B63]\]. Also other known isoenzymes of the group C (C1B, C1C, C2, and C3) have been reported to have signal peptides varying in length from 9 amino acids (C1C) to 29 amino acids (C3). By observing the alignment of all previously known and newly discovered isoenzymes (Additional file [2](#S2){ref-type="supplementary-material"}), existence of signal sequences also in other previously known and most of the new isoenzymes seemed very probable. According to the signal sequence prediction (SignalP) performed, all isoenzymes seem to have a signal sequence varying in length from 18 to 31 amino acids (Table [3](#T3){ref-type="table"}, Additional file [2](#S2){ref-type="supplementary-material"}). Isoenzyme C1C, previously reported to have a signal sequence of nine amino acids, was predicted to have - better corresponding to the sequences of the very closely related isoenzymes C1A and C1B - a signal sequence of 29 amino acids. In the case of unclear signal sequence prediction with more than one option for the length of the signal peptide, both forms were taken into consideration when planning the constructs for enzyme production in *P. pastoris*.

Conclusions
===========

To facilitate the possibilities for heterologous expression and isoenzyme characterization, we have elucidated the nucleotide sequences of 28 horseradish peroxidase isoenzymes by using the data obtained from *A. rusticana* 454 transcriptome sequence analysis with manual verification of PCR amplified genomic DNA. Although studies including transcriptome analysis of non-model species have become increasingly popular since the emergence of the NGS technologies, methods for the utilization of the 454 technology for the purpose of isoenzyme discovery in non-model plant species have not been established. In this project, transcriptome sequencing reads are further processed with alternative assemblies and manual sequence verification to determine the nucleotide sequences of all HRP isoenzymes. This study does not only contribute a set of transcripts, which can be used for marker development and genomic studies to understand agriculturally important traits in *A. rusticana*, but also provides valuable information of the peroxidase gene structure. Twenty-two of the verified isoenzymes have been produced in a form that was found active towards the tested substrates in *P. pastoris* utilizing a new *P. pastoris* expression platform \[[@B64]\], validating the success of the approach and providing first insights into the versatility of this large group of isoenzymes discovered.

Methods
=======

Plant specimens, RNA extraction and quality analysis
----------------------------------------------------

Wild horseradish (*A. rusticana*) roots were purchased from local farmers and grown in the laboratory to obtain fresh roots, sprouts, stems and leaves. Tissues were collected in aliquots, frozen in liquid nitrogen and stored at -80°C. Total RNA from all available plant parts was isolated using RNaqueous kit (Applied Biosystems/Ambion, Austin, TX, USA) according to the manufacturer's recommendations. Quality assessment to ensure RNA integrity was performed with Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

Normalized cDNA library construction and sequencing
---------------------------------------------------

Transcriptome sequence was obtained by a commercial service from LGC Genomics (Berlin, Germany). The methods used are roughly summarized as follows: mRNA was purified from total RNA using mRNA-ONLY™ Eukaryotic mRNA Isolation Kit (Epicentre, Madison, WI, USA). One μg of mRNA was used for first-strand cDNA synthesis and amplification according to the Mint-Universal cDNA Synthesis Kit user manual (Evrogen, Moscow, Russia), followed by a normalization reaction using the Trimmer Kit (Evrogen). Normalized material was re-amplified, digested (*Sfi*I), size-selected (\>800 bp, LMP agarose), purified (Qiagen, Hilden, Germany) and ligated to pDNR-lib vector (Clontech, Saint-Germain-en-Laye, France) using the Fast Ligation Kit (New England Biolabs, Ipswich, MA, USA). The desalted ligation was used to transform NEB10b competent cells (New England Biolabs).

Roughly a million clones were plated on LB (lysogeny broth) + chloramphenicol (Cm) plates, scraped off the plates and stored as glycerol stocks at -70°C. Plasmid DNA was prepared using standard methods (Qiagen, Hilden, Germany), and digested with *Sfi*I. cDNA inserts were gel-purified (LMP-Agarose/MinElute Gel Extraction Kit, Qiagen) and ligated to high-molecular-weight DNA using a proprietary *Sfi*I-linker.

Library generation for the 454 FLX sequencing was carried out according to standard protocols (Roche/454 Life Sciences, Branford, CT 06405, USA). In short, the concatenated inserts were sheared randomly by nebulization to fragments ranging in size from 400 bp to 900 bp. These fragments were end polished and the 454 A and B adaptors that are required for the emulsion PCR and sequencing were ligated to the ends of the fragments. The resulting fragment library was sequenced on half a picotiterplate on the GS FLX using the Roche/454 Titanium chemistry. Sequence data can be accessed via the EMBL-EBI European Nucleotide Archive under the study accession number PRJEB5793.

Assembly of the sequence reads to transcripts
---------------------------------------------

Raw reads produced by the pyrosequencing process were screened for the *Sfi*I-linker that was used for concatenation and the linker sequences were clipped from the reads. Poly A/T sequences were mostly (\~90%) removed with the linker. High-quality reads were selected using Newbler sequence filtering at default settings. The clipped, quality controlled reads were assembled into individual isotigs using the Roche/454 Newbler software (454 Life Sciences Corporation, version 2.5.3) with default settings (minimum read length 20, duplicate reads excluded, expected depth 0, seed step 12, seed length 16, seed count 1, minimum overlap length 40 bp, minimum overlap identity 90%, alignment identity score 2, alignment difference score -3).

Discovery of peroxidases in the assembled contigs
-------------------------------------------------

The PSSM (position-specific scoring matrix) corresponding to known horseradish peroxidases was obtained from NCBI\'s Conserved Domain Database (cd00693, CDD v3.01) \[[@B65]\]. It was used in a tblastn search (e-value cutoff 1e-5) \[[@B66]\] on the assembled contigs to yield a preliminary set of *HRP* candidate sequences. This set was refined by filtering sequences whose translation mapped back to a domain different to the original profile PSSM in an rpsblast classification of the entire CDD or had a bit score lower than the NCBI-specified threshold for a specific domain match. The read composition of the refined set of contigs was manually reviewed using SeqMan (DNASTAR, Madison, Wisconsin, USA). Isotigs where two apparently different variants were assembled into one contig by the assembler were split. Protein coding regions were extended using read information if two or more reads contained the same sequence.

Genome walking and manual verification of horseradish peroxidase sequences
--------------------------------------------------------------------------

The sequences of the identified peroxidase genes were manually verified by Sanger sequencing of PCR amplified genomic fragments using Phusion^TM^ High-Fidelity Polymerase (Finnzymes Oy, Espoo, Finland) and primers listed in Additional file [4](#S4){ref-type="supplementary-material"}. If no flanking regions were available for primer design, genome walking \[[@B67]\] was utilized to clarify and complete the sequences of the C-- and N-termini. Therefore, 2 μg aliquots of genomic DNA were singly digested with *Bsp*143I, *Hin*dIII, *Psu*I (Fermentas, St. Leon-Rot, Germany), *Bsa*WI (New England Biolabs GmbH, Frankfurt am Main, Germany) or *Xho*II (Promega GmbH, Madison, WI, USA) in order to get fragments of 1 - 5 kb size. The digestion was stopped by heat inactivation of the enzymes, the fragmented DNA was precipitated with ethanol and the pellet was dissolved in 30 μL of distilled water. An adaptor was created by annealing adaptor strand 1 (5\'-GTAATACGACTCACTATAGGGCACGCGTGGTCGACGGCCCGGGCTGGT-3\') either to adaptor strand 2.a (3\'-TCCCCGACCACTAG-5\') for *Bsp*143I/*Psu*I-/*Xho*II-digested DNA, 2.b (3\'-TCCCCGACCATTAA-5\') for *Bsa*WI-digested DNA or 2.c (3\'-TCCCCGACCATCGA-5\') for *Hin*dIII- digested DNA.

In the annealing reaction, adaptor strand 1 was mixed in 1:1 molar ratio with adaptor strand 2.a/2.b/2.c (i.e. 13.7 μL of 100-μM adaptor strand 1 + 4.0 μL of 100 μM adaptor strand 2) and heated to 95 °C for 5 min. To anneal the two strands to a functional adaptor molecule, the mixture was allowed to slowly cool down to room temperature. The three differently annealed adaptors (1 + 2a, 1 + 2b, 1 + 2c) were ligated for three hours at room temperature with T4 DNA Ligase (Fermentas) to the digested DNA fragments, considering the specific 5\' overhangs that have been created by the respective restriction enzymes. The ligation reaction was stopped by incubation for 5 min at 70°C and 70 μL of TE (Tris-EDTA) buffer was added. Two gene-specific primers and two adaptor primers were designed. The gene-specific primers were designed to bind approximately 100 bp from the end of the known sequence, considering that no restriction site of the restriction enzymes used for DNA digestion laid between the primer-binding site and the end of the known sequence. AdaptorPrimer1 (5\'-GTAATACGACTCACTATAGGGC-3\') and GeneSpecificPrimer1 were used as a primer pair for a first PCR with 1 μL of the DNA + adaptor ligation product as template DNA. One μL of the first PCR mix was used as template for a second PCR with AdaptorPrimer2 (5\'-ACTATAGGGCACGCGTGGT-3\') and GeneSpecificPrimer2. This second primer pair was designed to bind within the first PCR product. Both PCR steps were performed with an elongation time of 50 seconds.

A gene-specific DNA fragment as product from the second PCR was isolated from a preparative agarose gel, purified (SV DNA extraction kit, Promega) and sent to Sanger sequencing (LGC Genomics GmbH, Berlin, Germany), using AdaptorPrimer2 and the corresponding GeneSpecificPrimer2. If unspecific primer binding or nucleotide polymorphisms were suspected, the PCR products were cloned into the pJET 1.2blunt vector (GeneJet cloning kit, Promega), transformed to *E. coli* Top10 F\' and plasmids from single colonies were isolated for sequencing to ensure the read consisted of only one allele.

The resulting gene sequences were submitted to EMBL \[[@B68]\] under accession numbers HE963800-HE963825 (Table [3](#T3){ref-type="table"}).

Codon usage, GC content, isoelectric point, signal sequence prediction, disulfide bridge prediction and phylogenetic analyses of the horseradish peroxidase isoenzymes
----------------------------------------------------------------------------------------------------------------------------------------------------------------------

Codon usages and GC contents of the HRP isoenzymes were analyzed using CAIcal \[[@B69],[@B70]\] and Mega5 \[[@B71],[@B72]\]. The sequences of the *A. thaliana* Class III peroxidases were downloaded from TAIR \[[@B73]\]. Sequences were aligned with ClustalW2 \[[@B74],[@B75]\]. The theoretical isoelectric points (pI) were calculated with ExPASy Compute pI/Mw tool \[[@B76]\]. Disulfide bridges were predicted with EDBCP tool \[[@B77]-[@B79]\]. Phylogenetic analyses were performed with CLC Main Workbench 6.6.2 (CLC bio, Aarhus Denmark) \[[@B80]\] and Mega5. The pylogentetic tree was generated with the \"Create Tree\" function of CLC Main Workbench using the UPGMA algorithm with a bootstrap analysis of 100 replicates, based on an alignment of the HRP amino acid sequences using the \"Create Alignment\" function with the following settings: Gap open cost: 10.0, Gap extension cost: 1.0, End gap cost: as any other, Alignment: Very accurate (slow). Signal sequences were predicted using SignalP 3.0 \[[@B81],[@B82]\].

Gene synthesis and heterologous expression in *Pichia pastoris*
---------------------------------------------------------------

The codon usages of 14 isoenzymes were optimized for the expression in *P. pastoris* using a novel algorithm (DNA2.0, Menlo Park, CA, USA, Mellitzer *et al.* manuscript in preparation). Further twelve isoenzymes including allelic variants were optimized using GeneDesigner 1.1.4.1 (DNA2.0) in accordance to the *P. pastoris* codon usage described by Abad *et al.*\[[@B83]\]. Signal sequence variants were generated by PCR amplification and all *HRP* genes were cloned into the shuttle vector pPpT4_alpha_S of a newly generated open source expression platform \[[@B64]\]. The vector pPpT4_alpha_S is a basic low-copy (1--5 copies/genome), zeocin™ resistance based expression vector for efficient secretory expression of heterologous proteins. Sanger sequencing of the plasmids verified successful cloning into the right frame. The linearized expression cassettes were transformed into *P. pastoris* wild-type CBS7435 based mut^s^ strain using standard protocols \[[@B84]\], and selected on zeocin™-containing plates. From each gene, 88 clones were picked to 96-well deep-well plates for cultivation and high-throughput screening of peroxidase activity. Two of the well expressing clones of each isoenzyme were streaked out to single colonies. Four single colonies of each clone were used for re-screening to estimate the reproducibility of the results. All media compositions and cultivation protocols used in this study were as previously described by \[[@B85]\]. Minimal media BMD1% (buffered minimal media with 1% dextrose) was supplemented with 5 mM ferrous sulfate heptahydrate (Sigma-Aldrich Handels Gmbh, Vienna, Austria) to ensure sufficient iron supply for heme biosynthesis.

Peroxidase assays
-----------------

ABTS (2,2'-azino-bis(3-ethylbenzthiazoline-6-sulfonic acid), TMB (3,3',5'5-tetramethyl benzidine), pyrogallol (1,2,3-trihydroxybenzene) and guaiacol (2-methoxyphenol) assays were used to detect peroxidase activity essentially as described in \[[@B3],[@B86],[@B87]\]. ABTS assays were performed in 50 mM sodium acetate buffer pH 4.5 with 1 mM ABTS and 0.0026% (v/v) H~2~O~2.~ For the TMB stock solution, TMB was dissolved in DMSO to a concentration of 4.16 mM. For the assay solution, TMB stock solution and 30% (v/v) H~2~O~2~ were diluted with 20 mM citrate buffer pH 5.5 to final concentrations of 0.416 mM and 0.006% (v/v), respectively. Guaiacol assays were performed in 10 mM sodium phosphate buffer pH 7.0 with 5 mM guaiacol and 0.0009% (v/v) H~2~O~2~. For the pyrogallol assay solution, pyrogallol (Sigma-Aldrich Handels Gmbh, Vienna, Austria) was dissolved in 10 mM potassium phosphate buffer pH6.0 containing 0.027% (v/v) H~2~O~2~ to a concentration of 45 mM. For all assays, 15 μl cultivation supernatant was mixed with 140 μl of the assay solution in a flat-bottom 96-well microtiterplate (Greiner Bio-One GmbH, Frickenhausen, Germany). The reaction kinetics were followed with Spectramax Plus^384^ spectrophotometer and SoftMax® Pro software (Molecular Devices, LLC) for 3--5 min at wavelengths 405 nm (ABTS), 650 nm (TMB), 470 nm (guaiacol) and 420 nm (pyrogallol). Enzyme activity was calculated using only time points fitting to linear increase of the absorbance (ΔmAU min^-1^).

Availability of supporting data
-------------------------------

*A. rusticana* transcriptome sequencing data is available via EMBL-EBI\'s European Nucleotide Archive (ENA) under the study accession number PRJEB5793 (<http://www.ebi.ac.uk/ena/data/view/PRJEB5793>). Nucleotide sequences of the novel identified HRPs have been deposited into ENA as well: accession numbers HE963800-HE963825 (<http://www.ebi.ac.uk/ena/data/view/HE963800-HE963825>). All other supporting data are included as additional files with this manuscript.

Abbreviations
=============

ABTS: 2,2'-azino-bis(3-ethylbenzthiazoline-6-sulfonic acid); AU: Absorption unit; BLAST: Basic local alignment search tool; BLASTX: Search protein databases using a translated nucleotide query; BMD: Buffered minimal media with dextrose; bp: Base pair; CDD: Conserved domains database; cDNA: Complementary DNA; DMSO: Dimethyl sulfoxide; EST: Expressed sequence tag; gDNA: Genomic DNA; HRP: Horseradish peroxidase; LB: Lysogenous broth; mRNA: messenger RNA; NGS: New generation sequencing; ORF: Open reading frame; PSSM: Position-specific scoring matrix; RIN: RNA integrity number; RSCU: Relative synonymous codon usage; SNV: Single nucleotide variation; tBLASTn: search translated nucleotide databases with a protein query; TE: tris-EDTA; TMB: 3,3',5'5-tetramethyl benzidine.

Competing interests
===================

AG, FWK, and LN report a patent application (No. 11185833.8 - 2403), owned by Graz University of Technology, which relates to the novel recombinant horseradish peroxidase isoenzymes with improved technological properties described in this work and the authors do not have any objection on the patent right of authorship and ownership. All other authors declare that they have no competing interests.

Authors' contributions
======================

AG, GGT and LN conceived the study and designed the experiments. LN collected the samples, isolated RNA, manually verified the read composition of the peroxidase contigs, and participated in the bioinformatic analysis. MS and GGT developed a pipeline to discover peroxidase isoenzymes in the contigs and participated in the bioinformatic analysis. LN and FWK verified the genomic sequences, planned the synthetic genes and performed protein expression in *Pichia pastoris*. AG initiated the project and contributed to the evaluation of the results. LN and GGT drafted the manuscript and all other authors revised it and approved the final version.

Supplementary Material
======================

###### Additional file 1

**Heatmaps of of the changes in the relative synonymous codon usages (ΔRSCU) of A) all the HRP isoenzymes verified in this study and B) the known*A. thaliana*peroxidases.** Each column represents one codon indicated along the bottom, each row one isoenzyme marked to the right side of the row. Isoenzymes are clustered by their codon usage similarity. Green cells correspond to underrepresented codons, red cells to overrepresented codons. Missing codons are marked with a grey cell.

###### 

Click here for file

###### Additional file 2

Alignment of the amino acid sequences of the HRP isoenzymes.

###### 

Click here for file

###### Additional file 3

**BLASTP of respective full length HRP amino acid sequence against non-redundant protein sequences (nr) database with*A. thaliana*(taxid:3702) as organism.**

###### 

Click here for file

###### Additional file 4

Primers used in this study.

###### 

Click here for file

Acknowledgements
================

We would like to thank Christian Schmid and Christoph Fischer for expert technical assistance. This work was supported by the Austrian Science Fund 'Fonds zur Förderung der wissenschaftlichen Forschung' \[FWF Grant W901 DK Molecular Enzymology\] and in part by the Austrian Ministry of Science and Research GEN-AU project BIN \[FFG Grant 820962\]. The ACIB contribution was supported by the Austrian BMWFJ, BMVIT, SFG, Standortagentur Tirol and ZIT through the Austrian FFG-COMET- Funding Program \[FFG Grant 824186\]. Funding for open access charge: Fonds zur Förderung der wissenschaftlichen Forschung.
