Genomic characterisation and dissection of the onset of resistance to acetyl CoA carboxylase-inhibiting herbicides in a large collection of Digitaria insularis from Brazil

An in-depth genotypic characterisation of a diverse collection of Digitaria insularis was undertaken to explore the neutral genetic variation across the natural expansion range of this weed species in Brazil. With the exception of Minas Gerais, populations from all other states showed high estimates of expected heterozygosity (HE > 0.60) and genetic diversity. There was a lack of population structure based on geographic origin and a low population differentiation between populations across the landscape as evidenced by an average Fst value of 0.02. On combining haloxyfop [acetyl CoA carboxylase (ACCase)-inhibiting herbicide] efficacy data with neutral genetic variation, we found evidence of the presence of two scenarios of resistance evolution in this weed species. Whilst populations originating from north-eastern region demonstrated an active role of gene flow, populations from the mid-western region displayed multiple, independent resistance evolution as the major evolutionary mechanism. A target-site mutation (Trp2027Cys) in the ACCase gene, observed in less than 1% of resistant populations, could not explain the reduced sensitivity of 15% of the populations to haloxyfop. The genetic architecture of resistance to ACCase-inhibiting herbicides was dissected using a genome-wide association study (GWAS) approach. GWAS revealed the association of three SNPs with reduced sensitivity to haloxyfop and clethodim. In silico analysis of these SNPs revealed important non-target site genes belonging to families involved in herbicide detoxification, including UDPGT91C1 and GT2, and genes involved in the vacuolar sequestration-based degradation pathway. Exploration of five genomic prediction models revealed that the highest prediction power (≥ 0.80) was achieved with the models Bayes A and RKHS, incorporating SNPs with additive effects and epistatic interactions, respectively.

Methods
Plant material

A total of 205 D. insularis weed populations were used in the study. The populations were collected from the most relevant soybean- (and corn) growing regions across Brazil inhabiting crop fields and crop margins in the southern, mid-western, south-eastern, and north-eastern regions. Briefly, these 205 populations originated from seven states including Mato Grosso (77), Mato Grosso do Sul (10), Minas Gerais (19), Bahia (24), Goiás (20), Paraná (PR) and Rio Grande do Sul (18). The detailed information on the geographic origin and geographic coordinates are provided in the supplementary table (Table S1). Some collection sites (1/5th of the collection) had a history of multiple herbicide treatments on the crop including glyphosate, ACCase inhibitors such as haloxyfop or clethodim alone, or in mixtures and/or ALS inhibitors (Table S1). The populations were collected in 2020 and 2021 and were investigated for their sensitivity to ACCase inhibitors haloxyfop and clethodim at multiple dose rates by the Weed Research team at Brazil Resistance Management Laboratory in Uberlândia, Brazil, as part of resistance monitoring studies.

Herbicide resistance screening of D. insularis collection

The rates used in the population screening were initially determined in a pilot study using ten D. insularis-sensitive populations sampled in Brazil's urban areas with no history of herbicide application. Informative herbicide rates were determined as 7.8; 15.0, 27.0, and 64.0 g a.i ha-1 for haloxyfop and 27.0, 54.0, and 108.0 g a.i ha-1 for clethodim.

Seeds from all D. insularis populations were germinated and planted in 1 L pots filled with a commercial substrate to produce around 15 plants per pot. Each combination of population and herbicide dose was replicated 3 times (45 plants tested per treatment) in a completely randomized design. The herbicide treatments were sprayed on plants at the 4 leaf-stage, in a spray chamber equipped with flat fan nozzles calibrated to deliver 200 L ha-1 at 200 kPa pressure. Plant control was evaluated 21 days after treatment, using a visual scale of 0 to 100%; 0% represents healthy plants and the absence of symptoms, and 100% represents the death of the plant.

DNA extraction, Genotyping-by-sequencing (GBS) library preparation and sequencing

