Barley plants of Hordeum vulgare cv. Golden Promise were used. Grains were germinated in soil and grown at 22°C under a 16 h light/8 h darkness photoperiod for 7 days in Sanyo MLR-350-H chambers. Barley transgenic lines overexpressing or silencing the HvPap-1 gene (OE Pap1 and KD Pap1, respectively) were generated using Agrobacterium-mediated gene transfer and haploid technology in collaboration with Dr. Jochen Kumlehn’s group (IPK-Gatersleben, Germany) as described in Diaz-Mendoza et al. (2016). The homozygous transgenic barley lines for the inserted constructions were validated by molecular and biochemical characterization of HvPap-1, including copy number, and mRNA and protein expression patterns as described in Diaz-Mendoza et al. (2016). Transgenic barley lines were grown at the same conditions described above.

The flagellin (flg22) peptide (amino acid sequence QRLSTGSRINSAKDDAAGLQIA) was synthesized by AnaSpec laboratories at a purity level of ≥95% (v/v). Stocks were prepared by dissolving the peptide in H₂O at a concentration of 10 mM, stored at -20°C, and subsequently diluted in H₂O to a 5 μM concentration for experiments. Chitosan solutions were prepared by dissolving hydrolyzed chitosan purified from crab shells (Sigma–Aldrich) in 0.02% (w/v) acetic acid at a concentration of 100 μg/mL. All elicitor solutions were applied as a foliar spray over wild-type 7-day-old barley leaves. Three pots of two plants each were used per treatment and three independent experiments were performed. The solvent used to solubilize elicitors (H₂O or 0.02% (v/v) acetic acid, for fgl22 or chitosan, respectively) was used for mock treatments. Plants were further incubated at the same growth conditions described earlier. Barley leaves were monitored at different time points. Finally, leaves were harvested after 24 h of elicitor treatments. These samples were frozen into liquid nitrogen and stored at -80°C for further analysis.

As example of biotic stress mediated by pathogens, 7-day-old barley plants of wild-type and transgenic lines overexpressing or silencing the HvPap-1 gene were infected with the fungus M. oryzae. For infection assays, the fungal isolate used in this study was the M. oryzae wild-type strain Guy11 (Leung et al., 1988), kindly provided by Dr. Sesma, CBGP-UPM-INIA, Madrid. The growth, maintenance and media composition of M. oryzae were as previously described in Tucker et al. (2010). Infection assays were performed in whole plant leaves by spray inoculations using an airbrush nebulizer compressor as described in Sesma and Osbourn (2004). Seven plants were used per treatment and three independent experiments were performed. Each pot was sprayed with 1 mL suspension of 10⁵ conidia/mL counted in a Neubauer counting chamber, in 0.25% (v/v) gelatin, or just gelatin in the case of control treatments. Plants were covered by a plastic bag to avoid conidia dispersion and the same confinement system was applied to the controls. The plants were further incubated at 25°C, 65% RH, under a 16 h light/8 h darkness photoperiod. Barley leaves were monitored at different time points to score disease symptoms. Leaves were harvested after 3 and 7 days post inoculation (dpi). This material was imaged and scanned, or frozen in liquid nitrogen and stored at -80°C for further analysis.

A colony of the two-spotted spider mite T. urticae London strain (Acari: Tetranychidae), provided by Dr. Miodrag Grbic (UWO, Canada), was maintained on beans in a Sanyo MLR-350-H growth chamber at 25°C under a 16 h light/8 h darkness photoperiod. This colony was transferred to barley where it was maintained under the same conditions for more than 30 generations to ensure host adaptation To induce biotic stresses mediated by pest attack, 7-day-old barley plants, wild-type and transgenic lines overexpressing or silencing the HvPap-1 gene, were infested with 20 barley-adapted adults of T. urticae per plant. Seven plants were used per treatment and three independent experiments were performed. Barley plants were confined in independent pots with a plastic cylinder covered on top by nylon nets to avoid dispersion of mites and the same isolation system was applied to control plants. Plants were further incubated at 25°C under a 16 h light/8 h darkness photoperiod. Barley leaf damage was monitored at different time points after spider mite feeding. Leaves were harvested after 7 and 14 days of mite treatment. Samples were imaged and scanned, or frozen into liquid nitrogen and stored at -80°C for further analysis.

