The high molecular weight genomic DNA was isolated from muscle tissue of a single shrimp using QIAGEN genomic-tip 100/G kit (Qiagen, Germany). The size selection was done according to the protocol described under “Procedure and Checklist–20 Kb template preparation using BluePippin size selection system.” The sequencing libraries were prepared using SMRTbell Express Template Prep Kit 2.0 (Pacific Biosciences, United States) and assessed for quality and quantity using the Pippin pulse field inversion gel electrophoresis system (Sage Science, United States). The sequencing was performed on Pacific Biosciences Sequel system using magnetic bead loading and 600-min movies. The Pacbio subreads of 5 Kb and longer were used to generate primary contig-level assembly.

For short-reads sequence data, 15 paired-end sequencing libraries (5 each with 350, 550, and 650 bp insert size) were prepared using Illumina Truseq Nano DNA Library Prep Kit (Illumina, United States). The PCR enriched libraries were sequenced on Illumina NextSeq500 using 2 × 150 bp chemistry in paired-end mode. The raw reads were trimmed for adapters and poor-quality bases/reads using a sliding window algorithm in paired-end mode as implemented in Trimmomatic V0.36 (Bolger et al., 2014). The trimmed reads equal to or longer than 75 bases only were further used for polishing (correction of base errors and indels) the contigs in primary genome assembly.

For HiC data, the library was prepared using Arima Hi-C kit (Arima Genomics, United States). Tissue crosslinking and proximity ligation was performed following Arima Hi-C animal tissue protocol. The Illumina sequencing library was prepared from proximally ligated DNA using Swift Accel-NGS 2S plus kit (Swift Biosciences, United States) following manufacture’s guidelines. The library was sequenced on Illumina NovaSeq6000 platform with 150 bp paired-end mode. The HiC reads were used for scaffolding the contigs in primary assembly.

Total RNA was extracted from the gill, hepatopancreas, muscle, pleopod and heart tissues using TRIzol method and utilized for cDNA synthesis. The cDNA sequencing libraries were prepared following the protocol of “Sure select strand-specific RNA library prep for Illumina multiplexed sequencing” (Agilent Technologies, United States). Sequencing was performed on Illumina NextSeq500 with 2 × 150 bp paired-end chemistry. Similar procedures were followed to generate pooled-RNAseq datasets for gill, hepatopancreas and muscle tissues, wherein, each sample was derived by pooling cDNA from nine different shrimp (three each from Chennai, Kanyakumari, and Puri Coast of India). The Pacific Biosciences Iso-Sequencing data for gills, hepatopancreas, muscle and pooled larvae of P. indicus as described previously (Katneni et al., 2020) along with the RNAseq data were used for genome annotation and the pooled-RNAseq data were used for identification of SNPs in candidate transcripts.

The Pacbio subreads with a minimum length of 5 Kb at 73× coverage were used in WTDBG2.5 (Ruan and Li, 2020) to generate a de novo contig-level assembly. To improve the quality of the genome, the contigs were polished in two steps, using Pacbio subreads in first step and Illumina short-reads in the second. In first step, polishing of contigs was performed with the Arrow algorithm of variantCaller tool¹ using Pacbio subreads. In the second step, POLCA tool (Zimin and Salzberg, 2020) was used with Illumina short-reads to correct error bases and indels in contigs. The polished contigs were then scaffolded with 3D-DNA pipeline (Dudchenko et al., 2017) using HiC reads to generate the final draft assembly.

The quality of the final assembly was evaluated by aligning the Pacbio subreads and Illumina short-reads on to the genome scaffolds using bwa-mem v0.7.15-r1140 (Li and Durbin, 2009) and bowtie2 v2.3.4.3 (Langmead and Salzberg, 2012), respectively. Assembly evaluation based on alignment statistics was also performed with Illumina RNAseq reads and Pacbio IsoSeq transcripts using Hisat2 v2.2.0 (Kim et al., 2019) and GMAP v2020-06-30 (Wu and Watanabe, 2005), respectively. In addition, the genome was assessed for completeness by benchmarking against the arthropoda_odb10 (September 10, 2020) dataset of BUSCO single-copy orthologs (Seppey et al., 2019). Further, the assembly of P. indicus was compared with other shrimp genomes for assessment of sequence contiguity based on un-gapped contig length distribution and number of gaps in the finished genomes.

Homology-based annotation of repeat elements was performed with Penaeidae subset (Taxonomy ID:6685) of RepBase library using the RMBlast search of RepeatMasker (Jurka et al., 2005)² module implemented in OmicsBox v1.3.11 (Bioinformatics, 2019). For utility in gene prediction, the genome was soft masked for repeat regions with the same procedure excluding the low complexity and simple repeats.

