Switchgrass (Panicum virgatum L.) genotypes A4, AP13, and VS16 were grown in a greenhouse under approximately a 14-h day/10-h night photoperiod at 20–30°C in Norman, Oklahoma. Plants were watered once a week with deionized water and fertilized every 3 months with ∼15-g Osmocote 19-6-12 (Scotts). Single ramets from each genotype were propagated in 5-L pots containing a soil mixture of peat and topsoil (1:1) until producing multiple tillers. A total of nine plants, three biological replicate clones of each genotype, were grown. Harvesting took place between April and July.

Switchgrass plant developmental stages were determined as previously described (Moore et al., 1991). Single tillers for the vegetative 3 (V3), elongation 4 (E4), and reproductive 3 (R3) stages were collected from each biological replicate. All harvesting took place between 9 and 11 a.m. local time. Each tiller was harvested by cutting the tiller 3–4 cm above soil level, just beneath the first node, N1 (Figure 1). The sample was immediately chopped into small pieces with heavy-duty scissors into 50-ml conical tubes, followed by quick-freezing in liquid nitrogen. For each plant, an additional E4 tiller was further dissected into three segments. These consisted of S1 (lower, more mature) including the stem from just below the first node 1 (N1) above the soil to just below the second node and including the entire bottom leaf sheath and blade; S2 (middle, intermediate maturity) including the stem just below the second node to the stem just below the third internode and the associated leaf; and S3 (upper, least mature) including the stem just below the third node to the stem just below the fourth internode and the associated leaf uppermost node 4. Node 4 and all distal leaf and stem material were discarded. Samples were homogenized into fine powder in liquid nitrogen using mortars and pestles and stored at −70°C.

Alcohol insoluble cell wall residue (AIR) preparation and subsequent cell wall analysis were carried out essentially as described (Bartley et al., 2013a). Briefly, ∼100 mg of fresh frozen, ground tissue was treated 4–5 times with a 95% (v/v) ethanol wash for 30 min at 100°C in a thermomixer and centrifugation at 10,000 g for 10 min, with the supernatant removed after each wash. AIR was further washed three times with 70% (v/v) ethanol and then vacuum dried. The AIR was destarched using amylase (0.3 U/10 mg AIR) followed by amyloglucosidase (0.33 U/10 mg AIR) and pullulanase (0.04 U/10 mg AIR).

Cell wall carbohydrates were measured as previously described (Bartley et al., 2013a). Briefly, 2-mg destarched AIR was treated with 2 M of trifluoroacetic acid (TFA) at 120°C for 1 h. Monosaccharides were quantified relative to standards using high-performance anion exchange chromatography (HPAEC) with pulsed amperometric detection on an Dionex ICS-3000 system equipped with an electrochemical detector and a 4 × 250 mm CarboPac PA20 column, following the method described previously. The remaining glucose in the pellet was measured with the anthrone method.

Acetyl bromide soluble lignin content was used as a proxy for lignin content (Fukushima and Hatfield, 2004). For this analysis, 2 mg of destarched AIR was treated with 100 μl of freshly made acetyl bromide solution (25% v/v acetyl bromide in glacial acetic acid) for 3 h in a screw-cap tube at 50°C with continuous 1,050-rpm shaking and occasional gentle vortexing every 15 min for the last hour. Quantification in a 96-well plate and lignin content estimation were as described (Bartley et al., 2013a).

Destarched AIR (3 mg) was saponified using 500 μl of 2 M of NaOH with continuous mixing at 300 rpm in a thermomixer at 25°C for 24 h. The solution was then neutralized with 100 μl of concentrated hydrochloric acid (HCl) and extracted with 300 μl of ethyl acetate. Samples were dried and dissolved in 50% (v/v) methanol for HCA analysis using a Dionex Ultimate 3,000 high-performance liquid chromatography (HPLC) system. Details for the HPLC protocol and run parameters were as previously described, except that detection was at 305 nm (Bartley et al., 2013a).