Magnaporthe oryzae and Tetranychus urticae lesions observed on barley leaves were scanned using a scanner hp scanjet (HP Scanjet 5590 Digital Flatbed Scanner) and foliar damage on transformed and non-transformed lines after treatments was analyzed. Damaged leaf surface area was measured using the Fiji-ImageJ software (Schindelin et al., 2012) and was quantified as mm² of injured area. Seven replicates per line were analyzed and three independent experiments were performed.

For real-time RT-qPCR studies, total RNA was extracted from frozen barley leaves by the phenol/chloroform method, followed by precipitation with 8 M LiCl (Oñate-Sánchez and Vicente-Carbajosa, 2008) and digested with DNase (Promega). cDNAs were synthesized from 2 μg of RNA using the RevertAid H Minus First Strand cDNA Synthesis Kit (Thermo Scientific) following the manufacturer’s instructions. RT-qPCR analyses were performed for triplicated samples by means of a CFX96 Real-time system (BioRad) using SYBR Green (Roche) as a detection system. Protease mRNA quantification was expressed as 2^∧-dCt (Livak and Schmittgen, 2001) and normalized to barley cyclophilin (HvCycl gene) mRNA levels (Hua et al., 2015). Expression levels of the M. oryzae small subunit of ribosomal RNA (Mo28S-rRNA) (Marcel et al., 2010) and T. urticae Ribosomal Protein 49 (TuRP49) (Morales et al., 2016) were quantified as 2^∧-dCt by subtracting the dC_t value of the barley cyclophilin from the dC_t value of the fungal or mite probes to normalize for the amount of barley tissue present in each sample. Fold change values were expressed using the ddC_t method (Livak and Schmittgen, 2001) relative to the normalized expression of the same genes in non-treated plants. The primers used are shown in Supplementary Table 1. Efficiency of the primers was determined based on the slope of a standard curve and was between 95 and 105% for all primer pairs tested.

Total protein was extracted from treated and control barley plants by grinding leaf tissue in liquid nitrogen before the addition of 500 μl of extraction buffer (150 mM NaCl, 50 mM sodium phosphate, pH 6 and 2 mM EDTA). After centrifugation at 15,600 × g for 10 min at 4°C, the supernatant was used for protein quantification according to the method of Bradford (1976), with the bovine serum albumin as a standard. After separation on SDS-polyacrylamide gels (12–15%, w/v) according to Laemmli (1970), proteins were electro-transferred onto nitrocellulose membrane (GE Healthcare) and blocked in PBS (137 mM NaCl, 27 mM KCl, 10 mM Na₂HPO₄, 2 mM KH₂PO₄ pH 7.4) antisera buffer containing 5% (w/v) powdered skim milk, for 1 h. Immunoblotting was performed as described in Diaz-Mendoza et al. (2016), with anti-peptide polyclonal antibodies specifically selected against each protease whose sequences are indicate in Supplementary Table 2. All protease antibodies were produced in rabbits by Pineda Antibody Services. Polyclonal antibody against the Large Subunit of Rubisco (anti-LSR) was supplied by Agrisera. Optimal dilutions of primary antibodies were used according to Diaz-Mendoza et al. (2016).

Total RNA was isolated from barley leaves as described above. Using poly-T oligo-attached magnetic beads, mRNAs were purified from the total RNA. Then, the mRNAs were fragmented and cDNA was synthesized using random hexamer-primers, DNA polymerase I and RNase H. The double-stranded cDNAs were purified with magnetic beads and ligated to adaptors for Illumina sequencing. The quality and quantity of the library was verified using an Agilent 2100 Bioanalyzer and an ABI StepOnePlus Real-Time PCR system, respectively. The cDNA libraries were sequenced using the Illumina HiSeq2000 platform by the Beijing Genomics Institute (BGI). More than 10 M single-end reads were obtained for each sample (three biological replicates).