Structural annotation of protein-coding regions in P. indicus genome was carried out by combining the gene models obtained from short-read RNAseq data, long-read IsoSeq data and proteins from related species (Supplementary Table 1), with ab initio gene predictions. The ab initio gene models were predicted on masked genome using Augustus v3.3.3 (Stanke et al., 2006) and GeneMark-ES v4.59 (Lomsadze et al., 2005) self training module. While predicting gene models in Augustus, hints generated by aligning IsoSeq data on to the genome using GMAP v2020-06-30 (Wu and Watanabe, 2005) were also given as input. The PASA v2.4.1 (Haas et al., 2003) was used to generate valid gene structures from IsoSeq data which utilizes near perfect alignments made by GMAP v2020-06-30 and BLAT v36 (Kent, 2002; Wu and Watanabe, 2005) to the genome. Similarly RNAseq data was aligned in a splice aware manner to genome using Hisat2 v2.2.0 (Kim et al., 2019) and gene models were generated using StringTie v2.1.4 (Pertea et al., 2015), from which likely coding sequences were identified using TransDecoder v5.5.0.³ Proteins from the related species were aligned to the genome using GenomeThreader v1.7.3 (Gremme, 2012) to derive valid gene models. All these ab initio and evidence based predicted gene models were combined as a weighted matrix in Evidence Modeler (Haas et al., 2008) to obtain the final set of non-redundant consensus gene models. Homology based annotation was performed using blastx against non-redundant protein database of Genbank and UniProt database to obtain functional description of the genes. The Interproscan and EggNOG Mapper module of OmicsBox v1.3.11 were used to obtain annotations of protein domains and orthology based annotations respectively (Bioinformatics, 2019). The annotations were merged and the final gene ontology annotations were obtained using OmicsBox v1.3.11. Pathway maps for the annotated genes were generated by mapping against the KEGG database (Kanehisa and Goto, 2000).

Protein-coding gene sets of 17 species in the phyllum, Arthropoda (Supplementary Table 2) including the P. indicus were subjected to gene family analysis using OrthoMCL v2.0.9 (Fischer et al., 2011). For this analysis, other Metazoan species belonging to the phylla, Mollusca, Chordata, and Echinodermata were also included. The datasets downloaded from NCBI⁴ were filtered based on length (minimum 50 amino acids) and selection of the longest isoform. Good protein list from 21 species was subjected to an all-versus-all search using blastp (Altschul et al., 1990) and then single copy orthologous gene (SCOG) sequences were extracted from orthomcl groups. MUSCLE v3.8.1551 (Edgar, 2004) was used to generate multiple sequence alignments which were then trimmed using “trimAl” v1.4 (Capella-Gutiérrez et al., 2009) tool. The FASconCAT v1.04 (Kück and Meusemann, 2010) was used to concatenate the alignments, ProtTest v3.0 (Darriba et al., 2011) was used to find the best evolutionary model and RAxML v8.2.12 (Stamatakis, 2014) was used to build Maximum Likelihood tree. Visualization of the tree was done using FigTree v1.4.4 (Rambaut, 2009).

Pooled-individual RNAseq datasets were trimmed with Trimmomatic v0.39 (Bolger et al., 2014) to remove poor quality reads and bases. Good quality reads were aligned using TopHat v2.1.1 (Trapnell et al., 2012) on to the gene sequences indexed with bowtie2 v2.3.4.3 (Langmead and Salzberg, 2012). The generated bam file was sorted using SAMtools v1.2 (Li et al., 2009) and then processed in bcftools-1.3.1 (Li, 2011) to generate variant call format (VCF) file. Those SNPs with a raw read depth of ≥20 at SNP site, at least 10 reads each supporting the reference and alternative alleles and phred quality scores of ≥100 were extracted from VCF file as good quality SNPs. The non-synonymous coding SNPs were analyzed for the possible functional impact on the proteins harboring them, using a standalone version of the PANTHER Coding SNP Analysis tool, PANTHER PSEP v1.01 (Tang and Thomas, 2016). The tool analyses the SNP by first identifying the Panther family of the protein using blastp (Altschul et al., 1990), followed by retrieving the ancestor sequence of the family and tracing the query protein though evolution, then reporting how long the SNP position was conserved in Millions of years. Finally based on the predicted values the results are classified as probably benign, possibly damaging or probably damaging if the values are <180, ≥180, or ≥380, respectively.

