The Central Inland Fisheries Research Institute (ICAR) ethical committee at Barrackpore, Kolkata 700120, West Bengal, India, authorized this work. The Committee on Animal Ethics approved this experiment. (Approval date: 10/09/2019; approval code: CIFRI/EC/2019/62) The Institute’s defined protocols were followed when handling the fish used in the experiment.

Healthy Labeo rohita (weight 70 ± 10g) were obtained from the ICAR-Central Inland Fisheries Research Institute hatchery in Barrackpore, Kolkata 700120, West Bengal, India, for the in-vivo gene expression experiments. The fish were kept in a 500-liter tank. Throughout the experiment, the water temperature ranged from 25 to 30 degrees Celsius, and its pH was 7.5–7.8 with continuous aeration and a 20% water change every alternate day. The fish were given commercial floating feed twice a day and allowed to acclimate for two weeks before the start of the experiment.

The entire experiment was conducted in two steps, one after the other. (1) The first phase involved cloning the LrSTAT1a gene and analysing its expression using immunostimulants (Poly I: C – a synthetic analogue of ds RNA, Sigma, St. Louis, MO, USA). (2) The second experiment was an in-silico investigation of the predicted protein of LrSTAT1a transcript.

Twenty- one fish was given intraperitoneal (IP) injections of 0.1 mL of poly I: C (Sigma, St. Louis, MO, USA) (1 mg/mL) PBS solution as stimulants, while another set of fish was kept as a control (Each day we collected the tissue samples from 6 nos. of fish where three nos. fish treated and 3nos. fish from control). All of the fish were grown in separate water tanks with appropriate care. The experiment was run for 14 days. Following fish anaesthesia and sacrifice, 20–50 mg of various tissues, including the brain, spleen, liver, intestine, heart, kidney, gills, and muscle, were separately collected at various time points (days 1, 2, 3, 4, 5, 7, and 14), and they were immediately stored in RNAlater (Amibion, USA) at -80°C for the LrSTAT1a expression study.

Total RNA was extracted from L. rohita samples using Trizol (Invitrogen, USA). cDNA template was synthesized from the total RNA using GoScript™ Reverse Transcription System (Promega, USA). The full-length cDNA of LrSTAT1a was obtained via RACE-PCR (rapid amplification of cDNA ends) using the SMARTer™ RACE cDNA Amplification Kit (Clonetech, USA). Specific primers ( Table 1 ) were designed according to the putative coding sequences obtained from the NCBI database. To the core cDNA fragments of LrSTAT1a primers, LrSTAT1aFP1 forward and LrSTAT1aRP1 reverse were used for the primary PCR. The PCR amplification was done under the following conditions: initial denaturation at 95°C for 5min followed by 40 cycles of 95°C/45sec, 56°C for 30sec and 72°C for 45sec and final extension at 72°C for 7min. To clone the 5′ fragments, nested PCR was performed using four sets of specific primers (LrSTAT1aFP2, LrSTAT1aRP2; LrSTATa1FP3, LrSTAT1aRP3; LrSTAT1aFP4, LrSTAT1aRP4; LrSTAT1aFP5, LrSTAT1aRP5). Similarly, primer pairs were used to produce the t3′ fragments. The PCR segments were sequenced by Barcode Bioscience Pvt. Ltd. (Bangalore, India) after being ligated to the pMD19-T vector.

According to Das et al. (2019) and Panda et al. (2023), the LrSTAT1a gene expression analysis using qRT-PCR was conducted using the Light Cycler^®480 detection system (Roche, Germany). A total of 0.5µl of 25 nanomolar for each forward and reverse primer, 2x asymmetrical cyanine dye, 5µl of SYBR Green, 1µl of master mixture (Roche Germany), 1mg/µl of cDNA and 3µl of nuclease-free water was used to make the qPCR reaction mixture. To amplify the 146 bp LrSTAT1a fragment, PCR reactions were set up with triplicate samples and 40 cycles of 95°C/10 sec, 55°C/10 sec, and 72°C/10 sec. Primary denaturation was then carried out for 10 min at 95°C. Using β-actin as a standard gene, the 2^-ΔΔCT method was used to determine the relative expression of the LrSTAT1a gene. Following PCR, the products were run on one-point-five percent (1.5%) agarose gel electrophoresis to verify the desired amplified product’s quality and gene expression. The housekeeping gene β-actin was used in this method to normalize fold change compared to the control gene. The average value added along with the estimated standard error was how the qRT-PCR data were expressed. The efficiency of the ladder (E) and the Ct deviation of the analysed amplicons with the corresponding reference gene were used to determine the relative expression (R) of the LrSTAT1a transcript.