Raw reads in fastq format were firstly filtered and reads with adaptor sequences and low quality reads were removed. The gene and genome sequences of H. vulgare retrieved from the PGSB/MIPS PlantsDB website¹ (Nussbaumer et al., 2013) were used as the reference databases (International Barley Genome Sequencing Consortium, 2012). Two different approaches were performed: (i) all the clean reads were mapped to the reference genome using SOAPAligner/SOAP2 (Li et al., 2009). The transcripts abundance was normalized by the RPKM (reads per kilobase of exon per million reads) algorithm (Mortazavi et al., 2008). Differentially expressed genes (DEGs) between groups were obtained using the NOISeq method (Tarazona et al., 2011) with a log2Ratio (fold change) higher than 1 and a probability of differential expression higher than 0.8; and (ii) all the clean reads were pseudoaligned to the reference High Confidence genes using Kallisto (Bray et al., 2016). The transcript abundance was quantified as TPM (transcripts per million) and 100 bootstrap samples were performed. DEGs between groups were obtained using the Wald test of the Sleuth method (Pimentel et al., 2016) with a b ratio (bias) higher than 1 and a q-value (false positive probability) lower than 0.001. Venn diagrams were created by the Venny 2.1 utility² (Oliveros, 2007–2015) and heatmaps were performed by the shinyHeatmaply application³. Gene enrichment analyses were performed with the Fischer’s exact test using topGO package in R⁴ and the GO file retrieved from the PGSB/MIPS PlantsDB website.

For statistical analysis, the expression of individual C1A CysProt genes in WT background was compared between control and elicitor/biotic stimuli samples. Foliar damage and individual C1A CysProt expression data were compared among WT and transformed lines in control conditions or after biotic treatment. All analyses were performed by one-way ANOVA, followed by t-student or Student Newman-Keuls (SNK) tests using the statistic software R Project (v.3.1.2) package. In figures where pairwise samples were compared, one or two asterisks indicated significant differences (t-test, P < 0.01 and P < 0.001, respectively). In figures where more than two samples were compared, different letters indicated significant differences (SNK test, P < 0.01).

To check the transcriptional responses of C1A CysProt against biotic stresses, three cathepsin L-like (HvPap-4, -6 and -16), one cathepsin F-like (HvPap-1), one cathepsin H-like (HvPap-12) and one cathepsin B-like (HvPap-19) genes were selected, all of them previously studied under severe senescence processes induced by continuous darkness and nitrogen deprivation (Velasco-Arroyo et al., 2016).

For elicitor treatments, flg22, a 22-amino acid sequence of the conserved N-terminal part of bacterial flagellin, was used. This peptide is known to activate plant defense mechanisms. The effect of chitosan, which is a structural element in the exoskeleton of crustaceans and cell walls of fungi, on protease expression was also tested. A molecular and biochemical characterization of selected C1A CysProt members was performed including mRNA quantification (Figures 1A, B) and protein accumulation (Figure 1C) in response to elicitor treatments. The results, expressed as mRNA levels normalized to the constitutively active barley cyclophilin gene, revealed that after 24 h of flg22 treatment the cathepsin F, HvPap-1 gene, was significantly induced in treated leaf samples (Figure 1A). After 24 h of chitosan treatment the most abundant transcripts were HvPap-1 and HvPap-19, encoding cathepsin F- and B-like, respectively (Figure 1B). Besides, differences in expression were detected after flg22 treatment in the case of HvPap-16 and after chitosan treatment for HvPap-6, both proteases belonging to cathepsin-L group. These results were companied by immunoblot assays using antibodies against specific peptides of HvPap-1, HvPap-6, HvPap-16 and HvPap-19 CysProt (Figure 1C). In some cases the results point out a direct correlation between transcript and protein accumulation patterns, as for HvPap-1 after chitosan treatment. Bands increased their signal in HvPap-19 not only after chitosan treatment but also after flg22 elicitation. A weak induction was observed for the cathepsin L-like proteins HvPap-6 and -16 in chitosan treated plants, while no remarkable differences were detected for flg22 treatments. Elicitors did not alter the amount of Rubisco in barley leaves.