ED pretreatment conditions were analyzed using a high-throughput assay system, with minor changes, as previously described (Santoro et al., 2010). Specifically, 2 mg of R3 greenhouse-grown ground total biomass from A4 and AP13 was pretreated with 150 μl of aqueous solution at different NaOH concentrations (0, 0.62, 1.5, 3.0, 6.2, 9.4, 12.5, 15.0, 30, and 62 mM) for 1 h at 90°C with continuous shaking. After cooling and adjusting the pH to 4.8, we carried out ED for 20 h at 50°C with a mixture of enzymes (Accellerase 1,000, Genencor, Rochester, NY) in 30 mM of citrate buffer (pH 4.5) plus 0.01% sodium azide.

For ED assays with developmental samples, destarched AIR (2 mg) was suspended in 100 μl of 100 mM citrate buffer (pH 5.0) at 30°C for 1 h and then pretreated with 1.5 mM NaOH at 100°C for 1 h with continuous shaking. After being cooled to room temperature and neutralized, the slurry was incubated with a 1:1:1 cocktail of the Novozymes enzymes cellulose (NS50013, diluted 1:10), β-glucosidase (NS50010, diluted 1:100), and xylanase and other minor enzymes (NS22002, diluted 1:10) at 50°C with continuous shaking. The released reducing sugars were quantified by 3,5-dinitrosalicylic (DNS) assay as reported (Bartley et al., 2013a).

Caffeic acid 3-O-methyltransferase (COMT) locus IDs and sequences were obtained from the literature and through BLASTP searches of genome databases. Sequences were screened for the presence of the O-methyltransferase domain (PF0089.18) using HMMER with the online search tool (Potter et al., 2018). Sequences lacking this domain were excluded. We obtained Arabidopsis COMT (AtCOMT1, AT5G54160) from TAIR and related sequence from the NCBI GenBank (NP_200227). We identified seven COMT sequences from the P. virgatum v4.1 draft AP13 genome sequence available through Phytozome (DOE-JGI¹; locus ID, protein name: Pavir.6NG060500.1.p, PvCOMT1a; Pavir.6KG070300.1.p, PvCOMT1b; Pavir.2NG564200.1.p, PvCOMT2a; Pavir.2KG513400.1.p, PvCOMT2b; Pavir.2NG567300.1.p, PvCOMT2c; Pavir.1KG549300.1.p, PvCOMT3a; and Pavir.1NG551300.1.p, PvCOMT3b). PvCOMT1 and PvCOMT2 were named and characterized previously (Wu et al., 2019); however, COMT1b and COMT2c lack PF0089.18 in this version of the annotation, and we excluded them from further analysis. Other COMT protein sequences were found with an orthologous groups search in Phytozome using AtCOMT, in the NCBI GenBank database (Benson et al., 2011), and from the literature on genetically and biochemically characterized COMTs. These are as follows: for Amborella trichopoda (AmTr_v1.0_scaffold00001.509, AmTrCOMT1), for Arabidopsis (AT3G5310, AtCOMTlike11), for Brachypodium distachyon (Bradi1g14870, BdCOMT1; Bradi3g16530, BdCOMT4, or BdCOMT6; Bradi2g02380, BdCOMT2; Bradi2g02390, BdCOMT3) (Wu et al., 2013), for safflower (Carthamus tinctorius: BAG71895, CtCAldOMT1), for barley (Hordeum vulgare: HORVU7Hr1G082280.1, HvCOMT1; HORVU3Hr1G116770.1, HvCOMT2; HORVU6Hr1G000040.1, HvCOMT3) (Daly et al., 2019), for alfalfa (Medicago sativa: AAB46623, MsCOMT1) (Zubieta et al., 2002), for tobacco (Nicotiana tabacum: X74452.1, NtCOMT1; X71430.1, NtCOMT2) (Maury et al., 1999), for rice (Oryza sativa: LOC_Os08g06100.1, OsCAldOMT1 or OsCOMT1; LOC_Os02g57760, OsCOMT2) (Koshiba et al., 2013), for poplar (Populus trichocarpa: ACC63886.1, PtCALdOMT1) (Wang et al., 2015), for Setaria italica (Seita.6G055900.1, SiCOMT1); for Sorghum bicolor (Sorbic.007G047300.1.p, SbCOMT1 or BMR12; Sorbic.004G351400.1.p, SbCOMT2) (Green et al., 2014), and for maize (Zea mays: GRMZM5G814904_P01, ZmCOMT or ZmBMR3; GRMZM2GO82007, ZmCOMT2).

