ISU Comparative Genomics and Genome Size Evolution

SUMMARY

We recently initiated a large genome resequencing project (led by Josh Udall, with Daojun Yuan and Thiru Ramaraj) funded by the National Science Foundation Plant Genome Program.  In this project, we are addressing biological questions of broad relevance to polyploid plant genomes, while at the same time generating detailed information on genomic and phylogeographic patterns of diversity in wild and domesticated cotton, advancing our capacity to utilize an expanded gene pool for cotton improvement.

The two cultivated allopolyploid cottons, Gossypium hirsutum and G. barbadense, represent a remarkable example of parallel independent domestication, both involving dramatic morphological transformations under selection from wild perennial plants to annualized row crops. We combine deep resequencing of 643 newly sampled accessions spanning the wild-todomesticated continuum of both species, and their allopolyploid relatives, with existing data to resolve species relationships and elucidate multiple aspects of their parallel domestication. We confirm that wild G. hirsutum and G. barbadense were initially domesticated in the Yucatan Peninsula and NW South America, respectively, and subsequently spread under domestication over 4,000 - 8,000 years to encompass most of the American tropics. We present a robust phylogenomic analysis of infraspecific relationships in each species, quantify genetic diversity in both, and describe genetic bottlenecks associated with domestication and subsequent diffusion. As these species became sympatric over the last several millennia, pervasive genome-wide bidirectional introgression occurred, often with striking asymmetries involving the two coresident genomes of these allopolyploids. Diversity scans revealed genomic regions and genes unknowingly targeted during domestication and additional subgenomic asymmetries. These analyses provide a comprehensive depiction of the origin, divergence, and adaptation of cotton, and serve as a rich resource for cotton improvement.

OBJECTIVES

We have used whole genome sequencing to generate an average of 21.5X coverage for a total of 633 accessions of the 2.4 Gb cotton genome. We are using these data of polyploid cotton to address fundamental questions for polyploid plant genomes. For example, how is genetic diversity partitioned between the two subgenomes in an allopolyploid?  How have humans modified these two subgenomes, and are the effects of directional selection the same in both genomes? Have genes in the two genomes contributed equally to advanced, modern fiber phenotypes?

PARTICIPANTS

--info from: https://www.eeob.iastate.edu/faculty/wendel/projects/comparative-genomics-and-genome-size-evolution

PROTOCOLS

Plant Sampling and resequencing

643 accessions selected from the U.S. Cotton Germplasm Collection and our own collections with a focus on sampling the most diverse and broadly representative collection possible for exploring the specific questions about phylogenetic relationships and apportionment of diversity in wild and domesticated gene pools. Sampling included 442 accessions of G. hirsutum (AD1; 210 domesticated, 232 accessions spanning the wild‐to‐landrace continuum), 182 accessions of G. barbadense (AD2; 81 domesticated, 101 accessions spanning the wild‐to‐landrace continuum), and 19 other tetraploid samples (including accessions of G. tomentosum (AD3), G. mustelinum (AD4), G. darwinii (AD5), G. ekmanianum (AD6), and G. stephensii (AD7). Supporting Information; including Germplasm Resources Information Network (GRIN) accession numbers, where available). Seeds from each accession were planted in the greenhouse or field of Brigham Young University (BYU, Provo, Utah), Iowa State University (ISU, Ames, Iowa), and the Southern Plains Agricultural Research Center of USDA‐ARS (College Station, TX). Young leaves were collected and shipped to BYU and extracted using the Cetyl Trimethyl Ammonium Bromide (CTAB) method. High‐quality DNA was shipped to BGI or the BYU Sequence Center (BYUSC) for Illumina sequencing on either the X‐Ten or HiSeq 2500, respectively. Libraries were constructed using the Illumina PCR‐free method, as the PCR‐free library method has: 1) no GC content limit; 2) improved coverage across high and low GC regions; 3) no PCR‐induced bias prior to cluster generation; 4) increased maximum library insert size for increased complexity; and 5) low read‐duplication ratio relative to PCR‐based methods. All the reads were deposited in NCBI SRA with project accession no. PRJNA414461. An additional 789 samples from previously reported resequencing were also downloaded from the SRA for possible inclusion; however, many were later dropped due to insufficient sequencing depth (see below).