A qRT-PCR was carried out to determine the relative expression level of LrSTAT1a in various tissues after Poly I: C induction. Fish injected with PBS were used as the negative control, whereas fish treated with Poly I: C was used as the positive control. The brain and kidney were chosen because they are considered the target organs for research on viral contagion.

The statistical analysis tool SPSS was used to examine all the data (version 16.0 SPSS, USA). The differences in significance between the various treatment groups were compared at the P≤ 0.05 percent significant level using a one-way ANOVA followed by Duncan’s multiple range tests.

According to Das et al. (2019), RACE was performed using the GeneRacer ™ method (Invitrogen, USA) to retrieve the complete LrSTAT1a cDNA sequence. RACE involves RNA ligation, mRNA capping to modify the 5´ end of the mRNA, and precise RNA oligo selection using the T4 RNA ligase enzyme. CIP was combined with total RNA to remove 5’-phosphates. Additional truncation of mRNA and non-mRNA occurred without altering the genetic composition of capped and full-length mRNA. TAP was applied to capped (dephosphorylated) RNA to remove the 5’cap configuration from the mRNA. Using GeneRacer™ RNA oligo, this procedure guaranteed the existence of 5’-phosphate, which is a necessary precursor for ligation. The 26-base oligo GeneRacer™ 3′ primer had the sequence 5′-GCTGTCAACGATACGCTACGTAA CG-3′ with a Tm of 78°C. A 26-base oligo with the sequence 5′-GGACACTGACATGGACTGAAGGA GTAGTA-3′ with a Tm of 78°C was the GeneRacer™ 5′ nested primer.

There were multiple phases involved in the RACE procedure using a kit from Invitrogen, USA. The phenol-chloroform technique was used to precipitate and purify the RNA before enzymatic treatment. The RNA was decapped on the first day, and the ligated RNA was stored in ethanol for precipitation overnight. The PCR product was cloned using TOPO^® on the second day, followed by reverse transcriptase PCR and sequencing of plasmid DNA containing either the 3′ or 5′ PCR product.

Primers for the 5’ and 3’ RACE PCR, however, were created for the corresponding end sequences of the first cloned 360 bp (LrSTAT1a-FP1 and LrSTAT1a-RP1) and 146 bp cDNA fragments prior to performing RACE PCR ( Table 1 ) methods.