For each COMT sequence, the O-methyltransferase domain was identified using a MOTIF search (KEGG, 2014). Multiple sequence alignment of the O-methyltransferase domains was performed in ClustalW (Thompson et al., 1994). Phylogenetic reconstruction was conducted using the maximum likelihood model in MEGA X with the LG amino acid substitution matrix and a discrete Gamma distribution (parameter = 1.1485) with five rate categories and 1,000 bootstrap replications (Kumar et al., 2018).

Primers for quantitative reverse transcription-PCR (qRT-PCR) of phenylpropanoid pathway genes were obtained from the literature (Shen et al., 2012, 2013). Primers for COMT gene family members were designed using default parameters of Primer Express software (version 3.0) (Applied Biosystems, Foster City, CA, United States) and selected based on low conservation among homologs. Sequences are listed in Supplementary Table 1.

Total RNA was isolated from 100 mg of frozen ground tissue using the RNeasy Plant Mini Kit (Qiagen, Germantown, MD, United States) following manufacturer’s instructions. RNA integrity was checked on a 1% (w/v) agarose gel and from the 260- to 280-nm absorbance ratio determined with a Take3 plate on a SynergyHT reader (BioTek, Winooski, VT, United States). For first-strand cDNA synthesis, 1 μg of total RNA was subjected to Turbo-DNase (Ambion, Austin, TX, United States) treatment as described in the manufacturer’s protocol followed by cDNA synthesis using the primer d(T)₂₀VN primer (Sigma, St. Louis, MO, United States) and SuperScript III reverse transcriptase (Life Technologies, Grand Island, NY, United States). Reactions were performed as described previously (Lin et al., 2016) using a Bio-Rad CFX96 thermocycler (Bio-Rad, Hercules, CA, United States) in a white 96-well plate with optical sealing film (Bio-Rad) in 10-μl final volume with 1 μM each of gene specific primers, 2 μl of cDNA, and 5 μl of 2 × SsoFAST EvaGreen Mastermix (Bio-Rad). The following protocol was used for all qPCR reactions: 95°C for 30 s, 40 cycles of 95°C for 1 s, and 60°C for 5 s, followed by default dissociation curves conducted by heating from 60 to 95°C to ensure a single target amplicon. LinRegPCR software was used to estimate PCR efficiency (Ramakers et al., 2003), and relative expression of each gene was determined from efficiency-adjusted ΔCq values. Microsoft Excel was used to calculate means and standard errors of three biological replicates, each measured with three technical replicates.

Analysis of variance (ANOVA), Tukey’s range test, correlation analyses, and principal component analysis (PCA) were carried out in R (version 2.15.2) (R Development Core Team, 2010). ANOVA was used to detect significant differences of cell wall composition and morphological characteristics among the three genotypes, developmental growth stages, and tiller segments. A Tukey’s range test was performed after ANOVA to detect significant differences at the P < 0.05 level, which are indicated by different letters in the figures and tables. The R package “stats” with the function “princomp” was used for PCAs.

Gini correlation analysis was employed to identify relationships among cell wall traits and between cell wall traits and GE. Gini correlation analysis is a hybrid parametric and non-parametric correlation method that uses both rank and value information (Ma and Wang, 2012). It is more tolerant to outliers, less dependent on sample size and data distribution, and better at detecting non-linear relationships than other correlation methods (Ma and Wang, 2012). Gini correlation coefficients (GCCs) and corresponding P-values were calculated using the “rsgcc” R package with 2,000 permutation tests (Ma and Wang, 2012). False discovery rate q-values were derived with the “q-value” R package (Storey, 2002; Storey and Tibshirani, 2003). For the correlations among cell wall parameters, we compared the GCC results with Spearman correlation coefficients (SCC) and Pearson correlation coefficients (PCCs) produced with the “psych” R package (Supplementary Table 2; Revelle, 2018). Most significant correlations were detected by all three methods, with similar coefficients, P, and q-values. PCC gave the fewest significant correlations, but all correlations that were significant via GCC were significant either with PCC or with SCC. The correlation network generated using the GCC values in the cell wall parameters dataset with q < 0.01 was visualized using Cytoscape version 2.8.3 (Shannon et al., 2003). A significance cutoff of q < 0.05 was applied for the correlations between GE and cell wall change values, whereas a cutoff of q < 0.01 was used for correlations between GE and cumulative cell wall parameters.

