Phenylpropanoid Biosynthesis Gene Expression Precedes Lignin Accumulation During Shoot Development in Lowland and Upland Switchgrass Genotypes

Efficient conversion of lignocellulosic biomass into biofuels is influenced by biomass composition and structure. Lignin and other cell wall phenylpropanoids, such as para-coumaric acid (pCA) and ferulic acid (FA), reduce cell wall sugar accessibility and hamper biochemical fuel production. Toward identifying the timing and key parameters of cell wall recalcitrance across different switchgrass genotypes, this study measured cell wall composition and lignin biosynthesis gene expression in three switchgrass genotypes, A4 and AP13, representing the lowland ecotype, and VS16, representing the upland ecotype, at three developmental stages [Vegetative 3 (V3), Elongation 4 (E4), and Reproductive 3 (R3)] and three segments (S1–S3) of the E4 stage under greenhouse conditions. A decrease in cell wall digestibility and an increase in phenylpropanoids occur across development. Compared with AP13 and A4, VS16 has significantly less lignin and greater cell wall digestibility at the V3 and E4 stages; however, differences among genotypes diminish by the R3 stage. Gini correlation analysis across all genotypes revealed that lignin and pCA, but also pectin monosaccharide components, show the greatest negative correlations with digestibility. Lignin and pCA accumulation is delayed compared with expression of phenylpropanoid biosynthesis genes, while FA accumulation coincides with expression of these genes. The different cell wall component accumulation profiles and gene expression correlations may have implications for system biology approaches to identify additional gene products with cell wall component synthesis and regulation functions.


INTRODUCTION
Due to increased transportation energy usage and urgency to reduce fossil fuel use, demand for advanced fuels is predicted to increase to 79.5 billion liters by 2022 (US CRS Report, 2009), about 7% of annual petroleum utilization (EIA, 2019). Biofuels from lignocellulosic biomass, such as the leaves and stems of perennial grasses, hold promise to sustainably fulfill a significant fraction of the alternative fuel requirement with low greenhouse gas emissions (Schmer et al., 2008;Gelfand et al., 2013). Biochemical processes are now being deployed that convert polysaccharides, typically cellulose, from plant cell walls into alcohol fuels (Youngs and Somerville, 2012;Torres et al., 2016). However, lignin and hydroxycinnamic acids (HCAs) covalently crosslink cell walls and reduce saccharification efficiency during biomass enzymatic digestibility (ED) (Sattler and Funnell-Harris, 2013). Altering expression of single genes, especially those from the phenylpropanoid biosynthesis pathway, which synthesizes lignin and HCAs, improves biomass processing efficiency (Baxter et al., 2014;Li et al., 2018). Still, questions remain as to which manipulations are optimal and how genetic diversity can be harnessed to achieve simultaneous biomass composition and yield improvements. One approach to address this is to associate transcriptomes with biomass properties at harvest, but designing such studies requires an initial understanding of the relationships between gene expression (GE) and composition across development and genotypes. This study aims to reduce this knowledge gap for switchgrass.
Among potential dedicated bioenergy grasses, switchgrass (Panicum virgatum L.) is a front-runner species for lignocellulosic feedstock production in the United States (Bouton, 2007;Casler et al., 2011;Bartley et al., 2013b). Switchgrass is a C4, warm-season perennial that produces high annual biomass yield (typically ≥12 Mg/ha) and exhibits broad environmental adaptation (Lowry et al., 2019). This primarily outcrossing species consists of upland and lowland ecotypes and possesses high genetic diversity (Zalapa et al., 2011;Lu et al., 2013). Lowland genotypes are typically tetraploid, whereas uplands are octaploid or tetraploid (Lu et al., 2013). Most switchgrass cultivars have only undergone a few rounds of selection, and more genetic diversity exists within a cultivar than among cultivars (Cortese et al., 2010). That said, there are established characteristics that typify the ecotypes. Relative to upland cultivars under the same conditions, lowlands tend to cease growth later and have longer, thicker stems, contributing to lowlands typically accumulating greater biomass than uplands (Lowry et al., 2014). On the other hand, uplands exhibit greater drought and cold tolerance than lowlands (Stroup et al., 2003;Ayyappan et al., 2017).
The relationships between GE and cell wall properties across development and among ecotypes and genotypes remain underexplored. Generally, as plants mature, secondary cell wall formation and lignification occur; as a result, mature tissue contains a higher proportion of lignin and is less digestible . Previous cell wall and GE analyses revealed variations among developmental stages and internodes from lowland switchgrass of the Alamo cultivar (Mann et al., 2009;Shen et al., 2009;Escamilla-Treviño et al., 2010;Shen et al., 2013). Cell wall digestibility, a key output of biomass composition, varies across switchgrass development due to changes in cell wall components, with strong negative correlations between digestibility and total lignin, lignin monomers, and HCAs in Alamo switchgrass (Shen et al., 2009;Hu et al., 2010). However, whether developmentally associated cell wall changes are consistent among genotypes has not to our knowledge been examined. When genotypic variation effects on switchgrass cell wall composition has been examined among cultivars, including across ecotypes, these studies have focused on a single stage (Lemus et al., 2002;Hu et al., 2010).
This work reports the cell wall composition, ED, and lignin biosynthesis GE of a series of developmentally matched samples from A4, AP13, and VS16 switchgrass genotypes. We find that the measured cell wall parameters, especially phenylpropanoid content, vary across development and in many cases among genotypes. For example, VS16 is more digestible than A4 and AP13 at earlier developmental stages but not at reproduction. Generally, expression of phenylpropanoid biosynthesis genes precedes lignin accumulation, suggesting that early GE may be an indicator of cell wall properties later in development.