Since both elicitors induced the expression of some C1A CysProt, the next approach was to evaluate the effects of a pathogen, the fungus M. oryzae, and an herbivore, the spider mite T. urticae, in treated barley leaves. Firstly, the phenotype of damage was observed and the disease rating and leaf damage produced in the plant was recorded by scanning leaves at different days after infection/infestation (Supplementary Figure 1). The same CysProt studied after elicitor treatments were analyzed under M. oryzae and T. urticae attack. Results showed that at 3 dpi of M. oryzae the cathepsin F-like HvPap-1 was up-regulated (Figure 2A). After 7 days of M. oryzae infection, HvPap-1 maintained the highest mRNA levels; meanwhile significant overexpression of other protease genes such as HvPap-19 and HvPap-6, cathepsins B- and L-like respectively, was also observed (Figure 2B). The M. oryzae infected and non-infected barley leaves were also analyzed by immunoblotting in order to detect whether the CysProt accumulated after 3 and 7 days of fungus infection (Figure 2C). As previously reported Velasco-Arroyo et al. (2016), the protein profile of some CysProt showed two bands of different size corresponding to the immature protein and the mature processed form (without the pro-peptide). Both bands of HvPap-1 increased their signal in damaged leaf samples at 7 days of M. oryzae infection. A similar induction pattern was observed for the active form of the HvPap-6 and slightly for the -16 cathepsin L-like proteins, which were accumulated after fungus infection, particularly at 7 days. In addition, an induction pattern was observed for the active form of the cathepsin B-like HvPap-19, which was already accumulated at 3 days of treatment. No differences were detected for Rubisco levels in barley leaves.

On the other hand, at 7 days of T. urticae infestation (Figure 2D), no appreciable differences were shown in the expression of selected C1A CysProt. However, at 14 d of mite attack, the CysProt genes HvPap-1 and HvPap-19 were significantly up-regulated in infested leaves (Figure 2E). The protein profiles of CysProt HvPap-1, HvPap-6, HvPap-16 and HvPap-19 were analyzed by immunoblotting assays using total protein from T. urticae infested and non-infested barley leaves (Figure 2F). The lower (active mature protein) and the higher (inactive form) bands of these CysProt were observed in HvPap-1 and -6. HvPap-1 active form strongly increased after 14 days of mite infestation, whilst HvPap-6 revealed an increment of the inactive form after infestation in both time points, 7 and 14 days. HvPap-19 active form increased at 14 days of mite attack, while no differences of level of HvPap-16 were observed, in control or infested leaves at the two observed time points. Besides, 7 days assayed barley leaves showed a higher quantity of Rubisco protein than 14 days plants. And at this time, a slightly reduction of Rubisco in infested leaves could be appreciated.

HvPap-1 CysProt may have a relevant role in plant responses to biotic stresses, since it was the most up-regulated protease after all tested biotic treatments. Therefore, we used homozygous barley plants overexpressing or silencing the HvPap-1 gene to carry out in vivo experiments to test the resistance or susceptibility of modified plants toward the selected barley’s pathogen and herbivore. The implication of the HvPap-1 CysProt in the response to M. oryzae was first analyzed by comparing transgenic and non-transgenic plants after fungus infection (Supplementary Figure 2A). As expected, leaves from all infected plants presented more damage than those grown under control conditions (Figure 3A). After 3 days of infection, leaves from HvPap-1 knock-down lines (KD Pap1, 1130 and 1175) showed more damage than HvPap-1 overexpression lines (OE Pap1, 919 and 937) or wild-type (WT) barley plants. After 7 days of M. oryzae attack, the highest susceptibility of KD Pap1 plants was apparent, whilst the OE Pap1 lines presented fewer symptoms of damage and seemed less susceptible to the fungus than WT plants. These observations were corroborated by analyzing the leaf area damaged after infection and by comparing total damaged area in WT and transgenic infected lines (Figure 3B). The WT plants showed significantly lower damaged foliar area than KD Pap1 plants and greater damage than OE Pap1 transgenic lines at both time points, 3 and 7 days.