Among measured cell wall components, phenylpropanoids displayed the most significant differences among genotypes and across development. Acetyl bromide soluble lignin significantly varied among growth stages and tiller segments (ANOVA, P < 0.05). At the whole tiller level, plants at mature R3 stage contain significantly more lignin than immature V3 and E4 stages in all three genotypes [Figure 2A, Tukey’s honestly significant difference (HSD), P < 0.05]. A4 and AP13 have more lignin than VS16 at the V3 and E4 stages but not at the R3 stage. Thus, lignin varies among genotypes in young V3 and E4 stages but converges by the R3 stage. The oldest segment, S1_E4, contains significantly more lignin than the younger segment, S3_E4 segments of A4 and AP13 but not of VS16 (Figure 2B, Tukey’s HSD, P < 0.05), another indication of the lignin development delay in VS16 relative to A4 and AP13. The HCAs, para-coumaric acid (pCA) and ferulic acid (FA), also differ across genotypes and stages (Figure 3). Consistent with lignin being the major destination for pCA in the cell wall, pCA tends to increase from young to old plants. The amount of pCA at the R3 stage is significantly higher (Tukey’s HSD, P < 0.05) than that at V3 and E4 stages in VS16 and A4 (Figure 3A). At each stage, VS16 has less pCA than A4 and AP13 and less FA than A4. Across E4 segments, the FA concentration does not vary significantly; in the Alamos, increases in pCA with maturity are suggested, but not statistically significant (Figure 3B).

Cellulose and most matrix and pectin monosaccharides are invariant or do not exhibit clear trends across development or among genotypes (Tables 1, 2). Xylose and arabinose decline significantly from the V3 to E4 stage in AP13. TFA-soluble glucose, which derives from mixed linkage glucan and amorphous cellulose, declines with development in AP13 and A4, but not in VS16. Compared with trace amounts of fucose in some other grass species (Lin et al., 2016), switchgrass possesses a higher amount of fucose, ranging from 2.9 to 11 μg/mg, which increases with development. Very few statistically significant differences were detected among segments.

Cell wall residues from the three switchgrass genotypes were further characterized for cellulase ED using pretreatment conditions optimized to show differences among genotypes. To select pretreatment conditions, we evaluated the effect of varying NaOH pretreatment concentrations on yields of lignocellulose-derived sugars by ED of biomass from mature R3 tillers of A4 and AP13. Supplementary Figure 1 shows differences in reducing sugar yields (Δ) between AP13 and A4. At and below 3.0 mM NaOH pretreatment, the ΔAP13-A4 sugar yield was 0.014 mg/g AIR with significant P-values (P ≤ 0.03); however, at higher NaOH pretreatment concentrations (6.2, 9.4, 12.5, 15, 30, and 62 mM), the difference in sugar yields was not detectable and non-significant (P > 0.05). The maximum difference in ED between A4 and AP13 was observed with 1.5 mM NaOH pretreatment at 100°C (Supplementary Figure 1). Under these conditions, differences likely reflect structural properties, such as crosslinking, that prevent enzyme access to cellulose, as opposed to compositional differences, such as variation in the amount of cellulose. For all three genotypes, R3 tillers showed lower ED than V3 and E4 tillers (Figure 4A). Compared with A4 and AP13, VS16 showed significantly higher saccharification efficiency at the V3 and E4 stages (P < 0.05). However, as with lignin content, the difference in saccharification efficiency diminished by the R3 stage. Our assay possessed too much noise to statistically distinguish the E4 segments, but the trend was decreased digestibility in the more mature S1_E4 segment compared with the others for the Alamo genotypes (Figure 4B).

