Molecular characterization and expression profiling of transformer 2 and fruitless-like homologs in the black tiger shrimp, Penaeus monodon

Nucleotide sequences and deduced amino acid sequences analyses

Nucleotide sequences of P. monodon candidate sex determining genes PmOvtra 2, Pmfru-1 and Pmfru-2 were obtained from our transcriptome data prepared from P. monodon central nervous tissue (unpublished data). The nucleotide and deduced amino acid sequences of the putative sex determining genes were analyzed with Expasy Translate tool (Gasteiger et al., 2003) (https://web.expasy.org/translate/) for open reading frame (ORF), Compute pI/MW tool (https://web.expasy.org/compute_pi/) for molecular weight (MW) and isoelectric point (pI) and Simple Modular Architecture Research Tool GENOMES (Letunic, Khedkar & Bork, 2021) (http://smart.embl-heidelberg.de/smart/set_mode.cgi?GENOMIC=1) for functional protein domains. Homologous amino acid sequences in other crustaceans and arthropods (Tables 1 and 2) were obtained from NCBI GenBank database with literature and BLAST searches (https://blast.ncbi.nlm.nih.gov/). Selected species included in the analyses were those in evolutionary lineages leading to P. monodon and with putative orthologs of Pmtra 2, PmOvtra 2, Pmfru-1 and Pmfru-2 proteins, and their homologs/orthologs copies were either functionally characterized in other species (literature searches) or with highest sequence similarity to the corresponding homologs of interest in P. monodon (BLAST searches). Multiple sequence alignments of the homologous amino acid sequences of Tra 2 and Fru gene sets were separately performed with MAFFT version 7 web service (Katoh, Rozewicki & Yamada, 2019) with a default setting (https://mafft.cbrc.jp/), the best substitution model for each protein sequence set was obtained with ProtTest 3 (version 3.4.2 (Darriba et al., 2011)), and the phylogenetic trees were performed with RAxML (version 8.2.12; (Stamatakis, 2014)) with 1,000 bootstrapping. Tra 2 phylogenetic tree was constructed with a GAMMA model of rate heterogeneity and the WAG substitution model with empirical amino acid frequencies, but Fru phylogenetic tree was constructed with a GAMMA model of rate heterogeneity and the VT substitution model with empirical amino acid frequencies.

Animal and tissue sampling

Male and female P. monodon broodstocks were obtained from Shrimp Genetic Improvement Center (SGIC), Suratthani, Thailand. All shrimps were tested for specific pathogen infections, i.e., White spot syndrome virus (WSSV), Yellow head virus (YHV), Taura syndrome virus (TSV), Monodon baculovirus (MBV), P. monodon densovirus (PmDNV), Infectious hypodermal and hematopoietic necrosis virus (IHHNV), Decapod iridescent virus 1 (DIV1) and Vibrio parahaemolyticus. Shrimp were acclimatized in 1000 L artificial seawater (30–35 ppt) in fiberglass tanks with constant aeration. The water temperature was maintained at 28 ± 2 °C. Shrimp were fed with commercial pellets, three times a day. Shrimp were euthanized by ice water immersion before tissue collection. A variety of tissues, including eyestalk, nervous tissue, gill, hepatopancreas, heart, gastrointestinal tract, testis and ovary were collected from female and male shrimp. The tissue samples were immediately kept in Tri-Reagent^® (MRC, USA) and stored at −80 °C until used. In addition, various developmental stages of larvae, including nauplii, protozoea, mysis, and postlarvae, were obtained from SGIC. Samples were collected, washed with Phosphate buffered saline (PBS), and preserved in 1 ml of Tri-Reagent^® solution and stored at −80 °C until used. All animals experiment was approved by the Animal Ethics Committee, Faculty of Science, Mahidol University (MUSC63-019-527).

RNA extraction and cDNA synthesis

Total RNA was extracted from the tissue samples using Tri-Reagent^® according to the manufacturer’s instruction. Briefly, the tissue samples were homogenized in Tri-Reagent^® solution and incubated at room temperature for 5 min. Homogeneous mixture was subjected to chloroform extraction and the total RNA was precipitated with isopropanol. The RNA pellet was washed by 80% and 100% ethanol and finally RNA sample was suspended in DEPC treated sterile distilled water. Concentration of RNA was determined by Nanodrop spectrophotometer (ThermoScientific, USA). For cDNA synthesis, DNaseI was used to remove genomic DNA in the samples. One microgram of the total RNA was treated with 1 U DNaseI (Thermo Fisher Scientific, Denmark), and the reaction was performed at 37 °C for 1 h. The reaction was terminated by EDTA following the manufacturer’s protocol. The extracted RNA was added to reverse transcription reaction which consisted of Random hexamer (Invitrogen, USA) and SuperScript^® III Reverse Transcriptase (Invitrogen, USA). The reaction was performed according to the manufacturer’s instructions. All cDNA samples were stored at −80 °C until used.

Reverse transcription polymerase chain reaction (RT-PCR)

Shrimp tissues, including eyestalk, nervous tissue, gill, hepatopancreas, heart, gastrointestinal tract, testis and ovary were collected from six healthy adult shrimp. Pooled samples were prepared for RT-PCR analysis. Expression patterns of PmOvtra 2, Pmfru-1 and Pmfru-2 in various tissues of adult males and females were investigated by RT-PCR. The specific primers were designed and shown in Table 3, and 16s rRNA was used as housekeeping gene. Reverse transcription PCR master-mix included cDNA template, 10×PCR buffer, 10 µM of forward primer and reverse primer, 10 mM dNTP, 2.5 U of Taq DNA polymerase (Invitrogen, USA) and DEPC treated sterile distilled water. The amplifications were carried out using temperature cycling conditions as follows: 1 cycle of 94 °C for 5 mins; 35 cycles of 94 °C for 30s, 55 °C for 30 s and 72 °C for 45 s; 1 cycle of 72 °C for 10 min. The PCR products were resolved by 1.5% agarose gel electrophoresis.

Quantitative real-time PCR (qPCR)

Developmental stages of larvae, including nauplii, protozoea, mysis, and postlarvae were collected and approximately 30 individuals in the same stage were pooled as one sample. Ovaries and testes were collected from five adult females and five adult males, and each was run as individual sample. Three replicates were prepared for each sample. Quantitative real-time PCR was performed in 96-well plates in a final volume of 20 µl. Each reaction contained 2 µl of cDNA, 10 µl of KAPA SYBR^®FAST qPCR Master-Mix (2X) (KAPA Biosystems), 0.4 µl of specific forward and reverse primers (10 µM/ µl) (List of specific primers are shown in the Table 3), and 7.2 µl of DEPC treated sterile distilled water. The reaction was carried out in the 7500 Real Time PCR system (Applied Biosystems) with the following qPCR conditions: initial denaturing at 95 °C for 3 min; 40 cycles of 95 °C for 3 s, 60 °C for 30 s, and 72 °C for 30 s and finally, 1 cycle of 95 °C for 15 s, 60 °C for 1 min, 95 °C for 30 s, 60 °C for 15 s. The production of specific products was confirmed by performing a melting curve analysis of the samples. The reaction without cDNA was used as the negative control. The same qPCR profile was performed for 16s rRNA amplification which served as a reference gene. To confirm that 16s rRNA was qualified as a reference gene used, the stability of 16s rRNA, EF1- α, and β-actin, was evaluated in the samples. The result suggested that 16s rRNA would be one of the versatile reference genes for this assay (Supplement 3 and Supplement 4). Data were analyzed with 7500 software v.2.3 (Applied Biosystems). The baseline was set automatically by the software to maintain consistency. The comparative CT method (2^−ΔΔCT method) was performed for analyze the expression levels of PmOvtra 2, Pmfru-1 and Pmfru-2 (Livak & Schmittgen, 2001).

Statistical analysis

Data represented the mean ± standard deviation (S.D.). Statistical analysis was performed using Graph-Pad Prism 5 software (Motulsky, 2007). Statistical significance between two groups (testis and ovary) were determined using unpaired t-test. Significant difference was considered at P < 0.01. Statistical significance between groups of larvae developmental stages was determined using Tukey’s multiple comparison tests and one-way ANOVA. Significant difference was considered at p < 0.05.

Molecular identification of PmOvtra 2 and Pmfru

The full-length cDNA sequence of P. monodon ovary associated transformer 2 (PmOvtra 2) and two isoforms of P. monodon fruitless-like gene (Pmfru-1 and Pmfru-2) were obtained from our transcriptome database. PmOvtra 2 (GenBank Acc. No. MT543028) included 1,774 nt. with 110 bp 5′-untranslated region (UTR), a 744 bp open reading frame (ORF) and 920 bp 3′-UTR (Supplement 1). The putative amino acid sequences analyses of PmOvtra 2 showed that PmOvtra 2 was comprised of 247 amino acids with predicted molecular weight of 28.41 kDa and the theoretical isoelectric point was 11.34. PmOvtra 2 contained a predicted RRM domain and two RS regions that were highly conserved among Tra 2 proteins in several species. The predicted RRM domain of PmOvtra 2 was located at 101–180 aa (Fig. 1).

The full-length cDNA sequences of Pmfru-1 (GenBank Acc. No. MT497519) and Pmfru-2 (GenBank Acc. No. MT503286) were 1,306 and 1,858 bp, respectively. Pmfru-1 consisted of a 101 bp 5′-UTR, a 1,182 bp ORF and a 23 bp 3′-UTR. The cDNA sequence of Pmfru-2 consisted of a 92 bp 5′-UTR, a 1,437 bp ORF and a 329 bp 3′-UTR (Supplement 2). The putative amino acid sequences of Pmfru-1 comprised of 393 amino acids with predicted molecular weight of 43.93 kDa, whereas the Pmfru-2 comprised of 478 amino acids with predicted molecular weight of 52.41 kDa. The theoretical isoelectric points of Pmfru-1 and Pmfru-2 were 5.38 and 5.56, respectively. Translated amino acid sequences of both Pmfru-1 and Pmfru-2 contained BTB domain with zinc finger domain which were highly conserved in Fru proteins in several insect and crustacean species. The BTB domain and zinc finger domain of Pmfru-1 were located at 31-126 aa and 339-389 aa, respectively, while that of Pmfru-2 were at 31-126 aa and and 373-395 / 401-424 aa (Fig. 2).

Multiple sequence alignments and phylogenic analysis of PmOvtra 2, Pmfru-1 and Pmfru-2

Due to evolutionary complexity of Tra 2 and Fruitless protein families, we focused on analyzing amino acid sequences of Pmtra 2, PmOvtra 2, Pmfru-1 and Pmfru-2 in P. monodon and of their homologs/orthologs (see Tables 1 and 2) in other species. Other putative paralogous proteins annotated in a draft genome sequence of P. monodon were not included because in this study we sought to infer evolutionary scenarios of Pmtra 2, PmOvtra 2, Pmfru-1 and Pmfru-2 but not of all paralogs in these two protein families (see Materials and Methods).

Tra 2 phylogenetic tree from maximum likelihood analyses of Tra 2 amino acid sequences (Fig. 3) revealed two distinct sister subclades of PmOvtra 2 and Pmtra 2 with their corresponding homologous sequences in the other crustaceans, and these two subclades containing PmOvtra 2 and Pmtra 2 formed a clade that was separated from D. melanogaster homolog (bootstrap re-sampling value 100%). The PmOvtra 2 and Pmtra 2 topology in the Tra 2 phylogenetic tree (Fig. 3) suggests that they were paralogous and duplicated in the last common ancestor of crustaceans. We named the two sister subclades containing PmOvtra 2 and Pmtra 2 paralogs as the PmOvtra 2 and Pmtra 2 clades, respectively (Fig. 3). Based on Tra 2 multiple amino acid sequence alignment analyses, PmOvtra 2 had higher percentage amino acid identities to the other previously reported Tra 2 orthologs in crustaceans within the PmOvtra 2 clade (98.68%, 89.47%, 86.18% and 86.18% to P. chinensis, C. quadricarinatus, M. rosenbergii and M. nipponense, respectively) than to the previously characterized paralog Pmtra 2 (62.83%). Both PmOvtra 2 and Pmtra 2 had a high percentage amino acid identity (98.68% and 99.32%, respectively) to their putative ortholog in the P. vannamei draft genome. Similarly, for the conserved RRM domain (80 residues), PmOvtra 2 had a higher percentage identity to its orthologs in P. chinensis (98.75%), P. vannamei (98.75%) and C. quadricarinatus (92.50%) than to the paralog Pmtra 2 (70%). Note that similar percentage amino acid identities to the well-characterized insect D. melanogaster ortholog (264 residues) were observed for both PmOvtra 2 (51.97%) and Pmtra 2 (48.64%).

A maximum likelihood phylogenetic tree (Fig. 4) of Pmfru-1 and Pmfru-2 protein homologs (Table 2) placed both Pmfru-1 and Pmfru-2 with a copy in P. vannamei (bootstrap re-sampling value 100%) and in turn placed this clade of three proteins as a sister clade of another clade containing previously characterized Fru proteins in E. sinensis and insects, e.g., Fru in D. melanogaster (bootstrap re-sampling value 85%). Such a placement of Pmfru-1 and Pmfru-2 in the phylogenetic tree suggested that the genomic locus (see below) of Pmfru-1 and Pmfru-2 was orthologous to that of a putative Fru copy in P. vannamei and had an evolutionary relationship with other previously characterized fruitless in E. sinensis and insects. Based on Fru multiple amino acid sequence alignment analyses, Pmfru-1 and Pmfru-2 had 99.31% amino acid identity to putative P. vannamei Fru copy, but 26.3% to E. sinensis Fru copy, the only characterized copy in crustacean so far. Interestingly, Pmfru-1 and Pmfru-2 had 38.1–43.6% amino acid identities to those well-characterized copies in insects. The BTB domain sequences of Pmfru-1 and Pmfru-2 were identical (100% identity due to the identical nucleotide sequences in the region; see Supplementary alignment 2), and they were also identical to that of P. vannamei (Fig. 2). They had 46.80%, 43.61% and 41.48% identity with that of N. vitripennis, D. melanogaster and A. aegypti, respectively, but they had 27.65% to that of E. sinensis Fru copy. Contrary to the BTB domain, the ZF domain sequences of Pmfru-1 and Pmfru-2 shared only 26.31% identity. Pmfru-1 ZF domain sequences had a lower sequence identity to those of D. melanogaster, A. aegypti and N. vitripennis than did Pmfru-2 (28.94 vs. 31.57%, 28.94 vs. 31.57% and 23.68 vs. 31.57%, respectively). It should be noted that ZF domains were not predicted for the protein sequences of P. vannamei and E. sinensis.

Gene organization and alternative splicing of Pmfru-1 and Pmfru-2

Since nucleotide sequences of Pmfru-1 and Pmfru-2 share identical sequences between 11-1109 nt of Pmfru-1 and 2-1100 nt of Pmfru-2 (see Supplementary alignment 2), we sought to determine whether Pmfru-1 and Pmfru-2 are results of alternative splicing events by searching the nucleotide sequences against a draft genomic sequence of P. monodon. Both Pmfru-1 and Pmfru-2 were aligned by BLASTN to only a single scaffold sequence (Accession number NC_051415.1 or PmonScaffold_30) and formed seven interval regions as putative seven exons of 87, 95, 594, 243, 90, 761 and 197 bp in length, respectively (Fig. 5). Five out of the seven exons showed the identical sequences between sequences of Pmfru-1, Pmfru-2 and genomic scaffold sequences, whereas the other two (exons 3 and 6) each had one mismatched base-pair (≥ 99.7% identity) to the genomic scaffold sequence. At these exon-intron junctions, splice dinucleotide sequences of GT donor sites and AG acceptor sites are observed. The result suggested that Pmfru-1 and Pmfru-2 were likely produced by constitutive splicing events of exons 1 to 5 and alternative splicing events of exon 6 to Pmfru-2 mature mRNA and exon 7 to Pmfru-1 mature mRNA (i.e., mutually exclusive exons; Fig. 5). The start codons (ATG) of Pmfru-1 and Pmfru-2 were both located in exon 2, whereas the stop codons (TGA) of Pmfru-1 and Pmfru-2 were located in exon 7 and exon 6, respectively.

Temporal expression profiles of PmOvtra 2, Pmfru-1 and Pmfru-2

The temporal expressions of PmOvtra 2, Pmfru-1 and Pmfru-2 were determined by qPCR. The specific primers were designed and shown in Table 3. The relative expression levels of PmOvtra 2, Pmfru-1 and Pmfru-2 were determined in different developmental stages of larvae and postlarvae, including nauplius, protozoea, mysis, and postlarvae stage 1, 5, 10, and 15 (PL1, PL5, PL10, and PL15). PmOvtra 2, Pmfru-1 and Pmfru-2 mRNA expression levels were normalized to the 16s rRNA transcript level. The results showed that PmOvtra 2 was detected in all the developmental stages of P. monodon. Notably, PmOvtra 2 expression was relatively low in nauplius and protozoea stages, and then gradually increased from mysis to PL1. The PmOvtra 2 expression level reached the highest level in PL1 but it abruptly decreased from PL5 to PL15 (Fig. 6A).

The Pmfru-1 and Pmfru-2 transcripts were detected in all larval stages of P. monodon. The Pmfru-1 transcripts showed low expression levels at nauplii stage and slightly increased from protozoea to PL5. The Pmfru-1 expression level reached the highest level at PL1 and gradually decreased from PL10 to PL15 (Fig. 6B). The Pmfru-2 expression level maintained a low level from nauplius to mysis and then gradually increased at the postlarval stages. The Pmfru-2 expression level reached high level in PL1 and then gradually decreased at PL5. The Pmfru-2 transcripts also increased at PL10 and reached the highest levels at PL15 (Fig. 6C).

Tissue distribution of PmOvtra 2, Pmfru-1 and Pmfru-2 in adult P. monodon

Expression of PmOvtra 2, Pmfru-1, and Pmfru-2 transcripts were determined by RT-PCR in several tissues, including eyestalks, nervous tissues, gill, hepatopancreas, heart, gastrointestinal tract, and gonads. Expression of 16s rRNA, showing the band at 152 bp, was utilized as internal control. The expression of PmOvtra 2, which exhibited the specific band at 449 bp, showed in nervous tissues, gill, heart, gastrointestinal tract, testis and ovary (Fig. 7). Apparently, the expression of PmOvtra 2 was high in ovary compared with other tissues.

Expression of Pmfru-1 exhibited the specific band at 262 bp. Reverse transcription PCR results showed that Pmfru-1 was expressed in several tissue, i.e., eyestalk, nervous tissues, gill, hepatopancreas, heart, gastrointestinal tract, testis and ovary (Fig. 7). Notably, the band of Pmfru-1 expression was intense in gonadal tissues, both testis and ovary. It was suggested that the gonad tissue might be a major organ for Pmfru-1 mRNA expression. Regarding the expression of Pmfru-2, its expression exhibited the specific band at 345 bp. The result showed the expression of Pmfru-2 in eyestalks, nervous tissues, gill, hepatopancreas, heart, gastrointestinal tract, and gonad, with high preference for ovary (Fig. 7).

Since PmOvtra 2, Pmfru-1 and Pmfru-2 were highly expressed in gonadal tissues as shown by RT-PCR, the expression patterns of these genes in testes and ovaries were further quantified by qPCR to see whether their expressions were different between sexes of P. monodon. The expression patterns of PmOvtra 2, Pmfru-1, and Pmfru-2 in testes and ovaries were normalized to the 16s rRNA transcript level. Relative expressions of PmOvtra 2, Pmfru-1, and Pmfru-2 in the ovary were 26.37 ±10.64, 71.27 ±42.92, and 60.99 ±29.45, respectively. However, the relative expressions of these transcripts in the testis were much lower, i.e., 1.20 ±0.41, 0.99 ±0.26 and 1.40 ±0.68, respectively. The result indicated significantly higher expression of PmOvtra 2, Pmfru-1, and Pmfru-2 in the ovary compared with the testis (p <0.01) (Fig. 8).