extfieldsofstudy:
Biology
Computer Science

text:
CoCo: RNA-seq read assignment correction for nested genes and multimapped reads

Abstract Motivation Next-generation sequencing techniques revolutionized the study of RNA expression by permitting whole transcriptome analysis. However, sequencing reads generated from nested and multi-copy genes are often either misassigned or discarded, which greatly reduces both quantification accuracy and gene coverage. Results Here we present count corrector (CoCo), a read assignment pipeline that takes into account the multitude of overlapping and repetitive genes in the transcriptome of higher eukaryotes. CoCo uses a modified annotation file that highlights nested genes and proportionally distributes multimapped reads between repeated sequences. CoCo salvages over 15% of discarded aligned RNA-seq reads and significantly changes the abundance estimates for both coding and non-coding RNA as validated by PCR and bedgraph comparisons. Availability and implementation The CoCo software is an open source package written in Python and available from http://gitlabscottgroup.med.usherbrooke.ca/scott-group/coco. Supplementary information Supplementary data are available at Bioinformatics online.

Introduction
Detection and quantification of RNA transcripts are critical steps to understand the mechanism of gene expression and its impact on cell function. Traditionally, transcript abundance has been evaluated using techniques that target one known RNA sequence at a time, as in the case of quantitative RT-PCR. More recently, the development of RNA-sequencing (RNA-seq) techniques revolutionized transcriptome analysis by providing the tools necessary to study, at least in theory, all RNA transcripts simultaneously. Diverse library preparation protocols exist, the most commonly used ones focussing on particular classes of RNA through enrichment steps. Such strategies include polyA enrichment, non-ribosomal RNA (rRNA) enrichment (e.g. rRNA depletion), small RNA enrichment and enrichment for RNAs bound to specific factors (Conesa et al., 2016;Hrdlickova et al., 2017;O'Neil et al., 2013). All these protocols detect certain levels of non-coding RNA. Strategies that employ rRNA depletion are the best approaches to detect a wide range of different classes of both coding and non-coding RNAs including lncRNAs (long noncoding RNAs), small nuclear RNAs (snRNAs) and 7SL RNA (Boivin et al., 2018;Lai et al., 2016). However, no matter how the sequencing library is created, the capacity to correctly quantify RNA abundance ultimately depends on the correct assignment of the sequencing reads, including those generated by non-coding RNA.
Accurate RNA quantification is not easy to achieve since it depends on the quality of the transcriptome annotation used as a reference, the complexity of the target RNA sequence and its genomic context. RNA quantification is particularly difficult in the case of small non-coding RNA (<200 nt) which are often produced from multiple genes and/or nested in other genes (Boivin et al., 2018;Luo and Li, 2007;Mohammed et al., 2014;Weber, 2006). According to Ensembl annotations (Yates et al., 2016), in human, 2596 genes, including miRNA, snRNA, snoRNA (small nucleolar RNAs) and tRNA (transfer RNAs) are located in introns and many overlap exons (Supplementary Table S1). In total, 1838 protein-coding genes, lncRNA and pseudogenes host or overlap smaller non-coding RNA (Supplementary Table S2). Although some non-coding RNA are expressed from independent promoters like tRNA (Paule and White, 2000) many others do not have independent promoters and are at least partially linked to the expression of their host gene (Boivin et al., 2018;Filipowicz and Pogacic, 2002;Matera et al., 2007). Regardless of the mode of expression, correct identification of reads from nested and multi-copy genes is essential for the accurate quantification of both coding and non-coding RNA.
Abundance estimation from RNA-seq data generally requires a pipeline that aligns reads to a reference genome and then assigns them to annotated genes (Conesa et al., 2016). In higher eukaryotes, mapping RNA reads to a genome requires a gapped alignment to avoid the DNA intronic sequences, which is accurately performed by standard software like the splice-aware aligners STAR (Dobin and Gingeras, 2016) and HISAT (Kim et al., 2015). Once aligned, the sequencing reads must be assigned to genes, a task made challenging by overlapping genes and multimapped reads. Indeed, many available tools and commonly used settings only quantify correctly genes generating uniquely mapped reads, and reads mapping to more than one locus are typically discarded (Fig. 1). Reads originating from non-coding RNAs overlapping exons (mainly retained introns) of longer protein-coding RNA or lncRNA are either wrongly assigned to their host gene or labelled as ambiguous and discarded (Fig. 1). In addition, reads originating from duplicated genes are often discarded by default and many such genes, that we term multimapped genes, are under-represented in RNA-seq assigned counts.
To increase the overall accuracy of RNA quantification and monitor the expression pattern of overlapping and repetitive genes, we developed a count corrector (CoCo) pipeline that rescues and correctly assigns otherwise ambiguous sequencing reads. CoCo employs the read assignment function featureCounts from Subread (Liao et al., 2013), providing it with a modified annotation file that highlights nested genes, while ensuring appropriate distribution of multimapped reads. CoCo salvages and/or reassigns over 15% of aligned RNA-seq reads, significantly changing the abundance estimates for several classes of RNA and providing insight into the expression dynamics of often ignored classes of repetitive and overlapping genes. We investigated the performance of the main available read assignment pipelines, comparing their accuracy to digital PCR abundance values and bedgraph estimates, showing that CoCo performs best using these metrics while requiring reasonable runtime and memory.

Sequencing datasets
To ensure representative abundance of transcripts from nested and multimapped genes, most of which are highly structured RNAs, we chose thermostable Group II intron reverse transcriptase sequencing (TGIRT-seq). TGIRT displays high processivity and fidelity, providing an accurate view of the ribodepleted transcriptome (Boivin et al., 2018;Nottingham et al., 2016;Qin et al., 2016). The samples considered (GEO series GSE99065) are GSM2631741, GSM2631742 (fragmented) and GSM2631743, GSM2631744 (not fragmented). Although fragmented datasets enable the study of the whole transcriptome, non-fragmented datasets emulate sizeselecting short RNAs (Boivin et al., 2018). To evaluate the read assignment performance of CoCo on RNA-seq datasets generated using different library preparation methodologies, the following SRA datasets were also considered: SRR191548, SRR191549 (small RNA-seq; Farazi et al., 2011), SRR5312294, SRR5312295 (hydro- Fig. 1. CoCo read correction scheme for nested and multimapped genes. (A) Representation of a standard gene annotation used for depicting a genetic locus containing one host gene and three nested genes. The dashed lines indicate introns while the dark blue boxes indicate exons. (B) Representation of the gene annotation produced using the correct_annotation module of CoCo showing a gap in the retained intron over the first nested gene. (C) Examples of potential read pairs overlapping the different features and multimapped read pairs. (D) Comparison of the read pair assignment using standard and CoCo pipelines, for each of the read pairs illustrated in (C). The reads that are differentially assigned by CoCo are highlighted in red. (E) Comparison of the read count estimates by the standard and the CoCo pipelines, based on the assignments listed in (D). (F) Flow chart of the CoCo pipeline. Pre-processing and alignment steps are shown before the cor-rect_count module application. The correct_count module then assigns reads with Subread's featureCounts using the gapped CoCo annotation (built with the cor-rect_annotation module). Read pairs resulting in multiple alignments are considered separately and distributed proportionally to the uniquely assigned read pairs. (Color version of this figure is available at Bioinformatics online.) tRNAseq; Gogakos et al., 2017) and SRR2012208, SRR2012209 (standard viral RT ribodepleted RNA-seq; Choy et al., 2015).

Read alignment
FastQC was used to check Fastq files for quality. Reads were trimmed using Cutadapt (Martin, 2011) (with parameters -g GATC GTCGGACTGTAGAACTCTGAACGTGTAGATCTCGGTGGTCG CCGTATCATT -a AGATCGGAAGAGCACACGTCTGAACTCCA GTCACATCACGATCTCGTATGCCGTCTTCTGCTTG -minimumlength 2) and Trimmomatic (Bolger et al., 2014) (with TRAILING: 30) to remove adaptors and portions of reads of low quality, respectively. Read pairs were then aligned to the human genome build hg38 using an annotation file from Ensembl (described below) with the splice-aware RNA-seq aligner STAR (Dobin and Gingeras, 2016) using -alignIntronMax 1250000, all other parameters at default values. Reads not aligned using STAR were aligned once again using Bowtie v2 (Langmead and Salzberg, 2012), which performs well for the alignment of shorter reads. Bowtie parameters were the following: -local, -p 24, -q, -I 13. Read pairs successfully aligned by STAR/Bowtie were merged into a BAM file and separated into two groups: those that align to one genomic position and those that align to more than one position (Fig. 1F).

Read assignment and corrections
The correct_count module of CoCo employs featureCounts (Liao et al., 2013) to assign aligned read pairs to their corresponding genes using the given annotation file as reference (Fig. 1). Parameter values for the correct_count module are described in Supplementary Table  S3. CoCo uses a corrected annotation produced from Ensembl supplemented annotation by removing regions corresponding to nested genes from longer host genes (described below). Read pairs aligning to only one genomic location increase the alignment count by one for the genomic feature encoded at that location. By default, correct_count distributes the reads between their assigned genes according to their number of assigned uniquely mapped reads. If no singly mapped reads exist for any of the assigned genes, the reads are distributed uniformly. Total counts assigned to a gene correspond to the sum of the uniquely mapped read pairs and the proportion of multimapped read pairs assigned to the gene. The correct_count module also corrects for over-attribution of reads to nested genes by subtracting the length normalized background counts attributed to the feature in which the nested gene is encoded ( Supplementary Fig. S1).

Annotation supplementation
An annotation file in gene transfer format (.gtf) was obtained from Ensembl (Yates et al., 2016) (hg38, v87). The annotation file was supplemented with 628 additional tRNA from GtRNAdb (Chan and Lowe, 2016) and with 20 snoRNA from Refseq (O'Leary et al., 2016) that were missing from Ensembl annotations. 63 gene annotations were removed from the gtf file because they overlap another gene of same biotype over >90% of their length reciprocally (keeping both annotations would result in reads aligning to them being labelled as ambiguous). These 63 redundant genes consist of 17 snoRNAs, 43 miRNAs, 2 long intergenic non-coding RNAs (lincRNAs) and 1 antisense RNA. Details are given in Supplementary Material S1.

CoCo's annotation correction
The correct_annotation module of CoCo produces a modified annotation file in which any exon position overlapping a snoRNA, small Caja-body RNA (scaRNA), snRNA, tRNA or miRNA is removed, resulting in a gapped annotation file (Fig. 1B). To do so, correct_annotation builds a list of gene coordinates corresponding to the above biotypes and then seeks overlaps with exons from all transcripts of genes from other biotypes using Bedtools intersect (Quinlan and Hall, 2010). Exons overlapping a nested gene have their overlapping positions removed, effectively splitting the exon in two if the nested gene is fully within the exon, or truncating the exon if the nested gene overlaps its end. From this, two new .gtf annotation files are made, one containing all the genes and transcripts, the other one including only the portion of the genes hosting a nested gene (respectively referred as 'CoCo Full gtf' and 'Introns gtf' in Fig. 1 and Supplementary Fig. S1).

Conversion from counts to transcripts per million
The counts obtained from CoCo's correct_count module were normalized by the length of the main transcript of the gene and by the read count to give transcripts per million (TPM) as described further in (Boivin et al., 2018). The final output of the correct_count module holds raw read counts, counts per million (CPM) and TPM values for all annotated genes. Supplementary Material S2 provides average fragmented TGIRT-seq abundance values in TPM, estimated by CoCo and other read assignment tools considered.

Gene biotype pooling
Ensembl biotypes protein_coding, pseudogene and long_noncoding were pooled as recommended (http://ensembl.org/help/faq? id¼468). The group 'other' corresponds to all other biotypes not listed

Comparison to other read assignment tools
CoCo was compared with the following read assignment tools: featureCounts (Liao et al., 2013), HTseq-count (Anders et al., 2015), RSEM (Li et al., 2010), Cufflinks (Trapnell et al., 2012) and STAR (Dobin and Gingeras, 2016). The parameter values used for each tool are indicated in Supplementary Table S3. Throughout the text, featureCounts with typical parameter values is used as the standard read assignment pipeline.

Polymerase chain reaction
Digital PCR was performed by the Université de Sherbrooke RNomics Platform (http://rnomics.med.usherbrooke.ca/). Droplet Digital PCR (ddPCR) reactions were prepared using 10 ll of 2Â QX200 ddPCR EvaGreen Supermix (Bio-Rad), 10 ng (3 ml) cDNA, 100 nM final (2 ml) primer pair solutions and 5ul molecular grade sterile water (Wisent) for a 20 ll total reaction. Each reaction mix was converted to droplets with the QX200 droplet generator (Bio-Rad). Droplet-partitioned samples were then transferred to a 96-well plate, sealed and cycled in a C1000 deep well Thermocycler (Bio-Rad) under the following cycling protocol: 95 C for 5 min (DNA polymerase activation), followed by 50 cycles of 95 C for 30 s (denaturation), 59 C for 1 min (annealing) and 72 C for 30 s (extension) followed by post-cycling steps of 4 C for 5 min and 90 C for 5 min (Signal stabilization) and an infinite 12 C hold for reading using the QX200 reader (Bio-Rad). The concentration reported is in copies/ll of the final 1Â ddPCR reaction (using QuantaSoft software from Bio-Rad). All primer sequences are available in Supplementary  Table S4. 3 Results

The CoCo correction for nested genes and multimapped reads
Standard gene quantification programmes assign reads according to the amount of overlap between the read and the feature being quantified. As a consequence, reads mapping with the same number of matches to a host gene and to a small non-coding RNA gene nested within its intron/exon are often considered ambiguous (e.g. Fig. 1C and D, read pair A). In addition, reads from nested genes that exceed the annotations [which, in many cases, are not accurate (Deschamps-Francoeur et al., 2014;Kishore et al., 2013)], by even only one nucleotide, are typically automatically assigned to the host gene (e.g. Fig. 1C and D, read pair C). To address these problems, we have developed the CoCo package, which consists of three main modules: (i) the correct_annotation module which generates gapped annotation files in which the regions of the host gene transcript features overlapping with nested genes are precisely removed (Fig. 1B), (ii) the cor-rect_count module which recuperates the reads associated with nested and multimapped genes using the modified annotation ( Fig. 1D and E) and (iii) the correct_bedgraph annotation which produces accurate representations of paired-end reads ( Supplementary Fig. S2).
To test the quantification accuracy of the CoCo pipeline, we examined its capacity to correctly assign and quantify sequencing reads using diverse RNA-seq datasets, and compared its quantification to those of the main read assignment pipelines available. Most comparisons we report involve samples generated using a library preparation protocol that accurately detects both coding and noncoding RNA classes, thanks to the use of the TGIRT reverse transcriptase and ribodepletion (Boivin et al., 2018). Two of the samples were fragmented which enables the quantification of both long and short RNAs (e.g. nested genes and their host gene) while the other two samples were not fragmented, providing a deeper view of short RNAs, which is comparable to size-selection RNA-seq. In addition to TGIRT-seq datasets, we also investigate the performance of CoCo in small RNA-seq (Farazi et al., 2011), hydro-tRNAseq (Gogakos et al., 2017) and standard viral RT ribodepleted sequencing (Choy et al., 2015) datasets. Following standard pre-processing steps, reads were mapped to the hg38 human genome using STAR, resulting in an average alignment rate of 90%. The 10% unaligned reads using STAR are mostly alignments deemed 'too short' by STAR and thus a second alignment step using Bowtie was employed, resulting in an overall average alignment rate of 99% (Fig. 1F). Aligned reads were then assigned to annotated genes using the correct_count module of the CoCo pipeline and the annotation files provided by CoCo's correc-t_annotation module. The correct_count module not only reattributes reads from host genes to nested genes but also corrects the reassignment considering the background read counts of the feature in which the gene is nested ( Supplementary Fig. S1), as discussed below.

Impact of CoCo on the quantification of nested genes and comparison to other read assignment tools
To evaluate the impact of CoCo on RNA detection and quantification, we compared the aligned read assignments obtained using the CoCo pipeline to values obtained using five currently available and commonly used read assignment pipelines, described in Section 2 and in Supplementary Table S3. The pre-processing and alignment steps (described in Section 2 and in Fig. 1F) were the same for all read assignment tools, except for RSEM which requires an alignment to the transcriptome instead of the genome (so STAR but not bowtie was used, with the additional parameter -quantMode TranscriptomeSAM). As indicated in Figure 2 and Supplementary  Figures S3-S11, the CoCo pipeline performs well in the quantification of nested genes, in TGIRT-seq fragmented and non-fragmented datasets, as well as in the small RNA-seq, hydro-tRNAseq and standard viral RT ribodepleted sequencing datasets, always obtaining raw read counts in close agreement with those estimated from the bedgraph, for all nested genes considered. In contrast, HTSeq-count, featureCounts with standard parameter values and STAR fail to detect many nested genes, in all datasets considered. Cufflinks was not considered for this comparison because it does not provide raw read counts that can be compared with bedgraph counts. Examples of poorly quantified nested genes include snoRNA (Fig. 2, Supplementary Figs S3, S5, S9 and S11), snRNA ( Supplementary Fig.  S4), tRNA (Supplementary Figs S6, S8 and S10) and a sca RNA ( Supplementary Fig. S7), which all overlap, partially or completely, retained introns or exons of their host gene. Manual inspection of the alignment file and visual inspection of the bedgraphs confirm that all these nested genes have corresponding aligned reads and should thus be detected and quantified as expressed. RSEM and featureCounts with optimized parameter values perform better than HTSeq-count, STAR and featureCounts with standard parameter values, but still deviate more from bedgraph estimates than does CoCo (see e.g. Supplementary Figs S3, S4 and S6-S11). In the presence of ambiguous reads (Fig. 1, read pair A), featureCounts_optimized splits the read counts equally between the nested gene and the host gene. For example, in the case of the PTCH2 gene locus ( Supplementary  Fig. S4), a visual inspection of the bedgraph indicates that the great majority if not all read counts should be attributed to the nested gene RNU5E-6P while featureCounts_optimized splits the read counts equally between the two genes for the fragmented datasets. Interestingly, in this case, CoCo not only corrects the abundance values of the nested gene but also adjusts the counts given to the host gene, reassigning 89% of the counts (313 read pairs out of 350) from the host gene to the snRNA. Based on these results, we conclude that CoCo is capable of rescuing and properly assigning reads from the nested genes lost to several read assignment tools and generates quantification values that are most consistent with the bedgraph profiles.

Background correction for nested genes
In order to accurately quantify read counts from nested genes and distinguish their reads from the host gene reads, CoCo carries out a background correction. This is achieved by subtracting the average read count of the feature in which the nested gene is encoded from the read count assigned to the nested gene. This correction does not significantly change the final read counts of genes nested in host gene features (most often introns) that are expressed at low levels ( Supplementary Fig. S12). However, this correction has a major impact when the host gene feature encoding the nested gene is highly expressed. For example, in the case of the CH507-513H4.1 locus which hosts miRNAs miR-3648 and miR-3687, the reads were originally attributed to the miRNA despite the absence of corresponding peaks in the bedgraph. This inappropriate assignment is no longer observed following background correction ( Supplementary  Fig. S13).

Quantification of transcripts from duplicated genes
A large number of genes, in particular those producing non-coding RNA, exist in more than one copy. The RNA produced from most of these different forms of repeated features cannot be quantified by standard read assignment modules due to their identical or near identical sequences that result in multimapping of sequencing reads (e.g. Fig. 1C, read pair K). The CoCo pipeline recuperates these otherwise lost reads by distributing the counts between all genes assigned according to the distribution of their uniquely mapped read pairs when possible. For small non-coding RNAs such as snoRNAs, uniquely mapped reads typically originate from flanking genomic sequences that are seldom included in the transcripts. Longer genes are more likely to have longer proportions of unique sequences. For example, as indicated in Supplementary Figure S14, the annotated mature forms of SNORD103A and SNORD103B are 100% identical. As a result, 96% (2314/2401) of all read pairs aligning to these snoRNAs map equally well to both SNORD103A and SNORD103B and would thus be discarded without considering multimapped reads. Considering the distribution of the 87 uniquely aligned read pairs (to non-identical genomic flanking sequences for SNORD103A and SNORD103B), CoCo distributed the 2314 multimapped read pairs proportionally to the uniquely mapped read pairs ( Supplementary  Fig. S14B). As a consequence, the sum of the read counts to the SNORD103 family increased by 25Â passing from 87 (considering only uniquely mapped reads) to 2401 read pairs (using the CoCo's correct_count module). A similar example from the small RNA-seq datasets, the paralogues MIR101-1 and MIR101-2, is shown in Supplementary Figure S15. Only MIR101-1 has uniquely aligned read pairs, resulting in all the multimapped read pairs assigned to MIR101-1 and none to MIR101-2. In one of the datasets, this results in 5495 read pairs instead of 160 uniquely mapped reads, representing a 34Â increase. Therefore, by using the uniquely mapped reads as a guide, CoCo rescues and re-distributes the multimapped reads to provide a realistic read distribution. Overall, the CoCo multimapped correction module increased the estimated abundance of 1443 multimapped genes, most of which are tRNAs and snoRNAs, by more than 2-fold in the TGIRT-seq fragmented datasets.

Experimental validation of CoCo-based quantification for nested and repeated genes
To experimentally validate the accuracy of the CoCo-based quantification, we chose nine overlapping or repeated genes and examined their abundance using ddPCR, comparing these values to sequencing abundance estimates. These nine genes include seven snoRNAs overlapping with protein-coding genes and two repeated gene families (RN7SK and RN7SL2). In the case of the multimapped genes RN7SK and RN7SL2, their abundance was estimated to be respectively two and four times higher, (Fig. 3) using the CoCo pipeline compared with featureCounts_standard estimates. When the ddPCR was compared with all read assignment tools considered (Fig. 3 and Supplementary Fig. S16), as for the bedgraph comparisons, HTSeq, STAR and featureCounts using standard parameter values did not detect the nested genes while CoCo, RSEM and featureCounts with optimized parameter values agreed with ddPCR values, obtaining Pearson correlation values above 0.99.

CoCo increases the number of sequencing reads considered and enhances the detection of non-coding RNA
The capacity to detect transcripts from repeated and overlapping genes not only provides new information about the expression dynamics of these genes but may also alter the overall distribution of read counts. Accordingly, we examined the impact of using CoCo on the overall distribution of transcript estimates in the human genome. Comparison of the read assignments obtained using the CoCo pipeline or the standard featureCounts pipeline indicates that while 65.8% of the read pairs in fragmented datasets were assigned identically by both strategies and 18.6% were assigned by neither, >15% were only assigned by CoCo (Supplementary Fig. S17). Most of the read pairs correctly assigned by the standard pipeline originate from protein-coding genes (74%) and tRNAs (12%) (Supplementary Fig.  S17E, left panel). Read pairs that were not assigned by either pipeline originate mostly (88%) from currently unannotated genomics regions ( Supplementary Fig. S18), as is the case for example for an intronic region in the AKAP6 gene in which >32 000 read pairs align but no feature is annotated (Supplementary Fig. S19).
Most of the read pairs only assigned by CoCo (15.7%) originate from multimapped genes, representing 12.9% of all aligned read pairs ( Supplementary Fig. S17A). As expected, the rescued multimapped read pairs originate mainly from non-coding RNA including tRNAs (37%), 7SL (27%), snRNAs (14%) (Supplementary Fig. S17C). However, 13% of the reads aligned to protein-coding genes suggesting that there is a substantial number of protein-coding genes that contain repeated sequence (Supplementary Fig. S17C). A small proportion of the read pairs only assigned by CoCo (2.4%) were aligned to overlapping genes that are labelled as ambiguous by standard pipelines ( Supplementary Fig. S17A, right panel). Most of the overlapping gene counts originated from non-coding RNA including snoRNAs (40%), snRNAs (37%) and 7SK (18%) (Supplementary Fig. S17D). In addition to the 15.3% of rescued ambiguous and multimapped read pairs, CoCo also reassigned a small proportion of reads (0.37%) that are misassigned by the standard pipeline ( Fig. 1C and D, read pair C and Supplementary Fig. S20). In all cases, misassignments by the standard pipeline result in erroneous association of reads to the host gene instead of the nested gene. SnoRNAs represent 94% of such reassignments ( Supplementary Fig. S20). Together, these observations indicate that CoCo greatly increases the percentage of usable read counts and reduces the number of misassignments.

CoCo provides a more accurate depiction of the transcriptome landscape
As a consequence of the CoCo correction, the overall distribution of the human transcriptome was modified considerably. The proportion of all read counts attributed to protein-coding genes in fragmented datasets was reduced by 12%, representing 74% of all assigned reads using the standard pipeline but only 62% using CoCo (compare Supplementary Fig. S17E, left and right panels). In contrast, noncoding RNA including tRNAs and snRNA gain 4 and 2% of the total read counts, respectively. The most dramatic change in total read count distribution was observed with the highly redundant (multimapped) gene coding for the signal recognition particle RNA 7SL, which was increased from 2% using the standard pipeline to 6% using CoCo. Analyses of transcript abundance distribution, which take into consideration the transcript length, indicate that the biggest impact of CoCo is in adjusting the proportion of coding to noncoding RNA and increasing representation of RNA families with a strong prevalence of repeated genes. As indicated in Supplementary Figure S17F, CoCo estimates of transcript abundance reduced the proportion of protein-coding genes by 5% while tripling the proportion of the 7SL and increasing the proportion of snoRNA and snRNA by 2% each, when compared with the standard pipeline.
The impact of CoCo on transcript quantification is most visible in cases where RNA abundance changes from completely undetected using standard pipelines to abundantly detected with the CoCo pipeline (e.g. Figs 2, 3, Supplementary Figs S3, S4 and S6-S11). The de novo detection of these genes gives a new view of an otherwise uncharted portion of the human transcriptome. Accordingly, we counted the number and biotype distribution of genes only detected using CoCo as compared with a standard read assignment pipeline to better understand their origin and contribution to the human transcriptome, using a fragmented dataset to consider the full RNA landscape. As shown in Figure 4A, the most affected RNA family is the snoRNAs with 22% of genes (137 out of 628 expressed snoRNAs) detected only with the CoCo pipeline. The distribution of the corrected abundance of these 137 'invisible' C/D and H/ACA box snoRNA genes is displayed in Supplementary Figure S21. More than 50% of these snoRNAs (80/ 137) have a corrected abundance above 100 TPM, and are thus amongst the most abundant RNAs in the cell (indeed, only 1.8% of all expressed genes have cellular transcripts with abundance >100 TPM). Ten such snoRNAs are even detected with an abundance as high as >1000 TPM (e.g. SNORD26, SNORD14C). These data clearly indicate that failure to detect snoRNAs using standard assignment methods is not restricted to rare or lowly expressed RNA. On the contrary, it is the highly expressed genes that are often missed, most likely due to the fact that highly expressed non-coding RNA tend to originate from more than one gene or to be nested in highly transcribed genes. In addition to snoRNAs, a smaller but still significant proportion of each of the other main classes of non-coding RNA is not detected using standard techniques including >10% of scaRNA and as many as 27 snRNA genes, 27 tRNA genes and 224 lncRNAs genes (Fig. 4A).
A less obvious but equally important difference between the CoCo and standard pipelines lies in the accuracy of read assignments, especially to duplicated genes. For example, while both CoCo and the standard pipeline agree on tRNA being the most abundant transcript biotype ( Supplementary Fig. S17F), this class of RNA features the largest number of genes with count corrections by CoCo. Indeed, 30% of tRNA genes more than doubled their abundance after using CoCo (Fig. 4B). Similarly, a great number of protein-coding genes (>50%) also increased in abundance, but unlike tRNA, the difference in read pair count was mostly modest with an average of 5% change in read counts (Fig. 4B). We therefore conclude that the correct counting and assignment of sequencing reads plays an important role not only in increasing the number of detected genes but in establishing the relative abundance of RNA within the transcriptome.

Runtime and memory requirements
CoCo can be run using multiple threads thus improving its running performance. The runtime and memory required to use the different read assignment tools considered in this study is shown in Table 1.

Discussion
Increasing sequencing depth and coverage is considered the obvious target for improving the quality of transcriptome analysis. However, in this study, we show that the precision, quality and depth of transcriptome analysis can be greatly enhanced by carefully choosing and tuning read assignment tools using existing sequencing data. Modification of the read assignment pipeline for the analysis of a TGIRT-seq fragmented dataset, which provides an accurate view of the whole transcriptome (Boivin et al., 2018) increased the proportion of assigned reads by 15% and modified the number of assigned reads for 50% of all expressed genes ( Fig. 4 and Supplementary Fig.  S17). Indeed, more than 750 additional transcripts were detected simply by correcting read association for nested genes (Fig. 4A), while multimapping affects the quantification of 15 121 genes. The problem with overlapping genes was solved by introducing gaps corresponding to nested gene positions in the reference annotation files (Figs 1, 2 and Supplementary Figs S3-S11). The gaps ensure that read pairs are not automatically assigned to the host gene even if they slightly exceed the often inaccurate annotation of nested genes (Deschamps-Francoeur et al., 2014;Kishore et al., 2013). CoCo is not the only tool that can correctly address the quantification of nested genes, although CoCo provides abundance values that are most consistent with bedgraph estimates regardless of the library preparation methodology used (examples shown in Fig. 2 and Supplementary Figs S3-S11) and performs well when compared with ddPCR (Fig. 3). Multimapped reads were dealt with by distributing them proportionally to uniquely mapped reads, as first introduced by MuMRescue (Faulkner et al., 2008) and ERANGE (Mortazavi et al., 2008). The corrections proposed in our study can be applied as a supplement to any read-assignment tool as a feature. By applying the correction tool CoCo, most sequence analysis pipelines will benefit from increased transcriptome coverage and more accurate transcript quantification.
Five different tools (one of which was run using two different settings) were chosen for comparison to the CoCo pipeline because they are widely used and easily implemented standalone pipelines. These tools were compared considering typical examples of nested and multimapped genes in datasets covering five library preparation methodologies (TGIRT-seq both fragmented and non-fragmented, small RNA-seq, hydro-tRNAseq and standard viral RT ribodepleted RNAseq) to investigate their wide applicability. HTSeq-count, STAR, Cufflinks and featureCounts using standard parameter values generally performed poorly to quantify nested and multimapped genes, both according to comparisons to bedgraphs and ddPCR quantification. These four tools failed to detect many nested genes and did not accurately quantify many multimapped genes. In contrast, RSEM, featureCounts with optimized parameter values and most often CoCo performed much better at quantifying nested and multimapped genes, both according to bedgraph and ddPCR comparisons, displaying good capacity for both small and long RNA-seq, fragmented or not, single and paired-end sequencing, regardless of the strandedness and type of libraries (Fig. 2, Supplementary Figs. S3-S11, S14 and S15). RSEM's greatest drawback is its long runtime and high memory usage while featureCounts with optimized parameter values does not assign ambiguous reads in an ideal manner as it simply splits the assignment counts equally for all features overlapping these reads. In addition, RSEM and featureCounts with optimized parameter values do not perform well when quantifying non-fragmented datasets, generally attributing too many reads to host genes. They are thus not recommended for size-selection datasets often used to quantify miRNAs. In Cufflinks and HTSeq require sorting before the read assignment can occur. The values given include sorting before the tool is run. general, if any of these tools are chosen to quantify RNA-seq datasets, we highly recommend the parameter values used in this study (see Supplementary Table S4) as using different parameter values can significantly change the quantification and thus the accuracy (compare for example the results obtained using featureCounts_standard and featureCounts_optimized).
Accurately quantifying multimapped reads has been investigated in several studies. featureCounts annotates such reads, making it easy for CoCo's correct_count module to distribute them according to uniquely mapped reads. featureCounts itself provides an option to deal with multimapped reads, by splitting them equally between all members of a multimapped group (using the -M, -fraction options). However, a significant proportion of these members are likely not expressed, particularly in the large families of non-coding RNAs with tens or even hundreds of copies, making CoCo's strategy (the proportional distribution of multimapped reads) more realistic. Evaluating the accuracy of multimapped read assignment is difficult since unlike mRNA, most non-coding RNA lack external unique sequences that could be used to differentiate between genes with shared sequences. In building our pipeline for the quantification of nested and multimapped genes, we opted for a computationally quick and simple solution that enables the evaluation of the overall abundance of RNA generated from genes with shared sequence. While we cannot guarantee the accuracy of the abundance for each individual gene from a repeated family, the overall abundance of the RNA generated from these repeated families is accurately determined by CoCo as validated experimentally by ddPCR analysis and gene-specific analyses (Fig. 3 and Boivin et al., 2018).
One of the most surprising observations in the application of CoCo is how many genomic regions do not have proper annotation and how this mis-or lacking annotation may affect the overall interpretation of transcript distribution. Indeed, while CoCo was able to recuperate the ambiguous reads due to gene overlap with small noncoding RNAs and multimapping, about 18.6% of all read pairs remained unassigned. Examination of these read pairs shows that the great majority align to unannotated regions in the genome ( Supplementary Fig. S18), indicating that improving the annotations is now becoming essential to increase the coverage and the accuracy of RNA-seq quantification.
The read assignment corrections shown here are essential steps to study the human transcriptome. Most of the multimapped and overlapping reads originate from highly expressed genes and thus any changes in their read assignment will significantly affect the overall distribution of the human transcriptome. In addition, accurate quantification of nested and multimapped genes not only enhances the detection of these types of RNA but also corrects the quantification of their protein-coding host genes. Thus, any sequence analysis pipeline that ignores or fails to accurately detect nested and multimapped genes would likely result in significant changes in read count distribution and ultimately in incorrect expression estimates, for large proportions of the transcriptome. Indeed, as sequencing depth increases and the capacity to simultaneously detect both coding and non-coding RNA improves, read assignment tools like CoCo will become essential for any sequencing analysis pipeline.


extfieldsofstudy: Biology Medicine

text:
Genetic and cytological analyses reveal the recombination landscape of a partially differentiated plant sex chromosome in kiwifruit

Background Angiosperm sex chromosomes, where present, are generally recently evolved. The key step in initiating the development of sex chromosomes from autosomes is the establishment of a sex-determining locus within a region of non-recombination. To better understand early sex chromosome evolution, it is important to determine the process by which recombination is suppressed around the sex determining genes. We have used the dioecious angiosperm kiwifruit Actinidia chinensis var. chinensis, which has an active-Y sex chromosome system, to study recombination rates around the sex locus, to better understand key events in the development of sex chromosomes. Results We have confirmed the sex-determining region (SDR) in A. chinensis var. chinensis, using a combination of high density genetic mapping and fluorescent in situ hybridisation (FISH) of Bacterial Artificial Chromosomes (BACs) linked to the sex markers onto pachytene chromosomes. The SDR is a subtelomeric non-recombining region adjacent to the nucleolar organiser region (NOR). A region of restricted recombination of around 6 Mbp in size in both male and female maps spans the SDR and covers around a third of chromosome 25. Conclusions As recombination is suppressed over a similar region between X chromosomes and between and X and Y chromosomes, we propose that recombination is suppressed in this region because of the proximity of the NOR and the centromere, with both the NOR and centromere suppressing recombination, and this predates suppressed recombination due to differences between X and Y chromosomes. Such regions of suppressed recombination in the genome provide an opportunity for the evolution of sex chromosomes, if a sex-determining locus develops there or translocates into this region. Electronic supplementary material The online version of this article (10.1186/s12870-019-1766-2) contains supplementary material, which is available to authorized users.

Background
Sex chromosomes have evolved independently in a number of plants [1,2]. They have evolved from autosomes through the development of a sex-determining locus in a region of non-recombination or the development of non-recombination around the sex determining region [3], which blocks the exchange of genetic material between the sex chromosomes and enables them to begin to diverge from each other [4]. The combination of different selection pressures in each sex and the lack of recombination results in a range of evolutionary processes including Muller's ratchet, background selection and genetic hitchhiking [5,6].
Angiosperm sex chromosomes are much younger than mammalian sex chromosomes, which are estimated to have evolved approximately 167 million years ago (Mya) [7]. Mammalian Y chromosomes are generally severely degraded compared with their X chromosome counterparts, and contain few genes. This is thought to be a consequence of the lack of recombination between X and Y chromosomes in the male-specific region. The ancestor of the angiosperms was a hermaphrodite 140-200 Mya, and dioecy and sex chromosomes have evolved independently in different angiosperm families [3]. The sex chromosomes in Silene latifolia are estimated to have evolved 5-10 Mya [8] and those in papaya are also thought to have recently evolved through two chromosomal inversions occurring approximately 7.0 and 1.9 Mya [9]. The study of plant sex chromosomes can provide important insights into the evolution of sex chromosomes in general, especially the early stages [10,11].
A key event in the formation of sex chromosomes in organisms that contain more than one sex determining gene is the suppression of recombination around these genes, which allows the multiple genes to be inherited as a single unit so that only male and female progeny are produced. Questions remain as to how this occurs [12,13]. Studies of sex chromosomes in the early stages of evolution, such as those found in dioecious angiosperms, are needed to answer these questions.
The genus Actinidia (kiwifruit) contains more than 50 species of long-lived perennials [14]. All known Actinidia species are functionally dioecious, with male and female flowers on different plants [15]. Actinidia have a basic chromosome number of x = 29, many of the species are polyploid, ranging from diploid to octoploid including A. chinensis var. chinensis, which has diploid and tetraploid races [16]. Female plants bear flowers that appear complete; however, they produce only empty pollen grains, while flowers of male plants produce fertile pollen, but pistil growth is suppressed before style elongation or ovule initiation. Sex determination in Actinidia appears to be monofactorial at all ploidy levels, as male and female plants normally occur in a 1:1 ratio [17]. However, there are populations in which the inheritance of a Y chromosome is not 1:1 in the progeny, for example in populations where the male parent produces unreduced gametes [18]. Likewise, following the colchicine doubling of diploid males of A. chinensis var. chinensis to produce a tetraploid male (XXYY), male and female progeny are produced in a 5:1 ratio (Wu, unpublished). Rare 'inconstant' diploid male A. chinensis var. deliciosa plants that occasionally produce fruit have been identified; these fruiting males produce male and female progeny in a 3:1 ratio when selfed [19]. This indicates that in Actinidia there is a strict genetic control of sex expression by an active Y system, with diploids having homogametic females (XX) and heterogametic males (XY). Although there are no obvious dimorphic sex chromosomes in Actinidia, it has been suggested that in A. chinensis var. chinensis, the nucleolar organiser region (NOR)-containing chromosomes are the sex chromosomes [20]. This suggestion was based on observations of chromosome pairing during pachytene, where the satellite regions of the NOR-containing chromosomes paired in female plants, but did not pair in male plants. Genetic mapping of diploid A. chinensis var. chinensis also suggests that the sex-determining region (SDR) is located in the subtelomeric region of a single chromosome pair, linkage group 17 [21], which corresponds to the subtelomeric region of chromosome 25 [22]. It is not known how many genes are involved in sex determination in Actinidia.
As all extant Actinidia species are dioecious, regardless of the ploidy or species, it is very likely that their common ancestor was dioecious and contained sex chromosomes. The oldest fossil of the Actinidiaceae is suggested to be from 89 Mya [23] and the fossil record of Actinidia in China extends into the early Miocene or 20-26 Mya [24]; the dating of a backbone phylogeny of Actinidia suggested that the initial separation of the clades occurred c. 26.9 Mya [25]. The two other extant genera in the Actinidiaceae are Saurauia and Clematoclethra. Saurauia can be monoecious or functionally dioecious and Clematoclethra is monoecious [26], theoretically extending the upper limit of the age of the ancestral sex chromosomes in Actinidia into the Cretaceous, although it is likely that they have evolved much more recently than this. Compared with the sex chromosomes in mammals, sex chromosomes in Actinidia are in the early stages of their evolution.
To improve the identification and characterisation of the SDR in diploid A. chinensis var. chinensis, we have used high-density genetic mapping, microsatellite (SSR) mapping and genomic sequence data to investigate the recombination rates around the sex-determining locus. We have also identified hybridisation sites on pachytene chromosomes of Bacterial Artificial Chromosome (BAC) probes that we have found to be linked with the SDR.

Recombination rates along pseudochromosome 25
Recombination rates along Actinidia sex chromosomes (chromosome 25) as determined by recombinant mapping are shown in Fig. 1. The high density genetic map based on an interspecific F1 cross between Actinidia rufa 'MT570001' and A. chinensis var. chinensis 'Guihai No4' (Mapping population I) contains 210 markers along chromosome 25 [22], and exhibits suppressed recombination along the first 5-6 Mb of chromosome 25 in the A. chinensis var. chinensis (paternal) map, indicating minimal or no recombination between X and Y chromosomes in this region. The A. rufa (maternal) map also shows suppressed recombination in the terminal 5-6 Mb of chromosome 25, indicating minimal or no recombination between X chromosomes in this region.

Microsatellite mapping
Based on the inheritance of microsatellite marker alleles within the terminal 6 Mb of chromosome 25, the mapping population III progeny (an intraspecific diploid A. chinensis var. chinensis F1 mapping family from a cross between female parent X 1 X 2 and male parent X 3 Y 1 ) could be categorised into one of four groups: female progeny either X 1 X 3 or X 2 X 3 and male progeny either X 1 Y 1 or X 2 Y 1 . The markers were ordered by their position in the 'Hongyang' genome, and the allelic content of microsatellites for these groups of progeny are shown in Fig. 2. For 18 of 21 markers, almost all the microsatellite alleles obtained were as expected for inheritance of this region from their parents with no recombination. However, as can be seen in Fig. 2, there were some exceptions, for example MP121 had a b allele for marker Ke1251 from chromosome X 2 rather than an a allele from chromosome X 1 as expected from the surrounding markers. Similarly in MP118 for marker Ke1009 there was a b allele from chromosome X 3 rather than a c allele from the Y chromosome as expected from the surrounding markers.

Distribution of annotations and genetic variation on Pseudochromosome 25
Inspection of the distribution of gene models and repeat annotations showed that the gene density is lower and repeat density higher in the terminal 0-6 Mb region than in the remainder of chromosome 25. A very gene-poor, repeat-rich region was identified at 4.5-5 Mb (Fig. 3). Genomic sequence data was available for 14 A. chinensis genotypes and these were used to calculate nucleotide diversity. Nucleotide diversity was found to be lower in the 0-5 Mb region than in the remainder of the chromosome. Windowed fixation index (Fst) measures of differentiation between males and females showed peaks in the 3-8 Mb region.
To further assess meiotic recombination in A. chinensis var. chinensis sex chromosomes, we calculated kinship pairwise in 1-Mb windows between three full-sib males using whole genome sequence data. Each of the three full-sib males was expected to have received the same haplotype of the non-recombinant region of the Y chromosome from their paternal parent, and one of two possible X haplotypes from their maternal parent. Kinship between 0.46 and 0.48 was observed for the first 5 Mb of the pseudochromosome in all three full-sib-to-full-sib pairings, indicating highly similar to identical genotypes of this region in all three full-sib males (Fig. 4). Inheritance of the same X and Y haplotypes of this region reveals that they inherited the same X haplotype from their maternal parent (by chance, as  [20], the chromosome of origin of each allele is shown (X 1 (yellow), X 2 (dark blue), X 3 (light blue) or Y 1 (red) or undetermined (.) and these are grouped by chromosome pairings. This was a consequence of using informative microsatellite markers, e.g. for fully informative marker Ke511 the parental alleles were a(217)b(205) and c(215)d(219) giving a(217)c(215) and b(205)c(215) female progenies and a(217)d(219) and b(205)d(215) male progenies, where the four alleles were detected as length differences in base pairs of amplified fragments during capillary electrophoresis b The marker positions on chromosome 25 of the 'Hongyang' genome are shown, with non-segregating markers in blue and segregating markers in red each could theoretically have received either of the two haplotypes) and the same Y haplotype from their paternal parent through the entirety of the 5 Mb of the pseudochromosome. There was no recombination between the two X chromosomes or the X and Y chromosomes in this region during the meiosis that led to the production of these three full-sib males. For the rest of the chromosome, the windowed kinship coefficients ranged from zero or less (no haplotypes shared) to 0.47 (all haplotypes shared), consistent with recombination occurring throughout the rest of the chromosome.

Fluorescent in situ hybridisation
The observed hybridisation patterns were similar for all three BAC probes derived from the SDR of the female parent. In the female A. chinensis var. chinensis, CK51_05, a strong hybridisation signal was observed at the terminal end of one chromosome pair, as well as some signals at centromeres (Fig. 5). The major hybridisation site was associated with a secondary constriction, identified as the NOR. The NOR accounts for approximately 5% of the total chromosome length in A. chinensis var. chinensis [20].

Discussion
In our study we have employed several methods to characterise the recombination landscape of the SDR in A. chinensis var. chinensis, previously mapped genetically [21,22]. In vitro fluorescent hybridisation using probes specific for the SDR indicated that a single chromosome pair was involved and furthermore, an analysis of recombination rates along chromosomes indicated the location and extent of suppressed recombination on X and Y variants of chromosome 25.
Recombination rates between the X and Y chromosomes were found to be very low in the terminal region of chromosome 25 where the SDR is located. Suppression of recombination was observed between X and Y chromosomes over the first 5-6 Mb in both male parents of the three bi-parental genetic mapping families studied. In addition, kinship analysis of three full-sib males showed that they all inherited the same haplotype from their paternal parent for the first 5 Mb of the Y chromosome. This is consistent with an active-Y type system, as sex determination is thought to be controlled by at least two tightly linked sex-determining genes on the Y chromosome, and any recombination between the sex-determining genes in male (XY) plants would produce recombinant asexual or hermaphroditic offspring, rather than the unisexual male and female offspring observed in Actinidia.  The suppressed recombination between X and Y chromosomes is regarded as the first step in the formation of sex chromosomes [27]. Recombination rates between the two X chromosomes in A. chinensis var. chinensis 'Hort16A' and those in A. rufa 'MT570001' were suppressed over the first 5-6 Mb. This is consistent with the hypothesis, as was the inheritance pattern for microsatellite markers, where the alleles originating from the same X haplotype were inherited together over the terminal 6 Mb. The kinship analysis of the three full-sib males showed that each inherited the same haplotype from their maternal parent for the first 5 Mb of the X chromosome, further evidence of suppression of recombination in this region. The similarities in the extent of the recombination suppression between the two X chromosomes, as well as between X and Y chromosomes, suggests that the recombination suppression in this region is not due to differences such as inversions between the X and Y chromosomes alone.
The similar pattern of recombination in both X and Y chromosomes of A. chinensis var. chinensis, of which 30-35% is a terminal non-recombining region where the sex-determining genes are located, while the remaining 65-70% undergoes normal recombination, suggests a mechanism for the evolution of sex chromosomes in A. chinensis var. chinensis from autosomes. Recombination rates are uneven across a genome, with both hotspots of recombination and regions of suppressed recombination [28]. This suppression of recombination around sex-determining genes is a key feature of sex chromosomes, allowing the two chromosomes of a pair to evolve separately. A key question in addressing how sex chromosomes are initially established is "how is recombination initially suppressed?" [29][30][31]. In A. chinensis var. chinensis, it appears that the SDR has evolved in a region of suppressed recombination because of the proximity of the NOR and centromere. The close proximity of the NOR and centromere to the sex-determining genes is likely to have created a region of suppressed recombination, allowing recombination to be restricted between the sex-determining genes, and has allowed sex chromosomes to evolve. In A. chinensis var. chinensis there are thought to be at least two sex-determining genes [32] and these would include the gene responsible for suppression of ovary development (SuF, according to the terminology of Westergaard M [33]), as well as a gene responsible for pollen maturation (M) [3]. The Y-encoded gene responsible for the suppression of ovary formation in Actinidia males has recently been identified as a male-specific type-C cytokinin response regulator [34]. The evolution of sex chromosomes requires not only evolution of alleles for these genes, but also for them to be located in a region of suppressed recombination, to enable the production of unisexual male and female offspring only [35].
Reduced rates of recombination around the centromere have been shown in several animals as well as in plants and fungi [36]. The inhibition of recombination around the centromere is important for the even segregation of chromosomes during meiosis. It is achieved through the Ctf19/CCAN kinetochore complex. It has been suggested that the complex suppresses recombination in two ways: inhibiting the formation of double strand breaks in the region close to the centromere, and recruiting cohesion to the pericentromeric region, which promotes inter-sister chromatid repair rather than inter-homologue repair [37]. There is also evidence for reduced rates of recombination around the ribosomal DNA (rDNA) in the NOR. Reduced rates of meiotic recombination along the NOR-containing chromosomes are observed in Arabidopsis [38], while in yeast, meiotic recombination between non-sister chromatids is strongly suppressed around the rDNA genes [39], although mitotic and meiotic recombination between sister chromatids occur frequently [39]. Low recombination rates in the SDR of A. chinensis var. chinensis are likely to be due to the proximity of both the NOR and centromere to the SDR. NORs are also present on the sex chromosomes in other organisms, for example Drosophila [40], plants such as liverwort [41] and spinach [42] and these NORs could have also played a role in suppressing recombination around sex-determining genes, allowing the establishment of sex chromosomes in these organisms.
The hybridisation patterns of sex-linked BAC probes indicated that the major hybridisation site was at the terminal region of the chromosome pair that contains a secondary constriction. The hybridisation signal spanned the terminal region, encompassing both the secondary constriction and the centromere (Fig. 5), suggesting that the SDR spans this region and that the sex-determining genes lie between these two features. As this was the only secondary constriction observed in A. chinensis var. chinensis chromosomes, it is likely to be the NOR. The presence of DNA repeats in the BACs which occur both on the sex chromosomes and some of the autosomes is the likely explanation for the hybridisation signal observed on other chromosomes. DNA repeats isolated from sex chromosomes of other plants have been demonstrated to be present on autosomes [43]. The FISH hybridisation results support He and colleagues' [20] assertion that the NOR-containing chromosomes are the sex chromosomes (chromosome 25), and the telomeric position is consistent with genetic mapping results which show that the sex-determining locus is near the end of the linkage group corresponding to chromosome 25 [21,22].
Although recombination was suppressed in the SDR, there was still a low level of genetic exchange between X chromosomes and between X and Y chromosomes. This was observed in the analysis of SSR inheritance in the SDR when some markers showed the presence of alleles that were not expected based on the surrounding markers. Two main possibilities would give rise to this pattern of inheritance: a double crossover event or gene conversion. Since a double crossover event would be extremely unlikely in this region because of suppressed recombination, this pattern is consistent with non-crossover gene conversion. A non-crossover gene conversion is the non-reciprocal transfer of short DNA segments between homologs. In the SDR region of A. chinensis var. chinensis there is evidence for X-X, X-Y and Y-X type gene conversions. For example, there is evidence of an X-X gene conversion at Ke1180 in MP131, evidence of an X-Y gene conversion was seen at Ke1009 in MP118, and evidence of a Y-X gene conversion was seen at Ke1013 in MP088. Overall, when considering the segregation of the markers within the terminal 5 Mb of chromosome 25, i.e. all markers apart from Ke511, Ke216, Ke302 and Ke470, there was evidence for four X1-X2 gene conversions (four unexpected alleles out of 559), three X2-X1 gene conversions (three unexpected alleles out of 495), two X3-Y gene conversions (two unexpected alleles out of 695) and three Y-X3 gene conversions (three unexpected alleles out of 638). Non-crossover gene conversions affect a short region of the DNA. The tracts have been estimated to span 50-1000 bp [44], and gene conversions are one of the more frequent genetic exchanges between homologs in plant genomes [45]. Inter-chromosomal gene conversions have been identified in the sex chromosomes of other organisms, including humans, primates, birds and mosses [46][47][48][49]. A combination of a lack of recombination and the occurrence of gene conversion would account for our results from recombination mapping, kinship analysis and SSR inheritance.

Conclusions
The theory of sex chromosome evolution suggests that sex chromosomes were once an autosomal pair that evolved different morphology and gene content because they lost their ability to recombine [50]. The suppression of recombination around the SDR is a key step in this process. We have confirmed that in the sex chromosomes of A. chinensis var. chinensis, recombination is suppressed around the SDR between both X and Y and X and X chromosomes, suggesting that suppression of recombination in this region is not specific to differences between X and Y chromosomes. We propose that the proximity of the NOR and centromere to this region has created a region of suppressed recombination and that this, along with the evolution of male-sterile and female-sterile mutations, has enabled the sex chromosomes to evolve in Actinidia.

Plant material and DNA isolation
Three bi-parental mapping populations of Actinidia were used to study the recombination rates along chromosome 25. Mapping population I was an interspecific bi-parental mapping family A. rufa × A. chinensis var. chinensis (A, 'MT570001' × 'Guihai No4') [22]. Mapping population II was an intraspecific diploid A. chinensis var. chinensis family developed from a cross between 'Hort16A' and male parent P1. Mapping population III was also an intraspecific diploid A. chinensis var. chinensis family derived from a female from seed from Henan province, and the male parent from a seed accession from Guangxi province, China [21].
Seedlings from mapping population II were raised in tissue culture and 236 individuals selected for genotyping. Young expanded leaf tissue, weighing approximately 100 mg, was harvested and stored at − 80°C by snap freezing in liquid nitrogen. Genomic DNA extraction was performed using the CTAB method [51]. DNA quantification was carried out using Qubit™ fluorometric analysis.

Genotyping-by-sequencing (GBS), variant calling and selection of single nucleotide polymorphism (SNP) markers
The method for developing GBS libraries of population II followed [52], modified by the use of BamHI for the restriction digestion step. The libraries were individually amplified and successful preparation verified by analysis of an aliquot by agarose gel electrophoresis, before pooling the amplicons prior to sequencing [53]. Libraries prepared from 236 genotypes were sequenced over 5 lanes using Illumina™ HiSeq2000 in single-end mode, with each lane generating more than 200 million single-end 100 bp reads. Two plates of libraries (192 genotypes) were sequenced over 2 lanes and half a plate (44 genotypes) was sequenced over 1 lane. SNP calling was performed using the reference-guided TASSEL pipeline on the Red5 (version PS1.1.68.5 [54], an earlier version of the published genome [55]) and the 'Hongyang' reference genomes [56]. SNPs were filtered at the filtering criteria of 0.7 for coverage across all genotypes, gen-erating~44-50 K SNP sites from each genome across 29 pseudochromosomes and unassigned scaffolds in linkage group Chr 30. SNP markers were de-convoluted for each parent using the following criteria. First, the SNP markers were selected that were heterozygous (ab) in one parent and homozygous (aa) in the other parent and vice versa. This generated a set of markers which would theoretically segregate as 1:1 < ab×aa> (pseudo-testcross) and are unique to each parent; these were used for construction of the male and female genetic linkage maps individually.

Genetic map construction
Genetic maps using 'Hort16A' and P1 SNP markers were constructed using Joinmap3® (www.kyazma.nl). SNP marker data were processed in JoinMap using the 'CP' format for the population structure. Linkage groups were developed using default settings for grouping with modifications, including: a) the threshold range for Independence logarithm (base 10) of odds (LOD) started from a LOD score of 10 to 20 and, b) use of a regression mapping algorithm.

Segregation of microsatellite alleles within the SDR
The inheritance pattern of SSR markers within the SDR was investigated in an intraspecific A. chinensis var. chinensis mapping population previously described [20]. Twenty-one microsatellite markers were selected for analysis; eighteen of these amplify from within the terminal 6 Mb of chromosome 25, the remaining three amplify from the distal portion. The parents and eighty-seven progeny of this mapping population were screened with these markers in the same manner as previously described [57]. The sequences of these primers and annealing temperatures are given in Additional file 1: Table S1. Based on the pattern of segregation of alleles from fully informative, female-informative and male-informative markers, alleles were grouped into one of four groups depending on the chromosome from which they had originated i.e. group 1 originated from chromosome X 1 , group 2 originated from chromosome X 2 (inherited from the female parent), group 3 originated from X 3 , and group 4 from Y 1 (inherited from the male parent).

Whole-genome sequence alignment, variant calling and kinship analysis
Adaptors and low-quality or undetermined sequences of short-insert Illumina whole genome sequence were filtered and clipped using fastq-mcf (fastx toolkit, version 0.0.13) [58]. Paired reads at approximately 30x coverage were mapped to the 'Hongyang' genome [56] using bwa-mem 0.7.15 [59]. Variant calling was performed VCF files were phased using Beagle 4.0 and default values [61]. Fst, Tajimas Pi and kinship coefficients were generated for each 1-Mb window of pseudochromosome 25 using vcftools 0.1.14 [62] and kinship determined using the vcftools relatedness2 option [63]. Rates of linkage disequilibrium (LD) decay was determined in 1-Mb windows using PopLDDecay (https://github.com/BGI-shenzhen/ PopLDdecay). LD decay within each window was summarised as the median R 2 in the 1 kb-10 kb sub-interval.

Chromosome preparation
The female genotype of diploid A. chinensis var. chinensis used for this study, CK51_05, is the female parent of the mapping population III was used previously to construct a genetic map in A. chinensis var. chinensis [21]. Young flower buds at various stages were collected and placed immediately into 3:1 ethanol:acetic acid and stored at 4°C for at least one day. If buds were to be stored for more than two weeks, they were transferred to 70% v/v ethanol and stored at − 20°C until required.
Flower buds containing meiocytes at pachytene stage were identified by squashing an anther from each in FLP (formo:lacto:propiono) orcein and observing the meiotic stage. Once buds at pachytene had been identified, chromosome preparations were made following the method of Andras SC, Hartman TP, Marshall JA, Marchant R, Power JB, Cocking EC and Davey MR [64], using the remaining anthers in each bud and modified to use anthers rather than root tips. Anthers were hydrolysed in 1 M HCl at 37°C for 15-20 min and then digested for 80 min in the following enzyme mixture: 4% (w/v) Onozuka R10 cellulase (Merck 102,321), 4% (w/v) cellulase (Sigma C-9442), 2% (w/v) pectoylase (Sigma P-3026) and 1% (w/v) cytohelicase (Sigma C-8274) dissolved in 0.01 M citrate buffer pH 4.5. The drop techniques of Felsenstein J [65] and Henegariu O, Heerema NA, Lowe Wright L, Bray-Ward P, Ward DC and Vance GH [66] were subsequently used to prepare chromosome spreads.

Construction and screening of BAC library
The genomic DNA BAC library from A. chinensis var. chinensis CK51_05 was prepared by Bio S&T, Quebec, Canada and printed on 23 nylon filters in a 4 × 4 printing configuration. The genes ndhA (NADH-dehydrogenase subunit A) and cox2 (cytochrome c oxidase) were employed to estimate contamination from chloroplasts and mitochondria, respectively. Only 0.6% of the BAC library contained organellar DNA. From a sample of 309 BAC clones, we determined the average insert size to be 71.32 ± 46.15 kb and that 30% of the sampled BAC clones contained large inserts of 80-260 kb.
Polymerase chain reaction (PCR) probes developed from a non-polymorphic marker derived from the sex linked SmX marker [67], and two genetic markers flanking the sex locus (Ke225 and udkac096) [21] were used to identify the SDR of the female parent. Purified PCR probes were labelled non-radioactively with digoxigenin-11-dUTP (Roche Diagnostics), hybridized at 65°C overnight and detected as specified [68]. The corresponding BAC clones were isolated and the smallest clones for each of the three markers selected and labelled with biotin by nick-translation (Roche Diagnostics). These three clones were 47F17 -containing Ke225 (59.7 kb), 180D13 -containing SmX (49.3 kb), and 156B2 -containing udkac096 (85.2 kb).

Fluorescent in situ hybridisation (FISH)
The FISH procedure followed a previously published method [64], with some modifications. Briefly, probes (50-100 ng) were dissolved in 2X SSCP (0.3 M NaCl, 0.03 M sodium citrate, 0.04 M sodium dihydrogen phosphate, pH 6.5), 50% formamide and 10% dextran sulphate and denatured at 85°C for ten minutes. Thirty μL of this hybridisation mixture was added to each slide and the slide was covered with a plastic cover slip. Preparations were denatured for 5 min in 0.15 M NaOH in 70% ethanol and then dehydrated through an ice-cold ethanol series (70, 85, and 96% for 3 min each). Incubation was at 37°C for 24 h prior to a wash in 2X SSC (0.3 M NaCl, 0.03 M sodium citrate) at 42°C, followed by two stringent washes of 0.2X SSC for 15 min each at 42°C (Stringency~68%), and a final wash in detection buffer (0.1 M Tris, 0.15 M NaCl, pH 7.5) for 5 min at room temperature. The probe was detected using 50 ng/μL cy3-strepavidin conjugate (Sigma), 5% (w/v) bovine serum albumin (Sigma) in detection buffer for 60 min at 37°C. Slides were then washed twice in detection buffer containing 0.05% (v/v) Tween® 20 for 15 min and the chromosome preparations were counter-stained for 5 min in 1 mM 4′,6-diamidino-2-phenylindole, dihydrochloride (DAPI) (Sigma) in 1X PBS pH 7.4 and mounted in 40 μL of mounting solution (0.2% (v/v) 1, 4-diazabicyclo-[2.2.2]octane (DABCO) (Sigma), 50% (v/ v) glycerol in 1X PBS pH 7.4). Following storage at 4°C for 2-3 days before observation, chromosome spreads were observed with an Olympus Vanox AHT3 light microscope using epi-fluorescence, and images were captured with an RS Photometrics CoolSNAP digital camera. Two images were captured per cell, at excitation wavelengths of 358 nm and 550 nm. Entire images were then manipulated using Adobe® Photoshop Version 6.0. Chromosome measurements were made using the computer application MicroMeasure version 3.3 [69].