A PCA of the cell wall parameters indicated that developmentally older samples and the switchgrass genotypes are partially distinguishable (Figure 5A). Principal component 1 (PC1) was negatively influenced by ED but positively influenced by lignin, HCAs, and carbohydrates (Figure 5B). Principal component 2 (PC2) was negatively influenced by ED, HCAs, and carbohydrates including xylose, arabinose, and (non-crystalline cellulose) glucose but positively influenced by lignin and minor monosaccharides. Based on PC1 and PC2 (Figure 5A), A4 and VS16 had divergent cell wall properties, whereas AP13 was intermediate. Plant development also impacts cell wall composition and properties. Mature R3 and S1_E4 samples, shown by squares in Figure 5A, generally have higher values in both principal components than younger samples. As indicated by the spread in the PCA, the cell wall differences among genotypes are larger in younger samples (V3, E4, S3_E4, and S2_E4) than in older samples (R3 and S1).

To examine the influence of cell wall composition on saccharification efficiency, in particular, we analyzed the correlations among measured cell wall properties using the Gini correlation method (Figure 6). ED at different time points is strongly correlated. As expected, lignin and pCA negatively correlate with ED, as do minor sugars such as rhamnose, galactose, and fucose. Other cell wall carbohydrates, like cellulose and arabinoxylan, are not significantly correlated with ED possibly because of their low variation across samples. The correlations among cell wall components in this study are partially consistent with correlations based on a cell wall composition atlas of rice tissues from different stages (Lin et al., 2016). Glucuronoarabinoxylan components, including xylose, arabinose, and glucuronic acid, positively correlate with each other as in rice; however, there are also correlations not observed in rice but observed for switchgrass, like the positive correlation between pCA and lignin.

As lignin precursors and HCAs are synthesized by phenylpropanoid biosynthesis genes, we measured expression of selected switchgrass phenylpropanoid biosynthesis genes across the switchgrass genotypes, stages, and internode segments. Transcripts examined correspond to those described in the switchgrass literature, as follows (Shen et al., 2012, 2013): 4-coumarate:CoA ligase (4CL), coumarate 3’-hydroxylase (C3’H), cinnamate 4-hydroxylase (C4H), cinnamyl alcohol dehydeogenase (CAD), caffeoyl CoA 3-O-methyltransferase (CCoAOMT), cinnamoyl coenzyme A reductase (CCR1), ferulate 5-hydroxylase (F5H), and hydroxycinnamoyl-CoA/shikimate/quinate hydroxycinnamoyltransferase (HCT). Due to the importance of S:G ratios in cell wall recalcitrance, as partially determined by the action of COMT enzymes (Fu et al., 2011), we identified and separately measured transcript abundance from three COMT homologs, numbered 1 through 3, with some of our qRT-PCR primers recognizing multiple homeologs in the current version of the switchgrass genome (Supplementary Table 1). Phylogenetic reconstruction of the relationship of the three COMT gene product targets relative to COMTs in the literature indicates that PvCOMT1 and PvCOMT2 are orthologous to biochemically and genetically studied barley (Daly et al., 2019), sorghum (Green et al., 2014), rice (Koshiba et al., 2013), and Brachypodium (Wu et al., 2013) COMTs, which are also sister to the Arabidopsis COMT protein. On the other hand, PvCOMT3 is in a distinct clade that includes OsCOMT2 and is sister to NtCOMT2. In tobacco, NtCOMT2 transcript was induced in response to virus infection and did not vary with development (Pellegrini et al., 1993).

Quantitative real-time PCR measurements for the switchgrass phenylpropanoid biosynthesis genes showed that most genes have similar expression patterns in all three genotypes but vary across development (Figure 7 and Supplementary Table 3). PCA indicates that genotypes, stages, and E4 segments are not distinguished on the basis of the expression of these genes (Supplementary Figure 3), with the lowland Alamos, A4 and AP13, differing more from each other than from the upland VS16. Among the three developmental stages, most genes had the highest expression at the E4 stage, with geometric mean of relative expression being approximately threefold higher than in the R3 and twofold higher than the V3 stage. In the internode segments of A4 and AP13, expression of the COMTs increased from the young S3_E4 segment to the older S2_E4 and S1_E4 segments, but this pattern was not clear in VS16. VS16 instead shows an increase with segment maturity for 4CL and HCT. Among the segments across genotypes, the geometric mean across tested genes was the highest in the more mature S1_E4 segment, but besides those noted, patterns differed among genes and genotypes.

