Ten plants per line were grown in both experiments. Seeds for the 2009 experiment were sown into pots (22 cm diameter) in a phytotron maintained at 20/15°C day/night temperatures with 60% humidity. Group C Rhizobium inoculum, in the form of clustered clay granules (ALOSCA Technologies Pty Ltd.), was applied 7 days after sowing.

Seeds for the 2011 experiment were sown into hydrated Jiffy-9 peat pots (Jiffy Products Ltd., Norway) in a glasshouse. Three seeds were sown per pot, which were randomly thinned to a single seedling per pot 20 days later. Group C Rhizobium inoculum (Nodulaid, BASF Australia Ltd) was sprinkled onto peat pots 7 days after sowing and watered in gently. Peat pots were watered daily, while soluble fertilizer (Phostrogen, Bayer CropScience Ltd) was watered on weekly at a rate of 5 g/10 L. Seedlings were transplanted to the field into a moist, weed-free seedbed, described as a loam-dressed sand (Nichols et al., 2009), with pH (CaCl₂) 6.5. The trial site was hand-weeded and irrigated by overhead sprinklers as required. Other site details and management are the same as described in Nichols et al. (2009).

The edible plant parts (leaves and stems less than 10 cm long) were harvested before flowering for in vitro fermentability testing. The material was then freeze-dried in a Bench-top Freeze Dryer (VirTis, Germany) and ground (Glen Creston, Stanmore, England) to pass through a 1.0 mm screen. Material was stored at room temperature in sealed containers for further analysis.

The methanogenic potential was examined in an in vitro batch fermentation system commonly used to examine plants and their extracts with rumen fluid (Cardozo et al., 2005; Busquet et al., 2006; Bodas et al., 2008). We followed the protocol by Durmic et al. (2010). Briefly, 1 day before the experiment, 0.1 g of plant material or oaten chaff was weighed in bellco tubes and transferred to an anaerobic chamber (Coy anaerobic chamber; 80% N₂:10% CO₂:10% H₂) to expel oxygen from the tubes. On the experimental day, rumen fluid was collected from two fistulated sheep fed a general maintenance diet consisting of lupins and oaten chaff (1 kg oaten chaff + 250 g lupins + 25 g mineral mix) 2 h after feeding. This study was approved by and the donor sheep were cared for in accordance with the guidelines The University of Western Australia (UWA) Animal Ethics Committee (AUP number RA/3/100/1171). An electric vacuum pump fitted with a plastic tube was used to collect rumen fluid and then transported to the laboratory where it was strained through cheesecloth to remove large particles. After straining, the rumen fluid was pooled, transferred into the anaerobic chamber, and buffered to pH 7.2 using McDougall’s buffer (McDougall, 1948); 10 mL of this mix was dispensed into the prepared bellco tubes. Assay controls included a negative control (buffered rumen fluid only, NC) and a positive control (buffered rumen fluid + 0.1 g oaten chaff, PC). Inside the chamber, the tubes were sealed with rubber stoppers, crimped and then incubated for 24 h at 39°C, with constant shaking at 50 rpm, to mimic the rumen environment. At the end of the incubation period, tubes were placed in a water bath at 39°C, and the contents were transferred to an exetainer tube (Labco, United Kingdom) for methane analysis by gas chromatography. Methane concentration in the gas sample was determined using a concentration by gas chromatography (Bruker 450 GC, Bruker Technologies, Australia) with two packed columns and a Compass CDS data acquisition software (Bruker Technologies, Australia). The methanogenic potential of the plant was expressed as total methane produced per gram dry matter of substrate incubated (mL/g DM incubated) (Soliva et al., 2008).

A linear mixed model was used for the analysis of the two-phase experiments (Smith et al., 2006). Typically, in a multi-phase experiment, possible trends associated with residuals for each of the phases are modeled. In this study, spatial trends in the glasshouse/field and temporal trends in the laboratory were considered. The model accounted for block structure (biological replicates), spatial variation in the field, laboratory day variation, the possible correlation of ordered laboratory measurements within a day, the linear trend of ordered measurements within a day if present, and methane control values (fitted as a covariate).