To study the effect of the transgenic barley plants on the fungus performance the presence of M. oryzae was measured by quantifying mRNA levels corresponding to the small subunit of ribosomal RNA (Mo28S-rRNA) (Supplementary Figure 3A). The results showed that at 5 days of infection a significantly higher quantity of fungus mRNA was appreciated in KD Pap1 barley plants indicating that these plants seem more susceptible to M. oryzae in comparison to WT and OE Pap1 plants. Finally, all lines showed increased levels of fungus mRNA at 7 d of infection, although the OE Pap1 lines presented the lowest levels, remaining as the less susceptible to the attack.

To compare plant responses against a pathogen and an herbivore, WT and homozygous plants overexpressing or silencing the CysProt HvPap-1 were used to carry out in vivo experiments to test the susceptibility of modified plants upon T. urticae infestation (Supplementary Figure 2B). The participation of HvPap-1 in the response to T. urticae was first analyzed by comparing transgenic and non-transgenic plants after mite infestation (Figure 4). As expected, leaves from infested plants were more affected than those grown under non-infested conditions (Figure 4A). After 7 days of infection, the most damaged leaves were observed in the OE Pap1 (919 and 937) lines, followed by the WT, whereas the KD Pap1 (1130 and 1175) lines presented only slight damages. At 14 days of mite infestation the phenotypic observations were further manifested, the differences among lines were stronger, OE Pap1 lines were much more injured than KD Pap1 lines, and WT plants showed an intermediate damaged phenotype. Furthermore, spider mite feeding effects were analyzed by quantifying the injured leaf area (Figure 4B). These results corroborated the phenotypical observations. The OE Pap1 lines had significant larger damaged foliar area than WT plants, and the KD Pap1 lines showed significant lesser injury than WT plants.

To study whether the transgenic barley plants affected T. urticae performance, the presence of the mite at different time points of infestation was analyzed by quantifying T. urticae Ribosomal Protein 49 (TuRP49) mRNA levels (Supplementary Figure 3B). The results showed a noticeable increment of mite mRNA in the OE Pap1 lines followed by WT barley plants at 7 days of infestation. At 14 days of T. urticae feeding, all lines showed increased levels of mite mRNA, although the KD Pap1 lines presented the lowest levels.

The expression patterns of several CysProt were investigated in control and treated plants, either by M. oryzae or T. urticae, in WT and transformed lines. As shown in Figure 5A and Supplementary Figure 4, HvPap-1 transcripts increased in plants infected with M. oryzae independently of the transgene insertion, although transcript levels were lower in transgenic lines than in WT. The mRNA profile of other genes encoding CysProt was also analyzed in these transgenic plants. HvPap-19 and HvPap-12 were up-regulated in infected WT, OE Pap1 and KD Pap1 lines. HvPap-6 sharply increased in KD Pap1 and WT when compared to the controls without infection but it did moderately in OE Pap1. No differences among treatments or lines were found for HvPap-4, and the HvPap-16 gene was repressed in response to M. oryzae infection. HvPap-1, HvPap-19, HvPap-6 and HvPap-16 proteases were detected by immunoblot in protein extracts from control and M. oryzae infected leaves of transformed and non-transformed plants using specific antibodies (Figure 5B). The HvPap-1 protein increased not only in the overexpressing OE Pap1 lines in comparison with the WT, but also after the M. oryzae infection treatments. In contrast, HvPap-1 active form diminished in KD Pap1 lines. A slight increase of HvPap-19 and HvPap-6 proteins was also observed in leaves grown under infection. No alterations in the HvPap-16 protein levels were detected in transgenic lines compared to WT plants or comparing among treatments. No differences were detected in Rubisco levels in barley leaves neither lines nor treatments.

As shown in Figure 6A and Supplementary Figure 5, HvPap-1 transcripts increased in plants infested with T. urticae independently of the transgene insertion. The expression of the HvPap-19 gene was up-regulated in infested WT, OE Pap1 and KD Pap1 lines, but KD Pap1 showed the lowest levels of this protease before and after infestation. HvPap-12 was slightly induced after T. urticae attack and HvPap-6 was highly induced in KD Pap1 and OE Pap1 lines. No remarkable differences were found for HvPap-4 and HvPap-16 genes. Immunoblotting assays were performed using protein extracts from control and T. urticae infested leaves (Figure 6B and Supplementary Figure 6). The inactive and active HvPap-1 protein forms increased not only in the overexpressing OE Pap1 lines in comparison with WT and KD Pap1, but also after the T. urticae infestation treatments in all lines. A slight increase of HvPap-19 and the inactive form of HvPap-6 was also observed in leaves grown under infestation treatments. No remarkable alterations in the HvPap-16 protein levels were detected. Low differences were detected in Rubisco levels, which were slightly reduced in all lines after mite treatment.

