Pollen tubes shown in Figure 1 were visualized following silk fixation using a combination of methods performed by Lausser et al. (2010) and Zhang et al. (2012). Briefly, pollen tubes were allowed to germinate for a predetermined time before silk tissue was harvested and incubated in FAA (10:85:5 v/v/v formaldehyde: 95% ethanol: acetic acid) for 24 h at 4°C. Silks were then slowly rehydrated by washing in 75, 50, 25% ethanol for 3 min each. Samples were washed in 0.1 M potassium phosphate buffer (pH 8.0) before a 2 h incubation in 8 M NaOH. Silk tissue was stained in 0.1% aniline blue (Fisher Scientific, Pittsburgh, PA) dissolved in 0.1 potassium phosphate overnight at 4°C. Silks were placed on glass slides and pollen tubes were visualized at near UV excitation using a fluorescent microscope (Nikon Diaphot 300). At least 10 silks per plant were measured and three biological replicates were completed for each experiment. In transmitting tracts with more than one pollen tube germinated, measurements were recorded from the pollen grain closest to the ovule.

W22 (ga1) and Ga1-s (W22) (Maize Genetic Stock Center number 401D) plants were grown in the greenhouse at Iowa State University. Three biological replicates of unpollinated, emerged silks from each genotype were collected for RNA isolation. Following isolation, RNA was sent to the DNA facility at Iowa State University for single-end RNA sequencing on Illumina MiSeq with a read length of 50 bp.

For qRT-PCR analysis, W22 (ga1), Ga1-s (W22) and Hp301 were grown in the USDA greenhouse at Iowa State University. Five to seven biological replicates of unpollinated, emerged silks from each genotype were collected for RNA isolation by Qiagen RNeasy Plant kit.

Sequencing was performed at the Iowa State University on a MiSeq. Libraries were sequenced for 51 cycles. TopHat (version 2.0.3, Trapnell et al., 2009) was used to align reads to the B73 reference genome (version 3 ref). The program samtools (Li et al., 2009) was used to remove unreliably mapped reads. The resulting mapping files (bam) were imported into the statistical program R (R Core Team, 2016) using the Bioconductor package Rsamtools (Morgan et al., 2010). Mapped reads were submitted to NCBI SRA under BioProject PRJNA382364.

The R graphics program ggplot2 (Wickham, 2009) was used to compare sample replicates for technical reproducibility (data not shown). Samples of low quality were not used in further analyses. The Bioconductor package edgeR (Robinson and Smyth, 2007, 2008; Robinson et al., 2010; McCarthy et al., 2012) was used for single factor, pairwise comparisons to calculate normalization factors, estimate tagwise dispersion and determine differential expression (DEG) based in comparison to the reference genome B73_v3. 1452 DEG were identified at a cutoff of p > 0.05. (Table S3).

RNA-Seq reads were also assembled de novo through the Trinity package. RSEM was used to estimate abundance after mapping reads back to assembled transcripts using bowtie2. The raw counts for transcripts were then normalized (upper quartile) and differential gene expression (DEG) analysis was performed (three independent methods were used: edgeR, voom, DESeq2). EdgeR results are presented in this paper. ORFs were predicted from these transcripts and predicted peptides were generated. Annotation of the assembled transcripts was performed using Trinnotate which uses a curated dataset from swissprot to predict function as well as the pfam database, tmhmm, signalp, and gene ontology for additional annotation. Some sequences that were reported as differentially expressed most likely resulted from a sequence polymorphism between the two genotypes sequenced because a corresponding transcript with minor sequence differences and approximately the same read count was also found on this list. These transcripts were manually removed from the gene list shown in Table 1, but are still present in Table S1.