Flow cytometry analysis using propidium iodide stained hemocytes indicated a genome size of 2.47 Gb for P. indicus (Swathi et al., 2018). Early attempts made by us to assemble P. indicus genome with 145× coverage of Illumina short reads sequence data were unsuccessful (Supplementary Note). Various assemblers like SOAPdenovo2 (Luo et al., 2015), CLC Genomics Workbench v10.0.1 (CLCbio, Denmark), and Platanus v1.2.4 (Kajitani et al., 2014) produced primary assemblies with large number of contigs. A final assembly reduced with Redundans v0.14a (Pryszcz and Gabaldón, 2016) contained 358,878 contigs covering about a quarter of genome length (607.78 Mb) with N50 of 1698 bases. Though short-reads could not generate a quality genome, k-mer analysis on them indicated high repetitive nature of P. indicus genome thereby suggesting the necessity of long sequence reads to assemble repeat-rich P. indicus genome.

About 73× coverage of PacBio Sequel long reads (Supplementary Table 3) generated using DNA of a P. indicus female shrimp were processed in WTDBG2.5 (Ruan and Li, 2020) to generate primary contigs. These contigs were corrected for error bases/indels in POLCA tool (Zimin and Salzberg, 2020) using Illumina reads (Supplementary Table 4). The assembled genome was of 1.98 Gb length and consisted of 12,051 contigs with N50 of 1.4 Mb (Table 1). Among large crustacean genomes of more than 1 Gb size, only the genome of Eriocheir sinensis (Tang et al., 2020) has a better contig N50 length (3.16 Mb) than obtained for P. indicus. Compared to the other shrimp genomes generated for P. vannamei (Zhang X. et al., 2019; Yuan et al., 2021), P. monodon (Uengwetwanit et al., 2021), and P. chinensis (Yuan et al., 2021), the assembly obtained for P. indicus has 24–, 24–, and 18-fold improvement, respectively for contig N50 length.

The scaffolding of polished contigs with 104× coverage (1.7 billion reads/258 Gb) of HiC reads in 3D-DNA (Dudchenko et al., 2017) resulted in an assembly with 11,168 scaffolds. The final assembly was of 1.935 Gb length (∼78% of genome) with N50 and longest scaffold length of 34.4 and 51.57 Mb, respectively (Table 1). Final assembly obtained for P. indicus had 44 scaffolds that were longer than 5 Mb, which equaled the reported haploid chromosome number (Figure 1). These 44 longest scaffolds assumed as pseudochromosomes span about 1.57 Gb length and cover about 81.3% of the assembly (Supplementary Table 5). The genome length in P. indicus that is accounted by the pseudochromosomes alone is almost same as the total genome length presented for P. vannamei and P. chinensis (Yuan et al., 2021). The assembly obtained for P. indicus genome in the present study is the only shrimp genome that meets the reference standard of 1 Mb contig N50 for primary contigs and 10 Mb scaffold N50 length for final scaffolds (Figure 2A).

Benchmarking of the P. indicus genome with arthropoda_odb10 (September 10, 2020) BUSCO gene groups (n = 1013) indicated the presence of 77% of complete orthologs in the genome. About 3.5% of BUSCO orthologs were fragmented and the remaining 19.5% were missing. Against the same dataset, the P. vannamei genome (Zhang X. et al., 2019) had 76% complete, 4.1% fragmented, and 19.6% missing BUSCO orthologs. In comparison, the proportion of missing genes reduced to 14.4% if odb9 dataset (n = 1066) was used for benchmarking. On the other hand, the P. indicus genome when assessed with gVolante 1.2.1 (Nishimura et al., 2017) against CEG ortholog set in CEGMA pipeline (Parra et al., 2007), has a completeness score of 98.39%. We observed that the completeness scores changed with the orthologous gene set and the data release associated with it. We suggest exercising caution to draw conclusions while directly comparing quality of genomes based on different gene sets. Moreover, with the shrinkage of orthologous gene datasets due to increased number of sequenced genomes, we opine that their usage soon becomes debatable.

Read mapping statistics indicated high quality of the assembly with high percentage of reads aligning on to the genome (Supplementary Table 6). For DNA data, about 98.89% of Pacbio reads and 95.4% of Illumina reads could be aligned back to the genome scaffolds. For RNA-based data, about 94.29% of RNAseq reads and 99.25% of Pacbio IsoSeq transcripts could be mapped on to the genome. High quality of the P. indicus genome assembly presented in this study has been demonstrated with high mapping statistics and BUSCO completeness scores as comparable to other shrimp genomes.