After the full-length LrSTAT1a, CDS were generated from the sequence contigs and analysed through various bioinformatics analyses. LrSTAT1a amino acid sequences were inferred by using the ExPASy website (https://web.expasy.org/translate/). The L. rohita STAT1a gene’s open reading frames (ORFs) were predicted using the NCBI ORF finder, and the longest ORF was selected for further investigation. Using the ProtParam tool of ExPaSy, the physicochemical characteristics of the encoded STAT1a protein were ascertained. The molecular weight (MW) and isoelectric point (pi) of STAT1a protein were calculated by Expasy tool. The grand average of hydropathicity (GRAVY) was calculated by ProtParam of Expasy tool with protein sequence of STAT1a gene. The motifs of STAT1a protein were identified by motif tool of Genome Net. The promoter region of STAT1a gene was extracted from the NCBI database by taking 1000 bp of genomic DNA sequence upstream to the STAT1a gene. Several TATA boxes were identified by OProf of Expasy tool in predicted promoter region of STAT1a gene. Using the MEGA version 11 tool and the multiple sequence alignment of the STAT1a proteins, the neighbour-joining (NJ) tree was inferred ( Table 2 ; Figure 1 ). The bootstrap test with 1000 repeats was used to assess the reliability of the phylogenetic tree.

We used the Robetta web server (https://robetta.bakerlab.org/) to generate the initial model structure using RoseTTA Fold as the algorithm.

PROCHECK (https://saves.mbi.ucla.edu/) was used to validate the protein models. It helps create the Ramachandran plot and check for the presence of amino acid residues in both allowed and banned regions in a descriptive way.

The ligand (SH2) SDF (structure data file format) confirmation was obtained from PubChem (https://pubchem.ncbi.nlm.nih.gov/), and SMILES Translator (https://www2.chemie.uni-erlangen.de/services/translate/) was used to translate it to PDB format. It was then run via the energy reduction procedure of the Avogadro platform (https://avogadro.cc/). For each of these compounds, the MMFF94 force field and the steepest descents were used.

The virtual screening tool Python prescription (PyRx) 0.8 was used to perform molecular docking research using Autodocking vina (22). Positions on the X-, Y-, and Z-axis of the grid were first specified. The centre of the grid was positioned in the active site pocket. In PDBQT format, receptor and ligand files were analysed. The default docking algorithms followed the protocol. Individual docking experiments were performed for each ligand-protein combination (23). The findings were presented by increasing docking energies, The representative value was determined by taking the lowest binding energy of each cluster (24).

By using Desmond [D. E. Shaw Research Institute (https://www.deshawresearch.com/resources.html)], Molecular dynamic (MD) simulation of the protein complexes was done. We used the Protein Preparation Wizard to create receptor-ligand complexes. The protein-ligand complexes with the highest scores were simulated using Desmond. A 10 Å box wall distance was used to solvate each protein-ligand combination in an orthorhombic SPC (single point charge) water model. The system was neutralized by adding the correct number of counter ions and a salt concentration of 0.15M. Acetyl and methyl amide were used to cap the STAT1a receptor, while the titratable amino acid residues were left in their dominant condition at pH 7.0. For protein interactions, this system used the OPLS-4 force field. Following the removal of overlapping water molecules, the systems were neutralized using Na+ ions. After the complex system was relaxed with the sharpest downward energy decrease, it was subjected to progressive heating to 310 K while being constrained. The system was simulated with the Berendsen NVT ensemble and kept at 310 K to keep the heavy atoms on the solute. The MDS was conducted under isothermal isobaric ensemble conditions (300 K, 1 atm pressure, and 200 ps thermostat relaxation period). The pressure and temperature scales were regulated using the semi-isotropic Parrinello-Rahman barostat and the Nose-Hoover thermostat for 100 ns at 300 K and 1 atm, respectively. Every 50 picoseconds, step-by-step MD simulation images were captured. Post-MD analysis, which included intermolecular interaction studies, was conducted on the extracted MD trajectories. This analysis encompassed dynamical stability and flexibility.

The results of the BLAST analysis suggested the existence of the LrSTAT1a gene’s present protein sequence. The LrSTAT1a gene had its maximum point of uniqueness with a segment of Sinocyclocheilus rhinocerous STAT1 (92%) as followed by Sinocyclocheilus grahami (90%)>, Cyprinus carpio (90%)>, Carassius auratus (90%).

Constitutive expression of STAT1a was deliberated in the brain, spleen, liver, intestine, heart, kidney, gill, and muscle of L. rohita for fourteen days post-induction. By using the semi-quantitative PCR approach, it was possible to determine that the brain, muscle, kidney, and intestine had the highest levels of STAT1a expression on day two after Poly I: C induction ( Figure 2 ), which was confirmed in qPCR. The expression decreased after third day. While the mean STAT1a: β-actin ratio in the treated groups was after Poly I: C injection, the mean ratio on the second day was recorded to be approximately 12.19 and decreased to 1.40 between day 5 and day 14. Many of the untreated control groups showed no LrSTAT1a mRNA expression. However, based on mRNA configuration, statistical analysis showed a positive difference (p<0.05) between fish injected with poly I: C and control fish starting on day two and continuing until day four.

The expression of the LrSTAT1a gene in the brain tissues on day 2 (6.1fold) and then rapidly decreased on day 7 (2.7-fold change at lowest) and was nearly undetectable on day 14. LrSTAT1a expression exhibited a decreasing pattern on the spleen. The expression of LrSTAT1a in the liver tissues followed a similar trend. On day 2 (2) day 3 (1.9), and the heart tissue expression pattern peaked up to day 4, then it gradually declined until day 14 ( Figure 3 ).

Poly I: C encouraged STAT1a to express itself steadily and well; it peaked on day 2 and returned to normal on day 14. A semi-quantitative PCR was performed after poly I: C induction to determine the intensity of STAT1a expression in each of the vital organs (brain, spleen, liver, intestine, heart, kidney, gill, and muscle). The development of a quantitative real-time PCR is underway to assess the level of STAT1a expression in the key organs after Poly I: C induction. For the duration of the investigation, fish induced by poly I: C were taken into account, and fish injected with phosphate buffer saline (PBS) served as the baseline treatment control group. The results of the semi-quantitative real-time PCR reaction matched the expression of STAT1a transcript analysis by RT-PCR across all tissues. Fish that were given PBS as a control showed less STAT1a expression. On day 2, STAT1a expression was remarkably prevalent throughout all tissues, but maximal up-regulation of STAT1a expression was observed in the kidney, brain, and muscle, whereas the spleen, liver, heart, and gill also showed the greatest STAT1a expression, respectively ( Figure 4 ).

After the unknown ends of the cDNA were cloned using RACE and sequencing of rapid amplification of the cDNA ends result, the entire sequence of 3238 bp was acquired ( Figure 5 ). The molecular weight (MW) and isoelectric point (pi) of STAT1a protein calculated by the Expasy tool were 123144.46 Da and 8.32, respectively. The grand average of hydropathicity (GRAVY) value of STAT1a protein calculated by ProtParam was -0.101. Nine motifs were found in the STAT1a protein ( Figure 6 ). There were several TATA boxes (length 15) found in the predicted promoter site of STAT1a gene ( Figure 7 ). This was submitted to the NCBI GenBank with Accession Number OQ868190. An overlapping PCR fragment was spliced and combined to determine the whole sequence. Through ATG as starting codon at 145 base pairs and TGA as a termination signal at 2310 base pairs, an untranslated 5’ primary region of 42 base pairs and an untranslated 3’ region with 318 base pairs, STAT1a transcript exhibited 144 to 5’-UTR, an ORF possessing 2166 base pairs, the encoding with 721 amino acids protein of 928 base pairs, to 3’-UTR, which contains stop signals and a characteristic poly-adenylation signalling location (AATAAA) and In the Figure 8 showing multiple alignment of STAT1a amino acid sequence from different species (i.e., Homo-sapiens, fishes and birds).

The STAT1a protein is a member of the STAT family of proteins, which are transcription factors that play important roles in signal transduction and cellular communication. STAT1a is activated by phosphorylation by Janus kinases (JAKs) in response to various cytokines and growth factors. Once activated, STAT1a dimerizes and translocated to the nucleus, where it binds to specific DNA sequences to regulate gene expression. The 3D structure of the STAT1a protein ( Figure 9 ) is typically shown as a ribbon diagram, where each amino acid residue is represented by a coloured ribbon. The structure can be divided into several domains. The N-terminal domain typically contains the DNA binding site and is often located at the N-terminus (beginning) of the protein. In the ribbon diagram, it would likely be represented by a distinct section, potentially with a different colour than the neighbouring domains.SH2 domain: This domain is responsible for protein-protein interactions and is involved in the dimerization of STAT1a molecules. It is usually located near the N-terminus but not at the very beginning. Look for a compact, folded region in the ribbon structure. WHSH domain: This domain is also involved in protein-protein interactions and is thought to play a role in the regulation of STAT1a activity. It is often located between the SH2 domain and the C-terminal domain. C-terminal domain: This domain contains the phosphorylation site for JAKs and is responsible for the protein’s activation. It is typically located at the C-terminus (end) of the protein and might be another distinct section with a different colour in the ribbon diagram. It is important to note that the specific location and appearance of these domains can vary depending on the protein and the way it is visualized.

The Ramachandran plot is a valuable tool for protein structure analysis and validation. It can be used to identify errors in protein models and to assess the quality of protein structures. The allowed regions of the Ramachandran plot are shown in the image as the green and yellow shaded areas. These regions represent the combinations of φ and ψ angles that are sterically allowed, meaning that there are no steric clashes between atoms in the protein. The disallowed regions, which are not shaded, represent combinations of angles that are not possible due to steric clashes.

The Ramachandran plot ( Figure 10 ) of the model indicated that 91.8% of the residues were located in the core region, 6.6% in the allowed regions, 1.1% in the generously allowed regions, and 0.5% Residues in the disallowed region. Thus, the stereo chemically unstable areas contain ten residues. Thus, additional energy minimization and refinement were applied to these residues. 91.8% of the residues were discovered to be located in the core region after the final model was analysed. The overall shape of the projected structure was stabilized and optimized by the displacement of residues from the prohibited zone to the allowed region.

The binding energy between SH2 and STAT1a was -9.7 Kcal/mol. The amino acid residues SER107, GLN530, SER583, LYS584, MET103, and ALA106 form a typical H-bonding with bond lengths of 2.36 Å, 2.17 Å, 2.57 Å, 2.41 Å, 2.19 Å, and 2.47 Å, in that order. Hydrophobic interaction involved residues VAL50, VAL99, MET103, and ALA106 with bond lengths of 5.40 Å, 5.17 Å, 4.78 Å, and 5.40 Å, respectively. Two Halogen (Fluorine) interactions MET103 and ALA106 with bond lengths of 3.37 Å and 3.62 Å, and an Electrostatic interaction (ASP522) with a bond length of 4.98 Å were found by docked complex analysis. Figure 11 displays the binding location of STAT1a proteins, whereas Figures 12 , 13 depict the docked complex’s 3D and 2D interactions (STAT1a SH2 Domain).

While molecular docking is a quick and efficient method for determining the binding position of a ligand with the active site of a protein, it does not capture the conformational changes that occur during ligand-protein interactions. A molecular dynamics (MD) simulation was performed to evaluate conformational changes more accurately (25). One hundred nanosecond molecular dynamic simulations were executed to evaluate the flexibility, stability, and system stability of the bound ligand SH2. The molecular dynamic stability of each complex was investigated by computing the backbone root mean square deviation (RMSD) of the protein relative to its initial conformation. All systems reached equilibrium after 60 nanoseconds of simulation relative to the reference frame at time 0 nanoseconds, as shown in Figure 14 .

The left Y-axis of Figure 14 showed the protein RMSD, whereas the right Y-axis showed the ligand RMSD profile aligned on the protein backbone. All frames from the 100 ns trajectory are aligned on the reference frame backbone. Within the initial two nanoseconds of the trajectory, the variation in ligand RMSD was observed. Nevertheless, the modifications centered around 0.15 Å, suggesting that the complex had not experienced substantial alterations to its structure. In addition, fluctuations in ligand RMSD on the order of 4.5 Å that persist for a duration of 40 ns are initiated by the simulation at 5 ns. Until the conclusion of the 100-ns simulation, the ligand RMSD remained constant after 40 ns. Upon concluding the investigation, a slight divergence of 70 nanoseconds was identified. The fluctuation is considered insignificant since it falls within the permissible range of 1 to 10 Å.

Significant local changes in the RMSF ( Figure 15 ) pattern of C-atoms were observed in the STAT1a receptor systems following phytocompound binding, as compared to the other receptor systems. At the N-terminal end of each system, the amino acid residues implicated in interaction with the compounds exhibited significant changes (RMSF increases), indicating their involvement in ligand recognition. Amino acids varied between 1.5 Å and 10.5 Å in the RMSF of the STAT1a receptor system, but it steadied and remained consistent between 1.5 Å to 6.0 Å Figure 15 . Residues of the binding site of the docked protein and ligand are illustrated in the schematic interaction diagram in Figures 16 , 17 . The stability of the ligand-protein complex is shown to be affected by interactions involving hydrogen, hydrophobic groups, and ions. Ligand binding requires the formation of H-bonds. Hydrogen-bonding properties greatly influence drug selectivity, metabolization, and adsorption, making them crucial for drug design. The STAT1a complex interactions are shown in Figure 16 as normalized stacked bar charts over a 100 ns trajectory. Water bridges (TH518, THR519, ALA521, ASP522, LYS523, GLU526, GLU528, ALN529, GLN530, GLY531, PRO547, ASN548, LEU582, SER583, LYS584, GLU585, ARG586, ARG602, GLU605), ionic bond interactions (ASP522, GLU528, ALA529), hydrogen bonds (THR518, THR519, ASP522, LYS523, GLU528, GLN530, SER583, LYS584, GLU585), and hydrophobic interactions (VAL50, VAL99, MET103) are the many forms of these connections. Hydrogen bonds, hydrophobic contacts, and water bridges were the most prominent during the MD simulation, as seen in the bar chart. Figure 17 shows the interaction with the stat1a complex, namely molecular interactions with residues ASP522, SER583, which were constant during the 100 ns simulation.