Online Methods Plant Material

Online Methods

Plant Material.

A total of 1,027 F₂ plants derived from the original okra (NC05AZ21) x normal (NC11-2100) cross were used to identify recombinants for fine mapping the L-D₁ locus. A 538-member diversity panel was used for association mapping and genome wide association studies (GWAS). This panel consisted of the 384-member cotton diversity panel¹, 84 photoperiod insensitive lines and 91 photoperiod-sensitive accessions, which were obtained from USDA Cotton Germplasm Collection (Supplementary Table 1). Two sets of isolines were used in fine mapping and/or in gene expression and VIGS studies, a BC₈ set that included all four leaf shapes in the Stoneville 213 background² and a BC₃ pair of normal and okra in the Stoneville 7A background³.

Morphometric Analysis.

Four leaves from each of the accessions in the cotton diversity panel were sampled from a field in Central Crops Research Station, Clayton, North Carolina mid-August 2015. Leaves were arranged on a scanner (Epson Workforce DS-50000, Suwa, Japan) with a ruler and the abaxial side of the leaf was scanned. Files were named by the order they were scanned and appended with a field number corresponding to genotype information. In ImageJ⁴, the “Make Binary”, “Fill Holes”, “Open”, “Close-”, and “Image Inverter” functions were used to convert leaves to binary, polished objects. Individual binary leaves were then manually selected using the “Wand” tool and copy and pasted into individual files named by genotype. Binary leaf silhouettes were converted to chain code using the program SHAPE^{5, 6}. The nef code file from SHAPE (.nef file) was then imported into the Momocs package in R^{7, 8, 9, 10} using the NEF2COE function. Individual leaf contours in the Coe object were assigned phenotype factor levels of normal, sub-okra, okra, and super-okra and harmonics isolated for subsequent analyses. The PC.contrib() and pca() functions in Momocs were used to visualize eigenleaves and perform Principal Component Analysis (PCA) (respectively) on harmonics and the morpho.space() function was used to visualize the morphospace⁹. Linear Discriminant Analysis (LDA) on harmonic coefficients was performed using the lda function in conjunction with the MASS package¹¹. The predict function (stats package) and table function (base package) were used to reallocate leaves by their predicted phenotypic class. R package ggplot2¹² was used for all data visualizations unless indicated otherwise.

Cryo SEM of Okra and WT SEMs.

