Transcriptome data from three endemic Myrtaceae species from New Caledonia displaying contrasting responses to myrtle rust (Austropuccinia psidii)

The myrtle rust disease, caused by the fungus Austropuccinia psidii, infects a wide range of host species within the Myrtaceae family worldwide. Since its first report in 2013 in New Caledonia, it was found on various types of native environments where Myrtaceae are the dominant or codominant species, as well as in several commercial nurseries. It is now considered as a significant threat to ecosystems biodiversity and Myrtaceae-related economy. The use of predictive molecular markers for resistance against myrtle rust is currently the most cost-effective and ecological approach to control the disease. Such an approach for neo Caledonian endemic Myrtaceae species was not possible because of the lack of genomic resources. The recent advancement in new generation sequencing technologies accompanied with relevant bioinformatics tools now provide new research opportunity for work in non-model organism at the transcriptomic level. The present study focuses on transcriptome analysis on three Myrtaceae species endemic to New Caledonia (Arillastrum gummiferum, Syzygium longifolium and Tristaniopsis glauca) that display contrasting responses to the pathogen (non-infected vs infected). Differential gene expression (DGE) and variant calling analysis were conducted on each species. We combined a dual approach by using 1) the annotated reference genome of a related Myrtaceae species (Eucalyptus grandis) and 2) a de novo transcriptomes of each species.

The present study focuses on transcriptome analysis on three Myrtaceae species endemic to New Caledonia (Arillastrum gummiferum, Syzygium longifolium and Tristaniopsis glauca) that display contrasting responses to the pathogen (non-infected vs infected). Differential gene expression (DGE) and variant calling analysis were conducted on each species. We combined a dual approach by using 1) the annotated reference genome of a related Myrtaceae species (Eucalyptus grandis) and 2) a de novo transcriptomes of each species.
& 2019 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

Subject area
Genetics and Transcriptomics More specific subject area Transcriptomics of Myrtaceae species Type of data Table, figure How data was acquired Leaves of individual plant from three endemic Myrtaceae species from New Caledonia were sampled for total RNA extraction. Pairedend library were prepared and RNA-Sequencing was performed by the Illumina HiSeq™ 2500 system. The obtained data was subjected to 1) de novo transcriptome assembly, 2) alignment to a reference genome and de novo transcriptome assembly, 3) differential gene expression and 4) variant calling analysis.

Experimental features
For the RNA-Sequencing and transcriptome analysis, a total of 24 leaves samples from three host species have been collected: three infected and three non-infected individuals from A. gummiferum, two infected and five non-infected individuals from S. longifolium, four infected and one non-infected individuals from T. glauca from the nursery, four infected and two non-infected individuals from T. glauca from a natural population in a protected reserve. Data source location The nursery was located in Farino, South Province, New Caledonia (Long 165.772024:, Lat: -21.663800). The protected reserve was located in Bois du Sud, South Province, New Caledonia (Long 166.758640:, Lat: -22.169974) Data accessibility All raw data for Arillastrum gummiferum, Syzygium longifolium and Tristaniopsis glauca and the processed data (de novo transcriptome assemblies, transcriptome annotations and differential gene expression files) obtained in this study were deposited in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) with the Superseries accession number GSE106750 and the subseries accession numbers GSE106735, GSE106736, GSE106738, GSE106740, GSE106741, GSE106746, GSE106747 and GSE106749. The Superseries is available at https://www.ncbi.nlm.nih.gov/geo/ query/acc.cgi?acc ¼GSE106750 The VCF files can be downloaded at http://myrtaceae-omics.south green.fr/node/8 provided you are logged and have accepted the Terms of Use of these data. Related research article Soewarto

Value of the data
These are the first de novo transcriptomes of Arillastrum gummiferum, Syzygium longifolium and Tristaniopsis glauca.
The obtained transcriptomes data will be useful for further studies of the evolution of Myrtaceae or comparative genomics.
The data (de novo transcriptomes, reads assemblies, differential gene expression, SNP calling, etc…) will provide new valuable genetic resources for investigations of myrtle rust interactions and resistance-related pathways within Myrtaceae plant family.