The sequence contiguity assessed based on un-gapped contig lengths and the number of gaps in assembly is another quality metric that we used in this study to assess the quality of P. indicus in relation to other shrimp genomes. A distribution of un-gapped contig lengths (Figure 2B) indicated that P. monodon is the only shrimp genome other than the P. indicus that contained un-gapped contigs (n = 3) of over 1 Mb length. The P. indicus assembly has 346 un-gapped contigs of 1 Mb or higher length. Again, a plot of the number of gaps (Figure 2C) in the finished assembly also tables P. indicus genome over other shrimp genomes on assembly quality. Interestingly, the chromosome-scale assembly presented for P. monodon was found to have high gap number in the finished assembly. As observed in the Figure 2C, the intermittent elevations in P. monodon line plot indicate the presence of high number of gaps in some scaffolds. The coding sequences also were observed to contain N’s in P. monodon assembly while none existed in the P. indicus assembly (Figure 2D).

The repeat elements as derived on the basis of the number of bases masked, constituted 49.31% (954 Mb) of the assembled genome (Table 2). One of the prominent features was a high proportion of simple sequence repeats (SSR) which spanned 31.99% of the genome. The proportion of SSRs reported here for P. indicus was found to be the highest amongst all sequenced genomes in the animal kingdom. The role of SSRs in adaptive evolution was recently demonstrated for shrimp (Yuan et al., 2021). Other major repeat classes include LINEs (5.8%) and low complexity regions (4.57%). The satellites (1.36%), LTR elements (0.31%), DNA transposons (0.2%), and small RNA (0.07%) were the other minor repeat families observed in the P. indicus genome assembly.

Combined evidence from ab initio gene prediction, Illumina RNAseq data, Pacbio Iso-Sequencing data and protein sequences from related species, identified 28,720 protein-coding genes in the P. indicus genome (Figure 3). The predicted protein-coding gene number was higher than P. chinensis (26,343) and P. vannamei (25,596) genomes but lower when compared to P. monodon (30,038) genome. The mean exon and intron lengths were 259 and 2315 bp, respectively. The longest gene, exon and intron lengths were 98,168, 14,941, and 76,392 bp, respectively (Table 3). Overall, 81.79% of the predicted genes had evidence from RNAseq data or IsoSeq data or proteins from related species. Functional annotations yielded results for 98.36% of the predicted genes using Interproscan, non-redundant protein database of Genbank and the UniProt database. For majority of the genes, the P. vannamei was the top hit species showing homology (Supplementary Figure 1).

The gene family analyses with protein sequences of 21 species including P. indicus identified 148 single-copy orthologous genes amongst them. Out of 399,313 genes subjected to the analyses, 81.75% (326,455) were clustered into 35,611 orthogroups and 18.25% (72,858) were singletons (Supplementary Figure 2). About 1,504 orthogroups were shared by all the 21 species. In P. indicus, of the 9595 orthologous gene families, maximum number were found sharing with P. monodon (8387) followed by P. vannamei (8255) and P. chinensis (7928). We found that 6,722 orthologous gene families were shared among the four shrimp species and 1,987 gene families only among them. The phylogenetic tree generated with sequences of single-copy orthologous genes depicted three distinct clades representative of Chelicerates, Crustaceans, and Hexapods (Figure 4).

Exploiting the pooled-sample RNA sequencing approach, the study reports 15,554 coding SNPs present in 3,965 different protein-coding genes in P. indicus genome (Supplementary Table 7). Minimum raw-read depth of 20, minimum evidence of 10 reads for each of the alleles and minimum phred quality of 100 were the criteria followed for short-listing the good quality SNPs. About one-third of the genes (n = 1185) had only one SNP and 283 genes had ≥10 SNP positions (Supplementary Figure 3). Majority of the SNPs were transversions (n = 13,655) than transitions (n = 1,899). Of 15,554 only 2,571 SNPs situated in 1,262 unique gene sequences were non-synonymous in nature contributing to amino acid polymorphic sites in the resulting proteins. Among the non-synonymous SNPs, majority were observed to be transversions (n = 2038) rather transitions (533). Of the total SNPs, 28% of transition and 15% of transversion substitutions were observed to be non-synonymous in nature (Table 4). As observed, a transition substitution has more chances of becoming non-synonymous than transversion. The PANTHER PSEP v1.01 (Tang and Thomas, 2016) tool classified 76 of these non-synonymous SNPs into probably benign (19), possibly damaging (11), and probably damaging (46) as shown in Supplementary Table 8 and Supplementary Figure 4. The tool was unable to determine the score for the remaining non-synonymous SNPs either due to a mismatch of the amino acid at the mentioned position between the query sequence and the panther family sequence or the absence of the panther family to the given query sequence.