For comparisons of plastochron 2 (P2), vegetative shoot apices were hand-dissected from four-week old normal, okra, and super-okra BC₈ isolines² to expose the shoot apical meristem and the two most-recently initiated leaf primordia: plastochron 1 (P1) and plastochron 2 (P2). Apices were affixed to SEM stubs using cryo-glue, frozen in liquid nitrogen, and viewed using a Hitachi TM-1000 tabletop scanning electron microscope. Image contrast adjustment and scale addition were done in Fiji (http://fiji.sc/)..

Association Mapping.

The association mapping population consisted of the 384 member cotton diversity panel¹, to which an additional 84 photoperiod insensitive lines with okra leaf shape and 91 photoperiod-sensitive accessions with okra leaf shape were added. The diversity panel and photoperiod insensitive lines were grown under summer field conditions in Clayton, NC while the 91 photoperiod-sensitive accessions were grown in 10-inch single pots in the greenhouse under short-day conditions. All plants were phenotyped as described previously¹³.

A total of 47 STS markers were designed from the ten-gene candidate region (Fig. 2c). Three were polymorphic and run on the association mapping panel. Additionally, two SNPs were converted into a Kompetitive Allele Specific PCR (KASP) assay and analyzed on the population at the Eastern Regional Small Grains Genotyping Laboratory (Raleigh, NC, USA). Marker locations are summarized in Supplementary Fig. 1 and primer sequences are provided in Supplementary Table 10. DNA isolation, genotyping using SSRs and STS markers were done as described previously¹³. All the primers pairs were synthesized by Integrated DNA Technologies (Ames, Iowa).

To verify that the associations between leaf shape and the candidate genes tested were not due to population structure, 149 multiallelic simple sequence repeat (SSR) markers distributed throughout the genome were also analyzed on the 384 line diversity panel as well as 42 of the additional lines with okra leaf shape from the larger set tested above. Multiallelic SSR genotypes were converted to numeric allele content scores (0,1, or 2) using the “Expand Multiallelic Genotypes” option in JMP Genomics version 8 (SAS Institute, Cary, NC, USA). Principal components analysis was used to estimate population structure of the diversity panel. After initial analysis of population structure, it was clear that 17 wild accessions were distinct outliers along the first principal component axis (which accounted for 23% of the marker variation). In addition, only one line exhibited sub-okra phenotype and one line exhibited super-okra phenotype. The wild accessions and sub- and super-okra types were excluded from further association analyses, resulting in 54 SSR markers being monomorphic in the remaining sample of lines. These SSRs were excluded from further analysis and principal components analysis was performed again on the remaining sample of 404 lines and 95 SSR markers, which included 36 okra types and 368 normal leaf shape types.

All markers were then tested for association with the binary trait okra versus normal leaf shape using a logistic regression model in the PROC LOGISTIC procedure in SAS software version 9.4 (SAS Institute, Cary, NC). Population structure was controlled by including the first three principal components (explaining 9% of variation in marker profiles) as covariates in the model. To deal with complete and quasi-complete separation observed in some markers we used the FIRTH option available in SAS software which produces finite parameter estimates by means of penalized maximum likelihood estimation¹⁴.

In the initial scan, several candidate gene region markers and several SSRs had significant associations with leaf shape. Inspection of the marker data revealed that genotypes at significant SSRs were correlated with the candidate gene region marker genotypes. Therefore, we performed a second scan of all SSR markers in which the most significant candidate gene variant (GhLS-STS1) was included as an additional covariate in the model.

Semi-quantitative Expression Analysis.

Total RNA was collected from three field-grown plants each of six varieties: NC05AZ21 (okra), NC11-2100 (normal), LA213-63 (normal recurrent parent), LA213 Sea Island Leaf (sub-okra), LA213 okra, and LA213 super-okra. Samples were taken approximately 90 days after planting. Leaves were taken at the earliest possible time they could reliably be distinguished from the shoot apical meristem without the help of any equipment. At this time-point, leaves were ~30-50mm in length from tip to base.

RNA was isolated using the Spectrum™ Plant Total RNA Kit (Sigma-Aldrich, St. Louis, MO, USA) following the manufacturer’s instructions. Total RNA was converted to cDNA using the ImProm-II™ Reverse Transcription System (Promega, Madison, WI, USA) per the manufacturer’s instructions. cDNA was then used as template in 50µL PCR reactions and visualized on 3% HiRes agarose (GeneMate BioExpress, Kaysville, UT, USA). GLYCERALDEHYDE 3-PHOSPHATE DEHYDROGENASE (GAPDH) was used as the reference gene. Primers used in semi-quantitative expression analysis are provided in Supplementary Table 11.

Real-Time Quantitative PCR Analysis of GhLMI1-like Gene Expression.

cDNA from LA213 isolines in the preceding section was used in 25µL Real-Time qPCR reactions with Power SYBR® Green PCR Master Mix (Applied Biosystems, Foster City, CA, USA) and the ABI 7300 Real-Time PCR System (Applied Biosystems). The positive control/reference gene was UBI14¹⁵. Technical replicates were run in triplicate. Ct values were analyzed using the ΔΔCt method. Fold changes and their standard deviations were plotted using Microsoft Excel. Primers used in RT-qPCR are provided in Supplementary Table 11.

Sequencing of GhLMI1-like Genes.

Sanger sequencing was used to obtain the genomic DNA sequence of GhLMI1-D1a and GhLMI1-D1b in 20 different tetraploid Gossypium accessions listed in Supplementary Table 6. Sequencing was performed at the North Carolina State University Genomic Sciences Laboratory (Raleigh, NC). Genome-specific primers for PCR were developed by aligning homeologous sequences from G. raimondii and G. arboreum and targeting differences between the two donor diploid genomes (Supplementary Table 12, 13). Genome-specificity of the primers was confirmed by analyzing amplification in a panel of diploid species from both genomes (Supplementary Table 14). Nested PCR was necessary for the 5’ end of GhLMI1-D1a. A complete list of the primers used for Sanger Sequencing can be found in Supplementary Table 15. Sequencing results were analyzed and assembled using Sequencher 5.2.3 (Gene Codes, Ann Arbor, MI, USA) with additional alignments performed using Clustal Omega Multiple Sequence Alignment (http://www.ebi.ac.uk/Tools/msa/). Publicly available predictions for Gorai.002G244200 and Gorai.002G244000 (Phytozome, www.phytozome.net) were used to determine the exon/intron structure of GhLMI1-like genes with expansions to a stop codon when necessary. Protein translations were performed using the ExPASy Translate Tool (http://web.expasy.org/translate/). Renderings of GhLMI1-like genes were drawn with fancyGENE (http://bio.ieo.eu/fancygene/) and redrawn to scale in Microsoft PowerPoint.

Virus-Induced Gene Silencing of GhLMI1-D1b in LA213 Okra.

All enzymes used in the construction of the TRV silencing vectors were supplied by New England Biolabs (Ipswich, MA, USA). Kits with ‘Zymo’ in the name were supplied by Zymo Research (Irvine, CA, USA). To construct TRV2: GhLMI1-D1b, a 461-bp fragment of the GhLMI1-D1b gene was amplified from cDNA derived from LA213 Okra. This 461-bp silencing fragment of GhLMI1-D1b included the last 29 bp of the second exon, the entire 239 bp of the third exon, and the first 193 bp of the proposed 3’ UTR. This fragment, along with the pYL156 (TRV2) vector (obtained from the Arabidopsis Biological Resource Center) was then digested with the restriction enzyme Acc65I overnight according to the manufacturer’s instructions. Following the digestion, the restriction enzyme was inactivated by placing the reactions at 65^oC for 20 minutes. Digested vector DNA was then de-phosphorylated with Antarctic Phosphatase per the manufacturer’s instructions. Both digested vector and GhLMI1-D1b fragment were separated on a 0.8% agarose gel, excised, and purified using the Zymoclean Gel DNA Recovery Kit and ligated using T4 DNA ligase. When combined with the TRV1 vector of the bipartite TRV VIGS system, this treatment was named TRV:GhLMI1-D1b.

In addition to the TRV:GhLMI1-D1b experimental treatment, two negative controls were used as before⁶. TRV:Mock consisted of only an empty TRV2 vector, which does not support the spread of infection due to the absence of TRV1 and thereby controls for effects of the inoculation process. TRV:GFP, consisting of TRV1 plus TRV2 containing a silencing fragment for GFP, is capable of spreading throughout the plant. However, cotton lacks an endogenous GFP gene so that potential effects of infection only could be monitored. Additionally a TRV:CHLI treatment that blocks chlorophyll production was used as a visible marker to ensure that environmental conditions were suitable for VIGS and to time the phenotyping of LMI knockdowns.

Two of the VIGS control constructs, TRV:GFP and TRV:ChlI, were produced by digesting the TRV2 plasmid pYL156 with Acc65I. The ends of the digested vector were blunted using DNA Polymerase I, Large (Klenow) Fragment and dephosphorylated with Antarctic Phosphatase. A 499-bp GFP fragment flanked with StuI restriction sites was amplified from transgenic cotton carrying the mGFP5-ER transgene¹⁶ using the primers mGFPerF: 5PerF: 5e primers mGFPerF: 52002rom transgenic cotton carrying the I, Large (KlenowA TGG TTG TCT GGT AAA AG -3A. The GFP PCR products were desalted using the Zymo DNA Clean and Concentrator kit The cleaned PCR product was then digested with StuI. A 501-bp blunt-ended fragment of the cotton ChlI gene was digested out of pJRT.CLCrVA.009¹⁷ using MscI, gel purified, and extracted from the agarose gel using the Zymoclean Gel DNA Recovery Kit. Blunt-ended GFP and ChlI gene fragments were ligated into the TRV2 vector using T4 DNA ligase.

All vector constructs were transformed into DH10-beta competent cells (New England Biolabs) and plated on Luria Bertani agar plates containing 50µg/ml each of kanamycin and gentamicin. Transformants were screened for insert direction using the primers: TRV2MCSF: 5’- CTT AGA TTC TGT GAG TAA GGT TAC C -3’ and mGFPerR2 5’ ATT AGG CCT AGG TAA TGG TTG CT GGT AAA AG 3’ for TRVGFP; and TRV2MCSF: 5’- CTT AGA TTC TGT GAG TAA GGT TAC C -3’and GhChlIR: 5’ GCT TGG CCA ATC AAA CCG TGC TCT TT -3’ for TRVChlI. Positive clones were confirmed by sequencing.

Two-week old okra seedlings were agro-inoculated as described previously¹⁸. In each experimental replicate, five plants were inoculated per treatment with TRV:Mock, TRV:GFP, and TRV:GhLMI1-D1b. Additionally, at least two plants in each replication were inoculated with TRV:CHLI. Plants were grown under a 26/22°C day/night cycle. Starting at three weeks post-inoculation all plants were photographed weekly. At four weeks post-inoculation, leaves were collected randomly from three of the five plants in the TRV:GhLMI1-D1b, TRV:GFP, and TRV:Mock treatments for expression analysis as described above. Three experimental replicates of the VIGS experiments were performed with consistent results. All primers used in the VIGS experiment are listed in Supplementary Table 11.

Phylogenetic Analysis of LMI1-like Genes in Gossypium.

In addition to sequencing GhLMI1-like genes in tetraploid cotton, LMI1-D1b was sequenced in lobed D-genome diploid species G. thurberi (PIs 530766 and 530789) and G. trilobum (PI 530967). Sequences of LMI1-A1b and LMI1-D1b were pulled, from G. arboreum and G. raimondii respectively, were obtained from www.cottongen.org. Alignment and phylogenetic analysis were performed using Clustal Omega. Helix-turn-helix prediction was carried out using https://npsa-prabi.ibcp.fr/cgi-bin/primanal_hth.pl. Leucine zipper was predicted using http://2zip.molgen.mpg.de/.

GhLMI1-D1b cDNA isolated from the okra leaf mapping population parent NC05AZ21 was sequenced via Sanger sequencing as described previously. Analysis, assembly, and alignment of the sequence were as described previously for the genomic sequencing of the LMI1-like genes. Extraction of RNA and conversion to cDNA are described under the semi-quantitative expression analysis section. Alignment of cDNA sequence to gDNA sequence confirmed the predicted exon/intron structure.

RNA-seq Transcriptome Analysis.

Sample Harvesting and RNA Preparation for RNA-seq: Three biological replicates of plastochron 2 (P2), corresponding to leaf 8 were hand dissected directly into ice-cold acetone from normal and okra shoot apices. Ten individuals were pooled for each biological replicate. Acetone was removed from the samples and replaced with extraction buffer from the PicoPure RNA Isolation Kit (ThermoFisher Scientific, Waltham, MA, USA). RNA was isolated following the manufacturer’s protocol with the optional on-column DNAse treatment. The RNA integrity was assessed by running the samples on an Agilent RNA 6000 Pico Chip (Agilent Technologies, Santa Clara, CA, USA). The Clonetech SMARTer cDNA synthesis kit (Clonetech, Mountain View, CA, USA) was used to amplify 10ng of total RNA into polyA tail-enriched dscDNA. 150ng of dscDNA was fragmented for 17 minutes with Fragmentase (New England Biolabs) and processed into Illumina sequencing libraries using the NEBNext Ultra DNA Library Prep Kit for Illumina (New England Biolabs). Illumina libraries were quantified with a Nanodrop and pooled to a final concentration of 20nM. The pooled libraries were sequenced for single-end 100bp reads on an Illumina HiSeq 2500 at the Washington University in St. Louis School of Medicine Genome Technology Access Center (https://gtac.wustl.edu/).

Bioinformatic Processing of RNA-seq Data: Illumina adapters and low quality bases were trimmed using Trimmomatic¹⁹ with the following default parameters: LEADING: 3 TRAILING: 3 SLIDINGWINDOW: 4:15 MINLEN: 36. Trimmed reads were aligned to the Gossypium hirsutum AD1_NBI genome (https://www.cottongen.org/data/download/genome²⁰) to produce Sequence Alignment/Map (SAM) files using Spliced Transcripts Alignment to a Reference (STAR) software²¹ with the following alignment parameters: --outFilterMatchNmin 16 --outFilterMismatchNoverLmax .05. SAM files were converted into compressed, sorted, and indexed Binary Alignment/Map (BAM) files using Samtools –view, -sort, and –index commands, respectively (http://samtools.sourceforge.net/²²). Reads mapping to annotated genes for the AD1_NBI genome were extracted using BEDtools multicov (http://bedtools.readthedocs.org/en/latest/content/tools/multicov.html²³). Genes with less than 1 count per million across at least 3 samples were discarded from the analysis. Significantly differentially expressed genes (FDR adjusted p-value ≤ 0.05) were identified by performing a pairwise comparison between normal and okra P2 samples in edgeR version 3.0 (https://bioconductor.org/packages/3.0/bioc/html/edgeR.html^{24, 25}).

GO Enrichment: The R package TopGO²⁶ was used to test for GO category enrichment in differentially expressed genes between okra and normal P2 transcriptomes. The Fisher’s Exact Test (p-value ≤ 0.05) was used to identify significantly enriched GO categories in gene sets that are significantly differentially expressed between okra and normal P2 samples.

Co-localization of GhLMI1-D1b^Okra-GFP in Nicotiana benthamiana.

Fluorescent protein fusions of GhLMI1-D1b were generated by Gateway cloning (Life Technologies, Carlsbad, CA, USA). The GhLMI1-D1b^Okra coding sequence without the STOP codon was amplified from okra cDNA using primers GhLMI1-D1b-Okra-TOPO-F (5’ – CACCATGGATTGGGATGGCACCATTCGACCCTTT - 3’) and GhLMI1-D1b-Okra-STOP-R (5’ – GGGATAAGAAGGGAGTTGAA - 3’), and cloned into pENTR/D/TOPO (Life Technologies) to generate pENTR::GhLMI1-D1b-Okra-stop. LR recombination was carried out between pENTR::GhLMI1-D1b-Okra-stop and the following destination vectors: pGWB5²⁷, pGWB8²⁵, and pUBQ10-C-GFP. pUBQ10-C-GFP is a modified version of pUBC-GFP²⁸ but includes the full pUBQ10 promoter. The following constructs were obtained and confirmed by sequencing: 35S::GhLMI1-D1b-GFP and pUBQ10::GhLMI1-D1b-GFP. Both constructs were introduced by transient Agrobacterium transformation into Nicotiana benthamiana plants carrying a RFP-Histone2B marker²⁹ for co-localization analysis.

A Zeiss LSM 710 confocal microscope with a 40x water objective (1.1 N.A.) was used to image fluorescence protein fusions. The excitation/emission wavelengths during acquisition were 488nm/492–570nm for GFP and 561nm/588–696nm for RFP.


References
1. Tyagi, P. et al. Genetic diversity and population structure in the US Upland cotton (Gossypium hirsutum L.). Theor. Appl. Genet. 127, 283-295 (2014).
2. Kennedy, C.W., Smith, W.C. & Jones, J.E. Effect of early season square removal on three leaf types of cotton. Crop Sci. 26, 139-145 (1986).
3. Andries, J.A., Jones, J.E., Sloane, L.W. & Marshall, J.G. Effects of okra leaf shape on boll rot, yield, and other important characters of Upland cotton, Gossypium hirsutum L. Crop Sci. 9, 705-710 (1969).
4. Abràmoff, M.D., Magalhães, P.J., & Ram, S.J. Image processing with ImageJ. Biophotonics Int. 11(7), 36-42 (2004).
5. Iwata, H., Niikura S., Matsuura S., Takano Y. & Ukai Y. Evaluation of variation of root shape of Japanese radish (Raphanus sativus L.) based on image analysis using elliptic Fourier descriptors. J. Hered. 93, 384-385 (1998).
6. Iwata, H. & Ukai Y. SHAPE: A computer program package for quantitative evaluation of biological shapes based on elliptic Fourier descriptors. Euphytica 102, 143-149 (2002).
7. Claude, J. Morphometrics with R. Springer, New York (2008).
8. Claude, J. Log-shape ratios, Procrustes superimposition, elliptic Fourier analysis: three worked examples in R. Hystrix. 24, 94-102 (2013).
9. Bonhomme, V., Picq S., Gaucherel C. & Claude J. Momocs: Outline analysis using R. J. Stat. Softw. 56(13), 1-24 (2014).
10. R Development Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.
11. Venables, W.N. & Ripley, B.D. Modern Applied Statistics with S. Fourth Edition. Springer, New York (2002). ISBN 0-387-95457-0.
12. Wickham, H. ggplot2: elegant graphics for data analysis. Springer Science & Business Media (2009).
13. Andres, R.J., Bowman, D.T., Kaur, B. & Kuraparthy, V. Mapping and genomic targeting of the major leaf shape gene (L) in Upland cotton (Gossypium hirsutum L.). Theor. Appl. Genet. 127, 167-177 (2014).
14. Firth, D. Bias reduction of maximum likelihood estimates. Biometrika 80 (1): 27-38 (1993).

15. Artico, S., Nardeli, S.M., Brilhante, O., Grossi-de-Sa, M.F. & Alves-Ferreira, M. Identification and evaluation of new reference genes in Gossypium hirsutum for accurate normalization of real-time quantitative RT-PCR data. BMC Plant Biology 10, 49 (2010).

16. Sunilkumar, G., Mohr, L., Lopata-Finch, E., Emani, C. & Rathore, K.S. Developmental and tissue specific expression of CaMV 35S promoter in cotton as revealed by GFP. Plant Mol. Biol. 50, 463-479 (2002).
17. Tuttle, J.R., Idris, A.M., Brown, J.K., Haigler, C.H. & Robertson, D. Geminivirus-mediated gene silencing from cotton leaf crumple virus is enhanced by low temperature in cotton. Plant Physiol. 148, 41-50 (2008).
18. 16. Tuttle, J.R., Haigler, C.H. & Roberston, D. Method: low-cost delivery of the cotton leaf crumple virus-induced gene silencing system. Plant Methods 8, 27 (2012).
19. 17. Bolger, A.M., Lohse, M.& Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 30, 2114-2120 (2014).
20. 18. Zhang, T. et al. Sequencing of allotetraploid cotton (Gossypium hirsutum L. acc. TM-1) provides a resource for fiber improvement. Nat. Biotechnol. 33, 531-537 (2014).
21. Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15-21 (2013).
22. Li H. et al. The Sequence alignment/map (SAM) format and SAMtools. Bioinformatics 25, 2078-2079 (2009).
23. Quinlan, A.R. & Hall, I.M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 15, 841-842 (2010).
24. Robinson, M.D., McCarthy, D.J. & Smyth, G.K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139-140 (2010).
25. Zhou, X., Lindsay, H. & Robinson, M.D. Robustly detecting differential expression in RNA sequencing data using observation weights. Nucleic Acids Res. 42, e91 (2014).
26. Alexa, A. & Rahnenführer, J. Gene set enrichment analysis with topGO. (2016) https://www.bioconductor.org/packages/3.3/bioc/vignettes/topGO/inst/doc/topGO.pdf
27. Nakagawa, T. et al. Development of series of gateway binary vectors, pGWBs, for realizing efficient construction of fusion genes for plant transformation. J Biosci. Bioeng. 104, 34-41 (2007).
28. Grefen, C., Donald, N., Schumacher, K., Blatt, M.R. A Ubiquitin-10 promoter-based vector set for fluorescent protein tagging facilitates temporal stability and native protein distribution in transient and stable expression studies. Plant J. 64, 355-365 (2010).
29. Martin, K. et al. Transient expression in Nicotiana benthamiana fluorescent marker lines provides enhanced definition of protein localization, movement and interactions in planta. Plant J. 59, 150-162 (2009).