To understand how expression of phenylpropanoid biosynthesis GE influences cell wall phenolic compounds and ED (cell wall parameters, CW), we analyzed Gini correlations between these components with various models. The first model (N_GE vs. N_CW) hypothesizes a correlation between concurrent expression of phenylpropanoid biosynthesis genes and cell wall parameters in the same sample (N). For example, this correlates expression in V3 with cell wall features in V3, or expression in S1_E4 with cell wall features in the S1_E4 segment. Next, in two “delta” models, we compared GE at time N with the change (Δ) in a cell wall component between N and samples later in the developmental series (N_GE vs. Δ[(N + 1) - N_CW]) and (N_GE vs. Δ[(N + 2) - N_CW]). These models tested the hypothesis that expression earlier in development would lead to cell wall properties later in development. For example, the (N_GE vs. Δ[(N + 1) - N_CW]) model includes V3 expression correlations with E4 composition and S1_E4 GE with S3_E4 CW. Supplementary Table 4 contains the GCC and associated P and q-values for these N_GE vs. N_CW and N_GE vs. ΔN_CW models.

Significant correlations (q < 0.05) in the concurrent (N_GE vs. N_CW) and Δ delay (N_GE vs. Δ[(N + 1) - N_CW] and N_GE vs. Δ[(N + 2) - N_CW]) models differ qualitatively (Figure 8). In the N_GE vs. N_CW model (Figure 8A), phenylpropanoid biosynthesis genes show a mix of expected and unexpected correlations. Without exception, correlations between expression and FA, when significant, are positive. With lignin and pCA, there are a few expected positive correlations and a few unexpected negative correlations. Similarly, ED, which we expect to be negatively correlated with phenylpropanoid GE, shows mostly positive correlations. In contrast, in the N_GE vs. Δ[(N + 1) - N_CW] model (Figure 8B), a majority of phenylpropanoids synthesis genes, including C4H, 4CL, CCR1, F5H, and two of the COMTs, positively correlate with the change in lignin and pCA. A few negative correlations, and no unexpected positive correlations, are significant between GE and ED in this model. Finally, in the model with the larger delay (N_GE vs. Δ[(N + 2) - N_CW], Figure 8C), transcript abundance and the change in lignin abundance are positively correlated for all lignin biosynthesis genes, except the COMTs. Likewise, transcript abundance and ED are negatively correlated for all but the COMTs. In this model, pCA abundance positively correlates with CCR1 and CCoAOMT expression. In summary, the delay models produce more expected significant correlations for lignin, pCA, and ED, while the concurrent model shows more expected correlations between GE and FA.

We next asked how phenylpropanoid GE correlates with cumulative cell wall properties in the developmental series, again with various models that include delays and separating genotypes, stages, and segments. We particularly were interested in the question of which samples’ GE indicate cell wall properties at maturity. Table 3 shows the various Gini correlation models tested, sorted based on precision, i.e., the number of expected significant correlations divided by the total significant correlations. The input data are available in Supplementary Table 5. Expected correlations include positive correlations between phenylpropanoid GE and accumulated abundance of phenylpropanoid-derived cell wall components (lignin, pCA, and FA) and negative correlations with ED. While unexpected correlations are the opposite, such as significant negative correlations between GE and phenylpropanoid components.

Though limited by the small number of samples, genotypes, and transcripts examined, our result is that cumulative CW properties are generally highly correlated with GE at a recent preceding series sample. For example, the highest performing model was for the N_GE vs. (N + 1)_CW for developmental stages (excluding the segment data). This data subset revealed 20 expected correlations and no unexpected correlations (q < 0.01, Table 3 and Supplementary Figure 4). The next four models with high numbers of expected correlations and no unexpected correlations are for E4_GE to R3_CW and for E4_GE segments (S1, S2, and S3) to R3_CW (Table 3). In these models, expression of CCR, the COMTs, and F5H in the entire E4 tiller or in E4 segments shows significant correlations with lignin or ED in R3 samples (Supplementary Figure 4). In contrast to the delta CW analysis in Figure 8, the models that return the fewest expected correlations are those that correlate expression in V3 with later composition, e.g., V3_GE vs. R3_CW.