The first phase component in the model accounted for possible spatial variation in the glasshouse/field, in particular, the possible presence of linear row or column effects, local trends, and some extraneous variation if present (Gilmour et al., 1997; Stefanova et al., 2009). The second phase component included the blocking structure in the laboratory, possible linear trends, and a covariance structure linked to the order of the samples within the day. Additionally, two covariates were included, positive and negative controls for methane and gas, respectively, measured each day.

Random effects in each model were predicted using empirical best linear unbiased prediction (E-BLUP), and fixed effects were estimated using empirical best linear unbiased estimation (E-BLUE). The obtained prediction error variance matrix for the variety effects was used to calculate broad-sense heritability (Falconer, 1960) for each trial and each response.

All statistical analyses were performed using ASREML-R (Butler et al., 2009) and R software¹. The adjusted (predicted) mean values for methanogenic potential were used for quantitative trait loci (QTL) and genome-wide association study (GWAS) analyses.

The high-density single nucleotide polymorphism (SNP) linkage map was constructed using MultiPoint 3.3² as described in Hirakawa et al. (2016). QTL mapping was conducted using the high-density SNP linkage map anchored to the high-quality chromosome-level genome assembly (Kaur et al., 2017) by the integrated genetic analysis software QTL IciMapping v4.0³ (Meng et al., 2015). QTL with both additive and dominant genetic effects were screened using an inclusive composite interval mapping (ICIM) approach implemented in QTL IciMapping v4.0 (Zhang et al., 2008). Missing phenotypes were deleted using the ‘Deletion’ command in the software. The walking speed was set at 1 cM. A suitable probability for entering marker variables in stepwise regression was chosen so that the variation explained by the model approximated the trait heritability. The regression model was then used for background genetic variation control in the ICIM QTL mapping. The LOD was calculated using 1,000 permutations with the Type 1 error being 0.05, and significant QTLs were defined accordingly.

The population structure for GWAS was of two sub-populations: the first sub-population comprised 27 cultivars released in Southern Australia for grazing; while the second sub-population of 97 accessions was a core germplasm collection – a stratified sample of the world collection of T. subterraneum. The core collection was formed by a two-step selection process using ecological data and stratified proportional strategy followed by combined datasets and a maximizing strategy to best represent the entire collection (Ghamkhar et al., 2008).

Genomic DNA (gDNA) was extracted from a single plant from each of the 97 accessions of the subterranean clover core collection and the 27 diverse Australian cultivars, and sequenced. High-quality whole-genome resequencing (WGRS) data were generated for all 124 accessions and cultivars using paired-end sequencing libraries with an insert size of approximately 550 bp sequenced on an Illumina HiSeq 2000 sequencer. A total of 13.89 billion paired-end resequencing reads were generated with 5.27 billion reads uniquely mapped to the advanced genome assembly (Tsub_Refv2.0) (Kaur et al., 2017) using BWA (Li and Durbin, 2009). SNPs were identified using samtools and bcftools (Li et al., 2009; Li, 2011). The resulting SNPs were filtered by removing SNPs with at least one heterozygous allele, those with an MAF ≤ 5%, and those that were not present in at least one individual. Consecutive SNPs were merged using PLINK v1.9 (Purcell et al., 2007; Chang et al., 2015) into haplotype blocks if their r²-values were above 0.8. Linkage Disequilibium was visualized using Haploview v4.2 (Barrett et al., 2005).

The GWAS was conducted via a mixed linear model in TASSEL v5.2.30 (Bradbury et al., 2007) using four PCs and a kinship matrix as covariates as reported by TASSEL. Q-values were calculated using the R-package q-value (Storey et al., 2015). The Manhattan plot was drawn with qqman (Turner, unpublished).

Phenotypic information obtained from the cross 92S80 (cv. Woogenellup × cv. Daliak) and the diversity panel of 124 core collection lines and cultivars was associated with specific regions in the advanced assembly (Tsub_Refv2.0) (Kaur et al., 2017) using QTL and GWAS analyses, respectively, as described above to map significant associations. Each significant marker-trait association (MTA) resulting from the GWAS was checked for any overlaps with haplotype blocks with r²-values above 0.8. In which case, 25-bp upstream and downstream from the SNP were extracted from the reference and used as input for primer3 v2.3.7 (Koressaar and Remm, 2007; Untergasser et al., 2012) (settings: primer product size 250–500, primer optimum size 300, primer minimum temperature 55°C, optimal temperature 57°C, maximum temperature 60°C).