Four pinches of seeds were sown in a punnet of size 18 cm x 6 cm filled with peat keeping 3cm between pinches. For each population, two such punnets were sown. The punnets were watered manually and put on trolleys in a glass house with controlled conditions (day temperature, 240C; night temperature, 180C; light, 16 hrs; humidity -65%). The punnets were watered daily for three weeks. After three weeks, a 1cm x 2cm of leaf sample was cut from 25 individual plants for each population and pooled into a 14ml falcon tube. The falcon tubes were stored in a −80 °C freezer until further manipulation.

The leaf samples were dried in a freeze dryer for three consecutive days and nights by keeping the shelves at a contact temperature of 1.0 0C and the freezer at -60 0C. After freeze drying, the samples were shipped to LGC Genomics GmbH, Germany for DNA extraction, reduced representation library preparation, and sequencing.

Genomic DNA extraction was performed from the pooled samples using the sbeadex™ maxi plant kit (LGC) on KingFisher Flex (after the lysis step) followed by a spectrophotometric quantification step using Nano Drop 8000 (Thermo Fisher Scientific). Reduced representation library preparation was done by the standardized ddRAD protocol at LGC Genomics. Briefly, 100 ng of genomic DNA were digested with 2 units each of Apek I and Pst I enzymes (NEB) in 1 times NEB buffer 3.1 in 20μl volume for one hour at 37°C. The restriction enzymes were heat-inactivated by incubation at 75°C for 60 min. The detailed protocol for the ligation reaction, library purification, amplification, and normalization were performed according to the standardized ddRAD protocol at LGC Genomics, GmbH. The library was size selected on a LMP-Agarose gel, removing fragments smaller than 300 bp and those larger than 500 bp. Sequencing was done on an Illumina NovaSeq 6000 (150bp paired-end read).

Genotypes and SNP filtering