To obtain more information on the molecular basis responsible to the differential and opposite responses of transgenic HvPap-1 lines against M. oryzae and T. urticae stresses, a RNA-seq analysis was performed in 7 days old non-stressed plants. To test the robustness of the analysis, two different approaches were used to obtain DEGs between lines. Only a limited number of genes were detected as DEGs by either analytical methods, with around 70% of DEGs exclusively detected by the Kallisto/Sleuth or the SOAP/NOIseq methods (Supplementary Figure 7). When the expression of the transgenic lines was compared with the wild-type, only one DEG, HvPap-1, was up-regulated in overexpressing lines and down-regulated in silencing lines (Figure 7A). No other DEG was inversely regulated in both lines. On the contrary, a substantial number of DEGs were up or down-regulated when compared with the expression in wild type in both silencing and overexpressing lines. In the same way, only three DEGs were detected in the three OE Pap1/WT, KD Pap1/WT and OE Pap1/KD Pap1 analyses (Figure 7B) with around 60% of DEGs only detected in a unique comparison.

To further analyze the expression of the DEGs, a heatmap showing the differential values obtained by the two methods in all comparisons was performed. As expected, the results of both methods on the same comparison tend to show similar patterns (Figure 7C). Besides, clustering shows that only a reduced number of genes showed the OE > WT > KD or KD > WT > OE expression patterns. When we focused in the OE Pap1 vs KD Pap1 comparison, a robust relationship of DEGs with any specific biological process was not observed (Supplementary Table 3). Thus, the DEG lists were screened to detect genes potentially involved in the response to biotic stresses (Supplementary Tables 4, 5). From these analyses, several genes were selected as potential candidates to be involved in the different responses to biotic stresses in OE Pap1 and KD Pap1 plants. Two of these genes belong to the germin family, previously related to basal resistance against fungi (Zimmermann et al., 2006), and two genes were protease inhibitors, proteins broadly related to herbivore resistance (Schlüter et al., 2010), of the I12 Bowman–Birk and I13 Potato-I families. A comparative analysis of the expression data obtained in the RNA-seq experiment showed the highest expression of germin-like genes in the OE Pap1 line and the lowest in the KD Pap1 line (Supplementary Figure 8). On the contrary, the protease inhibitor belonging to the I13 family was more expressed in the KD Pap1 line and less expressed in the OE Pap1 line. Finally, the protease inhibitor from the I12 family exhibited lower expression in the OE Pap1 than in the KD Pap1 and WT plants.

The selected DEGs were analyzed by RT-qPCR in control or treated HvPap-1 silencing and overexpression barley plants as well as in WT (Figure 8 and Supplementary Figure 9). The variable expression patterns of these genes in non-treated control plants validate the data obtained by the RNA-seq analysis, showing the highest expression of germin-like genes in the OE Pap1 line and the lowest in the KD Pap1 line, whilst the protease inhibitors genes were less expressed in the OE Pap1 line. After biotic stress treatments, germin-like and protease inhibitor genes increased their levels of expression in all lines but with a different pattern depending on the biotic treatment and the plant genotype (Figure 8). Germin-like genes were remarkably more induced by M. oryzae in the KD Pap1 plants than in WT or OE Pap1 plants. Conversely, the induction of germin-like genes was greater in OE Pap1 plants after 14 days of T. urticae infestation. Protease inhibitor genes followed a similar induction pattern after 3 days of M. oryzae infection in all genotypes, but their expression showed a differentially reduced pattern in WT plants after 7 days of infection. In response to spider mite infestation, protease inhibitor genes were considerably more induced in OE Pap1 than in WT or KD Pap1 plants.