Equal amounts (0.5 ng) of total RNA were used in a one step RT/qPCR reaction with SuperScript® III Platinum® SYBR® Green One-Step qRT-PCR kit (Invitrogen #11736-059). 18S was used as an internal control and ZmPme3 primers (Table S4) were used to detect transcripts. ZmPme3 ΔCt levels were normalized to 18S RNA levels and ZmPme3 transcript levels in Ga1-s and Hp301 are expressed as relative to those of W22.

PCR primers were designed using the de novo assembled transcript sequence to amplify the corresponding genomic sequence in four pieces (Table S4). Using genomic DNA from the maize genetic stock center stock 401D which carries the genotype Ga1-s/Ga1-s in a W22 background, Hp301 which is a popcorn line known to be Ga1-s/Ga1-s(obtained from the North Central Region Plant Introduction Station), W22 (ga1/ga1)(obtained from the North Central Region Plant Introduction Station) and NC390/NC394(Ga1-m/Ga1-m)(provided by Dr. M. Goodman, NCSU), products were amplified and sequenced in both directions at the Iowa State University DNA Facility. NAM B73 × Hp301 seeds were obtained from the Maize Genetics Stock Center (maizecoop.cropsci.uiuc.edu/).

The NAM lines used to map the ga1 locus with recombination breakpoints in the region of Ga1-s (Bloom and Holland, 2011) were screened for the presence of the active PME gene sequence by PCR (Table 2). NAM B73 × Hp301 seeds were obtained from the Maize Genetics Stock Center (maizecoop.cropsci.uiuc.edu/). Two primer pairs (PME marker and PME_A, Table S4) which distinguish the active ZmPme3 from the inactive genes in W22, were used. The PCR results were compared to the published pollen exclusion phenotypes of these lines (Bloom and Holland, 2011).

To visualize the structure of the W22 PME repeat region, the 3 Mb portion of the Chromosome 4 genomic DNA sequence containing regions of PME homology was analyzed using YASS (Noé and Kucherov, 2005) by plotting the sequence against itself (Figure 3). Genomic sequences of each PME pseudo gene were extracted from the genomic sequence using the genome coordinates produced by YASS. The resulting PME pseudogene sequences were aligned using BioEdit and Clustal W. The genome coordinates of the aligned sequences was used to produce Figure 4.

Pollen protein extract was produced from 20 mg of mature pollen grains by mortar and pestle (Zhu et al., 2011). Pistil (pre- and post-pollination) protein extract was produced by grinding tissue in liquid nitrogen using a mortar and pestle into a fine powder that was extracted with trichloroacetic acid (Méchin et al., 2007). Pollen and pistil protein extracts were analyzed by the Protein Facility at Iowa State University as a fee-for-service. The isotopic label-free relative quantification method used was as follows; 100 μg of each protein extract was digested in solution using a trypsin/Lys-C protease mix. Samples were dried down and spiked with 250 fmol of peptide retention time calibration (PRTC) standards during reconstitution. PRTC standards contained an equimolar mixture of 15 known peptides to allow for the quantification of unknown peptides. Newly digested peptides were then separated by liquid chromatography before analysis by MS/MS using a Q Exactive Hybrid Quadrupole-Orbitrap Mass Spectrometer (Fisher Scientific, Waltham, MA). The resulting spectral data was compared to theoretical fragmentation patterns using the MASCOT search engine (Matrix Science, London, UK) to identify the most likely protein accessions by considering the highest-scoring peptide sequence for each input fragment and alignment of 2 or more high-scoring peptide sequences to a single accession. A trinity-predicted peptide database from silk transcript de novo assembly was used to make protein identifications.

To confirm pollen tube growth arrest (Lu et al., 2014) and slow growth rate (Zhang et al., 2012) characteristic of ga1 pollen on Ga1-s silks (Lu et al., 2014) in our genetic materials, pollen tubes from compatible and incompatible pollinations were visualized. In compatible pollinations, pollen tubes were straight and grew within the transmitting tract of the silk to the ovule. In incompatible pollinations of ga1 (W22) pollen on Ga1-s (W22) silks, the pollen tubes had a slower rate of growth (Zhang et al., 2012), were not straight and often left the transmitting tract while growing (Figure 1C).

To identify candidate genes for the Ga1-s allele, we sequenced transcripts from the silks of two near isogenic lines differing at the ga1 locus (Ga1-s/Ga1-s vs ga1/ga1) in the inbred background W22. The sequencing reads generated by this study were deposited in the National Center for Biotechnology Short Read Archive (NCBI SRA BioProject accession PRJNA382364). There were 65,780,803 reads from the ga1/ga1 sample and 65,261,178 reads from the Ga1-s/Ga1-s sample. Since genetic experiments suggest that ga1 may be a null allele, it could be missing from ga1 genotypes such as W22 and B73 and mapping reads to these genotypes would not be expected to identify transcripts from the Ga1-s gene. We therefore processed the transcript data by de novo assembly of the pooled RNA reads. There were 69,160 Trinity assembled transcripts and 61,891 “genes” identified. We then identified differentially expressed genes (DEGs) using edgeR with a FDR of < 0.001. At that cutoff, there are 1958 significant DEG, with 1005 of them upregulated in Ga1-s silks (Table 1, Table S1). We used Trinnotate and BLAST to identify the assembled transcripts where 706 transcripts had significant BLAST hits, leaving 1252 without identification.

Keeping with the hypothesis that the Ga1-s transcripts are missing from ga1 genotypes, we focused on genes that were upregulated in Ga1-s. The most significantly differentially expressed (based on FDR) of the genes upregulated in Ga1-s is a putative pectinesterase/pectinesterase inhibitor 38. Among the de novo assembled transcripts, there is only one Pme38-like transcript and the reads used to assemble it came overwhelmingly from Ga1-s samples (Supplemental File). A BLAST search was performed against the B73 genome (RefGen3) with the 40 most significant DEG that were up regulated in Ga1-s to determine potential genome localizations. Only the pectinesterase gene hit within the ga1 locus on chromosome 4. A hevein-like preproprotein resides just outside the locus.

Two genes identified in the B73 genome that lie within the 100 kb region of interest on chromosome 4 such as GRMZM2G027021 (GTP-binding protein hflX) and GRMZM2G039983 (Wave-dampened 2) were detected in the de novo assembly (Trinity_13706 and Trinity_17477). Both genes appear to be expressed in silks but not differentially. There do not appear to be genotype differences in transcript sequence for Trinity_17477 (Wave-dampened 2) despite SNPs for GRMZM2G039983 genomic sequence between Ga1-s, Ga1-m and ga1 (data not shown).

In order to confirm that the de novo assembled PME transcript is encoded by a gene present in Ga1-s plants and is not a bioinformatic artifact, we used PCR to amplify genomic DNA from Ga1-s(W22) and Hp301 (a popcorn line known to be Ga1-s/Ga1-s). We amplified a genomic DNA sequence consistent with the computationally derived transcript (Figure 2A). The same PCR reactions were done in ga1(W22) and while amplification products were produced, they were not identical in sequence to Ga1-s and Hp301 and did not encode functional PMEs. Thus, the hits from the BLAST search with the PME transcript against the B73 genome were likely to be sequences with homology to the PME within the chromosome 4 region of interest. The PME coding sequence from Ga1-s(W22) and Hp301 is 1578 bases long with a 99-bp intron and is predicted to encode transcripts identical to the sequence from the de novo assembly. We designated the gene encoding this transcript ZmPme3. According to the literature, the coding sequence of PME genes range from 1.2 to 3.0 kb long (Wang et al., 2013). We used Hmmer to search for Pfam families within the predicted protein sequence and found only the PME domain (PF01095) with a signal peptide.

We further confirmed the expression data by qRT-PCR. RNA was isolated from unpollinated, exposed silks from W22, Ga1-s(W22) and Hp301. Using two primer pairs for ZmPme3, relative expression levels were 17- to 97-fold higher in Ga1-s genotypes than W22 (Figure S1). ΔCt values for ZmPme3 from W22 samples were near no template control levels as would be expected if this genotype lacks the ZmPMme3 gene (Figure S1).

While we identified sequences homologous to ZmPme3 in the genomic region known to contain ga1 using BLAST, the genome position of ZmPme3 was not yet clear. To determine whether ZmPme3 maps to the ga1 locus, we mapped ZmPme3 in the Hp301 × B73 RIL population developed for the NAM project (McMullen et al., 2009). Based on previously published information about this population (Bloom and Holland, 2011) we screened the lines used to map the ga1 locus with recombination breakpoints in the region of Ga1-s for the presence of the active PME gene sequence by PCR. ZmPme3 was detected in the NAM lines that excluded pollen and was not detected in lines that did not exclude pollen (Bloom and Holland, 2011; Table 2). These data confirmed that ZmPme3 from Ga1-s is in the ga1 interval that Bloom and Holland identified.

To determine if the ZmPME3 protein showed an expression pattern similar to the ZmPme3 transcript, protein extracts from mature pollen and silk tissue pre- and post-pollination were subjected to mass spectrometry analysis. Six samples were analyzed, all in the W22 background; Ga1-s and ga1 mature pollen, Ga1-s and ga1 silks, ga1 silks pollinated with Ga1-s pollen (compatible) and Ga1-s silks pollinated with ga1 pollen (incompatible). Identification of the peptides generated was performed by matching resulting peptide masses to the predicted proteins from the de novo assembled transcripts of ZmPme3. Peptides from ZmPME3 (Trinity_DN41247) were detected in both Ga1-s unpollinated and pollinated silks at a coverage of 7.8% (Figure 2B). Peptides from ZmPME3 were detected in neither of the pollen samples nor the ga1 unpollinated and pollinated silks. We conclude that the ZmPME3 protein expression is consistent with the expression of the ZmPme3 transcript.

The Ga1-m allele lacks the female function of the ga1 locus, so the sequence of ZmPme3 in this genotype was of great interest. If ZmPme3 is responsible for the female function of Ga1-s, we would expect it to be non-functional in Ga1-m genotypes. Through PCR and sequencing of ZmPme3 in the Ga1-m/Ga1-m containing hybrid NC390/NC394, we found a 2 base insertion that causes a frameshift resulting in a truncated and likely inactive protein (Figure 2A). This observation is consistent with the hypothesis that ZmPme3 is involved in the female function of Ga1-s.

The transcript profiling with de novo assembly and proteomics experiments presented here show that the Ga1-s genotypes examined have a single active differentially expressed PME gene, however we observed BLAST hits to PME sequences in the ga1 genotypes as well. Our initial observation of PME-like sequences in a ga1 genotype was when RNA reads were mapped to B73_RefGen v3 (Schnable et al., 2009). The top 40 most significant (lowest FDR) upregulated genes are all annotated in AGPv3 as transposable elements, with 34 out of 40 on chromosome 4 and within a 2 Mb region containing the ga1 locus but outside of the smallest 100 kb mapped region identified by Liu et al. (2014). To characterize these gene models further, homology was examined by BLAST. We examined the 34 genes from the ga1 locus on chromosome 4 and all but one have pectinesterase/pectinesterase inhibitor 38, 17, or 46 as their top BLAST hit with the one remaining sequence having pectin methylesterase (PME) 63 as its top BLAST hit. The remaining six genes fall on other chromosomes and either have no significant BLAST hits, are related to kinesin (GRMZM2G327923) or nucleobase-ascorbate transporter 12 (GRMZM2G031728).

We compared the de novo assembled transcript to the B73 genomic sequence and the co-linear W22 sequence in the genetically mapped region of interest. In doing this, we found that this region contains a cluster of 58 PME-like sequences. In a 1.1 MB region of W22 there are 58 full length and partial PME genes (Figure 4; Table S2). All but three have pectinesterase/pectinesterase inhibitor 38/17 as their top BLAST hit. The remaining three have Pme63 as their top hit, but appear to be partial genes. By aligning the 58 PME sequences from W22 and the assembled transcript, there does appear to be an intron in the W22 sequences. Although similar to the differentially expressed transcript, none of the B73 or W22 sequences matched it exactly. Using the ExPASy.org translate tool, we found that none of these genes have a full-length open reading frame, even accounting for an intron, while the de novo assembled transcript does. In several cases, a simple one base difference introduces a stop codon in the W22 sequence. All PME-like sequences in the gene cluster are oriented in the same direction. The 1.1 MB region containing the gene cluster has a complex organization of nested repeats (Figure 3), however the PME-like sequences are dispersed throughout the region with no obvious association with a larger repeat unit. Alignment of the PME-like sequences according to their genome position (Figure 4) shows that partial PME sequences are interspersed with complete sequences. In some cases the sequence organization suggests that partial sequences are derived from insertions into complete PME sequences because partial sequences adjacent to and on either side of some PME-like sequences are contiguous.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.