Read Mapping and SNP Detection

Illumina reads (PE150) were trimmed and filtered using SOAPnuke v1.6.0,[73] and cleaned reads were mapped to the reference genome of TM1[74] using BWA‐mem v0.7.16.[75] Because G. hirsutum is a polyploid species, only uniquely mapped reads were retained, which allowed alleles and/or homoeologs to be distinguished.[76] The output of BWA‐mem was converted to a sorted Binary Alignment Map (BAM) using SAMtools v1.7.[77] Only those reads with high mapping quality against the reference genome (‐F 4 ‐q 30) were kept. Following mapping, PCR duplicates were removed from each accession with Picard v1.123 (https://github.com/broadinstitute/picard/). Then, the RealignerTargetCreator algorithm from GATK v3.5[72] was used to identify candidate insertion/deletion (InDel) regions, and IndelRealigner was subsequently used to correct local misalignments.[72] Indels identified by GATK were retained for further analysis and the realigned BAM files were used to call SNPs.

SNPs were identified by two independent SNP callers, GATK and FreeBayes v1.1.0,[78] and the intersection of these were considered high‐confidence SNPs. The enormous size of the dataset required that the genome be divided into 500kb segments for most analyses. GATK was run as per the Best Practices for genomic resequencing.[79] For FreeBayes, putative SNPs were required to have a minimum of 10X coverage containing 3 or more reads supporting that variant, a base quality of at least 20, and a mapping quality greater than 30 (‐m 30 ‐q 20 –min‐coverage 10 ‐C 3).  Consensus SNPs from both programs were generated in GATK using “SelectVariants”, and the quality value of each SNP was recalibrated (“BaseRecalibrator” from GATK) to minimize false‐positives. Only bi‐allelic SNP were retained. HaplotypeCaller from GATK was used to generate a genomic variant call format file (gVCF) for each sample, which was used for joint genotyping among all samples (via “GenotypeGVCFs” from GATK). Raw SNPs were filtered with “QD < 2.0 || FS > 60.0 || MQ <40.0 || MQRankSum < −12.5 || ReadPosRankSum < −8.0 ”, and then all 500kb segments were merged together to calculate the number of sites with missing data (per sample). Nearly 30% of individuals (408) had a missing site rate greater than 25%, and were subsequently removed from further analyses (Note S2 and Figure S3, Supporting Information); most of these samples were derived from previously generated, lower coverage data downloaded from the SRA. After removing samples with >25% missing sites, the SNP sites were filtered from the remaining 1024 samples based on the sample missing rate (>25%) and minor allele frequency (5%). Only chromosomal SNPs were retained for downstream analyses. Preliminary annotation of variant effects was performed using SnpEff v4.3.[80] All scripts are available from https://github.com/Wendellab/BYUReseq.

Plant Material and Polymorphism Detection

643 genetically representative accessions to sequence were selected, including 442 G. hirsutum, 182 G. barbadense, and 19 other tetraploid cottons. SNPs were identified by two independent SNP callers, GATK and FreeBayes v1.1.0, and the intersection of these were considered high‐confidence SNPs.

Population Genetic Analysis

PCA was performed with the smartpca program embedded in the EIGENSOFT package v7.2.0;[81, 82] parameters are available at https://github.com/Wendellab/BYUReseq. Maximum likelihood phylogenetic reconstruction was performed by RAxML v8.1.15 with the substitution model “GTRCAT” (‐m GTRCAT ‐p 12345 ‐b 1000 ‐T 4 ‐N 1000)[83] using only 4D SNPs, and results were visualized with the iTOL v5.[84] While SNP‐based phylogenies were not the only approach for phylogenetic reconstruction, it is noted that the interspecific SNP‐based phylogeny here recapitulates what has been previously demonstrated by classical phylogenetic analysis and by the STRUCTURE analyses here. The variant file of 4D sites was converted to structure file format using PGDSpider v2.0.9[85] with default parameters. STRUCTURE v2.3.4[86, 87] was used to infer population clusters. Ten independent runs were performed for each K between 2 and 10, with a burn‐in period of 100 000 iterations followed by 200 000 iterations of the Markov Chain Monte Carlo algorithm. Evanno's ΔΚ was applied to identify the most probable groups (K) that best fit the data using Structure Harvester v0.6.94.[88] The program CLUMPP v1.1.2[89] was used to align 30 iterations of the optimal K to generate a consensus structure, and the output was visualized using DISTRUCT v1.1.[90] All analyses are available at https://github.com/Wendellab/BYUReseq.

Genetic Diversity and Detection of Selection Sweeps

Nucleotide diversity (π) is a measure of genetic variation, which is defined as the average number of nucleotide differences per site between any two DNA sequences chosen randomly from the sample population.[91] The level of genetic diversity was measured with VCFtools v0.1.13[92] using 100‐kb windows sliding 20 kb. Fixation index (Fst) is a measure of population differentiation, genetic distance, based on genetic polymorphism data.[93] While typically used to describe population structure within species, recent methods have also used Fst to describe variation among species.[94-97] The Fst value was measured using VCFtools v0.1.13[92] 100‐kb window with sliding steps of 20 kb. Candidate domestication‐sweep windows were identified as the top 5% genomic regions exhibiting the greatest reduction in diversity (πL/πc) values and the top 5% of regions with the greatest Fst between landrace and cultivar.

Detection of Introgression

Introgression was evaluated using a previously published method,[34] which produces results congruent with ABBA‐BABA tests while also permitting imputation of introgressed regions. Based on the results of PCA, phylogenetic tree, and population structure, 36 accessions of wild G. hirsutum and 46 accessions of G. barbadense were used to infer SNPs that distinguish the two species. Because G. stephensii and G. ekmanianum were both very close relatives of G. hirsutum,[35, 36, 43] representatives of these species were also used to represent the ancestral SNPs of G. hirsutum. Mapped reads from each species were used in conjunction with interSNP,[98, 99] which generates a species‐specific SNP index (see Note S5, Supporting Information, for more detail). Reads were separated into “hirsutum” or “barbadense” bam files by PolyCat,[98] which categorizes reads based on the SNPs present, regardless of species identity. A BED file containing position data for segments of contiguous reads (minimum size 500 bp and 10x coverage) was generated from each bam file using eflen.[99] Subsequently, adjacent segments with ends closer than 30 kb were merged via awk (code available from https://github.com/Wendellab/BYUReseq) and only segments of putative introgressed sequence were retained for each species. Genes located in putative introgressed regions were identified via Bedtools v2.27.1.[100] Gene and segment output was tabulated and parsed in R v3.6.3. Code for analyses is available from https://github.com/Wendellab/BYUReseq.

Gene Expression Analysis

RNA‐seq reads for multiple, bulked tissues were downloaded from NCBI (Project ID: PRJNA490626).[24] Raw reads were cleaned by SOAPnuke v1.5.2[73] and subsequently aligned to the reference TM‐1 genome[74] using STAR v2.7.1a.[101] Quantification of gene expression was performed with Cufflinks version v2.2.1.[102] To detect tissue‐dominant or tissue‐specific expression, an enrichment test was performed with TissueEnrich.[101, 103] Gene expression in wild and cultivated G. hirsutum fiber was also compared using previously generated results.[104, 105] Expression heatmaps were drawn using the online tool ClustVis.[104] Gene in putative regions of selection was cross‐referenced with fiber expression for G. hirsutum, and GO enrichment categories were identified for all genes in regions of selection using the R package cluster Profiler.[106] Functional annotations were derived from the release on CottonGen

--info from: Yuan D, et al. Parallel and intertwining threads of domestication in allopolyploid cotton. 2021.

PULICATION

Yuan D, Grover CE, Hu G, Pan M, Miller ER, Conover JL, Hunt SP, Udall JA, Wendel JF. Parallel and intertwining threads of domestication in allopolyploid cotton. Adv. Sci. 2021 Mar 15; 2003634:1-17.

DATA SOURCE