Plant Material
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.

Sample Collection
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 Residue Preparation
Alcohol insoluble cell wall residue (AIR) preparation and subsequent cell wall analysis were carried out essentially as FIGURE 1 | Diagram of switchgrass stages and segments characterized in this experiment. N indicates nodes. Panicle was removed from the R3 tillers.
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).

Carbohydrate Composition
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.

Lignin Content
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,050rpm 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).

Hydroxycinnamic Acid Content
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).

Enzymatic Digestibility
ED pretreatment conditions were analyzed using a highthroughput 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.

Identification and Analysis of Caffeic Acid 3-O-Methyltransferase Sequences
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. 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, 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).

RNA Isolation and Quantitative Reverse Transcription-PCR
Primers for quantitative reverse transcription-PCR (qRT-PCR) of phenylpropanoid pathway genes were obtained from the literature (Shen et al., 2012. 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) 20 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.

Data Analysis and Statistics
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 FIGURE 2 | Lignin of switchgrass genotypes A4, AP13, and VS16. Samples are arranged from developmentally young tissues to old tissues. (A) Whole tillers at V3, E4, and R3 developmental stages and (B) S3, S2, and S1 tissue segments of the E4 tiller as shown in Figure 1. Error bars indicate standard error (N = 3). Means with unshared letters are likely to be significantly different [Tukey's honestly significant difference (HSD), P ≤ 0.05].
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.

RESULTS
Switchgrass development has been classified into vegetative (V), elongation (E), and reproductive (R) stages (Moore et al., 1991;Hardin et al., 2013). At each stage, different segments (phytomers) develop acropetally, with the upper internodes developing after the bottom segments. To address how changes in development among genotypes alter cell wall properties, we analyzed lignin, hydroxycinnamates, cell wall carbohydrates, cell wall ED, and phenylpropanoid biosynthesis GE of whole tillers at the V3, E4, and R3 stages and three tiller segments at E4 stage, from the bottom segment, S1 E4 , to the top segment, S3 E4 (Figure 1). These samples were studied for three switchgrass genotypes: two lowland, Alamo genotypes (A4 and AP13) and one upland, Summer genotype (VS16). Supplementary Table 3 contains the cell wall and GE data generated in this study.

Cell Wall Properties
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. TFAsoluble 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 lignocellulosederived 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 FIGURE 3 | Hydroxycinnamic acids (HCAs), ferulic acid (FA) and para-coumaric acid (pCA), of switchgrass genotypes A4, AP13, and VS16. Samples, diagramed as in Figure 1, are arranged from young to old stages/segments. (A) Whole tillers at V3, E4, and R3 developmental stages and (B) S3, S2, and S1 tissue segments of the E4 tiller. Error bars indicate standard error (n = 3). Means with unshared letters are likely to be significantly different [Tukey's honestly significant difference (HSD), P ≤ 0.05].
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).
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.

Phenylpropanoid Biosynthesis Gene Expression
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: 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 , sorghum , rice , and Brachypodium  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.

Correlations Between Gene Expression and Cell Wall Parameters
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 FIGURE 6 | Significant correlations among cell wall components and enzymatic digestibility (ED). Red and blue lines indicate positive and negative correlations, respectively. The thickness of lines is proportional to the absolute value of Gini correlation coefficient. Cell wall components are represented by circular nodes and color-coded based on the cell wall polymers from which they originate. Xyl, xylose; Ara, arabinose; Glc, glucose; Gal, galactose; Rhm, rhamnose; GalA, galacturonic acid; FA, ferulic acid; pCA, p-coumaric acid; ED, enzymatic digestibility. The glucose is non-crystalline cellulose, trifluoroacetic acid (TFA)-soluble glucose. GAX is glucuronoarabinoxylan. Only correlations with a q < 0.01 are shown. Correlation coefficients calculated with Gini, Pearson, and Spearman methods and corresponding P and q-values are provided in Supplementary Table 2. 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 FIGURE 7 | Relative gene expression of phenylpropanoid biosynthesis genes by quantitative RT-PCR. Cq values, relative to switchgrass Ubi10, were log 10 transformed. Color intensity indicates greater expression. 4CL, 4-coumarate:CoA ligase; C3'H, para-coumarate 3'-hydroxylase; C4H, cinnamate 4-hydroxylase; CAD, cinnamyl alcohol dehydrogenase; CCoAOMT, caffeoyl CoA O-methyltransferase; CCR, cinnamyl CoA reductase; COMT, caffeic acid 3-O-methyltransferase; F5H/Cald5H, ferulate 5-hydroxylase/ coniferaldehyde 5-hydroxylase; HCT, para-hydroxycinnamoyl-CoA: Quinate/shikimate p-hydroxycinnamoyltransferase.
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 .

DISCUSSION
This study examined the generality of switchgrass cell wall synthesis GE and cell wall composition among genotypes toward establishing methods for analyzing the determinants of cell wall composition and biorefining suitability across genotypes. Among the examined switchgrass genotypes at the V3 and E4 stages, lignin content, HCAs, and cell wall digestibility varied significantly, but the differences diminished at the R3 stage. The low abundance of lignin and HCAs in upland VS16 samples relative to lowland A4 and AP13 samples was associated with high ED. Due to superior digestibility, VS16 biomass at V3 and E4 stages has the potential to yield more biofuel than A4 and AP13, but VS16 loses this advantage by the R3 stage. A caveat of this work is that the greenhouse conditions might have had a greater influence on upland, VS16 (Casler et al., 2004). Still, the combined effect of developmental stage and genetic background should be considered when selecting cultivars with promising cell wall traits or when genetically modifying plant biomass for higher saccharification efficiency (Ashworth et al., 2017). For example, a plant with delayed lignin accumulation, such as that exhibited by VS16 under the tested conditions, may have a greater time to accumulate biomass before lignin deposition decreases saccharification efficiency.
In contrast to phenylpropanoids, we observed fewer and less pronounced differences in sugar composition among genotypes and stages. While lignification and hydroxycinnamates are determining factors for reduced saccharification and cell wall digestibility (Dien et al., 2006;McCann and Carpita, 2008;Buanafina, 2009), non-cellulosic polysaccharide components from xylan, pectins, arabinogalactan proteins, and xyloglucan also reduce digestibility (DeMartini et al., 2013;Chung et al., 2014;Li et al., 2019). Indeed, the negative correlations that we observe between ED and the pectic sugars, rhamnose and galactose, are consistent with these observations. Still, in our results, polysaccharide component abundance changed little across developmental stages and tiller segments. Thus, lignin remains a key genetic engineering target since it varies across development and among genotypes, while polysaccharide content appears relatively static and perhaps less tolerant of manipulation.
The relative timing of phenylpropanoid biosynthesis GE and cell wall properties provides targets for further study of cell wall-related GE strategies to control harvested biomass digestibility. Phenylpropanoid biosynthesis GE is the highest at the E4 stage, apparently leading to the approximately twofold increase in lignin from the E4 to R3 stage. This result was consistent with our overall analysis indicating that GE in prior stages (N GE vs. [(N + 1) -N CW ] and N GE vs.
[(N + 2) -N CW ] models) rather than concurrent expression (N GE vs. N CW model) captures expected positive correlations between phenylpropanoid biosynthesis transcript and a change in lignin abundance (Figure 8). Indeed, lignin deposition may be a slow process, in that it includes precursor synthesis, dehydrogenation, and polymerization (Raes et al., 2003). An additional explanation consistent with the data is that lignin undergoes turnover, i.e., breakdown and recycling, that is more rapid than accumulation during earlier developmental stages. Besides evidence for changes to the lignin-containing Casparian strip of roots (Vermeer et al., 2014), we are unable to find clear evidence of this in the literature, suggesting that additional research is required. Consistent with our results, transcriptomics of maize internode subsegments also showed that phenylpropanoid biosynthesis peaks prior to maximal lignin accumulation (Zhang et al., 2014). In our analysis, which includes different genotypes, the observation that correlations between GE and a decrease in ED are the most negative in the N GE vs. [(N + 2) -N CW model hints that high phenylpropanoid biosynthesis GE earlier in development may establish later digestibility recalcitrance. Furthermore, the delay models suggest that expression of all the examined phenylpropanoid biosynthesis genes correlates with lignin accumulation, albeit with different relationships with different genes. Potentially indicating an absence of functional conservation across species, this includes COMT3, for with the tobacco ortholog did not vary with development (Pellegrini et al., 1993).
The HCAs show different accumulation kinetics and GE correlations compared with lignin. Though FA is also synthesized by phenylpropanoid biosynthesis genes, it is mainly esterified to arabinoxylan (Bartley et al., 2013a). This different polymer destination may explain why FA positively correlates with phenylpropanoid biosynthesis genes in the concurrent (N GE vs. N CW ) model but not in the delay models. It may also indicate active turnover of FA and/or xylan in the wall (Franková and Fry, 2011). These patterns-early expression controlling lignin accumulation, but sustained expression controlling FA accumulation-may be used to insinuate regulators and other enzymes involved in these alternative processes. Likewise, the difference in correlations between phenylpropanoid biosynthesis GE with pCA and lignin abundance (Figure 8) hints at extra control of pCA synthesis, which is consistent with new evidence that hydroxycinnamoylated monolignols are synthesized by different enzymes than un-acylated monolignols (Takeda et al., 2018). This notion is supported by the absence of a correlation between C3'H and pCA abundance, since reduction of C3'H in rice does not alter abundance of cell wall monolignols esterified with pCA (Takeda et al., 2018).

CONCLUSION
In conclusion, this study advances our knowledge of the relationships between switchgrass developmental stage and genetic background with GE, cell wall composition, and digestibility. Under greenhouse growth conditions, lignin, HCAs, and ED differ significantly, and in a developmentally dependent manner, among genotypes. Though developmental cell wall profiles differed among genotypes, general relationships emerged. Correlations between cell wall phenolics and phenylpropanoid biosynthesis GE revealed that FA accumulation may precede, or be less stable than, pCA and lignin accumulation. Expression of phenylpropanoid biosynthesis genes in mid-development seems likely to be an indicator of cell wall properties, especially lignin, pCA, and ED, at harvest. Systems analysis with more genotypes and global expression analysis will be needed to confirm these conclusions. The recommendation from the analysis presented here would be that such studies should probe GE in whole or partial tillers at mid-development as an indicator of biomass properties at harvest.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
PS, FL, and LB conceived and designed the study and wrote the manuscript. PS and ST generated the data. All authors analyzed the data, edited the manuscript, and approved the final version.