Data
The transcriptomes were extracted from leaf samples of individuals of Myrtaceae species endemic to New Caledonia (Arillastrum gummiferum, Syzygium longifolium and Tristaniopsis glauca) naturally infected by the fungal pathogen Austropuccinia psidii (myrtle rust) and displaying contrasting responses to the infection. Due to facilities restrictions regarding the artificial inoculation of the pathogen in New Caledonia, the natural infection method has been chosen to rapidly screen the infection status of these species. A previous study based on natural myrtle rust infection concluded that A. gummiferum, S. longifolium and T. glauca were all susceptible to the disease and displayed variations in the disease incidence and severity [1].
A total of 24 leaves samples from the three species have been collected for the RNA-Sequencing on: three infected and three non-infected individuals from A. gummiferum, two infected and five noninfected individuals from S. longifolium, four infected and one non-infected individuals from T. glauca from a commercial nursery located in Farino (FAR), four infected and two non-infected individuals from T. glauca from the protected reserve of Bois du Sud (BDS). After sequencing, 48 RNA-Seq fastq files (paired end Illumina HiSeq) have been generated.
A total of four de novo transcriptomes have been implemented in this study. The individual displaying the most reads at sequencing was chosen as a representing of each species (A. gummiferum, S. longifolium, T. glauca) and location (FAR, BDS).
Height differential gene expression analysis (DGEA) were conducted: four by using each species's reads mapped on the annotated reference genome of Eucalyptus grandis and four by using the reads mapped on each species corresponding de novo transcriptome assembly.
Sixteen variant calling analysis where conducted for the four group of species using two types of reference sequences for the mapping (E. grandis genome and de novo transcriptome assemblies) and two types of variant calling methods (in-house SNP calling and GATK methods).

Plant material
Individuals from the three following species Arillastrum gummiferum, Syzygium longifolium and Tristaniopsis glauca, were sampled from a commercial nursery located in Farino, New Caledonia. The infection status of each individuals has been monitored at least one month before the sampling to  glauca have also been sampled from a natural population occurring in the protected reserve of Bois Du Sud, New Caledonia. To distinguish the individual plants from the two locations, the natural population from Bois Du Sud will be referred as T. glauca BDS and the plants originated from the nursery will be referred as T. glauca FAR. Thus, this study counts four group of species and population as followed: A. gummiferum, S. longifolium, T. glauca FAR and T. glauca BDS. When A. psidii successfully infects a host, disease symptoms can appear within 12 days and are visually characterized by the formation of pustules covered by yellow and powdery urediniospores [2]. The disease can infect various plants parts including actively growing leaves, shoots, fruits, flowers and buds [2,3]. Although natural infection was showed to be effective to infect at least once every host plants standing in the nursery during the previous monitoring [1], it does not constitute a reliable way to conclude on the resistant status of the individual plants that did not display any sign of symptoms. Thus, the present study will consider that all the individuals that showed signs of myrtle rust symptoms at the sampling time or during a previous monitoring will be classed as infected; and all the individuals that never showed sign of myrtle rust symptoms will be classed as non-infected.
All the samples have been harvested at the same time, in April the 28th, 2015 for the plants in the nursery and in May the 1st, 2015 for the plants in the protected reserve of Bois Du Sud. Each time the most recent leaves were always chosen for harvesting ( Fig. 1). All samples were frozen at the time of collection and then stored at À 80°C until total RNA was extracted. The sampling details are provided in Table 1.

Total RNA extraction
Total RNA was extracted using 3-10 g of fresh material from the cetyltrimethylammonium bromidebased protocol (CTAB) [4]. Briefly, 3-10 g of frozen leaves from each individual was ground in liquid nitrogen using a mortar and pestle. Then 3-10 mL of pre-heated extraction buffer (2% CTAB, 2% polyvinylpyrrolidone, PVP-40 (2% w/v), 2% β-mercaptoethanol, 100 mM Tris-HCl, 25 mM EDTA and 2 M NaCl) was added to the ground samples, and incubated at 65°C for 30 min. An equal volume of mixture of chloroform:isoamyl alcohol (24:1) was added and mixed immediately for 2 min using a vortex mixer. The samples were then centrifuged at 10,000g for 10 min. The upper aqueous phase was transferred to new tubes and 1/3 volume of 10 M LiCl was added. The samples were mixed and stored at 4°C overnight. The samples are then centrifuged at 18,000g for 20 min. The supernatant was removed and the pellet is washed with 75% ethanol and air-dried. The pellet was suspended in 30 μl of RNase-free water and 70 μl of SSTE buffer (1M NaCl, SDS (0.5% w/v), 10 mM Tris-HCl, 1 mM EDTA). An equal volume of acid:phenol:chloroform: isoamyl alcohol (25:24:1) is added to each sample, and vortexed. The samples were then centrifuged at 12,000 g for 10 min and the upper aqueous phase was transferred to new tubes with 2 volumes of cooled ethanol 100% and 1/10 volumes of NaAc (pH5.2). The samples are then mixed and incubated at À 20°C for 2 h. The samples were centrifuged at 18,000g for 20 min. The supernatant was removed and the pellet was washed three times with 75% ethanol before being air-dried. Finally, the pellet is resuspended in 30 μl of RNase-free water. The RNA extraction is followed by removal of DNA with the TURBO DNA-free ™ kit (Ambion) according to the manufacturer's instructions. RNA quantity and quality control was performed using a 2100 Bioanalyzer (Agilent Technologies).