Demultiplexing of all libraries for each sequencing lane was done using the Illumina bcl2fastq v2.20 software. Demultiplexing of library groups into samples was done according to their inline barcodes and verification of the restriction site. No mismatches or Ns were allowed in the inline barcodes, but Ns were allowed in the restriction site. Reads with final length < 20 bases were rejected and reads with 5’ ends not matching the restriction enzyme site were also discarded. The reads were quality trimmed at 3’-end to get a minimum average Phred quality score of 20 over a window of ten bases. The mapping of quality trimmed reads on the D. insularis reference genome v01.0 (available at Weedpedia, https://weedpedia.weedgenomics.org/) was done using BWA-MEM v0.7.12. One combined alignment file of all samples in the BAM format was used for variant discovery and genotyping of samples with Freebayes v1.2.0. Filtering of variants was done using the following GBS-specific rule set;

1.     The read count for a locus must exceed 8 reads

2.     Genotypes must have been observed in at least 66% of samples

3.     Minimum allele frequency across all samples must exceed 5%.

Genetic diversity and population differentiation

The genetic diversity indices expected heterozygosity (HE) and inbreeding coefficients (FIS) were calculated using the R packages ‘adegenet’ and ‘hierfstat’ (Jerome Goudet, 2005). The polymorphic information content (PIC) was calculated using an in-house R package. The two- and three-dimensional principal component analysis (PCA) was conducted using the R packages ‘stats’ and ‘rgl’. A Bayesian clustering approach implemented in the program STRUCTURE version 2.3.4 (Pritchard et al., 2000) was used to analyse population genetic structure by setting the replication number to 10,000 for the burn-in and Markov Chain Monte Carlo (MCMC) iterations each and using options of admixture model and correlated allele frequencies. The number of subpopulations i.e., K was set from 1 to 7 and three independent runs were performed for each K. The Structure Harvester (https://taylor0.biology.ucla.edu/structureHarvester/) was used to analyze the results from the STRUCTURE software, which constructs a deltaK vs K plot using the method of Evanno et al. (2005). The weighted neighbour joining (NJ) was constructed in DARwin 6.0. (Perrier and Jacquemoud-Collet, 2006).

Target site mutations in ACCase

Two primer pairs were used to amplify the ACCase gene sequence in D. insularis: FE35332 Forward (5′-ATGTCCACTCCTGAATTCCCA-3′), FE35333 Reverse (5′-CATTCTGAGGGAAGTATCAT-3′). PCR was performed in 25 μL reaction volume containing 5.0 µL of GoTaq Buffer, 0.5 µL of 10 mM dNTPs, 1.5 µL of 25 mM MgCl2, 0.5 µL of 10 µM of each forward and reverse primers, 0.2 µL of GoTaq G2 Hot Start Polymerase (Promega), and 14.8 µL of ultrapure nuclease-free water (Sigma). PCR cycling conditions were: one cycle at 95 °C for 2 min, 35 cycles at 94 °C for 30 s, 58 °C for 30 s, 72 °C for 90 s, and final extension at 72 °C for 10 min. PCR product was run on 1.0% agarose gel to verify the amplicon size of 1.5kb. The amplified samples were purified and sequenced in a Genetic Analyzer 3500 instrument (Applied Biosystems, Thermo Fisher) following the manufacturer’s instructions. Four individuals per population were sequenced using the original amplification primers and the following three internal sequencing forward primers: FE35334 Sequencing Forward 1 (5’- TGGGAGAGCAAAGCTTGGGGT-3’), FE35407 Sequencing Forward 2 (5’- GAAGTGCTGCTATTGCCAGTGC -3’) and FE35408 Sequencing Forward 3 (5’- GACCCACCAGACAGACCTGTTA -3’). The chromatograms were manually read using Bioedit version 7.2.5 software (Hall, 1999) to screen the seven known target-site mutations.

Genome-wide association study (GWAS) and in silico analysis of significant SNPs

The normality of the resistance data scores was checked in PAST3 program (Hammer et al., 2001). The resistance scores data of haloxyfop at a dose rate of 7.8 g a.i ha-1 and clethodim at a dose rate of 27.0 g a.i ha-1 showed normal and near-normal distribution, respectively, and hence were used in GWAS. Both the general linear model (GLM) and mixed linear model (MLM) were used for GWAS in TASSEL software ver 5.0 (Bradbury et al., 2007). The kinship matrix was calculated using VanRaden algorithm (VanRaden, 2008) in the GAPIT package 2.0 (Lipka et al., 2012). In the GLM, PCA was used as a fixed variate and in the MLM, PCA and kinship matrices were used as fixed and random variates, respectively. The threshold to declare significant marker-trait associations (MTA) was ≥10−3 (log10p) after applying a correction for a false discovery rate (FDR) at p < 0.05.

The VCF file was annotated with SnpEff version 5.1 using the IWGC Digitaria insularis reference genome and annotation v01.0. In addition, the in-silico analysis of the significant SNPs was conducted using nucleotide Basic Local Alignment Search Tool (BLAST) in the EnsemblPlants database (https://plants.ensembl.org/index.html). The EnsemblPlants database has cDNA/transcript sequences of more than 80 monocots and dicots and of model plants, which were used to find the homologies. The genes found in the overlapping region and within 1.0 Mb upstream and downstream of the matched regions were selected as candidate genes. To determine their molecular functions, the protein sequences of the candidate genes were downloaded from EnsemblPlants database and used in protein BLAST analysis in NCBI server (https://blast.ncbi.nlm.nih.gov/Blast.cgi) and their molecular functions were determined after ascertaining their homologies with known proteins in grasses. 

Genomic prediction models

Four parametric models (Ridge regression, Bayes A, Bayes B, Bayes C) and a non-parametric model (RKHS) were used and all these models are implemented in the ‘BGLR’ package (De Los Campos et al., 2022) in R. Ridge regression (RR) method considers common variance for all markers and shrinks the marker effects toward zero (Meuwissen et al., 2001). All Bayesian models do not consider the common variance of markers and incorporate additive genetic effects. The RKHS model uses a kernel function and captures non-additive effects (epistatic interactions). Either all available SNP markers were used or different sets of SNP markers were employed according to preliminary GWAS results. The SNP markers were ranked according to increasing p-values in GWAS analyses. The prediction accuracy of all models was calculated through “Pearson correlation coefficient” between observed and predicted values based on 100 iterations and 10-fold cross-validations.