Putative candidate genes were proposed for each significant MTA by extracting the genes upstream, downstream or overlapping with GWAS candidate SNPs, or genes located between the left and right markers of candidate QTLs.

The methanogenic potential measured in the IVFT was assessed using a linear mixed model fitting procedure (Butler et al., 2009). High correlation values for methane production were found among laboratory triplicates within the 2009 (r-values: 0.686, 0.672, and 0.759) and 2011 (r-values: 0.845, 0.793, and 0.860) measurements. Furthermore, the broad-sense heritability (Falconer, 1960) calculated for the 2009 trial (H² = 0.4) reflected the complexity of the assayed trait.

The phenotyping information for methanogenic potential (Supplementary Tables S1, S2) was associated with specific regions in the high-density SNP linkage map. A significant association [P < 0.0002 at a threshold logarithm of odds (LODs) of threshold 2.5 defined by permutation tests] for methanogenic potential, Q_MP_LG8, was identified on Linkage Group 8 (LG8) with the phenotypic data available (Table 1). Q_MP_LG8 was identified as a significant QTL accounting for 7.8% of the variation in methanogenic potential. Within the bi-parental cross, the ‘Daliak’ alleles had a positive additive effect (increased methanogenic potential) (Table 1).

WGRS data from the diversity panel were analyzed with phenotyping data collected from the GWAS method using TASSEL v5.2.30 for estimating MTAs (Supplementary Figure S1). The GWAS resulted in four markers, which reached suggestive P-values below 1e^-5 on chromosomes 2, 5, 7, and 8 (Figure 1 and Table 2). The SNP located on Chr 5 was located in a haplotype block containing six other SNPs with a total length of 27.63 Kbp.

Putative candidate genes were proposed for each significant MTA by extracting the genes upstream, downstream or overlapping with GWAS candidate SNPs, or genes located between the left and right markers of candidate QTLs.

The Q_MP_LG8 at 231.7 Kbp long contained 11 of the genes suggested as putative candidates for this trait. Of these, there was a DNA-damage-repair/toleration DRT100-like protein (TSub_g26967), a suppressor-of-white-APricot splicing regulator (TSub_g26969), a transcription factor bHLH122-like protein, two disease resistance proteins (TSub_g26971, TSub_g26972), two ribonuclease proteins (TSub_g26974, TSub_g26975), a flowering time control FY-like protein (TSub_g26973) and a far1-related sequence 10-like protein (TSub_g26976) in addition to two uncharacterized proteins (TSub_g26968 and TSub_g26977) (Table 1).

For the GWAS SNPs, one suggestive SNP was located at the first base of an esterase/lipase/thioesterase protein family gene (TSub_g24126 on Chr 7). On Chr 8, we identified an upstream tubulin beta-1 chain (TSub_g27846) and a downstream transmembrane protein (TSub_g27847) at a distance of 12,050/3,060 bp respectively, from the suggestive SNP. The SNP located on Chr 5 was mapped in a haplotype block with TSub_g18548 (plant peroxidase-like superfamily) upstream at 22,314 bp and TSub_g18549 (lectin receptor protein kinase) downstream at 1,150 bp (Figure 2 and Table 2). Being in the haplotype block, these two genes are the most stable and promising candidates for designing molecular markers to track haplotype variation for this complex trait. The haplotype block containing the MTA SNP on Chr 5 with a total length of 27,629 bp was used to design PCR-ready markers for MAS for this complex trait (Supplementary Table S3).

The 124 diverse germplasm panel entries for methanogenic potential in the field ranged from 22.4 to 33.7 mL CH₄/g DM. The two lowest values were ssp. subterraneum (22.4 mL CH₄/g DM) and ssp. brachycalycinum (23.3 mL CH₄/g DM), both collected from Italy (Supplementary Table S2). The diversity panel was analyzed using one-way ANOVA to test for differences in methanogenic potential between the means of the germplasm collection and the Australian cultivars, the means of the three subspecies, and accessions from different countries of origin but no differences were identified (Supplementary Figures S2A–C and Table S4). Similar analyses were run to identify any significant associations of methanogenic potential with the available passport data from collection sites (latitude, longitude, altitude, soil pH) and 19 BioClim climatic variables, as described in Ghamkhar et al. (2014), but no significant associations were found (Supplementary Figure S3 and Table S5).