Transcriptome analysis by RNA-seq

Three leaf fragments (~1 cm² each) of 8–10 plants per plot were collected at the IAL field assay, placed in liquid nitrogen, and stored at –80ºC. Samples from R5–R6 were used for RNA-seq. Total RNA was extracted from pulverized samples with RNAeasy (Qiagen). RNA quality and integrity were checked by absorbance (A) measurements (A260/A280 >1.8, A260/A230 >2.0) and electrophoresis. RNA was analyzed by BGI (San Jose, USA) by sequencing eight libraries. An average of 82 468 181 clean reads/sample with more than 95% of the bases having a Phred value>20 were reported.

Raw paired-end reads were first quality trimmed with Trimmomatic (version 0.36; Bolger et al., 2014) and then aligned to the Glycine max W82 genome, v4 (Schmutz et al., 2010; from Phytozome V13, Goodstein et al., 2012) using STAR (version 2.5.2b, Dobin et al., 2013) with a maximum intron length of 1200 bp. Using samtools (version 1.8; Li et al., 2009), only primary alignments with a minimum MAPQ of 3 were kept. Read quality before and after trimming was analyzed with FastQC (version 0.11.5; Andrews, 2010) and, together with mapping efficiency, was summarized with MultiQC (version 1.7; Ewels et al., 2016). Read counts on each gene were calculated with featureCounts (version 1.6.2; Liao et al., 2014) using the gene and exon annotation from Phytozome (v13; Goodstein et al., 2012). Differentially expressed genes (DEGs) were determined with DESeq2 (Love et al., 2014; R Core Team, 2018), filtering out genes with counts below 10 in all samples. This analysis pipeline was run with the aid of the Snakemake workflow engine (Köster and Rahmann, 2018). Gene ontology analysis was performed online with agriGO (v2; Tian et al., 2017).