cDNA library preparation and sequencing
Paired-end Illumina mRNA libraries were generated using the TruSeq RNA-Seq Sample Prep kit according to the manufacturer's protocol (Illumina Inc., San Diego, CA, USA). Briefly, Poly-A containing mRNA molecules were isolated using poly-T oligo-attached magnetic beads. The purified mRNA was then chemically fragmented. Reverse transcription of firstand second-strand of cDNA are performed and were followed by end repair. Single 'A' nucleotide is added to the 3 0 ends of each fragment before ligation of adapters. The purified cDNA templates were enriched by PCR to form libraries of 300 pb. Each indexed cDNA library was verified and quantified using a 2100 Bioanalyzer (Agilent Technologies). The final libraries were then quantified by qPCR with the KAPA Library Quantification Kit for Illumina Sequencing Platforms (Kapa Biosystems Ltd, SA) and normalized to 20 nM before being pooled. The quality of our samples was suitable for the Illumina HISAT 2500 requirement sequencing (RIN between 6.9 and 8.4). Details of RNA samples were supplied in Supplementary table 1. Sequencing was then conducted on a single lane of a flow cell on Illumina HiSeq ™ 2500 (Genotoul platform, INRA) as paired-ends reads of length 150 bp.

Raw reads cleaning
RNA-sequencing for the three species, produced a total of 668,165,595 raw reads. The quality of paired-end raw reads in fastq format was assessed using FastQC sofware (http://www.bioinformatics. babraham.ac.uk/projects/fastqc/) with the script arcad_hts_0_fastqc_in_chains.pl. Raw reads were then processed with Cutadapt [5] using the script arcad_hts_1_cutadapt_in_chain.pl and the TruSeq index sequences corresponding to the samples. We also used Cutadapt to improve mean reads quality by trimming the start and the end of each reads. We then filtered the reads on the basis of their mean quality score, keeping those with a mean quality higher than 30 using arcad_hts_2_Filter_Fast-q_On_Mean_Quality.pl. Reads with length inferior to 35 bp were discarded. Thereafter, single reads s fter quality-based trimming Adapter content of the raw data The same data after removing adapter Per base sequence content The same data after trimming ''Per base quality score'' of the raw data. The same data a has been performed. (i.e. those for which the mate pair was discarded in the previous steps) were separated from pairedreads using arcad_hts_3_synchronized_paired_fastq_end.pl. Each processing steps for reads cleaning were followed by a FastQC quality control of reads (arcad_hts_0_fastqc_in_chains.pl). Each preprocessing cleaning steps is illustrated in Fig. 2. After the cleaning stage, we kept 96% of the initial reads (644,707,222 reads). An example of the efficiency of quality control on our data is shown in Fig. 3. A summary of the RNA-seq raw and clean data is presented in Table 2 and detailed in Supplementary table 2.   [6]. The aligned data were passed to StringTie [7] for transcript assembly. A reference genome annotation file in GFF3 format (intron/exon positions) was provided to guide the transcripts assembly. The resulting binary files (.bam) were then filtered on quality using SAMtools view (with parameters -q 0 -F 4) to remove unmapped and multimapped reads [8]. Finally, the reads were sorted using SAMtools sort. Mapping statistics were verified using SAMtools flagstat [8] ( Table 3). We obtained mapping rates of approximately 68% for A. gummiferum, S. longigolium, T. glauca-FAR and T. glauca-BDS. And around 5% of the total reads mapped correspond to singletons (Table 4).

Calling
SNPs using E. grandis reference genome. Variant calling algorithms compare mapped reads to a reference genome and identify potential variants. The analysis pipeline we used is illustrated in Fig. 4. We followed the best practices recommended by the GATK (Genome Analysis Toolkit) pipeline [9]. Consistent with GATK's recommendations, mapped reads against E. grandis reference genome have been submitted to cleaning process before SNP calling step. We used Picard tools to remove PCR duplicates (MarkDuplicates) and reorder reads to match the contigs according to the reference genome (ReorderSam). Then, we used the GATK tool SplitNCigarReads (with parameters -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS), which splits reads into exon segment and hardclip any sequences overhanging into the intronic regions. The reads were realigned around INDELS and base quality values were recalibrated using RealignerTargetCreator and IndelRealigner GATK tools. Once the reads were pre-processed with Picard and GATK tools, variant calling was undertaken by two programs: Haplotype Caller (GATK) and an in-house script (Martin et al., in prep., Baurens et al., in prep.). Contrary to GATK, which uses statistics based on population genetics (which is not the case here), the in-house program only count the number of reads supporting each bases (A, T, G, C) at each covered sites for each accession. Based on this base count for each accession, a genotype was emitted based on a binomial test. We obtained between 3 and 5 million of SNPs with Haplotype Caller for each species, and the in-house script identified between 8 and 11 million of SNPs (Supplementary Table 3).  (Fig. 4). Because in-house SNP calling method implied to detect any difference between the reads and the reference (including false positive SNP), the initial outputs were too large and must be pre-filtered in order to allow comparison with those resulting from GATK SNP calling method. Firstly we applied an in-house pre-filter script (VcfPreFilter.1.0.py) on SNPs resulting from the in-house SNP calling method in order to remove SNPs due to sequencing errors. This pre-filter included the following parameters: a minimum site coverage (5 reads), a maximal site coverage (1000 reads) per SNP, a minimum allele frequency (0.01) and a minimum allele coverage (2 reads per allele). Then for all datasets (from GATK and in-house SNP calling methods) we used an in-house VCF filter script (VCF_Filter.py) to keep only the more robust SNPs according to the following criteria: a minimum coverage (10), a maximal coverage (1000), a minimum frequency (0.05), a minimum allele coverage (3) to keep a genotype; and no missing data were allowed to keep a variant site. Another in-house script was then used to remove monomorphic and multiallelic variants. And finally, as recommended by GATK, we used the Variant Filtration program to remove clustered SNP: in a window size of 10 pb we considered that 3 SNPS constitute a cluster. The final SNP set comprised the identified SNPs that had passed all filters. We compared the overlap among SNP positions obtained for each calling methods using vcf-compare from VCF tools [10] (Table 4 ).
With E. grandis genome as reference for mapping, we found around 30% of shared SNPs positions using Haplotype Caller, and around 60% of shared SNPs positions using In-house script ( Table 4).

De novo transcriptome assembly
We chose the individual with the highest read number after cleaning for each species for the de novo assembly as described by Sarah et al. [11]. De novo transcriptome assembly was carried out using Trinity software with 50 GB of memory (Inchworm, Chrysalis, and Butterfly modules) [12]. Briefly, overlapping k-mers are extracted from the cleaned paired-reads. Inchworm module assembled the reads into contigs. Next, Chrysalis module ranked Inchworm contigs into clusters and constructed complete de Bruijn graphs for each cluster. Finally, Butterfly module processed the individual graphs in parallel to reconstruct transcript sequences in a manner that reflects the original cDNA molecules. To avoid redundant transcripts, we kept the longest isoform for each "trinity gene". Assembly statistics (N50, contig length, GC content, etc) were computed by TrinityStats.pl embedded in Trinity (Table 5; Supplementary table 2). Transcriptome assembly for A. gummiferum resulted in 117 839 putative transcripts with an average contig length of 501 bp and N 50 of 1378 base pairs. Transcriptome assembly for S. longifolium yielded a total number of 89,782 putative transcripts with an average contig length of 530 bp in length and N 50 of 1406 base pairs. Transcriptome assembly for T. glauca-FAR resulted in 108,823 putative transcripts, with an average contig length of 547 bp and N 50 of 1396 base pairs. And finally, de novo transcriptome assembly for T. glauca-BDS resulted in 74,684 putative transcripts with an average contig length of 525 bp in length and N 50 of 1315 base pairs. However, Trinity de novo assembly resulted in larger number of transcripts than expected number of genes, likely because of alternative splicing. To avoid redundant transcripts, we kept the longest isoform for each "gene" identified by TRINITY (unigene). Overall size of filtered de novo assembly Table 5 Statistics of the de novo transcriptome assembly for each species using Trinity assembler.

Aligning to the de novo transcriptomes
We used the Burrows-Wheeler alignment tool (BWA-MEM) to map the cleaned reads from the 24 samples to the corresponding de novo transcriptome assemblies of each species [13]. The read aligner, BWA-MEM, aligns each mate of a paired-end read at the same time and produces SAM/BAM files containing the alignments. Samtools view from the SAMtools 1.3 package (http://www.htslib.org/doc/ samtools-1.1.html) was used to sort and index the BAM files by coordinate and remove multimapped reads (parameters -q 0 -F 4). Alignment statistics for the mapping to de novo transcriptome assembles are displayed in Table 3.

Calling SNPs using de novo transcriptome assemblies
The analysis pipeline we used is illustrated in Fig. 4. We processed the aligned reads to a cleaning process with the Galaxy [14] instance of the South Green platform http://galaxy.southgreen.fr/galaxy/. Based upon established GATK best practices, we used Realigner targer creator and IndelRealigner programs (LOD ¼ 5) to fix the misalignments due to the mapping process and perform a local realignment near the INDELS [15,16]. Variant calling was launched both with Unified Genotyper (GATK), using a minimum phred-scaled confidence threshold of 30 [6], and an in-house scripts (Martin et al., in prep., Baurens et al., in prep.). We obtained between 600 and 900 thousands of SNPs with Unified Genotyper for each species, and the in-house script identified between 9 and 12 millions of SNPs (Supplementary Table 4).

Variant filtering from mapping to the de novo transcriptomes
We applied the same filtering process as explained earlier with the variant filtering from mapping to E. grandis genome. We compared the overlap among SNP positions obtained for each calling methods using vcf-compare from VCFtools [10] (Table 6). With de novo transcriptomes as reference for mapping, we found around 2=3 of the SNPs identified by Unified Genotyper (GATK) and in-house methods sharing identical positions on A. gummiferum and S. longifolium ( Table 6). The proportion of identical SNPs positions is even bigger with T. glauca-BDS and T.glauca-FAR, up to 96% using the inhouse-script and 83% using Unified Genotyper (Table 6). Table 6 Overlapped and uniques SNPs called using two different calling methods (GATK and in-house script) from mapping using de novo transcriptomes.

Species
Methods

Differential gene expression analysis in infected versus non-infected plants
EdgeR (Bioconductor package) was used to identify Differentially Expressed Genes (DEGs) in the three species aligned with the two types of reference (E. grandis genome and de novo transcriptome) [17].
From the raw read counts, EdgeR normalized the size of the sample libraries and computed genewise tests for differences in the means between the groups of infected samples versus the group of non-infected samples. It outputted CPM (counts per million), log of fold change (logFC) between the two groups, along with the corresponding p-value and false discovery rate (FDR).
Prior to the normalization step, the genes were filtered. Only the genes whose sum of CPM values (calculated on all the samples) was greater than 1, and which were expressed in at least 2 samples, were kept. Differentially expressed genes were selected based on a fold change Z 2 (logFC Z1 and logFC r1) and an FDR-adjusted p value threshold of 0.05. An FDR of 0.05 implied that we were willing to accept that 5% of the differentially expressed genes were false positives.
The total number of identified differentially expressed genes, along with their up-/down-regulation, are summarized in Table 7 and displayed as a plotSmear in Supplementary Figure 1. In EdgeR, dispersion was estimated on a common dispersion basis for all species. The differential expression of genes was analyzed in infected individuals compared to non-infected ones. When using the E. grandis reference genome we showed that 12.69% of the total expressed genes were differentially expressed in A. gummiferum, 10.32% in S. longifolium, 2.45% in T. glauca-BDS and 1.75% in T.glauca-FAR. When using the de novo transcriptome of each species, around 5.6 % of the genes were differentially expressed in A. gummiferum, 8.46% in S. longifolium, 1.37% in T. glauca-BDS and 1.24% in T.glauca-FAR.
As gene expression differences existed between the two groups of individuals (non-infected/ infected), it should be expected that biological replicates of the same condition will cluster together. We used a multidimensional scaling (MDS) plot to see a spatial configuration of how similar or dissimilar the non-infected and infected individuals were according to the reference used (Supplementary figure 2). We observed clustering of individuals with the same phenotype for A. gummiferum When using the E. grandis genome as a reference for mapping, a cross-comparison of differentially expressed genes (DEGs) was presented as a Venn diagram illustrating the overlapped DEGs (Fig. 7). We found 33 over-expressed genes in infected individuals that overlapped between the three species A. gummiferum, S. longifolium and T. glauca (FAR/BDS) ( Table 8), but only 28 had an identified product. No under-expressed genes were found in common between the three species, suggesting that the resistance process toward A. psidii is potentially specific to each species. Of the 33 over-expressed genes in infected plants, we found several genes potentially involved in disease response processes (Table 8), such as LOC104438326 coding for a pathogenesis-related protein STH-2-like and various chitinase coding genes (LOC104415213, LOC104419011, LOC104456214, LOC104456215, LOC104456217, LOC104456219, LOC104456220, LOC104456221, LOC104456223). Chitinases are enzymes that hydrolyze the polymer chitin from many chitinolytic biotic aggressors (fungi, bacteria, viruses, viroids) [18]. They are considered as pathogenesis-related proteins playing a crucial role in resistance against pathogens [19]. We also identified gene LOC104445691 coding for a probable WRKY transcription factor 31 isoform X1. This transcription factor specifically interacts with cisacting elements of plant defense genes that are expressed in reaction to an elicitor. Elicitors are compounds of pathogen origin stimulating any type of plant defense.

Annotation
We ran FrameDP V1.2.2 software with default parameters [20], for the prediction of coding regions in the unigenes using the E. grandis protein database from Universal Protein Resource (UniProt-Swissprot). FrameDP V1.2.2 software automatically reverse-complement the sequences, however all the previous analysis were performed on the initial DNA strand orientation. Therefore, we changed the reverse-complemented orientation from FrameDP results to the initial genes orientation. The polypeptides sequences obtained were submitted to the Basic Local Alignment Search Tool (BLAST) to search against homologs in the UnitProt (Swiss-Prot and Trembl) databases. We identified 45,517 protein-coding genes in T. glauca-BDS, 50,173 protein-coding genes in S. longifolium, 57,366 proteincoding genes in T. glauca-FAR and 65,410 protein-coding gene in A. gummiferum.   6. Numbers of SNPs per calling methods and for each studied species using de novo transcriptome of each species as reference for mapping.

Genome browser
To make easier the analysis of the massive amounts of genetic information that was generated during this study, we used the JBrowse tool from GMOD project (Generic Model Organism Database project) [21]. JBrowse is a web-based genome browser, allowing interactively visualizing and exploring a large genomic dataset [22]. Each track consists of a particular type of sequence feature along a reference sequence, as showed in Fig. 8. We deployed five different JBrowses with the data from each species, taking first the E. grandis genome as a reference for alignment, then the de novo transcriptomes (A. gummiferum, S. longifolium, T. glauca-FAR and T. glauca-BDS). To build the JBrowses, we rely on a released workflow [23] and we imported the following files: reference sequence files (de novo transcriptomes (.faa) or E. grandis reference genome (.fna), alignments files (.bam), variants calling files (.vcf), annotations files (.gff3), and EdgeR output files containing the RPKM (expression level) of each locus for the samples (.tsv). To load the following data we used the following perl scripts provided by the JBrowse developers: prepare-refseqs.pl for fasta files and flatfile-to-json.pl for .gff3 and .tsv files. To load the others data types we had to edit two files tracks.conf and trackList.json, respectively displaying the dataset-specific names and configuration. We also used a JBrowse plugin called Multi-VariantViewer (multivariantviewerjbrowsepluging) to display the corresponding genotype for each SNP from each individual of this study (https://github.com/elsiklab/multivariantviewer). These genotypes are displayed in three different colors: cyan for heterozygotes, grey for homozygotes for the reference allele and deep blue for the homozygotes of the alternative allele. To load the MVV tracks (sample name, category, and genotype), we edited the trackList.json from each JBrowse. The JBrowses were integrated into the open source CMS Drupal [24], which is distributed under the terms of the GNE's Not Unix General Public License (http://myrtaceaeomics.southgreen.fr).